J Moulin
email: johan.moulin@u-psud.friman.shahosseini@u-psud.frfrancisco.alves@lgep.supelec.fr

I Shahosseini

F Alves

F Mazaleyrat
email: frederic.mazaleyrat@satie.ens-cachan.fr

Ultrasoft Finemet thin films for magneto-impedance microsensors

FeCuNbSiB thin films have been deposited using RF sputtering. Characterizations have shown that oxygen contamination and residual stress are mainly responsible for magnetic hardening. The sputtering and annealing conditions have been optimized and films with coercive field as low as 10 A m -1 (0.125 Oe) have been achieved. In addition, the influence of film thickness on the magnetic properties has been studied. Thus, magnetic field microsensors based on the magneto-impedance effect have been fabricated by stacking up Finemet/copper/Finemet films. The highest sensitivity (4000 V/T/A) is reached for 750 nm thick films. It is in the same range as cm-sized macroscopic devices realized using 20 µm thick ribbons.

Introduction

Since the discovery of the magneto-impedance (MI) effect in amorphous microwires [1] and films [2], magnetic field sensors based on this property have become competitive with other techniques like GMR and fluxgates. Two domains of application are emerging: Non Destructive Testing using eddy current for the detection of microdefects in conductive sheets [3] and the detection of magnetic nanoparticules for biomedical applications [4,5]. NDT involves the measurement of the stray field generated by eddy currents induced in the conductive sheet to be controlled. A magnetic cartography allows detecting the deformation of the stray field induced by the presence of a defect in the material. Classical techniques involve pick-up coils. However, for detection of deep defect (typ. 10 mm deep), the frequency of the eddy current (thus of the stray field) is low (some tens of Hz). In addition, the stray field is very weak. In this case, pick-up coils are not suitable and GMI sensor are good candidates. The MI effect finds its origin in the variation of the skin depth δ with an applied magnetic field Hext (through the magnetic permeability) while exciting the conductor using a high frequency current:

( )

e ext t ext H f (H ) ρ δ = π µ (1)
with f the excitation frequency, ρe the electrical resistivity and µt the transversal (or rotational in case of a wire) magnetic permeability.

The impedance of the conductor depends on its geometry:

( ) ( )

0 dc 1 J ka ka Z R 2 J ka = ( 2
)
for a wire [START_REF] Landau | Electrodynamics of continuous media[END_REF], with J0 and J1 the Bessel function, Rdc the DC resistance of the wire, a the wire radius and k = (1-j)/δ , or :

dc ka Z R coth(ka) 2 = ( 3
)
for a thin film [START_REF] Kittel | [END_REF], where a is the film thickness. For both geometries, the calculation of the derivative of the modulus of Z with the value a/δ shows that the variations of Z are the most sensitive to the one of a/δ for a/δ ≈ 2, i.e. a ≈ 2δ . It corresponds to an optimal excitation frequency for a given geometry. The highest MI ratios are obtained in amorphous Febased microwires, in which external and internal stresses are respectively circumferential and axial. This stress orientation corresponds to an optimal magnetic structure. However, collective fabrication, for instance in case of sensors arrays, is much easier while using thin films technique. It has been found that for a planar structure the MI variations can increase strongly by inserting a conductive film between the ferromagnetic layers [8]. In the case on large plates allowing neglecting border effect, the impedance on the trilayer structure can be expressed as:

( ) ( ) ( ) ( ) ( ) 1 1 2 2 1 1 2 2 1 1 2 2 coth coth coth coth DC k d k d Z R k d k d k d k d ζ ζ ζ + = + + (4)
( )

2 1 2 1 2 1 1 ² 1, 2 e e t i e e i j k i ρ δ ρ µ ζ ρ δ ρ δ - = = = =
with 2d1 and d2 the thicknesses of the ferromagnetic and conductive layers respectively [9]. In this multilayer structure, the optimal excitation frequency for obtaining (∆Z/Zsat)max is much lower than for a single film, as the conductivity of common ferromagnetic materials is about 100 times lower than that of a good conductor. Previous calculations using finite elements method on such a structure [10] have shown that the highest MI effect is obtained around 10 MHz for a 1 µm thick ferromagnetic film with a relative permeability of 4000 with a thickness ratio of 1/2/1 for the multilayer structure. In addition, if the magnetic domains are oriented perpendicularly to the applied magnetic field, the MI rises because of the increase in transverse permeability in the ferromagnetic layers. For that purpose, the ferromagnetic layers were annealed under transverse magnetic anisotropy field. The sensor characteristics, in particular hysteresis and sensitivity, are firmly dependent on the ferromagnetic layer softness. Figure 1 presents numerical calculations using finite elements for a G500 sensor configuration (see table 1 for dimensions). The MI increases quite linearly with the transverse permeability, which shows the importance of depositing a soft material. Several works have succeeded in integrating Finemet® alloys into thin film MI sensors [11][12][13][14]. Indeed, this material presents ultrasoft properties due to its microstructure: after a heat treatment of the amorphous precursor, ferromagnetic Fe-Si nanograins nucleate in a ferromagnetic matrix which remains amorphous. This ensures quasi-null magnetostrictive anisotropy by the means of negative and positive magnetoelastic effects balance between grains and matrix. Moreover, the magnetocristalline anisotropy in such nanocristalline materials is decreased following the Random Anisotropy Model proposed by Herzer [15]. The MI results obtained with such a material are promising, despite the fact that the developed sensors present lower MI ratios than those with FeCo-based alloys. However, MI ratio is not the only quality determinant factor for a sensor, variation of the MI with the magnetic field is more relevant. For this purpose, ultra-soft Finemet, which can be obtained by an optimization of the elaboration process, is a good candidate.

Sensors elaboration process

We developed a MI microsensor fabrication process compatible with thicknesses close to 1 µm. Each layer of the sensors was patterned using lift-off technique with a 3 µm thick AZ5214 photoresist (see figure 2). In order to obtain such thick resist, the rotation speed during spin coating was limited to 900 rpm. The lithography exposure dose and development time were also adapted to such a thickness. The materials were deposited using sputtering technique: Finemet® and copper for the ferromagnetic and conductive layers respectively. The Finemet has been deposited on a SiO2/Si substrate by RF sputtering in Ar plasma under the following conditions: residual pressure in the range of 10 -7 mbar, working pressure 40 mbar and RF power 250 W. The copper layer was realized by DC sputtering with a working pressure of 6.5 mbar and a power of 50 W. Before each sputtering, a 10 nm thick layer of titanium was deposited in order to increase adhesion. The sensors are constituted of rectangular stripes whose dimensions vary from 50 to 430 µm in width and from 1 to 5 mm in length. For our samples, we attempt to obtain a central copper layer totally shielded by ferromagnetic layers, i.e. the Finemet layers are wider than the copper track. However, the latter is longer and is ended by two connection pads for bonding.

Finemet elaboration process

Sputtering conditions

In order to obtain the lowest coercivity, the magnetic properties of the Finemet® films have been optimized in terms of elaboration and heat treatment. Former study has showed that the coercivity is mainly influenced by the oxygen contamination and the residual stresses in the films [16]. The former is correlated to the deposition conditions. In practice, the deposition time has been minimized in order to limit oxygen contamination. This implies a short distance between the target and the sample, as well as a high argon pressure. In addition, a liquid nitrogen trap is needed. Fig. 3 shows that the coercive field follows quite linearly the content in oxygen in the film. The reader can refer to [17] for details on the measurement of the composition. The stress in the films has two origins: the atomic disorder induced by the deposition and the differential thermal expansion between the substrate and the film (difference between room and process temperature). The measurement of the wafer curvature using a laser scanning yields an estimation of the internal stress of a 2" continuous film as a function of temperature (heating and cooling rates are about 2.5°C/min). Figure 4 presents the results for different final temperatures: 200, 300, 400 and 500°C. The stress induced by the deposition is high (around 550 MPa for a 2" continuous film) due to particular sputtering con- ditions and does not depend on the film thickness (from 200 nm to 2 µm). It can be reduced by annealing only. MOKE characterizations on large samples (several cm 2 ) have revealed planar isotropic behavior and gave an estimation of the coercive field of 4.6 kA m -1 , what is consistent with large stress in the film. In case of non-continuous films, it is supposed that the residual stress due to sputtering is drastically reduced. The films have been patterned by mean of shadow masking for magnetic characterization using AGFM or VSM (2 mm × 2 mm squares). The coercive field (measured using AGFM or VSM) is in that case lowered down to 60 A m -1 . This result shows that the patterned film is subjected to lower stress than its continuous counterpart. It is however difficult to estimate the stress in the pattern, except by indirect method. Considering that the stress is the only responsible for the increasing of the coercive field, one can estimate that the stress is reduced down to several MPa after patterning. As this state, the material is subjected to secondary process conditions, like the initial temperature of the substrate (sample is not cooled in our equipment). Indeed, the differential thermal expansion between the substrate and the film during the cooling after the deposition is responsible for stress:

