Two-dimensional short-range interacting attractive and repulsive Fermi gases at zero temperature
aa r X i v : . [ c ond - m a t . qu a n t - g a s ] M a r EPJ manuscript No. (will be inserted by the editor)
Two-dimensional short-range interactingattractive and repulsive Fermi gases at zerotemperature.
Gianluca Bertaina a Institute of Theoretical Physics, Ecole Polytechnique F´ed´erale de Lausanne EPFL, CH-1015Lausanne, Switzerland
Abstract.
We study a two-dimensional two-component Fermi gas withattractive or repulsive short-range interactions at zero temperature. Weuse Diffusion Monte Carlo with Fixed Node approximation in order tocalculate the energy per particle and the opposite-spin pair distributionfunctions. We show the relevance of beyond mean field effects and verifythe consistency of our approach by using Tan’s Contact relations.
Low dimensional configurations of degenerate Fermi and Bose gases are being objectof many experimental and theoretical studies [1] as a means of simulating stronglycorrelated systems and as interesting systems per se, from a fundamental point ofview.Many important areas of investigation on Fermi gases are being extended totwo-dimensional (2D) configurations, which have already been pursued in the three-dimensional (3D) case, such as the BCS-BEC crossover in a superfluid gas with reso-nantly enhanced interactions and the possible onset of itinerant ferromagnetism in agas with repulsive interactions. Other interesting phenomena that should be promi-nent in the 2D case are the Fulde-Ferrell-Larkin-Ovchinnikov superfluid state, theultracold gases analogue of the quantum Hall effect in presence of non-abelian gaugefields and the long range correlations induced by dipolar interactions.In particular, 2D ultracold Fermi gases have received a lot of attention in recentyears, having been realized using highly anisotropic pancake-shaped potentials. Den-sity profiles of the clouds have been measured using in situ imaging [2,3]; the singleparticle spectral function has been measured by means of rf spectroscopy [4], thepresence of many-body polaronic states and their transition to molecular states havebeen characterized [5].On the theoretical side, the solution of the BCS equations for the 2D attractiveFermi gas has been first investigated in [6] and later in [7]. The perturbative anal-ysis of the 2D repulsive Fermi gas has been performed in [17]. Recent experimentsare in a regime where beyond mean-field contributions become relevant. In such arespect, many results are now available especially for the highly imbalanced case [9],both for the ground-state and for the so-called Upper Branch (UB). The 2D case is a e-mail: [email protected] Will be inserted by the editor particularly challenging from the point of view of theory, being a marginal case infield theories, when the leading dependence of the observables on the coupling is nonalgebraic. Recently we have obtained the first determination using Quantum MonteCarlo methods of the equation of state at T = 0 of a balanced homogeneous 2D Fermigas in the BCS-BEC crossover [8]. Such an equation of state has been positively com-pared to experiments using the local density approximation [3], it has been used todiscuss recent results on collective modes in a 2D Fermi gas [10] and to judge thetemperature dependence of the Contact in 2D trapped gases [11].In this article we report new results concerning the equation of state of the repul-sive gas and its density-density correlation function. Moreover we emphasize the needof a consistent choice for the coupling constants and the coefficients of the beyondmean-field contributions to energy, both in the strongly interacting molecular regimeand in the weakly interacting Fermi liquid regime. In Section 2 we introduce themodel potentials and we discuss the trial nodal surfaces used in the Quantum MonteCarlo method; in Section 3 we present the equation of state of the weakly attractiveor repulsive gases, we show results for the Upper Branch of the attractive gas andwe discuss the equation of state of the composite bosons in the molecular regime;finally in Section 4 we discuss the extraction of the Contact from the density-densitycorrelation function and we show results for the Contact of the repulsive gas. We consider a homogeneous two-component Fermi gas in 2D described by the Hamil-tonian H = − ~ m N ↑ X i =1 ∇ i + N ↓ X i ′ =1 ∇ i ′ + X i,i ′ V ( r ii ′ ) , (1)where m denotes the mass of the particles, i, j, ... and i ′ , j ′ , ... label, respectively,spin-up and spin-down particles and N ↑ = N ↓ = N/ N being the total number ofatoms. We model the interspecies interatomic interactions using three different typesof model potentials: an attractive square-well (SW) potential V ( r ) = − V for r < R ( V > V ( r ) = 0 otherwise; a repulsive soft-disk (SD) potential V ( r ) = V for r < R ( V > V ( r ) = 0 otherwise; and a hard-disk (HD) potential V ( r ) = ∞ for r < R and V ( r ) = 0 otherwise. Due to the logarithmic dependence on energy ofthe phase shifts in 2D, different definitions of the scattering length have been used [12];we define the scattering length a D = R in the case of the HD potential, so thatfor the SD potential one gets a D = R e − I ( κ ) /κ I ( κ ) and for the SW potential a D = R e J ( κ ) /κ J ( κ ) , where J and I are the Bessel and modified Besselfunctions of first kind and κ = p V mR / ~ , κ = p V mR / ~ . In order to ensurethe diluteness of the attractive gas we use nR = 10 − , where n is the gas numberdensity. The Fermi wave vector is defined as k F = √ πn , and provides the energyscale ε F = ~ k F / m . In order to assure universality in the repulsive case we checkthat the results for the HD potential and for the SD potential with different valuesof R are compatible.For the attractive SW potential in 2D the scattering length is always non negativeand diverges at null depth and at the zeros of J , corresponding to the appearance ofnew two-body bound states in the well. For a SW potential therefore a bound stateis always present, no matter how small the attraction is. The shallow dimers havesize of order a D and their binding energy is given by ε b = − ~ / ( ma D e γ ), where γ ≃ .
577 is Euler-Mascheroni’s constant. A different and widely used definition ofthe 2D scattering length is b = a D e γ /
2, such that ε b = − ~ /mb , analogously to ill be inserted by the editor 3 ´ ´ ´ V (cid:144) ¶ F a D (cid:144) R (cid:144) k F R (a) V R (cid:144) Ñ a D (cid:144) R a D = R (b) Fig. 1. (a) Scattering length of the SW potential, as a function of the depth V /ε F , in theregime used in the simulations. The dashed line corresponds to 1 /k F a D = 1. (b) Scatteringlength of the SD potential, as a function of the barrier V . The dotted line corresponds toone of the values of R used in the simulations. the 3D case. The dependence of a D on the depth V in the region where the wellsupports only one bound state is shown in Fig. 1(a); the dashed line indicates thevalue of scattering length at which the bare molecules have a size comparable to themean interparticle distance 1 /k F . The region k F a D ≫ k F a D ≪ k F a D ∼
1; this can be seen at the two-bodylevel by considering the low-energy scattering amplitude f ( k ) = 2 π/ [log(2 /ka D e γ ) + iπ/
2] [12], which is enhanced for k ∼ /a D (with logarithmic accuracy).For the repulsive SD potential the scattering length is always positive and smallerthan R , like in 3D (see Fig. 1(b)). Although the real physical potentials always havean attractive part, the purely repulsive models that we use can be useful in describingthe cases in which the scattering length is small and positive compared to the meaninterparticle distance, when it is possible to prepare metastable repulsive gas-likestates without a significant production of molecules [5]. We use the fixed-node diffusion Monte Carlo (FN-DMC) method. This numericaltechnique solves the many-body Schroedinger equation by an imaginary time pro-jection of an initial guess of the wavefunction. Provided that the initial guess has afinite overlap with the true ground-state, this method provides the exact energy of asystems of bosons, with a well controllable statistical error. For fermions, FN-DMCyields an upper bound for the ground-state energy of the gas, resulting from an ansatzfor the nodal surface of the many-body wave function that is kept fixed during thecalculation (see Refs. [14]).The fixed-node condition is enforced using an initial and guiding trial functionthat we choose of the standard form ψ T ( R ) = Φ S ( R ) Φ A ( R ), namely the productof a purely symmetric and a purely antisymmetric term. Φ A satisfies the fermionicantisymmetry condition and determines the nodal surface of ψ T , while Φ S is a positivefunction of the particle coordinates and is symmetric in the exchange of particles withequal spin (Jastrow function).Two opposite regimes are described by the Φ A component. The deep attractiveBEC regime, where the opposite-spin fermions are expected to pair into a condensate Will be inserted by the editor of dimers, can be described by an antisymmetrized product Φ A ( R ) = A (cid:0) φ ( r ′ ) φ ( r ′ ) ...φ ( r N ↑ N ↓ ) (cid:1) of pairwise orbitals φ corresponding to the two-body bound state of the potential V ( r ). This wavefunction has been proposed by Leggett as a projection of the grand-canonical BCS function to a state with a finite number of particles, later extended tothe polarized case [13] and extensively used in 3D Quantum Monte Carlo simulations[14]. The weakly interacting regime, where a Fermi liquid description is expected tobe valid, can be instead described by a typical Jastrow-Slater (JS) function with Φ A ( R ) = D ↑ ( N ↑ ) D ↓ ( N ↓ ), namely the product of the plane-wave Slater determinantsfor spin-up and spin-down particles. This description is expected to hold both in theweakly repulsive branch and in the attractive BCS regime of a weakly interacting gaswhere the effect of pairing on the ground-state energy is negligible.The symmetric part is chosen of the Jastrow form Φ S ( R ) = J ↑↓ ( R ) J ↑↑ ( R ) J ↓↓ ( R ).The diluteness of the gas allows us to consider just two-body Jastrow functions, with J ↑↓ ( R ) = Q i,i ′ f ↑↓ ( r i,i ′ ), J ↑↑ ( R ) = Q i 2, both in the Jastrow factors and in themolecular orbitals in the BCS wavefunction; this is obtained by smoothly matchingeach of the previously described functions to a sum of exponentials f e ( r ) = c + c [exp( − µ ( L − r )) + exp( − µr )] at some healing length ¯ R , with µ and ¯ R parameters tobe optimized. The PBC also select a specific set of compatible plane-waves to be usedin the JS wavefunction. It is known that finite-size effects of the JS wavefunction canbe strongly depressed by considering numbers of fermions which provide a maximallysymmetric Fermi surface (closed shells): we choose N/ N/ ill be inserted by the editor 5 -0.020.000.020.040.060.080.100.120.14 -0.4 -0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 0.4 E / N - m f [ E F G ] g eq. (2)fitSW JS-DMCSW UB-VMCHD JS-DMC 0.00.20.40.60.81.01.21.41.61.82.0 -0.4 -0.2 0.0 0.2 0.4 E / N [ E F G ] Fig. 2. Energy per particle of a weakly interacting 2D Fermi gas (inset), with the mean-field contribution subtracted (main figure). The attractive regime corresponds to g < 0, therepulsive regime to g > 0. The results from the simulations with the SW and HD potentialsare compared to the Fermi liquid theory result (2). The dashed line indicates the quadraticfit to the SW JS-DMC data. Empty symbols refer to calculations with the SD potential,with R /a D = 2 , 4, only shown for two representative values of g . The Fermi Liquid Theory (FL) prediction for the zero temperature equation of stateof a weakly short-range interacting 2D gas [17] is the following: E F L /N = E FG (cid:2) g + (3 − g (cid:3) , (2)where E F G = ~ k F / m = ε F / g is the coupling constant. The peculiar logarithmic dependence of the2D scattering amplitude on the available kinetic energy leads to an arbitrary de-pendence of the 2 D coupling constant g = 1 / log( E A /E K ) on a reference energy E K , the E A parameter being relative to the specific choice of the potential. Whenfinite-range effects are negligible E A = 4 ~ / ( ma D e γ ), which is equal to | ε b | in caseof an attractive potential. Since we consider a weakly interacting Fermi liquid, wemake the appropriate choice E K = 2 ε F . Therefore we use the coupling constant g = − / log( na D c ) = − / k F b ) with c = πe γ / 2. An important remark is inorder: while the coefficient of the first order term is independent of c , the reportedvalue of the coefficient of g depends on the chosen c , since a modification of c givesrise to higher order contributions to the equation of state.In Fig. 2 we show the FN-DMC results with the JS nodal surface. Negative cou-plings g < Will be inserted by the editor main figure we show the energy with the mean-field result E F G (1 + 2 g ) subtracted,in order to emphasize beyond mean-field contributions. Although the three series ofdata obviously confirm the mean-field prediction, interesting different behaviors arefound when checking the residual contributions to energy. In the g < A of the second order term obtaining A = 0 . − . A ≈ − c = 2 π (as has beendone in [8]). On the opposite side g > 0, the Upper Branch (UB-VMC) results areapproximately close to the second order fit done for the g < g . The discrepancy is due to the lack of optimization of the UB wavefunction, which isproblematic since there is not a variational principle for this state. The JS results forthe SD potential with R = 2 a D are compatible with those with the HD potential.We also perform simulations with shallower SD potentials ( R = 4 a D ), which resultin departures from universality at g ≈ . 25. Agreement between eq. (2) and the HDJS-DMC results is found up to g ≈ . 1, which corresponds to na ≈ − . The be-yond mean-field contributions in this paramagnetic phase soon start to decrease. In arelated system a ferromagnetic transition has been excluded [18]; further investigationconcerning the strongly repulsive regime is however beyond the scope of this paper.The intriguing difference in the beyond mean-field terms, between the attractive caseand the purely repulsive case, would also deserve further investigation, hopefully withan analytical treatment. We now pass to discuss the main result published in [8], namely the characterization ofthe composite bosons equation of state in the 2D BCS-BEC crossover. The 2D mean-field BCS equations can be analytically solved [6,7] along the BCS-BEC crossover, interms of the interaction coefficient x = | ε b | / ε F . For the BCS order parameter oneobtains ∆ = 2 ε F √ x , while for the energy per particle one obtains E/N = E F G (1 − x )and for the chemical potential µ = ε F (1 − x ). For small binding energy one recoversthe non interacting limit; when x ∼ 1, the chemical potential becomes zero andthen negative, so that the role of the dimers becomes more important; for very strongbinding the chemical potential of the fermions is equivalent to half the binding energyof a molecule, so the system is made of non interacting bosons. Although very usefulfor providing a global self-consistent picture and for setting a stringent variationalupper bound to the energy per particle, the BCS solution fails in various aspects. Inthe BCS regime it neglects the Hartree-Fock contributions to energy (discussed in theprevious Section), which are dominant, since the gap is small. In the BEC regime itmisses the correct interaction energy between the bosons. In general it is not able toreproduce the logarithmic dependence of energy on the density, which is typical of2D.In Fig. 3 we show the FN-DMC results (from [8]) for the equation of state of ashort-range attractively interacting 2D gas as a function of the interaction parame-ter η = log( k F a D ) in units of E F G , with ε b / η . 1, while the JSfunction is more favorable for larger values of η . In the deep BEC regime, besides themolecular contribution, the remaining fraction of energy corresponds to the interac-tion energy of the bosonic dimers.In the BEC regime the FN-DMC results are fitted with the equation of state ofa 2D gas of composite bosons. The functional form is derived from the analysis ofBeane [19], in the framework of quantum field theory; it also corresponds to previousanalytical and DMC studies (see [19] and references therein). Following Beane weintroduce the running coupling g d ( λ ) in terms of the scattering length of the bosons ill be inserted by the editor 7 E / N + | ε b | / [ E F G ] η =log(k F a ) DMC JSDMC BCSfitmean fieldComposite bosonsleading order Fig. 3. Energy per particle of an attractive 2D Fermi gas along the BCS-BEC crossoverwith the bare binding energy of the molecules subtracted. The BEC regime corresponds to η → −∞ , the BCS regime corresponds to η → + ∞ , while the resonant regime correspondsto η ∼ 1. The various curves are described in the text. Reprint of the data in [8]. a d and the particle density of the bosons n d : g d ( λ ) = − / log ( n d λπ e γ a d ), where λ is an arbitrary dimensionless cutoff parameter, which is present due to the truncationin the perturbative expansion; it is the analogous of the c parameter discussed inSubsection 3.1. In the following m d is the mass of the bosons. Up to the second orderin the running coupling, one can express the energy density in the following way E = 2 π ~ n d m d g d ( λ ) " g d ( λ ) (cid:16) log ( πg d ( λ )) − log λπ + 12 (cid:17) . (3)It is evident from the above expression that fixing the scattering length a d whichappears in the definition of g d ( λ ) is not sufficient for determining the energy densityif one does not declare its choice for λ , that is the form of the coupling constant.Notice again that the coefficient of the second order term does depend on the choiceof λ . A convenient choice for simplifying the expression is to set λ = e − γ /π , so thatwe can introduce the coupling g d = − / log ( n d a d ) and we obtain E = 2 π ~ n d m d g d " g d (cid:16) log ( πg d ) + 2 γ + 12 (cid:17) . (4)Now let us consider the case when the bosons are dimers, consisting of two pairedfermions with mass m = m d / n = 2 n d . There must exist aregime where the binding is so tight that the equation of state of such compositebosons is the same as in the case of point-like bosons, with the simple replacement a d → αa D . Let us therefore introduce η = log ( k F a D ) in the previous expres-sion, so that the composite bosons coupling turns to be g d = − / log ( nα a D / 2) =1 / (log 4 π − η − α ). In such a situation the energy per fermion can be writtenin the following way: EN F = − | ε b | ε F g d " g d (cid:16) log ( πg d ) + 2 γ + 12 (cid:17) , (5) Will be inserted by the editor g ↑↓ (r) r/a JS SWshort-range (a) JS HDshort-range (b) Fig. 4. (a) g ↑↓ correlation function for the attractive gas at η = 2 . 15 ( g = − . g ↑↓ correlation function for the repulsiveHD gas at g = 0 . where −| ε b | / a d = 0 . a D , in agreement with the four-bodycalculation in Ref. [12]. In order to exemplify the importance of the second orderexpansion (5), in Fig. 3 we also show the leading linear contribution alone, with thesame choice of the coupling constant. Without the second order term it would beimpossible to determine a d with accuracy.To our knowledge a theoretical analysis for the running coupling of a 2D fermionicnon-relativistic fluid, analogous to the one in Beane [19] for the bosonic counterpart,is still lacking. It would be highly useful for interpreting the results of Subsection3.1, putting the comparison of the Quantum Monte Carlo data and the Fermi liquidequation of state on firmer grounds.The solid curve in Fig. 3 is a global fit of the data; the function to be fitted isa piece-wise defined function of η , the matching point being a fitting parameter andthe two matched functions being the ratio of second order polynomials of η such thatthe global function is continuous at the matching point up to the second derivativeand the extreme regimes coincide with the known perturbative results. The Contact parameter C is a property of short-range interacting gases, measuringthe amount of pairing between the particles and relating a large number of observableswith each other [20]. For example it can be obtained from the derivative of the equa-tion of state with respect to the coupling parameter C = (2 πm/ ~ ) d ( nE/N ) /d (log k F a D )or from the short-range behavior of the antiparallel pair distribution function g ↑↓ ( r ) → r → C/k F log ( r/a D ) [21] (where r → R < r ≪ /k F , R being the range ofthe potential).We calculate the g ↑↓ correlation function along the crossover using FN-DMC andextract the Contact by means of the following more refined formula, which betterdescribes the behavior at small r (see Fig. 4(a)): g ↑↓ ( r ) → r → Ck F " − log (cid:18) ra D (cid:19) (cid:18) ra D e γ (cid:19) ! + (cid:18) ra D e γ (cid:19) . (6)The additional terms permit a more precise determination of the Contact in the deepBEC regime, where the peak in g ↑↓ at r → a D is quite smallwith respect to the mean interparticle distance. ill be inserted by the editor 9 C [ k F ] gEq. of stateg ↑↓ HD JS Fig. 5. Contact parameter in the weakly repulsive regime. The solid line corresponds to thecalculation from the derivative of the equation of state (2). The triangles correspond to theextraction from the short-range behavior of g ↑↓ for the HD gas. The results for the BCS-BEC crossover have already been presented in Fig. 4of Ref. [8], where they have also been compared to the Contact extracted form thederivative of the global fit to the energy. The overall agreement between the two deter-minations of C is a useful check of the accuracy of the trial wavefunctions used in theFN-DMC approach. Small deviations in the region log( k F a ) ∼ C from the g ↑↓ of therepulsive HD gas. The short-range details are model dependent, nevertheless it is stillpossible to find an intermediate region which is universal. In this case the leading orderformula for g ↑↓ is used, since eq. (6) is valid only for r ≪ a D , which corresponds tothe non universal region for the repulsive gas. In Fig. 5 we compare C extracted from g ↑↓ to its analytical expression obtained from eq. (2): C/k F = [1 + (3 + 4 log 2) g ] g .Similarly to Fig. 2 deviations appear around g = 0 . To conclude we have detailed the functional forms needed for accurately extract usefulinformation from the FN-DMC simulations of 2D Fermi gases. Both for the energyin the weakly interacting case and in the BEC regime and for the g ↑↓ correlationfunction the knowledge of the next-to-leading order correction is crucial in order toavoid ambiguity in the measured properties. We have also shown new results on theequation of state and the Contact of the repulsive gas in the weakly interacting regime. I want to acknowledge the support of the University of Trento and the INO-CNR BECCenter in Trento, where a large part of the results have been obtained, and the valuablediscussions with S. Giorgini.0 Will be inserted by the editor References 1. S. Giorgini, L.P. Pitaevskii and S. Stringari, Rev. Mod. Phys. , (2008) p. 1215; I. Bloch,J. Dalibard and W. Zwerger, Rev. Mod. Phys. , (2008) p. 885.2. K. Martiyanov, V. Makhalov and A. Turlapov, Phys. Rev. Lett. , (2010) p. 030404.3. A. Orel et al. , New J. Phys. , (2011) p. 113032.4. B. Fr¨ohlich et al. , Phys. Rev. Lett. , (2011) p. 105301; M. Feld, et al. , Nature ,(2011) p. 75; A. Sommer et al. , Phys. Rev. Lett. , (2012) p. 045302.5. M. Koschorreck et al. , Nature , (2012) p. 619.6. K. Miyake, Prog. Theor. Phys. , (1983) p. 1794.7. M. Randeria, J.-M. Duan and L.-Y. Shieh, Phys. Rev. Lett. , (1989) p. 981; Phys. Rev.B , (1990) p. 327.8. G. Bertaina and S. Giorgini, Phys. Rev. Lett. 106, (2011) p. 110403.9. M. Parish, Phys. Rev. A , (2011) p. 051603, S. Z¨ollner, G. Bruun, and C. Pethick,Phys. Rev. A , (2011) p. 021603, M. Klawunn and A. Recati, Phys. Rev. A , (2011)p. 033607, R. Schmidt et al. , Phys. Rev. A. , (2012) p. 021602.10. E. Vogt et al. , Phys. Rev. Lett. (2012) p. 070404, J. Hofmann, Phys. Rev. Lett. , (2012) p. 185303, E. Taylor and M. Randeria, Phys. Rev. Lett. (2012) p. 135301.11. Fr¨ohlich, B. et al. , Phys. Rev. Lett. (2012) p. 130403.12. D.S. Petrov, M.A. Baranov and G.V. Shlyapnikov, Phys. Rev. A , (2003) p. 031601(R);D.S. Petrov and G.V. Shlyapnikov, Phys. Rev. A. , (2001) p. 012706.13. A.J. Leggett, Modern Trends In The Theory Of Condensed Matter, Lecture Notes InPhysics Vol. 115 , (Springer-Verlag, Berlin 1980) p. 13, J.P. Bouchaud, A. Georges and C.Lhuillier, J. Phys. (Paris) , (1988) p. 553.14. P.J. Reynolds et al. , J. Chem. Phys. , (1982) p. 5593; G.E. Astrakharchik et al. , Phys.Rev. Lett. , (2004) p. 200404; S.-Y. Chang et al. , Phys. Rev. A , (2004) p. 043602;G.E. Astrakharchik et al. , Phys. Rev. Lett. , (2005) p. 230405.15. Pilati, S. et al. , Phys. Rev. Lett. (2010) p. 030405; Chang, S.-Y. et al. Proc. Natl.Acad. Sci. (2011) p. 51.16. Ceperley, D. M., Phys. Rev. B (1978) p. 3126; Ceperley, D. M. and Alder, B. J. ,Phys. Rev. B (1987) p. 2092.17. P. Bloom, Phys. Rev. B , (1975) p. 125; J.R. Engelbrecht, M. Randeria and L. Zhang,Phys. Rev. B , (1992) p. 10135; J.R. Engelbrecht and M. Randeria, Phys. Rev. B ,(1992) p. 12419.18. Drummond, N. et al. , Phys. Rev. B (2011) p. 195429.19. G.E. Astrakharchik et al. , Phys. Rev. A , (2009) p. 051602(R); C. Mora and Y. Castin,Phys. Rev. Lett. , (2009) p. 180404; S.R. Beane, Phys. Rev. A , (2010) p. 063610.20. S. Tan, Ann. Phys. , (2008) 2952 and ibid , (2008) p. 2971.21. F. Werner and Y. Castin, Phys. Rev. A.86