Gap solitons in elongated geometries: the one-dimensional Gross-Pitaevskii equation and beyond
aa r X i v : . [ c ond - m a t . qu a n t - g a s ] M a y Phys. Rev. A , 053610 (2011) Gap solitons in elongated geometries: the one-dimensional Gross-Pitaevskii equationand beyond
A. Mu˜noz Mateo ∗ and V. Delgado † Departamento de F´ısica Fundamental II, Facultad de F´ısica,Universidad de La Laguna, 38206 La Laguna, Tenerife, Spain
Boris A. Malomed ‡ Department of Physica Electronics, School of Electrical Engineering,Faculty of Engineering, Tel Aviv University, Tel Aviv 69978, Israel (Dated: 15 February 2011)We report results of a systematic analysis of matter-wave gap solitons (GSs) in three-dimensionalself-repulsive Bose-Einstein condensates (BECs) loaded into a combination of a cigar-shaped trapand axial optical-lattice (OL) potential. Basic cases of the strong, intermediate, and weak radial(transverse) confinement are considered, as well as settings with shallow and deep OL potentials.Only in the case of the shallow lattice combined with tight radial confinement, which actuallyhas little relevance to realistic experimental conditions, does the usual one-dimensional (1D) cubicGross-Pitaevskii equation (GPE) furnish a sufficiently accurate description of GSs. However, theeffective 1D equation with the nonpolynomial nonlinearity, derived in Ref. [Phys. Rev. A ,013617 (2008)], provides for quite an accurate approximation for the GSs in all cases, includingthe situation with weak transverse confinement, when the soliton’s shape includes a considerablecontribution from higher-order transverse modes, in addition to the usual ground-state wave functionof the respective harmonic oscillator. Both fundamental GSs and their multipeak bound states areconsidered. The stability is analyzed by means of systematic simulations. It is concluded thatalmost all the fundamental GSs are stable, while their bound states may be stable if the underlyingOL potential is deep enough. PACS numbers: 03.75.Lm, 05.45.Yv
I. INTRODUCTION
Matter-wave gap solitons (GSs) are localized modesthat can be created in elongated Bose-Einstein conden-sates (BECs) loaded into one-dimensional (1D) optical-lattice (OL) potentials, with the intrinsic nonlinearityinduced by repulsive interactions between atoms. GSshave been the topic of a large number of original works.Results produced by these studies were summarized inseveral reviews [1–5]. Within the framework of the mean-field approximation, the description of matter-wave pat-terns is based on the Gross-Pitaevskii equation (GPE)for the macroscopic wave function ψ [6]: i ~ ∂ψ∂t = (cid:18) − ~ M ∇ + V ⊥ ( r ⊥ ) + V z ( z ) + gN | ψ | (cid:19) ψ, (1)which has proven to be very successful in reproducingexperimental results for the zero-temperature BEC (forinstance, for multiple vortices, as shown in detail in Refs.[7, 8]). In this equation, M is the atomic mass, the normof the wave function ψ is unity, N is the number of atoms,and g = 4 π ~ a/M is the interaction strength with a beingthe s -wave scattering length. In this work, we consider ∗ Electronic address: [email protected] † Electronic address: [email protected] ‡ Electronic address: [email protected] only repulsive interactions (i.e., a > V ⊥ ( r ⊥ ) = (1 / M ω ⊥ r ⊥ is the radial-confinement poten-tial and V z ( z ) is the axial potential, which may includethe axial harmonic trap and the 1D OL, with depth V and period d : V z ( z ) = (1 / M ω z z + V sin ( πz/d ) . (2)The energy scale in the underlying (no-interaction) lin-ear problem is set by the OL’s recoil energy, E R = ~ π / (2 M d ).Most commonly, theoretical studies of GSs in elon-gated geometries are carried out in terms of the 1D GPE[1, 2, 9–13]: i ~ ∂φ∂t = − ~ M ∂ φ∂z + V z ( z ) φ + g N | φ | φ, (3)where g = 2 a ~ ω ⊥ . This is an effective evolution equa-tion for the axial dynamics —described by the 1D wavefunction φ ( z, t )— which can be derived from the fullthree-dimensional (3D) GPE after averaging out the ra-dial degrees of freedom under the assumption that theradial confinement is so tight that the transverse dynam-ics are frozen to zero-point oscillations. These condi-tions, however, are not easy to realize using typical ex-perimental parameters [2], and when such conditions arenot met the transverse excitations can no longer be ne-glected, making an essentially 3D analysis necessary. Inthis work, we aim to investigate different physically rele-vant regimes and capture 3D effects in the generation andstability of matter-wave GSs in elongated BECs. Thisanalysis should make it possible to determine the rangeof validity of the 1D GPE for the description of the GSsunder typical experimental conditions, as well as charac-teristic features of the fundamental and higher-order GSsin realistic situations.In principle, when the transverse excitations are rel-evant, Eq. (3) fails and one has to resort to the full3D GPE (1), as recently done in Ref. [14], or, alterna-tively, use extended 1D models that, within the frame-work of certain assumptions, take into account effects ofhigher-order radial modes on the axial dynamics of thecondensate, leading to effective 1D equations with thecubic-quintic [15] or nonpolynomial [16, 17] nonlinear-ities. In this work, we will consider both the full 3DGPE and the effective 1D model with a nonpolynomialnonlinearity, which was derived, for the case of the self-repulsive nonlinearity, in Ref. [17]. The latter one, whichrepresents a simple generalization of the usual 1D GPE,that reduces to Eq. (3) in the appropriate limit, hasdemonstrated an excellent quantitative agreement withexperimental observations [18–21] (chiefly, for delocalizeddark solitons, which are natural patterns in the case ofthe self-repulsion; however, the comparison was not re-ported before for localized GS modes). We will demon-strate that, while the range of applicability of the 1DGPE (3) is severely limited in realistic situations, theabove-mentioned generalization gives a good descriptionof stationary matter-wave GSs in virtually all cases ofpractical interest.The paper is organized as follows. The model is formu-lated in Section II, where we also recapitulate the deriva-tion of the effective nonpolynomial 1D equation followingthe lines of Ref. [17], as this derivation is essential forthe presentation of the results for the GSs. The mainfindings are collected in Sections III, where families ofGS solutions are reported in several physically relevantregimes for strong, intermediate, and weak transverseconfinement and shallow or deep OL potential. Bothfundamental GSs and their multipeak bound states areconsidered. Only in the case of tight confinement com-bined with a shallow OL does the ordinary 1D GPE pro-vide for a sufficiently accurate description of the GSs incomparison with results of full 3D computations. On theother hand, the effective nonpolynomial 1D equation pro-vides for good accuracy in all cases; for the fundamentalsolitons and their bound states alike. The stability ofthe GSs is studied in Section IV by means of systematicsimulations of the evolution of perturbed solitons. Theconclusion is that the fundamental GSs are stable (ex-cept in a narrow region close to the upper edge of theband gaps), even in the case of strong deviation fromthe usual 1D description. Multisoliton bound states arestable if the OL potential is deep enough. Conclusionsfollowing from results of this work, including applicabilitylimits for the mean-field approximation and 1D approxi-mations, are formulated in Section V. II. THE MODEL
The model that was developed in Ref. [17] for a BECin the absence of OL potentials resorted to the adiabaticapproximation, to neglect correlations between the trans-verse and axial motions. This approximation assumesthat the axial density varies slowly enough to allow thetransverse wave function to follow these slow variations.Because the OL imposes a new spatial scale, which maybe more restrictive, it is necessary to find out if the adi-abatic approximation remains valid in the presence ofthe OL. To this end, we will now briefly recapitulate thederivation of the effective 1D equation based on this ap-proximation.The starting point is the ansatz based on the factorized3D wave function ψ ( r , t ) = ϕ ( r ⊥ ; n ( z, t )) φ ( z, t ) , (4)with n ( z, t ) being the local condensate density per unitlength characterizing the axial configuration: n ( z, t ) ≡ N Z d r ⊥ | ψ ( r ⊥ , z, t ) | = N | φ ( z, t ) | . (5)To derive Eq. (5), we have assumed that the transversewave function ϕ is normalized to unity. Next, the sub-stitution of Eq. (4) into the 3D equation (1) leads to (cid:18) i ~ ∂φ∂t + ~ M ∂ φ∂z − V z φ (cid:19) ϕ = (cid:18) − ~ M ∇ ⊥ ϕ − ~ M ∂ ϕ∂z + V ⊥ ϕ + gn | ϕ | ϕ (cid:19) φ − ~ M ∂ϕ∂z ∂φ∂z . (6)Because of the very different time scales of the axial andradial motions, it is natural to assume that the slow axialdynamics may be accurately described by averaging Eq.(6) over the fast (radial) degrees of freedom. Doing this,one obtains i ~ ∂φ∂t = − ~ M ∂ φ∂z + V z φ + µ ⊥ φ − ~ M (cid:18)Z d r ⊥ ϕ ∗ ∂ϕ∂z (cid:19) ∂φ∂z , (7) µ ⊥ ( n ) ≡ Z d r ⊥ ϕ ∗ (cid:18) − ~ M ∇ ⊥ − ~ M ∂ ∂z + V ⊥ + gn | ϕ | (cid:19) ϕ. (8)Since n ( z, t ) enters the last term of Eq. (8) merely as anexternal parameter, it is clear that, whenever the axialkinetic energy associated with the transverse wave func-tion may be neglected; namely, K z [ ϕ ] ≡ Z d r ⊥ ϕ ∗ (cid:18) − ~ M ∂ ∂z (cid:19) ϕ ≪ Z d r ⊥ ϕ ∗ (cid:18) − ~ M ∇ ⊥ + V ⊥ ( r ⊥ ) + gn | ϕ | (cid:19) ϕ, (9)the radial dynamics decouple and µ ⊥ ( n ) can be deter-mined without the knowledge of the axial wave function.Actually, the right-hand side of Eq. (9) is the chemi-cal potential of an axially homogeneous condensate char-acterized by the density per unit length n and wavefunction ϕ . For a sufficiently small value of the lineardensity ( an ≪ ϕ ( r ⊥ ) is given by the Gaussian wave function of theground state of the radial harmonic oscillator, with thecorresponding chemical potential being ~ ω ⊥ (1 + 2 an ).As the linear density increases, more radial modes be-come excited and, in general, the ground state of thecorresponding homogeneous condensate involves manyharmonic-oscillator modes.In previous works, it was shown using different ap-proaches that, also in this case, an analytical solutioncan be constructed [17]. In particular, by using a varia-tional approach based on the direct minimization of thechemical-potential functional, it was shown that, for any(dimensionless) linear density an , a very accurate esti-mate for the chemical potential of the ground state isgiven by ~ ω ⊥ √ an [22]. This expression can beeasily extended to incorporate the case in which the con-densate contains an axisymmetric vortex [17] (see alsoRef. [23], where the intrinsic vorticity was includedinto the derivation of the 1D nonpolynomial nonlinearSchr¨odinger equation in the case of self-attraction). How-ever, in this work we restrict the consideration to zerovorticity. Using Eq. (9), we see that a sufficient condi-tion for the second derivative in z appearing in Eq. (8)to be negligible is K z [ ϕ ] ≪ ~ ω ⊥ √ an . (10)Taking into account an estimate, K z [ ϕ ] ∼ ~ / (2 M ∆ z ),where ∆ z is the characteristic length scale in the axialdirection, Eq. (10) can be rewritten as ~ M ∆ z ≪ ~ ω ⊥ √ an . (11)When this condition holds, µ ⊥ ( n ) coincides, to a goodapproximation, with the transverse local chemical poten-tial of the stationary radial GPE: (cid:18) − ~ M ∇ ⊥ + V ⊥ ( r ⊥ ) + gn | ϕ | (cid:19) ϕ = µ ⊥ ( n ) ϕ, (12)and is given by µ ⊥ ( n ) = ~ ω ⊥ √ an . (13)Finally, substituting Eqs. (13) and (5) into Eq. (7)and taking into account that ϕ is a real normalized wavefunction (which implies the vanishing of the integral inthe last term), we arrive at the following effective 1Dequation to govern the slow axial dynamics of the con-densate: i ~ ∂φ∂t = − ~ M ∂ φ∂z + V z ( z ) φ + ~ ω ⊥ p aN | φ | φ. (14) We note that the contribution from interatomic in-teractions enters the above equation through term ~ ω ⊥ ( p aN | φ | − φ , which vanishes for N → ~ ω ⊥ , which isirrelevant for the dynamics but simplifies the form of theequation and makes the global chemical potential thatfollows from this effective 1D equation exactly coincidingwith the corresponding 3D result.The above derivation demonstrates that the validity ofEq. (14) relies on two conditions:(i) The typical time scale ∆ t of the axial motion mustbe much larger than the typical time scale of the radialmotion ( ∼ ω − ⊥ ).(ii) The axial kinetic energy K z [ ϕ ] associated with thetransverse wave function must be negligible.The former requirement is necessary to allow the ra-dial wave function to adiabatically follow the axial dy-namics. While the specific temporal scale ∆ t dependson the particular initial conditions, typically, in the pres-ence of an OL, the fulfillment of this condition gets moredifficult as the lattice period d decreases, or the lineardensity an increases. In particular, a sufficient condi-tion for the validity of the adiabatic approximation is: d ≫ a ⊥ , an ≪
1, with a ⊥ ≡ p ~ / ( M ω ⊥ ) being the ra-dial harmonic-oscillator length. Nevertheless, such a con-straint, which is hard to satisfy in realistic situations, isnot a necessary condition. Actually, for stationary states∆ t → ∞ and condition (i) always holds. In the presentwork we are interested in this case, which is the mostrelevant one to seek for matter-wave GSs.A sufficient condition for the fulfillment of condition(ii) is given by the inequality (11). In the presence of anOL, the characteristic length scale ∆ z is typically on theorder of the lattice period d , hence Eq. (11) becomes E R ≪ ~ ω ⊥ √ an , (15)where E R is the corresponding recoil energy. We willdemonstrate that this condition is well satisfied for sta-tionary GSs in most cases of interest.It is clear that, when the linear density is small enough( an ≪ E R ~ ω ⊥ , V E R , ω z ω ⊥ , N aa ⊥ . (16)As said above, the recoil energy E R sets the energy scaleof the underlying linear problem, while the nonlinear cou-pling constant N a/a ⊥ determines the order of magnitudeof the mean-field interaction energy in units of the radialquantum ~ ω ⊥ . Note that the dimensionless linear density an is proportional to N a/a ⊥ .In this work, we assume the axial confinement to be soweak that the corresponding harmonic-oscillator length, a z ≡ p ~ / ( M ω z ), is much larger than the lattice period d . Under these conditions, it is safe to neglect the mod-ulation induced in the condensate density by the axialharmonic trap and set ω z = 0, so that we are left with auniform 1D lattice potential acting along the axial direc-tion. III. MATTER-WAVE GAP SOLITONS INDIFFERENT PHYSICAL REGIMES
The stationary solutions of the 3D GPE are wave func-tions of the form ψ ( r , t ) = ψ ( r ) exp( − iµt ), with ψ obeying the time-independent GPE: µψ = (cid:18) − ~ M ∇ + V ⊥ ( r ⊥ ) + V z ( z ) + gN | ψ | (cid:19) ψ , (17)where µ is the chemical potential. When Eq. (14) isapplicable, one can instead generate the stationary axialwave function φ ( z ) by solving the effective 1D equation µφ = (cid:18) − ~ M ∂ ∂z + V z ( z ) + ~ ω ⊥ p aN | φ | (cid:19) φ , (18)which, in the limit of an ≪
1, reduces to the time-independent version of the usual 1D GPE (3): µφ = (cid:18) − ~ M ∂ ∂z + V z ( z ) + g N | φ | + ~ ω ⊥ (cid:19) φ . (19)Matter-wave GSs are characterized by a chemical po-tential lying in a band gap of the energy spectrum of theunderlying linear system. To construct solutions for real-istic 3D matter-wave GSs in the effectively 1D geometry,we looked for numerical solutions to Eqs. (17) and (18) indifferent physically relevant regimes, and compared theobtained results with those produced by the usual 1DGPE (19) to determine to what extent 3D contributionsare relevant. In particular, the comparison allows us toestimate the accuracy and range of applicability of theeffective 1D equations (18) and (19).Experimentally relevant values for the OL period d range from 0 . . µ m [2], which implies that, for thecondensate of Rb, E R / (2 π ~ ) ranges from 3 . ω ⊥ / π . E R / ~ ω ⊥ & /
4. On the otherhand, for the condensate of Na one has E R / ~ ω ⊥ & Rb condensate is most relevant for the experimental re-alization of GSs in the quasi-1D setting [27, 28], in whatfollows below we use particular parameters of this con-densate to illustrate our results. Nevertheless, the resultsof the present work, which are always expressed in termsof relevant dimensionless parameters, are valid for anyBEC. Actually, this is a direct consequence of the scalingproperties of Eqs. (1), (3), and (14).
A. Tight radial confinement: E R / ~ ω ⊥ ≪ In this case, the quantum of radial excitations is muchgreater than the typical energy scale in the linear prob-lem. Note that, to realize this regime in a Rb conden-sate, even in an OL of period d ≃ . µ m, one needs touse a harmonic trap with radial frequency ω ⊥ / π & E R / ~ ω ⊥ =1 / e µ ≡ ( µ − ~ ω ⊥ ) /E R , (20)as a function of the dimensionless lattice depth, s ≡ V /E R . Since in this work we are interested in GSs withzero vorticity, this diagram, which represents the band-gap structure of the underlying linear problem, has beenobtained from the zero-vorticity solutions of the linearversion of Eq. (17). Regions I, II, and III correspondto the lowest finite band gaps, which separate shadedbands, where linear solutions exist. An important pointis that, within the region of interest in the parameterspace (for V /E R .
25 and up to the third band gap),this 3D diagram is indistinguishable from that obtainedusing the linear version of the effective 1D equation (18)[which, obviously, coincides with the linear version of thestationary 1D GPE (19)]. In fact, Eq. (18) leads to aband-gap diagram that differs from Fig. 1(a) solely inthe band marked by the arrow, which does not appear inthe 1D case. This extra band is, essentially, a replica ofthe lowest one, shifted up in energy by 2 ~ ω ⊥ /E R . Tak-ing into account that the energy spectrum of the radialharmonic oscillator is given by E = (2 n r + | m | + 1) ~ ω ⊥ , (21)where n r = 0 , , , . . . is the radial quantum number and m = 0 , ± , ± , . . . is the axial angular-momentum quan-tum number, it is evident that the additional up-shiftedband corresponds to the first-excited radial mode with m = 0. In the notation of Ref. [14], it corresponds toquantum numbers ( n = 1 , m = 0 , n r = 1), with n beingthe band index of the 1D axial problem. Thus, the ap-pearance of this band is a purely 3D effect that cannot FIG. 1: (Color online) (a) Band-gap structure of a noninter-acting 3D BEC with E R / ~ ω ⊥ = 1 /
10, as a function of thedimensionless lattice depth s ≡ V /E R . The representativecases s = 2 and 20, indicated by vertical dashed lines, areconsidered in more detail in (b) and (c), respectively, whichshow the dimensionless chemical potential e µ as a function ofthe quasimomentum q in the first Brillouin zone (in units of π/d ). The right panels display the location of the Rb gapsolitons considered in this work (points A–G). be accounted for by the above 1D models. Since the en-ergy shift increases as E R / ~ ω ⊥ decreases, it is clear that,within the parameter region of interest, Fig. 1(a) is uni-versal in the sense that it is valid (for both 1D and 3Dsystems) for any E R / ~ ω ⊥ ≪ , , E R / ~ ω ⊥ = 1 / V /E R = 2 or 20. These two representativecases correspond to the dashed vertical lines in Fig. 1(a)and are shown in more detail in Figs. 1(b) and 1(c),which display the corresponding band-gap structure asa function of the quasimomentum in the first Brillouinzone (in units of π/d ). Bold red lines in Figs. 1(b) and1(c) have been obtained from the linear version of theeffective 1D equation (18) while thin black lines corre-spond to results obtained by means of the linear versionof the 3D equation (17) that cannot be reproduced withthe 1D model. As mentioned above, the only appreciabledifference between the 1D and 3D results comes from the contribution of the first-excited radial mode [see the toppart of Fig. 1(c)]. The case of V /E R = 2, displayedin Fig. 1(b), corresponds to the condensate in relativelyshallow OL potentials. Under these circumstances, thelinear energy spectrum does not differ too much fromthat corresponding to the translationally invariant case( V z = 0), and the band-gap structure exhibits wide en-ergy bands separated by relatively small gaps, which aretangible only in the lowest part of the spectrum. On thecontrary, for V /E R = 20 [the tight-binding regime, seeFig. 1(c)] the condensate is trapped in the deep periodicpotential. In this case, the lowest part of the linear spec-trum is dominated by large gaps separating relativelynarrow energy bands. The panels on the right side ofFigs. 1(b) and 1(c) show the location, with respect tothe corresponding linear band-gap diagram, of the GSsthat will be considered below (points A–G). The horizon-tal axes in these panels indicate the number of atoms ineach GS for the Rb condensates (with the s -wave scat-tering length a = 5 .
29 nm) in the trap with radial fre-quency ω ⊥ / π = 2400 Hz. For other parameter values,the GS family is described by the same plots, if consid-ered in terms of the above-mentioned nonlinear couplingconstant, N a/a ⊥ [ a = 5 .
29 nm and ω ⊥ / π = 2400 Hzcorrespond to a/a ⊥ = 0 . a ′ and transverse confinement radius a ′⊥ , which is tantamount to the GS with particular val-ues of a , a ⊥ and N , is given by N ′ = aa ′⊥ a ′ a ⊥ N. (22)GS solutions have been obtained by numerically solv-ing the full 3D equation (17), as well as the effective1D equations (18) and (19). To this end, we have im-plemented a Newton continuation method, based on aLaguerre–Fourier spectral basis, that uses the chemicalpotential µ as a continuation parameter. To ensure theconvergence, methods of this kind require initiation of theiterative procedure with a sufficiently good initial guess.This means that one needs a good estimate for both theaxial and the radial parts of the wave function. Our com-putations used the following initial ansatz for the wavefunction: ψ ( r ⊥ , z ) = 1Γ a ⊥ √ π exp (cid:18) − r ⊥ a ⊥ (cid:19) φ ( z ) , (23)where the z -dependent condensate’s width, expressed inunits of a ⊥ , is Γ = (cid:0) aN | φ ( z ) | (cid:1) / , (24)and φ ( z ) = p k / k z ) is the axial wave function.The number of particles, N , and k are free adjustableparameters. Note that the radial part of the wave func-tion (23) is the same as the variational solution of Eq.(12), which was found in Ref. [17]. FIG. 2: (Color online) Atom density of the gap soliton corre-sponding to point B in Fig. 1(b), displayed as (a) an isosurfacetaken at 5% of the maximum density and (b) as a color mapalong a cutting plane containing the z axis. (c) Dimension-less axial density an obtained from the 3D wave function,as prescribed by Eq. (5) (open circles) along with the cor-responding predictions from the nonpolynomial 1D equation(18) (solid red line) and the ordinary 1D GPE (19) (dashedblue line).
1. Shallow optical lattice: V /E R = 2 Point A in Fig. 1(b) corresponds to a fundamental GSlocated near the bottom edge of the first band gap. Itlooks qualitatively similar to that shown below in Fig.3, but with the peak linear density an (0) ≃ .
04. Inthis case, the effective 1D equations (18) and (19) yieldresults which are indistinguishable, on the scale of the fig-ures, from those obtained from the full 3D equation (17).This is not surprising because, for these parameters, con-dition (15) is certainly satisfied, and inequality an ≪ e µ = 1 .
75, which implies µ = 1 . ~ ω ⊥ ≃ ~ ω ⊥ , hencethe radial wave function of the condensate should notdiffer too much from the ground state of the correspond-ing harmonic oscillator. These fundamental GSs, how-ever, can only accommodate 9 particles (for the Rbcondensate), which is too small to use the mean-field ap-proximation. Yet these solutions play an important roleas building blocks of higher-order GSs. Point B in Fig.
FIG. 3: (Color online) Same as Fig. 2 but for the gap solitoncorresponding to point C in Fig. 1(b). e µ = 1 .
75 and N = 35. Note that,for the relatively shallow OL ( V /E R = 2), interferenceeffects arising from the overlap between fundamental GSssitting in adjacent lattice cells play an important role.For this reason, the contrast between the three centralpeaks and minima separating them is rather low in Fig.2, and the total number of particles in the compoundsoliton does not coincide with the sum of the numbers ofits fundamental constituents. Extending this procedure,one can readily build a sequence of compound GSs withan increasing number of peaks.Panels (a) and (b) in Fig. 2 display the 3D density ofthe three-peak GSs as, respectively, an isosurface (takenat 5% of the maximum density) and a color map in thecross-section plane drawn through the z axis. From theseresults, which have been obtained by solving the full 3Dequation (17), we have also derived the axial linear den-sity n ( z ), defined as per Eq. (5). This density is shownin Fig. 2(c) by open circles, along with the correspond-ing results obtained from the effective 1D equations (18)and (19). The OL potential is also displayed by the dot-ted line (in arbitrary units). It is seen that the effective1D equation (18) (solid red lines) reproduces the 3D re-sults very accurately. The 1D GPE (19) (dashed bluelines) gives a slight discrepancy ( ≃ N, e µ ) plane, which sits nearthe top edge of the first band gap. It contains 28 particlesand has chemical potential e µ = 2 .
42, which correspondsto µ = 1 . ~ ω ⊥ . This quantity is again very similarto ~ ω ⊥ , so that one expects the 1D GPE to be a goodapproximation in this case. Figure 3 shows the 3D den-sity and the axial linear density an ( z ) of this GS. InFig. 3(c) one can see that the 1D GPE (19) (the dashedblue line) reproduces the 3D results obtained from thefull GPE (17) (shown by open circles) within a 5% devi-ation. The effective 1D equation (18) (whose results arerepresented by the solid red line) is again more accurateand reproduces the 3D results with an error < E R / ~ ω ⊥ ≪
1) and shallow OL ( V /E R ≤ an ≪ E R , an ~ ω ⊥ ≪ ~ ω ⊥ ), hence nohigher transverse modes are significantly excited. There-fore, the situation considered in this subsection may becategorized as the quasi-1D mean-field regime, in whichthe 1D GPE (19) accurately describes the matter-waveGSs, provided that condition N ≫
2. Deep optical lattice: V /E R = 20 In terms of the underlying OL, this is the case of thetight-binding regime, V ≫ E R . The nonlinear tight-binding regime is realized when the potential depth V ismuch larger than both the linear and nonlinear (mean-field) energies (i.e., V ≫ E R , an ~ ω ⊥ ). Under thesecircumstances, the condensate density is highly localizednear potential minima, making the overlap between den-sities associated with different wells negligible. Point Din the first band gap of Fig. 1(c) corresponds to a fun-damental GS with e µ = 6 .
79 and N = 18. Since its peakaxial density is an (0) ≃ .
2, and V /E R = 20 implies V = 2 ~ ω ⊥ , this fundamental GS belongs to the nonlin-ear tight-binding regime. Figure 4 shows a compound GSbuilt of three such fundamental solitons. It correspondsto point E in Fig. 1(c). Its chemical potential, e µ = 6 . FIG. 4: (Color online) Same as Fig. 2 but for the gap solitoncorresponding to point E in Fig. 1(c). is the same as that of the corresponding fundamentalGS, and its number of particles, N = 55, now coincideswith the sum of the number in its constituents (note thatwe always approximate N to the nearest integer). Ascan be seen from Fig. 4, this compound soliton exhibitsthree well-separated identical peaks (BEC droplets), eachone being practically indistinguishable from the above-mentioned fundamental GS. Figure 4(c) shows that theresults obtained from the effective 1D equation (18) (thesolid red line) agree very well with those produced by thefull 3D GPE (17) (open circles). In particular, the 1Dequation yields the particle number N = 56, which is veryclose to N = 55, as given by the 3D solution. Since N represents a measurement of the norm of the wave func-tion, it is clear that the error in the number of particlesreflects the error in the corresponding wave functions.The 1D GPE (19), corresponding to the dashed blue linein Fig. 4(c), gives N = 49, which implies an error ∼ e µ = 11 .
5, which implies µ = 2 . ~ ω ⊥ . At this valueof µ there is a significant probability of the excitationof higher levels in the radial-confinement potential,which is corroborated by the fact that the peak axialdensity of this GS is an (0) ≃ .
7. Because inequality an ≪ FIG. 5: (Color online) Same as Fig. 2 but for the gap solitoncorresponding to point G in Fig. 1(c). the ordinary 1D GPE (19) to be valid. In fact, ityields the number of particles N = 53, which impliesan error greater than 25%. Nevertheless, the effectivenonpolynomial 1D equation (18) remains valid and quiteaccurate. It produces N = 73, which corresponds to amaximum error of 2 . e µ = 17 . N = 182. Figure 5shows the 3D density and axial linear density an ( z ) ofthis BEC droplet. In this case, an (0) ≃ .
3, indicatinga large contribution of excited radial modes. As a conse-quence, the ordinary 1D GPE (19) [the dashed blue linein Fig. 5(c)] fails to reproduce the 3D results (open cir-cles). It predicts N = 119, which means the error exceeds34%. As seen from Fig. 5(c), the effective 1D equation(18) (the solid red line) remains accurate enough in thiscase, too. In particular, it yields N = 187, the respectiveerror being 2 . E R / ~ ω ⊥ ≪ FIG. 6: (Color online) (a) Band-gap structure of a noninter-acting 3D BEC with E R / ~ ω ⊥ = 1 /
4, as a function of thedimensionless lattice depth s ≡ V /E R . The representativecases s = 2 and 20, indicated by vertical dashed lines, areconsidered in more detail in (b) and (c), respectively, whichshow the dimensionless chemical potential e µ as a function ofthe quasimomentum q in the first Brillouin zone (in units of π/d ). The right panels display the location of the Rb gapsolitons considered in this work (points H–M). most favorable case for the applicability of the 1D ap-proximation, the usual 1D GPE (19) with the deep OLpotential cannot be reliably used beyond the first thirdof the first band gap. Since the validity of the meanfield treatment requires N ≫
1, the usual GPE cannotbe used close to the bottom edge of the band gap, ei-ther. It is thus clear that the range of validity of thisstandard 1D equation is limited. This conclusion be-comes even more severe if one takes into account thatthe tight-radial-confinement regime is not easy to real-ize using typical experimental parameters. Our simula-tions indicate, however, that the effective nonpolynomial1D equation (18) properly describes such essentially 3Dsituations, providing for an accurate description of thefundamental GSs in the entire band gaps.
B. Intermediate radial confinement: E R / ~ ω ⊥ ≃ / This regime is of particular interest because it corre-sponds to typical experimental parameters. It can berealized, for instance, with a Rb condensate in an OLof period d = 1 . µ m, radially confined by a harmonicpotential with ω ⊥ / π = 960 Hz. The particle numbers N of the GSs considered below are given for these physi-cal parameters, and, they can be converted into values of N for other situations by means of Eq. (22) (now, with a/a ⊥ = 0 . ~ ω ⊥ /E R . Since this quantity issmaller in the present case than it was before, the ef-fect of these contributions in the region of interest is nowstronger. While it is obvious that the 1D equations can-not account for the 3D effects originating from the shiftedbands, we will demonstrate that these equations may stillproduce an accurate description of common GSs (i.e.,those originating from the m = n r = 0 energy bands).Figures 6(b) and 6(c) display the band-gap structureas a function of the quasimomentum for lattice depths V /E R = 2 and 20, respectively. Thin black lines in thesefigures represent 3D results that cannot be reproducedby the 1D models, and points H–M in ( N, e µ ) plane markexamples of GSs that will be considered below.
1. Shallow optical lattice: V /E R = 2 Point H in Fig. 6(b) corresponds to a fundamental GSwith e µ = 2 .
19 and N = 50, which is similar to the one inFig. 3, but with the peak axial density an (0) ≃ .
2. Inthis case, the effective 1D equation (18) predicts N = 51,thus corresponding to a 2% error, while the 1D GPE(19) predicts N = 45 (an 11% error). Symmetric combi-nations of five such fundamental GSs generate the com-pound GS shown in Fig. 7, which contains 258 particlesand corresponds to point I in Fig. 6(b). The 3D den-sity of this soliton exhibits five weakly separated peaks,as shown in Figs. 7(a) and 7(b) by means of the isosur-face and the color map in the cross-section plane drawnthrough the z axis. As seen in Fig. 7(c), in this casethe effective 1D equation (18) yields an error of 1 . N = 262) while theuse the ordinary 1D GPE (19) generates an error of 12%( N = 226), showing that, already for this GS, the 3Deffects play an important role.Point J in Fig. 6(b) corresponds to a fundamental GScontaining 75 particles and located near the top edge ofthe first band gap. It is qualitatively similar to the oneshown in Fig. 3, but having e µ = 2 .
42 and an (0) ≃ . ≃
1% ( N = 76) and 12% ( N = 66), respectively. FIG. 7: (Color online) Atom density of the gap soliton corre-sponding to point I in Fig. 6(b), displayed as (a) an isosurfacetaken at 5% of the maximum density and (b) as a color mapalong a cutting plane containing the z axis. (c) Dimension-less axial density an obtained from the 3D wave function,as prescribed by Eq. (5) (open circles) along with the cor-responding predictions from the nonpolynomial 1D equation(18) (solid red line) and the ordinary 1D GPE (19) (dashedblue line).
2. Deep optical lattice: V /E R = 20 Points K, L, and M in Fig. 6(c) correspond to funda-mental GSs in the case of the tight-binding underlyinglinear structure. The first soliton, which contains 19 par-ticles, with e µ = 5 .
33 and axial density an (0) ≃ .
23, alsocorresponds to the nonlinear tight-binding regime and isthus highly localized at a single lattice site. This followsfrom the fact that, for V /E R = 20 and E R / ~ ω ⊥ = 1 / V / ~ ω ⊥ = 5, which is much greater than theabove-mentioned value of an (0). For this GS, locatednear the bottom edge of the first band gap, the ordinary1D GPE (19) yields an error > .
5% deviation.Point L in Fig. 6(c) indicates the position of a funda-mental GS near the top edge of the first band gap, with e µ = 11 . N = 243, and an (0) ≃ .
3. Since the inequal-ity an ≪ N = 133). Note0 FIG. 8: (Color online) Same as Fig. 7 but for the gap solitoncorresponding to point M in Fig. 6(c). that, even though point L is relatively close to a 3D en-ergy band [the thin black line in Fig. 6(c)], which cannotbe reproduced by the effective 1D equation (18), it can,however, reproduce this GS within a 4% error ( N = 253).The same is true for the fundamental GS displayed inFig. 8, which corresponds to point M in Fig. 6(c). Itrepresents a BEC droplet containing 696 particles, with e µ = 17 . an (0) ≃ .
4. Sincewe have V / ~ ω ⊥ < an (0) in the present case, the sys-tem is no longer in the nonlinear tight-binding regime.This is seen in Fig. 8, where the density is not stronglylocalized around the minimum of the potential well. Asseen in Fig. 8(c), the nonpolynomial 1D equation (18)yields results (the solid red line) that agree with the 3Dpicture produced by the full GPE (17) (open circles),within a 3% deviation [Eq. (18) yields N = 719 in thiscase]. On the contrary, the 1D GPE (19), which givesresults (the dashed blue line) with an error >
55% (ityields N = 297), clearly is not valid in this case.These findings demonstrate that, in the physicallyrelevant regime of the intermediate radial confinement( E R / ~ ω ⊥ & / an ( z ) and theparticle content N . For E R / ~ ω ⊥ = 1 / V /E R = 2, the fundamental GSs located in thefirst band gap, as predicted by the 1D GPE equation, fea-ture the error ≃ FIG. 9: (Color online) (a) Band-gap structure of a nonin-teracting 3D BEC with E R / ~ ω ⊥ = 1, as a function of thedimensionless lattice depth s ≡ V /E R . The representativecases s = 2 and 20, indicated by vertical dashed lines, areconsidered in more detail in (b) and (c), respectively, whichshow the dimensionless chemical potential e µ as a function ofthe quasimomentum q in the first Brillouin zone (in units of π/d ). The right panels display the location of the Rb gapsolitons considered in this work (points P–R).
GSs. For the deep OL ( V /E R = 20), the 1D GPE isvalid only in a small region near the bottom edge of thefirst band gap. In general, this equation is applicable onlywhere both conditions an ≪ N ≫ E R / ~ ω ⊥ increases, the relative strength ofthe radial confinement decreases, and the range of valid-ity of the 1D GPE becomes more and more narrow. Fromthe experimental viewpoint, the increase of E R / ~ ω ⊥ maybe relevant because, in this way, one can easily increasethe number of particles in the fundamental GSs. How-ever, the increase of E R / ~ ω ⊥ also implies a decrease inthe relative size of the radial-excitation energy quantum,which can compromise the stability of the solitons be-cause they can decay by exciting higher radial modes,even for condensates with a relatively small number ofparticles. We briefly consider this case below. The sta-bility properties of the GSs will be analyzed in Sec. IV.1 C. Weak radial confinement: E R / ~ ω ⊥ ≥ In this regime, the typical energy scale in the under-lying linear problem is sufficiently large to easily excitehigher transverse modes. As a representative example,we consider the case of E R / ~ ω ⊥ = 1, which can be re-alized in a Rb condensate in an OL with d = 1 . µ m and radial-confinement strength ω ⊥ / π = 240 Hz.The particle contents of the GSs considered below corre-spond to such condensates. As before, the correspond-ing values of N can be converted into those correspond-ing to other situations by means of Eq. (22) (this time, a/a ⊥ = 7 . × − ).Figure 9(a) displays the band-gap structure of the un-derlying 3D linear problem, to be compared with Fig.1(a) which shows the band-gap structure obtained fromthe linearization of 1D equations (18) or (19). As before,the latter equations cannot account for the 3D contribu-tions from the excited radial modes, which account forthe series of shifted bands separated by gaps of value2 ~ ω ⊥ /E R in Fig. 9(a). Figures 9(b) and 9(c) display thelinear band-gap structure as a function of the quasimo-mentum for V /E R = 2 and 20, respectively.Point P in Fig. 9(b) corresponds to a fundamental GSwith e µ = 1 .
59 and N = 74. Its 3D density plot, shownin Fig. 10, demonstrates that this fundamental soliton isspread over several lattice sites, which is a consequenceof the fact that its chemical potential is very close to thefirst linear energy band, where only extended solutionsof the stationary GPE exist. Figure 10(c) shows thatthe results obtained from the effective 1D equation (18)(the solid red line) coincide with the 3D results (opencircles) within 1% (it predicts N = 75). However, the1D GPE (19) (the dashed blue line) gives rise to an error ≃
10% (it predicts N = 66). We thus conclude that, inthe weak-radial-confinement regime ( E R / ~ ω ⊥ ≥ e µ = 2 . N = 400, and an (0) ≃ . N = 265), whilethe nonpolynomial 1D equation (18) limits the error to3 .
5% ( N = 414). Point R in Fig. 9(c) is an example of afundamental GS in a deep OL. It corresponds to a disk-shaped BEC droplet trapped in a single lattice cell [seeFig. 11(b)], which contains 2323 Rb atoms, and has e µ = 11 .
5. Its axial linear density is qualitatively similarto that shown in Fig. 8(c), but with a maximum value an (0) ≃
23. In this case, the 1D GPE predicts theparticle content N = 534, while the effective 1D equa-tion (18) yields N = 2437, which corresponds to an error ≃
5% in the norm of the wave function. Thus, the 1Dnonpolynomial equation provides for a sufficiently gooddescription of the stationary GSs even in the highly non-
FIG. 10: (Color online) Atom density of the gap soliton corre-sponding to point P in Fig. 9(b), displayed as (a) an isosurfacetaken at 5% of the maximum density and (b) as a color mapalong a cutting plane containing the z axis. (c) Dimension-less axial density an obtained from the 3D wave function,as prescribed by Eq. (5) (open circles) along with the cor-responding predictions from the nonpolynomial 1D equation(18) (solid red line) and the ordinary 1D GPE (19) (dashedblue line). linear ( an ≫
1) weak-radial-confinement regime.
IV. STABILITY ANALYSIS
We have investigated the stability of the GS solutionsby monitoring their long-time behavior after the applica-tion of a random perturbation [13, 29]. Specifically, wehave perturbed the corresponding wave functions with asmall-amplitude ( ∼ FIG. 11: (Color online) Evolution in time, after the applica-tion of a 2% white-noise perturbation, of the fundamental gapsolitons corresponding (a) to point P in Fig. 9(b) and (b) topoint R in Fig. 9(c). categorized as stable.We have found that both the full 3D GPE (1) andthe effective nonpolynomial 1D equation (14) lead to thesame conclusions regarding the stability of the GSs. How-ever, once an unstable soliton begins to decay, one cannotexpect, in general, the latter equation to reproduce thedynamical behavior correctly. The reason for this prob-lem is that, as already mentioned in Sec. II, in the pres-ence of the OL potential the above-mentioned adiabaticcondition (i) is hard to fulfill in time-dependent settings.Of course, the same is true for the 1D equation (3), whichis a limit form of Eq. (14) and, consequently, has a muchnarrower range of applicability. The results presented be-low have been obtained using the full 3D GPE (1). Foreach GS we have used at least two different basis setsand time steps to check that the results do not dependon details of the numerical procedure.Our simulations demonstrate that, in general, (1 , , t = 1 s. On the contrary,GSs in a shallow OL ( V /E R = 2), located close to thetop edge of a band gap, such as those corresponding topoints C, H, J, and Q in Figs. 1, 6, and 9, turn out tobe unstable. This instability manifests itself as a steadydecay of the norm of the GS through emission of radia- FIG. 12: (Color online) (a) Evolution in time, after the ap-plication of a 2% white-noise perturbation, of the compoundgap soliton corresponding to point B in Fig. 1(b). Panel (b)displays the long-time behavior, after the application of theperturbation, of a stable three-peak soliton (see text). tion, on a time scale that increases as one moves awayfrom the top edge of the band gap. In this regard, onefinds that the GS corresponding to point H decays muchslower than the other ones. As the lattice depth V /E R increases, in general, GSs become more stable and theinstability region shrinks. In particular, our simulationsindicate that GSs such as those corresponding to pointsF, G, L, M and R in Figs. 1, 6, and 9 (which are funda-mental GSs in the deep OL, with V /E R = 20, locatedclose to the top edge of the band gap) also remain stableup to t = 1 s. This is in good agreement with previousanalyses carried out in the context of deep optical latticesin terms of the ordinary 1D GPE (3) [9, 13].An important result is that (1 , ,
0) fundamental GSsremain stable even in the weak-radial-confinement regime( E R / ~ ω ⊥ ≃ an ( z, t ), whichhas been obtained from the 3D wave function ψ ( r ⊥ , z, t )by integrating out the radial dependence, as per Eq. (5).The left panel in this figure shows the 3D condensate den-3sity at t = 0 (i.e., just after the application of the pertur-bation) as an isosurface taken at 5% of the maximum den-sity, while the right panel represents the density at t = 1s. This GS has e µ = 1 .
59, which implies µ = 2 . ~ ω ⊥ .Note that, despite the fact that this chemical potentialis greater than the quantum ~ ω ⊥ , the GS remains stableup to 1 s. The same is true for the fundamental GS cor-responding to point R in Fig. 9(c), which correspondsto a disk-shaped BEC containing more than 2300 Rbatoms with chemical potential µ = 12 . ~ ω ⊥ ≫ ~ ω ⊥ . AsFig. 11(b) shows, this GS also remains stable.Our simulations also demonstrate that the compoundGS from Fig. 2, corresponding to point B in Fig. 1(b),is unstable. This is seen in Fig. 12(a), which shows howthe soliton decays on a time scale ≃
120 ms after theaddition of a 2% white-noise perturbation. This com-pound GS, which has e µ = 1 .
75 and N = 35, was built inthe shallow OL ( V /E R = 2) as the symmetric combina-tion of three fundamental constituents. We have foundthat, in such shallow lattices, GSs of this type are al-ways unstable. However, it is not difficult to find stablesolitons of this type in somewhat deeper lattices. An ex-ample is shown in Fig. 12(b), which represents a stablethree-peak soliton with e µ = 3 .
13 and N = 59, trappedin the OL with V /E R = 4 and E R / ~ ω ⊥ = 1 /
10. Sim-ilar dynamics is observed for the GS displayed in Fig.7, which corresponds to point I in Fig. 6(b). Thisis a five-peak compound generated in the shallow OL( V /E R = 2) in the regime of the intermediate radial-confinement strength ( E R / ~ ω ⊥ ≃ / V /E R = 4, E R / ~ ω ⊥ = 1 / e µ = 2 .
83, and N = 212. In general, GSs naturally be-come more stable as the lattice depth increases. Fordeep OLs ( V /E R = 20), the stability of compound GSsis essentially identical to that of their fundamental con-stituents. The GS in Fig. 4, which corresponds to pointE in Fig. 1(c), is an example of a stable three-peak soli-ton realized in a deep OL. V. CONCLUSIONS
To appraise the physical relevance of the results re-ported in this work, it is pertinent to recall that themean-field treatment of the GSs (gap solitons), based onthe GPE, is valid if N ≫ a n ≪
1, where a is the scatteringlength of the inter-atomic collisions and n is the atomicdensity. In shallow OLs (optical lattices) such conditionscan be easily met, though the former one imposes a seri-ous limitation on the range of applicability of the usual1D GPE (19), which fails as N increases. In deep OLs,where the tunneling between adjacent cells is stronglysuppressed, the above conditions must be fulfilled at eachsite of the lattice. Under these circumstances, GSs mayhave a rather high local density. A good estimate for the3D peak density n ( ) = N | ψ ( ) | in terms of the peakaxial density, n (0), can be obtained from the followingrelation [17]: n ( ) = n (0) πa ⊥ p an (0) . (25)Substituting the values of an (0) obtained in Section IIIfor different GSs into Eq. (25), one can easily verifythat condition a n ( ) ≪ a n ( ) ≃ − . For the GScorresponding to point R in Fig. 9(c), one finds a n ( ) ≃ × − .Despite the fact that applicability conditions for themean-field treatment can be readily met, it is muchharder to justify the validity of the usual 1D GPE (3),which requires both conditions an ≪ N ≫ Acknowledgments
V.D. acknowledges financial support from Ministeriode Ciencia e Innovaci´on through contract No. FIS2009-07890 (Spain). [1] V. A. Brazhnyi and V. V. Konotop, Mod. Phys. Lett. B , 627 (2004).[2] O. Morsch and M. Oberthaler, Rev. Mod. Phys. , 179(2006).[3] R. Carretero-Gonz´alez, D. J. Frantzeskakis, and P. G. Kevrekidis, Nonlinearity , R139 (2008).[4] Emergent Nonlinear Phenomena in Bose-Einstein Con-densates: Theory and Experiment , edited by P. G.Kevrekidis, D. J. Frantzeskakis, and R. Carretero-Gonz´alez (Springer, Heidelberg, Germany, 2007). [5] Y. V. Kartashov, B. A. Malomed, and L. Torner, Rev.Mod. Phys. , 247 (2011).[6] E. P. Gross, Nuovo Cimento , 454 (1961); J. Math.Phys. , 195 (1963); L. P. Pitaevskii, Zh. Eksp. Teor.Fiz. , 646 (1961) [Sov. Phys. JETP , 451 (1961)].[7] A. Mu˜noz Mateo and V. Delgado, Phys. Rev. Lett. ,180409 (2006).[8] J. A. M. Huhtam¨aki, M. M¨ott¨onen, T. Isoshima, V.Pietil¨a, and S. M. M. Virtanen, Phys. Rev. Lett. ,110406 (2006).[9] P. J. Y. Louis, E. A. Ostrovskaya, C. M. Savage, and Y.S. Kivshar, Phys. Rev. A , 013602 (2003).[10] F. Kh. Abdullaev and M. Salerno, Phys. Rev. A ,033617 (2005).[11] T. Mayteevarunyoo and B. A. Malomed, Phys. Rev. A , 033616 (2006).[12] Y. Zhang and B. Wu, Phys. Rev. Lett. , 093905(2009).[13] Y. Zhang, Z. Liang, and B. Wu, Phys. Rev. A , 063815(2009).[14] A. Mu˜noz Mateo, V. Delgado, and B. A. Malomed, Phys.Rev. A , 053606 (2010).[15] A. E. Muryshev, G. V. Shlyapnikov, W. Ertmer, K. Sen-gstock, and M. Lewenstein, Phys. Rev. Lett. , 110401(2002).[16] L. Salasnich, A. Parola, and L. Reatto, Phys. Rev. A ,043614 (2002); L. Salasnich, Laser Phys. , 198 (2002).[17] A. Mu˜noz Mateo and V. Delgado, Phys. Rev. A ,013617 (2008); , 063610 (2007); Ann. Phys. (NY) ,709 (2009).[18] S. Middelkamp, G. Theocharis, P. G. Kevrekidis, D. J. Frantzeskakis, and P. Schmelcher, Phys. Rev. A ,053618 (2010).[19] G. Theocharis, A. Weller, J. P. Ronzheimer, C. Gross, M.K. Oberthaler, P. G. Kevrekidis, and D. J. Frantzeskakis,Phys. Rev. A , 063604 (2010).[20] C. Wang, P. G. Kevrekidis, T. P. Horikis, D. J.Frantzeskakis, Phys. Lett. A, , 3863 (2010).[21] S. Middelkamp, J. J. Chang, C. Hamner, R. Carretero-Gonz´alez, P. G. Kevrekidis, V. Achilleos, D. J.Frantzeskakis, P. Schmelcher, and P. Engels, Phys. Lett.A , 642 (2011).[22] An independent derivation of this formula was reportedby F. Gerbier, Europhys. Lett. , 771 (2004).[23] L. Salasnich, B.A. Malomed, and F. Toigo, Phys. Rev. A , 063614 (2007).[24] A. Mu˜noz Mateo and V. Delgado, Phys. Rev. A ,065602 (2006).[25] S. Adhikari and B.A. Malomed, Europhys. Lett. ,50003 (2007); Physica D , 1402 (2009).[26] F. H. Mies, E. Tiesinga, and P. S. Julienne, Phys. Rev.A , 022721 (2000).[27] B. Eiermann, Th. Anker, M. Albiez, M. Taglieber, P.Treutlein, K.-P. Marzlin, and M. K. Oberthaler, Phys.Rev. Lett. , 230401 (2004).[28] J. Armijo, T. Jacqmin, K. Kheruntsyan, and I. Bou-choule, Phys. Rev. A , 021605(R) (2011).[29] J. Yang and Z. H. Musslimani, Opt. Lett.28