( ) ( ) ( ) exp 1 f therm w f f E T T σ µ µ α υ = - ∆ = ∆ - (5)
with µw, µf, Ef, νf and ΔT respectively the thermal expansion coefficients of the substrate and the film, Young's modulus and Poisson's ratio of the film and the temperature variation involved in the thermal expansion process. This stress influences strongly the magnetic properties of the films. The figure 5 shows that the coercive field of 2 mm x 2 mm squares deposited at the same residual vacuum varies quite linearly with the initial temperature. This phenomenon can be related to an increasing of the global anisotropy induced by magnetoelastic effect. The magnetoelastic anisotropy associated to an amorphous material is:

s s 3 E 2 = λ σ (6)
with λs the saturation magnetostriction coefficient and σ the global stress in the film. In case the thermal expansion only is responsible of stress variations, the magnetoelastic energy thus the coercive field varies linearly with the difference of temperature during sputtering and ambient temperature. Indeed, in the case of randomly oriented stress without anisotropy of form or crystal [18]:

s c S E H 0.64 J = (7)
with Js the saturation magnetization of the material. Finally, one obtains a linear relation between the coercive field and the difference between the initial temperature of the substrate and the room temperature:

( )

s c S H 0.96 T J λ = α ∆ (8)
what has been experimentally measured.

Figure 4 presents the results for different final temperatures: 200, 300, 400 and 500°C. From this work, the coefficient α in Eq. ( 8) has been estimated about 0.75 MPa/K while the temperature remains below 200°C. Indeed, in this range of temperature, the heating and cooling curves are superposed thus thermal expansion only is responsible for stress variations in the film. Using the magnetostriction coefficient and the saturation magnetization of such alloys (i.e. λs = 23 10 -6 and Js = 1.2 T), one can estimate the increasing of the coercive field due to thermal expansion in the range of 13.6 A m -1 /K. This is coherent with the experimental results that present a slope of 13 A m -1 /K for as-deposited films (see Fig. 5).

After an annealing at 300°C for 1h, the dependence of the coercive field on the initial substrate temperature is also linear, what reflects a similar behavior, except that the slope is reduced to 5 A m -1 /K. This result is associated to a reduction of λs due to a modification of the microstructure.

Annealing

The remaining stress has to be minimized in order to obtain softer films. This is one of the roles of postannealing.

The figure 6 presents the hysteresis loops of 500 nm films annealed at various temperatures for 1h. The coercive field decreases strongly after annealing at 200°C. Above this temperature, the coercive field decreases slightly, down to 10 A m -1 for annealing at 400°C. At 500°C, The coercive field drops to 15 A m -1 . Similar results have been obtained with a minimum coercive field at 280-330, 380 and 450°C for 50, 1000 and 2000 nm respectively [19][20][21]. Compare to these works, the characteristic temperatures of the film have been precisely defined using the measurement of the stress in the film as a function of the temperature, during heating and cooling. From Fig. 5, two particular temperatures can be extracted. The glassy transition temperature Tg corresponds to the limit with the supercooled liquid state in which the effect of the differential thermal expansion is reduced. Thus Tg corresponds to the first slope changing in the heating curve: 300°C for a 500 nm thick film. The cooling curve shows that the mechanical properties (Young's modulus and Poisson ratio) of the material are not changed while the final temperature remains below 400°C. The second slope breaking corresponds to an evolution of the microstructure leading to a new Young modulus value, i.e. crystallization. It is to notice that the crystallization temperature decreases with the film thickness [16], from 480 to 450 and 400 °C for 20, 2 and 0.5 µm respectively. These results are coherent with values measured using other technics in [19][20][21].

The magnetic behavior is partially correlated to these measurements. Indeed, the minimum coercive field is reached for an annealing temperature equal to the crystallization temperature Tx. At this state, stresses are released and nanograins of Fe-Si have nucleated but crystalline volume fraction remains small. The obtained nanostructure leads to an averaging out of the magnetocrystalline anisotropy, according to the Random Anisotropy Model [22]. Above Tx, the modification of the mechanical properties (especially Young's modulus) induces large residual stress due to differential thermal expansion during the cooling, thus to an increasing of the coercivity. However, the magnetic behavior of samples annealed at 200°C can not be explained using stress measurement of a continuous film. Indeed, the drastic decreasing of the coercive field is not correlated to an obvious modification of the mechanical properties of the film between room temperature and 200°C.

It is to notice that the relative permeability, as well as the saturation polarization are not affected by the annealing and remains constant with values of 1.2 T and 3400 respectively.

Field induced anisotropy

As reported in other works [23], the MI sensor sensitivity is increased after inducing an anisotropy perpendicular to the current direction in the ferromagnetic material. This has been performed by annealing in the magnetic field created by a 3 mm thick and 5 cm diameter magnetic circuit containing AlNiCo magnets. The magnetic field in the 7 mm × 14 mm air gap where samples were set out is in the range of 50 kA m -1 . The AlNiCo magnets allow working up to 500°C in the secondary vacuum of the sputtering chamber. Figure 7 depicts the hysteresis loops of 2 mm × 2 mm × 500 nm samples realized by shadow masking and characterized using AGFM. The field induced anisotropy leads low coercive field, as well as annealing without magnetic field: 10 A m -1 for zero and // field annealing and 15 A m -1 for ⊥ field annealing. The induced anisotropy can be estimated from loops measured in both the parallel and perpendicular directions of the annealing magnetic field. Its value is in the range of 150 J m -3 , much larger than that induced in bulk materials by this technique. In addition, the anisotropy field can be estimated to 4350, 300 and 700 A m -1 for zero, // and ⊥ field annealing respectively. The induced energy is correlated to a variation of the permeability: 2300, 3000 and 1300 for zero, // and ⊥ field annealing respectively.

After optimization of the ferromagnetic material properties, sensors with variable sizes and thicknesses have been processed as listed in table 1.

Sensors characteristics

MI measurements have been carried out using a voltage source and a lock-in amplifier in the 1-40 MHz frequency range. The DC magnetic field measured by the sensor was created by a set of Helmholtz coils, providing a 2 cm wide homogeneous field zone. The maximum magnetic field applied to the sensor is close to 7000 A m -1 .

The MI ratio variations are weak compared to other type of sensors, in the range of 10 %. Indeed, the resistance value is large (around 10 Ω) and varies little. By opposition, the inductance variations reach 400 % but the inductance value is low. These measurements highligth magnetic properties like anisotropy field, which is in the range of 700 A m -1 for sensors G500 and G750. The result is coherent with the value estimated from the hysteresis loop of the films annealed under transverse magnetic field. Papers related to MI effect generally compare the materials properties using the MI ratio only. However, this is not the most relevant factor to quantify a sensors quality. The most important figure is either the variation of the MI ratio with the magnetic field (in %/Oe) or the variation of the impedance with the magnetic field (in Ω/T or in V/T/A). Furthermore, measurements at high frequency are strongly sensitive to external perturbations, so the impedance value can be unreliable. To solve this problem, the sensors were characterized by measuring the output voltage as a function of the magnetic field, starting from the saturated state. Then the sensitivity was extracted by numerical derivation of the voltage by the magnetic field.

The characterizations have highlighted a strong dependence of MI properties on the shielding of the central copper track by the ferromagnetic material. Former investigation has shown that defects at the edge of the deposited layer appear while the sputtering direction is not vertical enough. This phenomenon is related to the deposition of the sputtered material on the side of the photoresist for thick film and high working pressure. The consequence of this defect on the properties depends on the involved layer and can be a lack of shielding for the central copper track or an over-di- mensioning of the sensor for the Finemet layers. The former reduces the effect of the transverse induced anisotropy, as the upper Finemet layer is separated in three magnetically independent parts, whereas the latter impedes down-sizing of the sensor and thus limits the spatial resolution. The sputtering conditions of copper have been optimized in order to limit over-deposition at the edge of the central conductive track. Thus the upper Finemet layer is continuous and provides magnetic shielding of the copper track. However, the defects at the edge of the finemet layer are not yet eliminated, as deposition conditions are driven by the minimization of the coercivity. Their influence on MI properties is theoretically weak, as calculated by digital simulation. However, these defects act as pinning centers for the magnetic domains movement and reduce the sensitivity of the device.

Sensors with variable dimensions and thicknesses have been realized (see table 1). Digital simulations predict, for a given excitation frequency, an increase in the sensitivity with the films thickness. Indeed, experimental characterizations show that the maximal sensitivity increases from 1900 V/T/A for the G500 sensor to 4100 V/T/A for the G750 sensor. Unfortunately, the size of defect above the copper track increases with the film thickness. This results in a low sensitivity for the sensor B1000. Fig. 8 presents the maximum sensitivity of the G750 sensor as a function of the frequency with a 1.7 mA excitation current. It reaches a maximum at 30 MHz. This optimal frequency is slightly larger than expected (in the range of 20 MHz for such thicknesses) what is supposed to be linked to an overestimation of µt in the calculation. Below this optimal frequency, the skin depth is larger than the film thickness thus the ferromagnetic film is poorly sensitive to skin effect. At the opposite, high frequency excitation current induces a skin effect located at the surface of the film and the magneto-impedance suffers from a magnetic hardening due to surface anisotropy. On the same figure, the results are compared to FeCoSiB/Cu/Fe-CoSiB [22] and FeSiB/Cu/FeSiB [23] MI microsensors whose sensitivity and excitation frequency are smaller. Results are also compared to the sensitivity of a macrosensor realized with Finemet and copper ribbons which thicknesses are 20 µm and 40 µm respectively. The total size of the sensor is 5 cm x 1 cm. It is noticeable that there is no drastic decrease of the sensitivity due to scale effect. Indeed, the structures of both the macro-sensor and the present devices are comparable. The maximal sensitivities are in the same range. However, the optimal frequency increases from 100 kHz to 30 MHz for the macro-and microsensors respectively, due to the thickness reduction.

Conclusion

This paper reports a technological process for depositing ultrasoft Finemet thin films using a common technique. In addition, it has been shown that this material can be integrated in MI microsensors. The sensitivity can reach 4100 V/T/A, which is the largest among all published similar structures. The results highlight the strong dependence of the sensor quality on the presence of a defect above the copper track which stems from the association of lift-off and non directive sputtering of thick films. In future works, strength will be put on the removal of this defect, for instance using a specific lift-off resist dedicated to thick film deposition. The optimal excitation frequency for the sensor G750 is 30 MHz, which remains high for standard electronic. Thicker films will allow decreasing this frequency to 10 MHz.

Figure 1 :

 1 Figure 1: Finite elements calculation of the magnetoimpedance as a function of the transverse permeability of the ferromagnetic layer. Dimensions correspond to G500.

Figure 3 :

 3 Figure 3: Influence of oxygen contamination on the coercive field

Figure 2 :Figure 5 :

 25 Figure 2: Fabrication process

Figure 4 :

 4 Figure 4: Internal stress of a 500 nm thick continuous Finemet film deposited on a 2" wafer during heating (red) then cooling (blue)

Figure 6 .

 6 Figure 6. Hysteresis loops of 500 nm Finemet film annealed for 1h at various temperature.

Figure 7 :

 7 Figure 7: Influence of field annealing at 300°C on hysteresis loops of 500 nm thick samples

Figure 8 :

 8 Figure 8: Dependence of the G750 sensor maximum sensitivity on the excitation frequency. Comparison with other amorphous material-based MI sensors

Table 1 .

 1 Dimensions of the sensors

		Film thickness	Finemet	Copper a
		(nm)	(µm)	(µm)
	G500 G750 B1000	500/500/500 750/750/750 1000/1000/1000	4000 × 200 4000 × 200 4000 × 200	4200 × 50 4200 × 50 4200 × 50

a electrical

pads are not taken into account

Acknowledgment

This work has been financially supported by the French Agence Nationale pour la Recherche, through the Project IMAGINE and by the Conseil Général de l'Essone through the use of the equipments of the Centrale de Technologie Universitaire IEF-MINERVE.


A V Kazeykina

Absence of traveling wave solutions of conductivity type for the Novikov-Veselov equation at zero energy

We prove that the Novikov-Veselov equation (an analog of KdV in dimension 2 + 1) at zero energy does not have sufficiently localized soliton solutions of conductivity type.

Introduction

In this note we are concerned with the Novikov-Veselov equation at zero energy

∂ t v = 4Re(4∂ 3 z v + ∂ z (vw)), ∂ z w = -3∂ z v, v = v, v = v(x, t), w = w(x, t), x = (x 1 , x 2 ) ∈ R 2 , t ∈ R, (1)
where

∂ t = ∂ ∂t , ∂ z = 1 2 ∂ ∂x 1 -i ∂ ∂x 2 , ∂ z = 1 2 ∂ ∂x 1 + i ∂ ∂x 2 .
Definition 1. A pair (v, w) is a sufficiently localized solution of equation ( 1) if

• v, w ∈ C(R 2 × R), v(•, t) ∈ C 3 (R 3 ), • |∂ j x v(x, t)| q(t) (1 + |x|) 2+ε , |j| 3, for some ε > 0, w(x, t) → 0, |x| → ∞, • (v, w) satisfies (1). Definition 2. A solution (v, w) of (1) is a soliton (a traveling wave) if v(x, t) = V (x -ct), c ∈ R 2 .
Equation ( 1) is an analog of the classic KdV equation. When v = v(x 1 , t), w = w(x 1 , t), then equation ( 1) is reduced to KdV. Besides, equation ( 1) is integrable via the scattering transform for the 2-dimensional Schrödinger equation

Lψ = 0, L = -∆ + v(x, t), ∆ = 4∂ z ∂ z , x ∈ R 2 .
(2)

Equation ( 1) is contained implicitly in [M] as an equation possessing the following representation

∂(L -E) ∂t = [L -E, A] + B(L -E), ( 3
)
where L is defined in (2), A and B are suitable differential operators of the third and zero order respectively and [•, •] denotes the commutator. In the explicit form equation (1) was written in [START_REF] Novikov | Finite-zone, two-dimensional, potential Schrödinger operators. Explicit formula and evolutions equations[END_REF], [START_REF] Novikov | Finite-zone, two-dimensional Schrödinger operators. Potential operators[END_REF], where it was also studied in the periodic setting. For the rapidly decaying potentials the studies of equation ( 1) and the scattering problem for (2) were carried out in [BLMP], [GN] [T], [LMS]. In [LMS] the relation with the Calderón conductivity problem was discussed in detail.

Definition 3. A potential v ∈ L p (R 2 ), 1 < p < 2, is of conductivity type if v = γ -1/2 ∆γ 1/2 for some real-valued positive γ ∈ L ∞ (R 2 ), such that γ δ 0 > 0 and ∇γ 1/2 ∈ L p (R 2 ).
The potentials of conductivity type arise naturally when the Calderón conductivity problem is studied in the setting of the boundary value problem for the 2-dimensional Schrödinger equation at zero energy (see [START_REF] Novikov | Multidimensional inverse spectral problem for the equation -∆ψ+(v(x)-Eu(x))ψ = 0. Funkt. Anal. i Pril[END_REF], [N], [LMS]); in addition, in [N] it was shown that for this type of potentials the scattering data for (2) are well-defined everywhere.

The main result of the present note consists in the following: there are no solitons of conductivity type for equation (1). The proof is based on the ideas proposed in [START_REF] Novikov | Absence of exponentially localized solitons for the Novikov-Veselov equation at positive energy[END_REF].

This work was fulfilled in the framework of research carried out under the supervision of R.G. Novikov.

2 Scattering data for the 2-dimensional Schrödinger equation at zero energy with a potential of conductivity type

Consider the Schrödinger equation ( 2) on the plane with the potential v(z), z =

x 1 + ix 2 , satisfying v(z) = v(z), v(z) ∈ L ∞ (C), |v(z)| < q(1 + |z|) -2-ε for some q > 0, ε > 0. ( 4
)
For k ∈ C we consider solutions ψ(z, k) of ( 2) having the following asymptotics

ψ(z, k) = e ikz µ(z, k), µ(z, k) = 1 + o(1), as |z| → ∞, (5)
i.e. Faddeev's exponentially growing solutions for the two-dimensional Schrödinger equation ( 2) at zero energy, see [F], [GN], [START_REF] Novikov | Multidimensional inverse spectral problem for the equation -∆ψ+(v(x)-Eu(x))ψ = 0. Funkt. Anal. i Pril[END_REF].

It was shown that if v satisfies (4) and is of conductivity type, then ∀k ∈ C\0 there exists a unique continuous solution of (2) satisfying (5) (see [N]). Thus the scattering data b for the potential v of conductivity type are well-defined and continuous:

b(k) = C e i(ky+ k ȳ) v(y)µ(y, k)dReydImy, k ∈ C\0. (6)
In addition (see [N]), the function µ(z, k) from ( 5) satisfies the following ∂-equation

∂µ(z, k) ∂ k = 1 4π k e -i(kz+ kz) b(k)µ(z, k), z ∈ C, k ∈ C\0 (7)
and the following limit properties:

µ(z, k) → 1, as |k| → ∞, (8)
µ(z, k) is bounded in the neighborhood of k = 0. (9)
The following lemma describes the scattering data corresponding to a shifted potential. As for the time dynamics of the scattering data, in [BLMP], [GN] it was shown that if the solution (v, w) of ( 1) exists and the scattering data for this solution are well-defined, then the time evolution of these scattering data is described as follows:

b(k, t) = e i(k 3 + k3 )t b(k, 0), k ∈ C\0, t ∈ R. ( 11
)
3 Absence of solitons of conductivity type Theorem 1. Let (v, w) be a sufficiently localized traveling wave solution of (1) of conductivity type. Then v ≡ 0, w ≡ 0.

Scheme of proof. From ( 10), ( 11), continuity of b(k) on C\0 and the fact that the functions k, k, k 3 , k3 , 1 are linearly independent in the neighborhood of any point, it follows that b ≡ 0. Equation ( 7) implies that in this case the function µ(z, k) is holomorphic on k, k ∈ C\0. Using properties (8) and ( 9) we apply Liouville theorem to obtain that µ ≡ 1. Then ψ(z, k) = e ikz and from (2) it follows that v ≡ 0.

Lemma 1 .

 1 Let v(z) be a potential satisfying (4) with the scattering data b(k). The scattering data b y (k) for the potential v y (z) = v(zy) are related to b(k) by the following formula b y (k) = e i(ky+ kȳ) b(k), k ∈ C\0, y ∈ C. (10) Proof. We note that ψ(zy, k) satisfies (2) with v y (z) and has the asymptotics ψ(zy, k) = e ik(z-y) (1+o(1)) as |z| → ∞. Thus ψ y (z, k) = e iky ψ(z -y, k) and µ y (z, k) = µ(z -y, k). Finally, we have b y (k) = C e i(kζ+ k ζ) v y (ζ)µ y (ζ, k)dReζdImζ = = C e i(kζ+ k ζ) v(ζy)µ(ζy, k)dReζdImζ = e i(ky+ kȳ) b(k).

CMAP, Ecole Polytechnique, Palaiseau, 91128, France; email: kazeykina@cmap.polytechnique.fr


Simone Pigolotti

Massimo Cencini

Fax +4535325425 Ph +4535325238

Coexistence and invasibility in a two-species competition model with habitat-preference

Keywords: Spatial models, Dispersal, Voter model, Heterogeneous habitat, Neutral Theory

The outcome of competition among species is influenced by the spatial distribution of species and effects such as demographic stochasticity, immigration fluxes, and the existence of preferred habitats. We introduce an individual-based model describing the competition of two species and incorporating all the above ingredients. We find that the presence of habitat preference -generating spatial niches -strongly stabilizes the coexistence of the two species. Eliminating habitat preference -neutral dynamics -the model generates patterns, such as distribution of population sizes, practically identical to those obtained in the presence of habitat preference, provided an higher immigration rate is considered. Notwithstanding the similarity in the population distribution, we show that invasibility properties depend on habitat preference in a non-trivial way. In particular, the neutral model results results more invasible or less invasible depending on whether the comparison is made at equal immigration rate or at equal distribution of population size, respectively. We discuss the relevance of these results for the interpretation of invasibility experiments and the species occupancy of preferred habitats.

Introduction

A central problem in community ecology is to understand the ecological forces leading to the observed patterns of coexistence or exclusion of competing species [START_REF] Ricklefs | Species diversity in ecological communities: historical and geographical perspectives[END_REF][START_REF] Brown | Spatial Variation in Abundance[END_REF]. This issue is important for understanding both simple communities made up of few species [START_REF] Chesson | Mechanisms of manteinance of species diversity[END_REF] and "biodiversity hotspots" with a large number of coexisting species [START_REF] Leigh | Why do some tropical forests have so many species of trees?[END_REF]. Historically, this problem has been approached at two distinct levels. On the one hand, focus has been put on the detailed mechanisms of interaction between species (e.g. intra-and inter-specific competitions) caused by their differentiation in exploiting resources, resulting in the concept of the ecological niche [START_REF] Chase | Ecological Niches: Linking Classical and Contemporary Approaches[END_REF]. For example, it has been shown how habitat heterogeneity [START_REF] Beckage | Seedling survival and growth of three forest tree species: the role of spatial heterogeneity[END_REF] or a tradeoff in dispersal range strategies [START_REF] Bolker | Spatial Moment Equations for Plant Competition: Understanding Spatial Strategies and the Advantages of Short Dispersal[END_REF]) may promote coexistence. An alternative explanation for the observed species richness and distribution is in terms of processes intrinsically due to chance, such as colonization, immigration and extinction [START_REF] Macarthur | The theory of island biogeography[END_REF], disregarding differences among species.

In recent years, the neutral theory of biodiversity [START_REF] Hubbell | Tree dispersion, abundance and diversity in a dry tropical forest[END_REF][START_REF] Bell | Neutral Macroecology[END_REF][START_REF] Hubbell | The Unified Neutral Theory of Biodiversity and Biogeografy[END_REF] considerably developed the latter approach by explicitly assuming equivalence among species at the individual level. The interest in the neutral theory has been triggered by its ability to successfully predict several biodiversity patterns observed in tropical forests, such as species-abundance distributions in different permanent sampling plots [START_REF] Bell | Neutral Macroecology[END_REF][START_REF] Hubbell | The Unified Neutral Theory of Biodiversity and Biogeografy[END_REF][START_REF] Volkov | Neutral theory and relative species abundance in ecology[END_REF] and species-area relations [START_REF] Durrett | Spatial Models for Species-Area Curves[END_REF][START_REF] Bell | Neutral Macroecology[END_REF][START_REF] Hubbell | The Unified Neutral Theory of Biodiversity and Biogeografy[END_REF]. Its success underlined the importance of stochasticity (ecological drift) and dispersal limitation in the assemblage of natural communities [START_REF] Chave | Neutral theory and community ecology[END_REF][START_REF] Alonso | The merits of neutral theory[END_REF], which are now recognized as key ingredient also in niche-based models [START_REF] Tilman | Niche tradeoffs, neutrality, and community structure: A stochastic theory of resource competition, invasion, and community assembly[END_REF]. However, niche-based models yield predictions for the biodiversity patterns which perform similarly to neutral ones when compared with data [START_REF] Chave | Comparing Classical Community Models: Theoretical Consequences for Patterns of Diversity[END_REF][START_REF] Mcgill | Strong and weak tests of macroecological theory[END_REF][START_REF] Mouquet | Community Patterns in SourceSink Metacommunities[END_REF][START_REF] Gilbert | Neutrality, niches, and dispersal in a temperate forest understory[END_REF][START_REF] Tilman | Niche tradeoffs, neutrality, and community structure: A stochastic theory of resource competition, invasion, and community assembly[END_REF]. This suggests that these patterns tend to average out the dependence on the details of the theory (see also the discussion in [START_REF] Pueyo | The maximum entropy formalism and the idiosyncratic theory of biodiversity[END_REF]) and thus cannot be used to discriminate the relative importance of niche-based and neutral forces. In this perspective, the study of dynamical properties such as invasibility can be a promising way to disentangle these effects [START_REF] Daleo | Biological invasions and the neutral theory[END_REF].

To understand the key differences between neutral and non-neutral competition, it is useful to consider models that can be continuously tuned from niche-based to neutral settings by varying some parameters [START_REF] Chave | Comparing Classical Community Models: Theoretical Consequences for Patterns of Diversity[END_REF][START_REF] Gravel | Reconciling niche and neutrality: the continuum hypothesis[END_REF][START_REF] Adler | A niche for neutrality[END_REF]. An obvious difficulty with this program comes from the unavoidable complexity of realistic niche-based models [START_REF] Chase | Towards a really unified theory for metacommunities[END_REF], where species are not equivalent and the environment is heterogeneous both in space and time. This suggests an approach whereby simplified models with few parameters are studied, for example by making some specific assumptions on how neutrality is violated.

In this paper we study the dynamics of two species A and B that compete for space. The model is devised in such a way that a single parameter controls the overlap between the niches occupied by the two species, from complete -neutral -to no overlap -two independent niches. The model is individual-based and incorporates the basic ingredients of neutral theory: coexistence results from immigration from a metacommunity, balancing demographic stochasticity, which alone would lead to extinction. Niches are introduced in this neutral scenario via preferential habitats: half of the sites in the ecosystem are favorable for the colonization of individuals of one species and the other half are favorable for the other species. The ecological advantage is realized through a biased "lottery" [START_REF] Chesson | Environmental variability promotes coexistence in lottery competitive systems[END_REF]. We consider a symmetric situation by choosing the same statistical bias, γ, for individuals of species A and B to colonize their respective preferred habitats. When γ = 0 (no habitat diversification) the model reduces to the voter model [START_REF] Holley | Ergodic Theorems for Weakly Interacting Infinite Systems and the Voter Model[END_REF]Cox and Griffeath, 1986), which is a prototype of neutral dynamics. Increasing γ, species acquire an advantage in colonizing some sites and a disadvantage in others. A very large γ eventually leads to segregation of the two species to their preferential habitats. Segregation will be complete when the choice of dispersal allows individuals to reach all their preference sites or incomplete in the presence of dispersal limitation.

The aim of this work is to use this simple model to understand the effect of habitat diversification on coexistence and dynamics of ecological communities. In particular, our concern will be on contrasting the effect of habitat diversification with the neutral model where no preferred habitat exists.

Model

We consider an individual based, spatially-explicit model of a community made of two competing species A and B with population N A and N B , respectively. The community lives in a patch made of N = L 2 sites on a square lattice of side L, on which we assume periodic boundary conditions. Each lattice site is occupied by a single individual of one of the two species. For the sake of simplicity, we assume that the patch is saturated, i.e. with no empty sites -each dead individual is immediately replaced, so that the total number of individual is constant, N A + N B = N. The latter hypothesis is commonly assumed for its convenience [START_REF] Hubbell | The Unified Neutral Theory of Biodiversity and Biogeografy[END_REF][START_REF] Chave | Comparing Classical Community Models: Theoretical Consequences for Patterns of Diversity[END_REF] and, strictly speaking, corresponds to considering infinite fecundity. However, a finite but reasonably high fecundity would lead to almost-saturated ecosystems with qualitatively similar dynamics [START_REF] Durrett | Spatial Models for Species-Area Curves[END_REF][START_REF] Chave | Comparing Classical Community Models: Theoretical Consequences for Patterns of Diversity[END_REF].

Our main interest here is to study the effect of habitat preference on competition. To this aim, we assign to each site a specificity: half of the sites are favorable (as below specified) to individuals of species A and the other half to individuals of species B. We denote such sites by a and b, respectively. The site specificity can have several different (often concomitant) ecological origins such as abundance of resources, predation pressure (see, e.g., the review by [START_REF] Amarasekare | Competitive coexistence in spatially structured environments: a synthesis[END_REF], and/or environmental conditions such as elevation, temperature, soil moisture or other parameters as in [START_REF] Zillio | The impact of neutrality, niche differentiation and species input on diversity and abundance distributions[END_REF] and as suggested by observations [START_REF] Beckage | Seedling survival and growth of three forest tree species: the role of spatial heterogeneity[END_REF]. The net effect of these different mechanisms is here assumed to increase by a factor γ the chance of individuals to colonize a preferred site. This is illustrated in the top panel of Fig. 1, to be compared with the bottom cartoon which shows the neutral model, without site specificity. Similar models have been proposed also in the context of heterogeneous catalysis [START_REF] Frachebourg | Heterogeneous Catalysis on a Disordered Surface[END_REF] and social dynamics [START_REF] Masuda | Heterogeneous Voter Models[END_REF].

For the sake of simplicity, site specificity is randomly assigned at the beginning and left unchanged during the dynamics. Of course, in natural ecosystems, spatial arrangement of sites with a certain specificity will usually be characterized by a certain degree of correlations, which will in general tend to enhance niche effects. In this respect, we expect that our choice will tend to underestimate the effect of habitat preference. Clearly, the model can be generalized by introducing asymmetries, i.e. different γ's for the two species or different fractions of advantageous sites. Here we shall limit our analysis to the simple symmetric case, so that no species has a net advantage and the degree of habitat preference is controlled by a unique parameter.

We also assume a continuous immigration in the patch of individuals A or B at rate ν. This inflow is necessary to avoid the drift to extinction of one of the two species.

For any given size L of the patch, which fixes the number of individuals N = L 2 , the model is controlled by two parameters only: the colonization advantage γ and the immigration rate ν. The elementary step of the dynamics is as follows i) a site is randomly chosen and the individual there residing is killed; ii) with probability (1ν), the individual is replaced by a copy of one of the four neighbors, chosen via a lottery which gives a competitive advantage (modeled as a weight γ) to individuals having that site as preferred habitat (see the sketch in Fig. 1 andeqs. (1) and (2) below). iii) with probability ν, the individual is replaced by an immigrant. For simplicity, we assume the two species being equipopulated at the metacommunity level, so that the probability of being replaced by an individual of one of the two species is 1/2, apart from the competitive advantage on the specific empty site.

In formulas, steps ii) and iii) can be expressed as follows. If the individuals

W a A (n A , n B ) = (1 -ν) (1+γ)n A (1+γ)n A +n B + ν 1+γ 2+γ W a B (n A , n B ) = (1 -ν) n B (1+γ)n A +n B + ν 1 2+γ , (1)
respectively, where n A and n B denote the number of individuals of species A and B in the neighborhood of the considered site. Similarly, if the individual is killed in a site of type b, we have

W b A (n A , n B ) = (1 -ν) n A n A +(1+γ)n B + ν 1 2+γ W b B (n A , n B ) = (1 -ν) (1+γ)n B n A +(1+γ)n B + ν 1+γ 2+γ .
The competitive lottery (1) and (2) used in step ii) represents a biased (nonneutral) generalization of that used in neutral models [START_REF] Hubbell | The Unified Neutral Theory of Biodiversity and Biogeografy[END_REF]. Indeed, for γ = 0, the neutral dynamics of the voter model with only two species is recovered [START_REF] Holley | Ergodic Theorems for Weakly Interacting Infinite Systems and the Voter Model[END_REF].

Notice that the model is set up in such a way that the fitness advantage to be on colonization: after having colonized a site, mortality and dispersal do not depend on being on a preference site. In other terms, the fitness advantage belongs to the seeds and not to the individuals themselves. In this interpretation, and at variance with [START_REF] Chesson | Environmental variability promotes coexistence in lottery competitive systems[END_REF], we excluded from the contribution to the (implicit) seed pool of the individual who died. Although the latter may be more realistic in modeling, e.g., perennial plants [START_REF] Lin | Demographic trade-offs in a neutral model explain death-rate-abundance-rank relationship[END_REF], we preferred the first to directly compare the neutral version of the model (γ = 0) originally proposed by [START_REF] Hubbell | Tree dispersion, abundance and diversity in a dry tropical forest[END_REF] with the non-neutral (γ > 0) variants. We do not expect big differences in the outcome of the two models, a part from a weak discrepancy in the percentage of occupancy of preferential habitats. We also note that this difference becomes less relevant when longer dispersal is introduced. Indeed, the model can be generalized to different form of dispersal mechanisms, e.g. with a given finite range. Previous investigations on similar models [START_REF] Rosindell | Species-area relationships from a spatially explicit neutral model in an infinite landscape[END_REF][START_REF] Zillio | The impact of neutrality, niche differentiation and species input on diversity and abundance distributions[END_REF][START_REF] Pigolotti | Speciation-rate dependence in species-area relationships[END_REF] have shown that different dispersal mechanisms yield qualitatively similar results, as far as they are finite-ranged and all species adopt the same dispersal strategy. Here, we consider the two limiting situations of nearest-neighbor and global dispersal.

In the global dispersal case, all individuals present in the patch can colonize the site left empty by the dead one. This limiting case strongly simplifies the simulations of the model since it makes the distance and spatial distribution of the sites irrelevant. The global dispersal model may be thought as a two islands model [START_REF] Wright | Evolution in Mendelian populations[END_REF], where each of the two islands has room for N/2 individuals and contains all the sites favorable to one of the two species.

The state of the system is then unequivocally identified by the numbers N Aa and N Bb of individuals of the two species occupying the respective sites of preference. All the other variables can be expressed in terms of N Aa and N Bb , for instance: the number of individuals of species A (resp. B) outside their preference sites is

N Ab = N/2 -N Bb (resp. N Ba = N/2 -N Aa ) and the total number of individuals of species A (resp. B) is N A = N/2 + N Aa -N Bb (resp. N A = N/2 + N Bb -N Aa ).
The evolution of the system is thus determined by the probabilities per elementary step that N Aa and N Bb increase or decrease by one unity:

W N Aa →N Aa +1 = 1 2 - N Aa N W a A (N A , N B ) W N Bb →N Bb +1 = 1 2 - N Bb N W b B (N A , N B ) W N Aa →N Aa -1 = N Aa N W a B (N A , N B ) W N Bb →N Bb -1 = N Bb N W b A (N A , N B ). ( 2
)
where W x Y is given by eqs. ( 1) and ( 2) with n A and n B replaced by N A and N B , respectively.

For γ = 0, the model is neutral and reduces to the ordinary voter model with only two species and immigration (the multi-species version was considered in [START_REF] Durrett | Spatial Models for Species-Area Curves[END_REF]) in the spatially explicit case, or to the Moran model with mutation [START_REF] Moran | Random processes in genetics[END_REF] in the global dispersal version. Conversely, when γ becomes very large, the competitive advantage gets so intense that individuals tend to localize in their preferred habitats, with an extremely low chance of colonizing the rest of the ecosystem, so that the two species do not compete anymore.

Results

Extinction times and fixation probabilities

In the absence of immigration (ν = 0) and for any choice of N, γ and dispersal range, persistent coexistence is not possible, since species cannot recover from a local extinction event. Ecological drift will ultimately drive one of the two species to extinction, and the dominance of the other will become stablereaching fixation as from population genetics terminology, see e.g. [START_REF] Gillespie | The causes of molecular evolution[END_REF]. Studying the properties of the dynamics toward extinction is however interesting and informative [START_REF] Chesson | Environmental variability promotes coexistence in lottery competitive systems[END_REF][START_REF] Chesson | The stabilizing effect of a random enviroment[END_REF]. In particular, the time needed for the extinction of one of the species is important to understand how crucially biodiversity depends on a steady immigration flux. Moreover, the way the extinction time and the probability of fixation of one of the two species depend on the deviation from neutrality allow to quantify to what extent habitat preference promotes coexistence.

We start describing a set of simulations in which the two species are equipopulated at the initial time, i.e. N A = N B = N/2, and individuals are randomly placed in the system. The dynamics is then followed until the extinction of one of the species. By averaging over many realizations of the dynamics (from 5×10 3 to 10 4 ), we can access the extinction time. Here we focus on the average time, T , though also the fluctuations have an ecological relevance [START_REF] Pigolotti | Species lifetimes for simple models of ecologies[END_REF]. The average extinction time T as a function of the community size N and for several values of γ is shown in Fig. 2 as obtained with the global dispersal model. As customary in models with overlapping generations, we measure T in generations, where the time unit corresponds to N successive time steps i)-iii) described in the model section. For γ = 0, we recover the result expected for the Moran model [START_REF] Moran | Random processes in genetics[END_REF] i.e. T ∝ N. While the presence of habitat preference γ > 0 leads to a dramatic increase of the average extinction time: for large enough populations, T becomes exponentially large in the community size N, i.e. where C depends on γ (as shown in the inset of Fig. 2a). The exponential dependence implies that the coexistence may be considered stable for large N [START_REF] Chesson | The stabilizing effect of a random enviroment[END_REF]. The spatially explicit version of the model (not shown here) presents the same qualitative features with small quantitative changes: in the neutral case γ = 0, logarithmic corrections are present T ∝ N log(N) [START_REF] Krapivsky | Kinetics of monomer-monomer surface catalytic reactions[END_REF]) and the exponential rates for γ > 0 are slightly different from those obtained with global dispersal (inset of Fig. 2a). In large communities, even a tiny habitat preference leads to a dramatic increase in the average extinction time, stabilizing the system on any realistic timescale.

To further confirm the stabilizing effect of habitat preference, we consider a different setting in which species A is present at initial time with a fraction of individuals x = N A /N . By averaging over 10 5 realizations for each x, we computed the probability that species A becomes fixated, P fix (x) (Fig. 2b).

In the neutral case, P fix (x) = x (Gillespie, 1994), while P fix (x) is closer to 1/2 when γ is increased, the effect being stronger the larger is the community size (see inset). Therefore, strong habitat preference tends to compensate any initial disproportion between two large populations, making equally likely the fixation of one of the two species. In particular, if A is the invading species, i.e. the less represented at the beginning (x 1/2), its chances of invading the systems greatly increase in the presence of habitat preference, in agreement with other models and observations [START_REF] Melbourne | Invasion in a heterogeneous world: resistance, coexistence or hostile takeover?[END_REF].

Regimes of coexistence

The presence of an immigration pressure allows species to recover from local extinction events, ensuring dynamical coexistence, as illustrated in the left panels of Fig. 3, where the evolution of the population N A is shown over 10 3 generations at varying the rate of immigration ν. The figures refer to the global dispersal model with a habitat-preference intensity γ = 0.3. We mention that qualitatively similar features are observed also in the nearestneighbor dispersal case and also when the model is neutral.

As expected from classical theories [START_REF] Macarthur | The theory of island biogeography[END_REF], increasing the immigration pressure enhances the degree of coexistence of the two species. At fixed N, if ν is small, one of the two species dominates for most of the time. Stochastic fluctuations, whose amplitude decreases at increasing N, lead to an alternation in the dominating species (top panel of Fig. 3). The turnover in the dominating species takes place on timescales rapidly growing with N (roughly of the order of the average extinction time). At intermediate values of ν, coexistence is possible, but episodic local extinctions can rule out a species from the system for several generations (middle panel). Finally, increasing ν even further, local extinctions cannot persist and coexistence becomes the rule: it is more and more likely to observe states in which the two species are roughly equipopulated (bottom panel). We quantify these three behaviors in terms of the functional shape of the probability P (N A ), averaged over time, of one of the two species to occupy a given fraction of the ecosystem (see, e.g. [START_REF] Loreau | Immigration and the manteinance of local species diversity[END_REF]. We recall that there is a complete symmetry between A and B so that P (N A ) = P (N -N A ) = P (N B ). In particular, with reference to the right panels of Fig. 3, from top to bottom we can distinguish the following classes of distributions corresponding to different regimes of coexistence, as determined by the number of maxima of P (N A ): Monodominance: P (N A ) is U-shaped, achieving its two maxima at N A = 0 and N A = N; Mixed : P (N A ) has three maxima at N A = 0 and N A = N as before plus N A = N/2, meaning alternation between states of monodominance and pure coexistence as defined below; Pure coexistence: P (N A ) is bell-shaped and has a single maximum at N A = N/2, so that most of the time the populations fluctuate around the equipopulated state.

We quantify the stabilizing effect of habitat preference by studying for which values of the parameters γ, ν, and N, the above defined regimes are observed.

Figure 4 shows, for different values of γ, how the monodominance (white), coexistence (gray) and mixed (patterned) regimes organize in the Nν plane for both global (left) and nearest-neighbor (right) dispersal.

As observed for the extinction times, nearest-neighbor and global dispersal models display qualitatively similar features. The only difference is that in the neutral case γ = 0 the global dispersal model never shows a mixed distribution; conversely, with nearest-neighbor dispersal, a tiny region of mixed phase exists also for γ = 0, compare Fig. 4a andd. For γ = 0 and global dispersal, the transition between monodominance and coexistence can be analytically obtained and occurs for ν > ν c (N) = 2/(2 + N) (see Appendix A).

The most surprising feature is that the critical immigration rate ν c (N) for observing coexistence seems to be independent of γ and of the dispersal. Indeed, we could not determine appreciable quantitative differences between the two extrema of global and nearest-neighbor dispersal. In all cases, the transition to coexistence occurs for values of the immigration rate well described by the neutral prediction, ν c (N ) ≈ 2/(2 + N). Unfortunately, we could not systematically explore larger values of γ or N as it requires huge statistics for distinguishing between the mixed and coexistence regimes. In fact at increasing γ and/or N the distribution becomes strongly peaked on N/2 with very low probabilities on the tails, so that the region where the differences between mixed and pure coexistence manifest is inaccessible (see Fig. 5).

Summarizing, the main effect of increasing γ (from top to bottom in Figure 4) is to stabilize the coexistence by increasing the portion of the Nν plane where the mixed regime is realized. An equivalent analysis could in principle be performed by plotting the transition lines in the Nγ and/or in the νγ plane. However, the difficulty in sampling the tails of the distribution strongly limits the range of values of γ and N which can be explored, so that these plots would not add much information to Figure 4).

Neutral vs Non-neutral: coexistence & invasibility

We now directly compare neutral and non-neutral dynamics by exploring whether the former (γ = 0) can reproduce/mimic the latter, e.g., in terms of generating similar patterns of coexistence such as the distribution P (N A ).

As both γ and ν promote coexistence, it is natural to expect that the lack of habitat preference in the neutral model could be compensated by a larger immigration rate ν to reproduce similar distributions. However, as we will see, the different stabilizing effects of immigration and habitat preference have interesting dynamical consequences.

In particular, we consider two attempts to reproduce non-neutral distributions with neutral ones for two different community sizes N = 100 (Fig. 5a) and N = 400 (Fig. 5b). In both cases, we proceed as follows. We fixed the habitat preference intensity, γ = 1 in this specific example. Then we simulated the neutral version of the model (γ = 0) and varied ν until the distribution looked as similar as possible to that obtained with habitat preference. In other words, we searched for the value of the immigration rate ν that compensated best the absence of habitat preference. As expected, this compensation is obtained by using a larger value of ν.

For N = 100, the agreement between the curves is very good in the central peak, meaning that the differences between the two models can be appreciated only when one of the two species is dominating. For N = 400 the distributions are almost indistinguishable. The fact that the non-neutral system of Fig 5a

is in the mixed phase is clear by looking at the distribution tails; however this feature could be difficult to detect in an experimental time series. In Fig. 5b, the tails have such a low probability that are essentially inaccessible even in long simulations, so that we cannot distinguish between mixed and coexistence phases.

For large N, both in the coexistence and mixed phase, the distributions around the central mode are well fitted by a Gaussian. In the neutral model, the limiting Gaussian distribution can be analytically derived (see Appendix B), showing that the variance in this case is proportional to N. Numerical simulations (not shown) suggest that variance in the non-neutral version of the model scales in the same way with N, with a prefactor rapidly decreasing at increasing γ reflecting the presence of niches.

Even if the difference between the models is essentially undetectable when looking at the distributions, it can play an important role when studying invasibility properties, i.e. by considering a situation in which one of the species is absent at the initial time (which means to directly test the tails of the distribution). As an example, for the same parameters of Fig. 5 we prepared the initial state with only species B present and we computed the average time it takes to reach perfect coexistence (N A = N B = N/2) for the first time. When N = 100 (a), this time is roughly equal to 46 generation for the model with habitat preference and 17 generation for the neutral one, while for N = 400 (b) the difference is even more pronounced (61 vs 21 generations, respectively).

Habitat preference statistics

We now investigate the role of preference sites in achieving coexistence and determining the spatial distribution of species. The natural quantity to look at is the average fraction of individuals of both species occupying their preference sites, P occ = (N Aa + N Bb )/N . This quantity is a suitable measure of the nonneutrality of the system: it is equal to 1/2 in the neutral case, as individuals have no reason to prefer one of the two site types. The average occupation increases with γ eventually reaching 1 when the distribution of sites strongly determines the spatial distribution of the two species.

The average occupation P occ as a function of γ is shown in Fig. 6 . For illustrative purposes, here, we fixed N = 100 and varied ν. As one can see the curves are weakly dependent on the choice of ν, especially when γ is large. The dependence on N is not shown but it is also rather weak. We remark that in the nearest neighbor case it is crucial to average over different realizations of the dynamics and of the distribution of advantageous sites: different spatial arrangements of the sites may be harder (or easier) to occupy for one of the species.

When the two species are exactly equipopulated and with global dispersal, one can easily derive from (2) that the average occupation must be equal to (1 + γ)/(2 + γ). This prediction is shown for comparison in Fig. 6a. Significant deviation are present for small γ and ν values, for which pure coexistence is not achieved and species are not equipopulated, i.e. the system is in the monodominance or mixed phase. However, it must be noticed that with nearestneighbor dispersal (Fig. 6b) even for large γ the prediction (1 + γ)/(2 + γ) is never realized due to dispersal limitation, which can prevent species from colonizing advantageous sites. In this case the average occupation does not reach 1 even for very large γ. The presence of patchiness in the distribution of preference sites is expected to enhance this effect: one could encounter situations in which a whole patch is unreachable because it is surrounded by patches being advantageous to the other species. On the other hand, a longer dispersal range could compensate for the presence of patches. In other words, the possibility for the occupation to reach one for large γ will depend on the ability of species to reach all the favorable patches, which is known to be a crucial feature in fragmented landscapes [START_REF] Keitt | Detecting Critical Scales in Fragmented Landscapes[END_REF].

We conclude this section by discussing whether the results of Fig. 6 are determined by the particular choice of the ecological advantage. To check the robustness of our results, we simulated the variant of the model in which the advantage gives rise to a lower mortality in the preference sites, instead of a larger colonization probability, obtaining very similar curves to those of Fig 6.

Discussion

We have introduced an individual-based model of competition between two species in the presence of habitat preference. As is typical in stochastic models of community assembly, the system drifts towards the extinction of one of the two species when immigration from the metacommunity is not taken into account [START_REF] Macarthur | The theory of island biogeography[END_REF][START_REF] Loreau | Immigration and the manteinance of local species diversity[END_REF]. Nevertheless, even without immigration, habitat preference has a strong stabilizing effect on the coexistence in large communities, leading to a dramatic increase of the average time for the extinction of a species, consistently with studies of niche models [START_REF] Tilman | Niche tradeoffs, neutrality, and community structure: A stochastic theory of resource competition, invasion, and community assembly[END_REF][START_REF] Adler | A niche for neutrality[END_REF][START_REF] Lin | Demographic trade-offs in a neutral model explain death-rate-abundance-rank relationship[END_REF]. In particular, [START_REF] Lin | Demographic trade-offs in a neutral model explain death-rate-abundance-rank relationship[END_REF], extending previous results [START_REF] Zhang | The effects of competitive asymmetry on the rate of competitive displacement: how robust is Hubbells community drift model?[END_REF][START_REF] Yu | Can high tree species richness be explained by Hubbells null model?[END_REF][START_REF] Wright | Plant diversity in tropical forests: a review of mechanisms of species coexistence[END_REF], introduced a model where both mortality and fecundity were varied among the species in such a way that the averaged life-history fitness of several species is the same, so to fulfill the weaker requirement of the equivalence of average fitness of the neutral theory [START_REF] Hubbell | Neutral theory and the evolution of ecological equivalence[END_REF]. The presence of birth/fecundity tradeoffs (globally equalizing the performances of different species) has a stabilizing effect, without affecting the ability of predicting biodiversity patterns.

The main message of the present study is that neutral and non-neutral coexistence may yield almost identical patterns such as the population distribution (see, e.g., Fig. 5b) while dynamic properties may be significantly different, as revealed by invasibility experiments [START_REF] Daleo | Biological invasions and the neutral theory[END_REF]. In this perspective, the advantage of the proposed model is to allow for quantitatively comparing niche and neutral competition in a simple setting and thus for providing an insight on the origin of such differences. It is thus worth discussing the predictions of the model for invasibility experiments.

The simplest setting is that of an isolated ecosystem (i.e. without migratory flux), in which a new species with a small population is introduced. This corresponds to Fig. 2b, where the probability of the invading species to take over grows with γ. The ecological interpretation is that niche-based assemblies are easier to invade than neutral ones if the niches are not saturated, simply because the availability of a free niche facilitates the invasion (see [START_REF] Melbourne | Invasion in a heterogeneous world: resistance, coexistence or hostile takeover?[END_REF] for a recent review on the role of heterogeneity on invasibility properties).

Another possibility is to consider the invasion of an ecosystem initially populated by a dominant species in the presence of a regular flux of immigrants species from the metacommunity. In this case, if we compare neutral and nonneutral dynamics at equal immigration rate there is hardly any difference with respect to what is observed in the absence of immigration.

Conversely, the response may be quite different if the comparison is made at equal realized distributions (see Fig. 5). The neutral model realizes patterns of coexistence similar to non-neutral ones by compensating with a larger immigration pressure. Under these circumstances, if some catastrophic event leads to the extinction of one of the species, the non-neutral dynamics takes longer than the neutral one to recover the equipopulated state. Therefore, non-neutral coexistence is more robust, but it is also harder to achieve. The presence of two habitats makes both the two species coexist and one species hard to invade; in other words, niche-based coexistence is history-dependent. An important practical consequence is that inferring immigration rates from the fit of the distributions with neutral models could result in overestimates, even in situations in which the quality of the fit is very good.

The above result agrees with the observation that large scale, species-rich patches are easier to invade [START_REF] Robinson | Invasibility of experimental habitat in a California winter annual grassland[END_REF][START_REF] Stohlgren | The rich get richer: patterns of plant invasions in the United States[END_REF], contrarily to the prediction of theories of competitive exclusion via niche par-titioning and data from small scale observations [START_REF] Tilman | Community invasibility, recruitment limitation, and grassland biodiversity[END_REF]. These contrasting evidences constitute the so-called "invasion paradox" [START_REF] Fridley | The invasion paradox: reconciling pattern and process in species invasions[END_REF]. Neutral theory predicts that the chances of establishment of an alien species depend on the frequency of introduction of new individuals only [START_REF] Daleo | Biological invasions and the neutral theory[END_REF]. On the other hand, in a non-neutral community this chance is strongly influenced by habitat diversity and availability of free niches. Our results suggest that, when comparing the invasibility properties of ecosystems, it is crucial to establish in which measure the diversity is maintained by the intensity of the immigration flux and/or by niche availability. The two-species framework considered here the stationary solutions have not the same complexity as those of multi-species models. However, previous studies [START_REF] Chave | Comparing Classical Community Models: Theoretical Consequences for Patterns of Diversity[END_REF][START_REF] Mcgill | Strong and weak tests of macroecological theory[END_REF][START_REF] Tilman | Niche tradeoffs, neutrality, and community structure: A stochastic theory of resource competition, invasion, and community assembly[END_REF]) also pointed out that species abundance distributions of neutral and niche-based multi-species models are very similar, provided suitable parameters are tuned and suggesting that at least the qualitative features of invasibility dynamics here identified should hold also in the multi-species case. We recall that in multi-species models speciation, as well as immigration, contributes to the introduction of new species, especially when large systems are considered. Therefore, the problem of inferring these rates from observed patterns becomes even more dramatic, since it is hard to compare these estimate with reliable measures of speciation rates. This could explain the high speciation rates predicted by neutral theory (see [START_REF] Ricklefs | A comment on Hubbell's zero-sum ecological drift model[END_REF] and also [START_REF] Hubbell | Modes of speciation and the lifespans of species under neutrality: a response to the comment of Robert E. Ricklefs[END_REF]; [START_REF] Haegeman | Neutral Models with Generalised Speciation[END_REF]).

Finally, the occupation probability provides a measure of the departure from the static equilibrium case in which each species is confined to its preference sites. It is known that models correlating species distributions with environmental and climatic features cannot fully determine the observed geographical range of species [START_REF] Svenning | Limited filling of the potential range in European tree species[END_REF]. In fact, a complete predictive power of species-distribution models would require the existence of a static equilibrium state [START_REF] Guisan | Predicting species distribution: offering more than simple habitat models[END_REF], an assumption that is directly challenged by the neutral theory and partially violated by niche-neutral models like the one presented in this paper. The departure from the equilibrium case is enhanced by dispersal limitation, as shown by the comparison between the nearest neighbor case and the global dispersal of Fig. 6. Indeed, this observable is the one showing the most significant dependence on the choice of the dispersal among those considered in this work. This difference will become more dramatic if one adopts the more realistic choice of a correlated environment, i.e. a distribution of preference sites being correlated in space.

As it has been studied for the case of fragmented landscapes [START_REF] Keitt | Detecting Critical Scales in Fragmented Landscapes[END_REF], the comparison of dispersal ranges and characteristic size of patches determines whether species will be able to reach all their preference sites or will be arrested by dispersal limitation.

Summarizing, we studied the effect of habitat preference on the stochastic competition between two species. In the absence of immigration, demographic stochasticity eventually leads to the extinction of one of the two species. However, habitat preference leads to a dramatic increase of the extinction time.

The drift to extinction is arrested by introducing an immigration rate. When the latter is large enough, the model generates a coexistence state in which both species occupy a significant fraction of the system. By compensating with a larger immigration rate the loss of specificity, the neutral case (γ = 0) can generate distributions which are essentially indistinguishable from those obtained in the presence of habitat preference (γ > 0). In spite of the pattern similarity, the dynamics of the two cases is very different. When only one species is present at the beginning of the simulation, the neutral variant of the model is much easier to invade than the non-neutral one. This latter result reproduces in a simplified framework what is observed in species-rich community models: neutral and non-neutral models reproduce similar patterns of species abundances, here of population distribution, but can be distinguished by looking at dynamical properties.

which must hold at stationarity since the process is one dimensional [START_REF] Gardiner | Handbook of Stochastic Methods: for Physics[END_REF]. The above expression can be used to determine the critical value of the immigration rate for the transition from monodominance to coexistence. It is enough to determine the value of ν for which P (N A ) passes from being a decreasing function of N A (we are assuming N A < N/2) to an increasing one.

Coexistence is thus obtained when P (N A + 1) > P (N A ), which using (A.1) and after some algebra can be recast as [(2 + N)ν -2](N -2N A -1) > 0 .

(A.3) Therefore, whenever ν > 2/(2 + N) the distribution is increasing up to N A = (N -1)/2 then decreasing, i.e. there is coexistence. When ν < 2/(2 + N) the opposite happens and there is monodominance. We conclude noticing that for ν = 2/(2+N), i.e. at the the transition between the two classes of distributions, the distribution becomes uniform, i.e. P (N A ) = 1/N . This is true only for the neutral model with global dispersal.

B Gaussian limit of the neutral model

In this Appendix we show that the stationary solution approaches a Gaussian in the large N limit, at least in the neutral, global dispersal case. We perform the calculation in the diffusion approximation; the procedure is similar to the classic results for the Moran model, see, e.g. [START_REF] Karlin | On a genetic model of Moran[END_REF], apart from small differences in how the immigration (mutation) mechanism is defined.

From the rates (A.1) the diffusion approximation yields the Fokker-Planck equation:

∂ t ρ(x, t) = -∂ x (a(x)ρ(x, t)) + 1 2 ∂ 2 x (b 2 (x)ρ(x, t)) (B.1)
where x = N A /N and

a(x) = ν(1 -2x) 2 ; b 2 (x) = 2x(1 -x)(1 -ν) + ν N . (B.2)
The stationary distribution is easily found to be

ρ s (x) ∝ [2x(1 -x)(1 -ν) + ν] Nν 2(1-ν) -1 2 . (B.3)
When N is large, an expansion around the maximum x 0 = 1/2 leads to the Gaussian distribution:

ρ st (x) ∼ 1 √ 2πσ 2 exp - (x -x 0 ) 2 2σ 2 (B.

Fig. 1 .

 1 Fig. 1. Sketch illustrating the model with (top) and without (bottom) habitat preference. Two representative configurations of 4 × 4 systems are shown, white squares are advantageous to A, gray ones to B, and on the right we sketch the lottery dynamics (eqs. (1 and (2)): the width of arrows represents the habitat preference intensity γ. Notice that in the bottom panel the width of arrows is insensitive to the site specificity is killed in a site of type a advantageous for individuals of species A, the probabilities of being replaced by an individual A or B are given by

TFig. 2 .

 2 Fig. 2. Results for the model with global dispersal without immigration (ν = 0): (a) Mean extinction times T as a function of N for different γ, as in label. The black lines indicate exponential fits of the form T = C exp(C (γ)N ) while the gray straight line indicate the neutral expectation T ∝ N . The inset shows C (γ) versus γ for both the global dispersal (GD) and nearest neighbor (NN) models. (b) Probability of fixation P fix of A vs the initial fraction x = N A /N for N = 40. Inset: P fix vs N holding fixed x = N A /N = 0.1.

Fig. 3 .

 3 Fig.3. Different regimes of coexistence in the presence of immigration in the global dispersal model with habitat preference: evolution of the population N A across 10 3 generation for a system of N = 100 individuals (left) and the corresponding butions P (N A ) (right). The three panels are obtained holding the habitat-preference intensity fixed γ = 0.3 and varying the immigration rate ν: (top) ν = 0.003 corresponding to monodominance; ν = 0.015 leading to a mixed phase; ν = 0.05 displaying coexistence (see text for details).

Fig. 4 .

 4 Fig. 4. Phases in the Nν parameter space. White regions denote monodominance, gray-filled regions pure coexistence, and dashed regions mixed phases. On the left (panels a-c) data obtained with the global dispersal model, while on the right (panels d-f) with the nearest-neighbor model. The three rows correspond to different values of γ = 0, 0.3, 1.0 (from top to bottom). Axis are in log scale. Notice that for γ = 1 and global dispersal (panel c) we could not show the curves for N > 240 due to lack of statistics.

Fig

  Fig. 5. Non-neutral vs neutral distribution of population P (N A ): (a) for N = 100 with γ = 1, ν = 0.005 (solid) and γ = 0, ν = 0.12 (symbols): (b) for N = 400 with γ = 1, ν = 0.001 (solid) and γ = 0, ν = 0.123 (symbols). In both cases, the y axis is in logarithmic scale.

Fig. 6 .

 6 Fig. 6. Average preference occupation P occ as a function of γ for N = 100 and several values of ν as in labels. (a) Results obtained with global dispersal, (b) nearest-neighbor dispersal. Data have been averages over 10 6 generations. Data obtained with nearest-neighbor dispersal are also averaged over 200 independent realization of the initial distribution of preference sites. The continuous curve represents the simple prediction P occ = (1 + γ)/(2 + γ) (see text).

  4) with σ 2 = (1 + ν)/[4(N + 1)ν -4]. This shows that the relative fluctuations around the maximum have the expected characteristic size 1/ √ N .

Acknowledgments

We thank P. Gerlee, K. H. Andersen, M. A. Muñoz and L. Bach for comments on a preliminary version of the manuscript.

A Condition for coexistence in the neutral model with global dispersal

For γ = 0, i.e. in the absence of habitat preference the model with global dispersal reduces to the Moran model with mutation [START_REF] Moran | Random processes in genetics[END_REF]. In this case the transition rates can be expressed in terms of a single population N A .

In particular, equations (2) reduce to the birth and death rates of population N A , i.e. W N A →N A ±1 = W ± (N A ) which, incorporating the effect of immigration, read [START_REF] Karlin | On a genetic model of Moran[END_REF]:

These rates can be used to calculate the equilibrium distribution P (N A ) through the detailed-balance relation