# Bose-Einstein condensates under a spatially-modulated transverse confinement

###### Abstract

We derive an effective nonpolynomial Schrödinger equation (NPSE) for self-repulsive or attractive BEC in the nearly-1D cigar-shaped trap, with the transverse confining frequency periodically modulated along the axial direction. Besides the usual linear cigar-shaped trap, where the periodic modulation emulates the action of an optical lattice (OL), the model may be also relevant to toroidal traps, where an ordinary OL cannot be created. For either sign of the nonlinearity, extended and localized states are found, in the numerical form (using both the effective NPSE and the full 3D Gross-Pitaevskii equation) and by means of the variational approximation (VA). The latter is applied to construct ground-state solitons and predict the collapse threshold in the case of self-attraction. It is shown that numerical solutions provided by the one-dimensional NPSE are always very close to full 3D solutions, and the VA yields quite reasonable results too. The transition from delocalized states to gap solitons, in the first finite bandgap of the linear spectrum, is examined in detail, for the repulsive and attractive nonlinearities alike.

###### pacs:

03.75.Ss,03.75.Hh,64.75.+g## I Introduction

It is well known that Bose-Einstein condensates (BECs) with weak attractive interactions between atoms can form stable solitons in “cigar-shaped” (nearly one-dimensional, 1D) traps. In these traps, the gas is strongly bound in the transverse plane, while being loosely confined along the longitudinal axis (). Using this configuration, stable bright solitons were created in the gas of Li atoms Strecker02 (); Khaykovich02 (). In the condensate of Rb atoms trapped in a similar configuration, stronger attraction between atoms leads to collapse and emergence of nearly 3D solitons in a post-collapse state Cornish ().

The strongly elongated (cigar-shaped) settings are described by effectively 1D equations which were derived, by dint of various approximations, from the full 3D Gross-Pitaevskii equation (GPE) PerezGarcia98 ()-denicola (). The derivation assumes an ansatz factorizing the 3D wave function into a product of a transverse mode and an axial (one-dimensional) wave function. As shown in Refs. sala1 (); sala2 (), the substitution of the factorized ansatz in the underlying cubic GPE leads, in the general situation, to a nonpolynomial Schrödinger equation (NPSE) for the axial wave function (or a system of coupled NPSEs for a two-component BEC we-first ()). In the case of weak nonlinearity, the NPSE can be expanded, which leads to a simplified 1D equation with a combination of cubic and quintic terms, the latter one being always attractive; if the cubic term is attractive too, the cubic-quintic equation gives rise to a family of stable solitons available in an exact analytical form Lev (). Unlike the ordinary 1D nonlinear Schrödinger (NLS) equation with the attractive cubic term, the NPSE, as well as its cubic-quintic truncation, may give rise to collapse, which reflects the occurrence of the collapse in the underlying cubic GPE in three dimensions sala1 (); sala2 (); Lev (). Despite the possibility of the collapse, solitons are stable in these models.

A relevant problem is to derive the NPSE for the cigar-shaped trap equipped with a periodic potential, which is created in the experiment as an optical lattice (OL), i.e., an interferences pattern, by a pair of counterpropagating laser beams illuminating the trap in the axial direction (the BEC dynamics in periodic potentials was recently reviewed in Ref. Markus ()). A 1D equation of the NPSE type including the OL potential was recently derived in Ref. we (). Using that equation, the influence of the periodic potential on the collapse threshold for axially localized states in the quasi-1D trap was investigated, both for single-peak solitons found in the semi-infinite gap of the periodic potential, and multi-peaked solitons found in finite bandgaps. The results were compared to direct numerical solutions of the full 3D GPE, showing good agreement.

In the experiment, the transverse potential which confines the atomic gas to the cigar-shaped configuration may be axially nonuniform, which corresponds to the corresponding trapping frequency being a function of the axial coordinate, (generally speaking, it may also depend on time). As proposed in Ref. denicola (), a specially designed nonuniformity (axial modulation) of may be used as an alternative tool for the control of dynamics of nearly 1D solitons, inducing an effective potential for them. In that work, the effective potential for a soliton with norm (scaled number of atoms) was found in the limit of weak nonlinearity and long-scale axial modulation,

(1) |

where is an effective nonlinearity coefficient, see Eq. (3) below (we will use normalization tantamount to setting ), is the coordinate of the soliton’s center, and integer is possible intrinsic vorticity of the quasi-1D soliton, which is defined below in Eq. (8) (ordinary solitons correspond to ). Expression (1) was derived by means of a variational approximation (VA) applied directly to the full 3D energy functional, cf. Eq. (9) below.

Of special interest is the case of periodic modulation of the transverse trapping frequency, with wavenumber , modulation depth , and amplitude :

(2) |

As suggested by Eq. (1), the periodic modulation may replace the OL potential. In particular, this setting may be especially relevant to a situation when the quasi-1D magnetic trap is not rectilinear, but rather circular (toroidal), which was realized in the experiment torus (). Indeed, the OL cannot be created in such a setting, but the induction of an effective periodic potential by means of the modulation of the transverse trapping frequency is quite feasible.

The objective of the present work is to derive an effective NPSE for the quasi-1D trap subject to the periodic modulation as per Eq. (2), and to investigate various self-trapped states in this geometry, both delocalized ones and solitons. The paper is organized as follows. The NPSE is derived, in a general form, in Sec. II. Then, the model with the repulsive nonlinearity is considered in Sec. III. At first, the analysis includes an additional parabolic trapping potential acting in the axial direction. Then, this potential is removed, and we analyze solutions demonstrating a transition from delocalized solutions to a gap soliton. The solutions predicted by the effective 1D equation are compared to their counterparts found from the full 3D GPE (bandgap structures generated by linearized versions of both equations are compared too). Section IV is dealing with the case of the attractive nonlinearity. The corresponding ground-state solitons are found by means of the VA and in a numerical form, which are also compared with results produced by the full 3D equations. Because the attractive nonlinearity may give rise to collapse, the collapse threshold is considered in detail too, by means of both the VA and numerical methods. In addition to the ground-state solitons, which reside in the semi-infinite gap induced by periodic modulation (2), we also explore solutions with the chemical potential belonging to the first finite gap, and demonstrate the transition from delocalized states to gap solitons in that case. The paper is concluded by Sec. V.

## Ii The nonpolynomial Schrödinger equation

Static BEC configurations can be derived from the normalized functional which determines the energy-per-atom in terms of order parameter (single-atom wave function) and includes a generic external axial potential :

(3) |

Here is the adimensional strength of the interaction between atoms, with the inter-atomic scattering length and the number of atoms in the condensate. In Eq. (3) lengths are measured in units of the characteristic transverse trapping (harmonic-oscillator) length, , where is the atomic mass, the energy is in units of , modulation wave number is in units of , and the wave function is subject to the ordinary normalization,

(4) |

The chemical potential corresponding to Eqs. (3) and (4) is

(5) |

Due to the above normalizations, and are not explicitly present in expressions (3) and (5).

An accurate investigation of the present setting can be performed by using the approach which was first developed for the GPE with the unmodulated transverse trap; after averaging the full 3D equation in the transverse plane, it leads to an effective one-dimensional NPSE sala1 (). The derivation of the NPSE starts with the factorization of the 3D wave function into a product of an arbitrary complex axial wave function, , and the ordinary transverse Gaussian ansatz with transverse width ,

(6) |

being subject to normalization condition

(7) |

Configurations including the above-mentioned intrinsic vorticity, characterized by positive integer , correspond to the following generalization of Eq. (6):

(8) |

where and are the polar coordinates in transverse plane . Nearly-1D solitons with intrinsic vorticity were studied, by means of methods similar to those used in the present work (although without deriving a closed-form NPSE for that case) in Ref. Luca (); for a limit case of a strongly elongated 3D trap, similar localized states were also studied in Ref. Dumitru (). In this work, we focus on the ordinary solitons, with .

Inserting expression (6) in Eq. (3), one obtains, neglecting the -derivative of , the following effective (1D) energy functional:

(9) |

The minimization of this functional with respect to , taking normalization (7) into regard, leads to the stationary NPSE,

(10) |

where the chemical potential appears as the Lagrange multiplier generated by the normalization condition.

An equation for the transverse width is obtained by the minimization of functional (9) with respect to :

(11) |

The substitution of this expression in Eq. (10) leads to a closed-form stationary NPSE for , although in a rather cumbersome form.

In the limit of zero scattering length, , Eqs. (10) and (11) reduce to the linear Schrödinger equation,

(12) |

with the effective axial potential,

(13) |

(note that, even in this limit, the full 3D GPE is not separable, due to the modulation imposed on the transverse trapping potential). If the nonlinear term is small, the NPSE may be approximated by the cubic NLS equation with the same effective axial potential and, in addition to that, with a periodically-modulated nonlinearity coefficient:

(14) |

Essentially the same result as given by Eq. (14) was obtained, in the simplest approximation, in Ref. denicola (), see Eqs. (1) and (2).

## Iii The model with repulsive nonlinearity

### iii.1 One-dimensional solutions

In the case of the repulsive inter-atomic interaction, i.e., positive and , the application of the Thomas-Fermi (TF) approximation, which neglects the second derivative, to Eqs. (10) and (11) yields an analytical expression for the normalized atomic density, :

(15) |

where an effective local chemical potential is

(16) |

Equations (15) and (16) generalize the result obtained in Ref. sala3 () for the repulsive BEC under the unmodulated transverse confinement (). In the limit case of strong nonlinearity, , Eq. (15) reduces to and in the opposite limit of , which implies , it takes the form .

To obtain accurate results, we solved the time-dependent variety of Eq. (10) (with replaced by and, accordingly, replaced by ), combined with Eq. (11) (without any change in the latter equation) numerically, by means of the finite-difference Crank-Nicholson method in imaginary time, following the approach elaborated in Ref. sala-numerics (). The explicit axial potential was chosen in the usual form, , where , and is the axial-confinement frequency. In this way, the profiles for , displayed by dashed lines in Fig. 1, have been obtained for different values of modulation depth and fixed nonlinearity strength, . The so obtained NPSE profiles are compared to those produced by the TF approximation (dotted lines in Fig. 1), see Eq. (15). In the absence of transverse modulation, (the upper panel in Fig. 1), the axial density profile is well approximated by the TF formula. At (the central and lower panels in Fig. 1), the TF approximation much overestimates the density contrast between points of local minima and maxima of the total axial potential.

### iii.2 Three-dimensional solutions

It is necessary to compare results produced by the NPSE to their counterparts found from a direct numerical solution of the full GPE in three dimensions. The latter equation is obtained by the minimization of energy functional (3). Figure 1 shows that the density profiles generated by the NPSE (dashed lines) are very close to ones obtained from the 3D equation (solid lines), unless is very large. It is noteworthy that the NPSE gives very accurate results for a model with non-separable potential.

One may expect that the effective axial periodic potential induced by the transverse modulation (without the inclusion of the axial parabolic trap) may support quasi-periodic Bloch states and gap solitons in the self-repulsive BEC, as in the case of the ordinary (direct) periodic potential GSprediction (); Markus (). To consider this possibility, a self-consistent numerical method was used, with periodic boundary conditions. We employed a spatial grid of points, covering the interval of , which corresponds to periods of the external modulation. To test the numerical scheme, we have checked that the lowest-energy state in the semi-infinite bandgap (i.e., the ground state of the system), produced by this method, is identical to that found above by the integration of the time-dependent NPSE in imaginary time.

In the upper panel of Fig. 2 we plot, as a function of nonlinearity strength , the first eigenvalues , as found by means of the above-mentioned numerical method from Eqs. (10) and (11) with and periodic boundary conditions. The lowest eigenstates belong to the first nonlinear band induced by the periodic modulation in the full NPSE, and the other states, which are well separated by a bandgap, form a second nonlinear band. For some of the intraband states, the numerical method does not converge to a single configuration, but rather oscillates between two or three configurations with very close eigenvalues. It is also possible that the nonlinear model may give rise to intraband states which have no counterparts in the linear limit. This challenging issue needs special treatment, and will be considered elsewhere. As concerns the nonlinear eigenstates for which our presently employed numerical method leads to convergence (including the gap soliton, see below), we have verified, by direct simulations of the time-dependent NPSE in real time, that they all are stable.

As a typical example, in the lower panel of Fig. 2 we plot the density profile of the 32th state (it has number in the set of numerically found states) for and three different values of . One can check that, in the linear approximation, this state lies at the top of the first Bloch band. With the increase of the nonlinearity strength , the energy of the 32th state grows, and, in doing so, it enters the first finite bandgap (as defined in the linear approximation) from below. Figure 2 demonstrates that, for , this state is still fully delocalized, being thus similar to a Bloch wave, while for it becomes localized, with a width much smaller than the length of the periodic box. Clearly, this solution may be identified as a gap soliton. At , the gap soliton compresses itself (see Fig. 2) into a still narrower state.

It is also necessary to check that the bandgap structure in the NPSE with the periodically modulated trapping potential in the linear limit () is not, by itself, an artifact of the approximation (reduction of the 3D equation to the 1D limit), but a true feature of the full 3D model. To this end, we solved the eigenvalue problem for the full linear GPE in three dimensions,

(17) |

where acts on and , and, as above, . The transverse part of solutions to Eq. (17) may be taken as an eigenstate of the corresponding 2D harmonic oscillator, with its quantum numbers and . Defining as the smallest reciprocal lattice vector and writing ( is an integer), we find that the solution can be written as

(18) |

with coefficients satisfying a linear eigenvalue problem,

(19) |

where is a wavenumber in the first Brillouin zone, and the matrix is

(20) |

We solved Eq. (19) by truncating the sum in expression (18) to . At first, we also truncated the summation to (i.e., only the contribution from the ground state of the transverse potential was taken into account), which yielded the result shown in the left panel of Fig. 3. Then we found more accurate solutions, by extending the truncated summation to (i.e., including the contribution from excited transverse states). In that case, Eq. (19) produces the dispersion relation shown in the right panel of Fig. 3.

These results were compared to the dispersion law as found directly from the NPSE, Eq. (10). We observe that the two dispersion laws are identical if the 3D analysis is confined to , while they are different if the effect of states with are considered. Nevertheless, Fig. 3 shows that the first band, the first gap, and half of the second band do not alter essentially, if the NPSE is replaced by the full 3D equation, which justifies the use of the NPSE approximation for values of the chemical potential up to the first half of the second band.

## Iv The model with the attractive nonlinearity

### iv.1 Numerical results

The model with the attractive inter-atomic interactions, i.e. , may be expected to generate bright solitons. We analyzed this possibility through the numerical solution of the NPSE equation, in the absence of the direct axial potential, . In Fig. 4, we plot axial density profiles of the thus found solitons for different values of modulation depth and fixed nonlinearity strength, .

For , the soliton’s profile, with a single maximum, may be well fitted by , which is an asymptotically exact solution (the usual NLS soliton) in the above-mentioned weakly nonlinear limit corresponding to , provided that varies on the entire real axis. On the other hand, if belongs to a finite interval, , with periodic boundary conditions, , it is known sala4 () that the NPSE with yields a spatially uniform ground state, , for sufficiently weak nonlinearity, ; the ground state develops a spatial structure at . As shown in Fig. 4, for nonzero the soliton profile features several local maxima and minima due to the action of the effective periodic potential. Thus, under such conditions, the Bose condensate self-traps into a multi-peaked soliton, which occupies several cells of the periodic modulation.

We also compared the results yielded by the NPSE with those found from the direct numerical solution of the full GPE in three dimensions. Figure 4 shows that the density profiles generated by the NPSE (dashed lines) practically coincide with the ones obtained from the 3D GPE (solid lines), unless is very large.

It is also interesting to analyze the behavior of the density profile of the bright soliton as a function of the wavenumber . In Fig. 5 we plot the density profile at and , for four values of . The figure shows that, at , the profile is strongly localized within one single site of the effective periodic axial potential: accordingly, we call the corresponding solution a single-site soliton (cf. Ref. we (), where the NPSE for the model with attraction and an explicit periodic axial potential was considered). At , delocalization of the bright soliton, which occupies more than one site, is observed. We call it a multi-site soliton (the distinction between the strongly and weakly localized solitons is also observed in the 1D cubic NLS equation with the self-attractive nonlinearity NLS ()).

### iv.2 Variational vs numerical results

From the numerical results presented above we see that the profiles of the localized solutions of Eq. (10) both in the single-site (strong attraction, large values of ) and in the multi-site soliton cases but with a weak transverse modulation (small values of ) are smooth. This observation suggests that one could achieve some analytical insight into the model with attraction by using a VA with the Gaussian ansatz (a review of the VA can be found in Ref. progress ()),

(21) |

where and are variational parameters accounting for the transverse and axial width of the BEC. Inserting this ansatz in Eqs. (3) and (5), with (and ), we obtain the respective expressions for the energy-per-atom functional and chemical potential,

(22) |

(23) |

Next, we minimize the energy, demanding , which yields the variational equations,

(24) |

(25) |

which can be easily solved numerically method (). The solutions provide for a minimum of the energy only if the curvature of the energy surface, , is positive, i.e.,

(26) |

As concerns the dynamical stability of the solitons against small perturbations, it may be, first of all, estimated by means of the VK criterion VK (). According to it, a necessary stability condition is (in the present notation), if the soliton family is described by dependence (note that the VK criterion does not apply to gap solitons in the model with the repulsive nonlinearity, therefore it was not used in the previous section).

In Fig. 6, we plot axial length and transverse width versus interaction strength for and four different fixed values of modulation parameter . The figure shows that the soliton in the self-attractive BEC exists up to a critical value of the nonlinearity strength, . At , the 3D collapse of the nearly-1D soliton occurs, which is a well-known result in the case of PerezGarcia98 (); sala1 (); sala2 (). Note that, for , the VA predicts , which is somewhat higher than obtained from the numerical solution of the full GPE in three dimensions gammal (); sala2 ().

The axial length of the soliton, , diverges as drops to zero, while the transverse width approaches , actually becoming equal to the above-mentioned harmonic-oscillator length, . On the other hand, as approaches , both and remain finite and smaller than .

New results are presented in Fig. 6 for (recall previous works were only dealing with the case of ). At small , the figure shows only a slight distortion of the curves. A qualitative change is observed at , when there appear two stable branches for both and , the curves for displaying a clear gap. The lower branches of the and dependences exists only in a finite interval, which we denote as

(27) |

while the upper branches extend up to , in interval , with . Physically, the lower branches (with smaller values of axial length ) correspond to single-site solutions, where the self-attractive BEC is strongly localized – essentially, within a single cell of the modulation structure. The upper branches of and correspond instead to the weakly localized solutions (with a larger axial length), which occupy several cells (sites).

Analysis of expressions (22) and (23) demonstrates that, in interval (27), the multi-site and single-site solution may assume the role of ground state (the one corresponding to the lowest energy). However, this analysis also suggests that both families are dynamically stable, as they always meet the VK criterion, . Direct numerical simulations (not shown in detail here) have confirmed this conjecture.

For , the numerical solution of Eqs. (24) and (25) demonstrates that the lower border of existence interval (27) for the single-site soliton, , vanishes, but this is an artifact of the VA, which occurs in other contexts too Salerno (). It is explained by the above-mentioned inadequacy of the Gaussian ansatz in the limit of weak nonlinearity, i.e., for widely spread small-amplitude solitons featuring a multi-peaked shape. In this situation, one may, in principle, apply a more sophisticated ansatz, combining the Gaussian and periodic functions, such as ; however, the generalized ansatz results in a cumbersome algebra Arik (), therefore we do not pursue such an approach here.

In the upper panel of Fig. 7, we plot the density profile, , of the soliton for different values of the self-attraction strength, , and fixed modulation parameters, and . As seen in the figure, the increase of may strongly compress the soliton in the axial direction, making the secondary maxima very small. In this case, the condensate actually self-traps into a single-peak soliton, which occupies only one cell of the modulation structure.

Contrary to the numerical solution of the NPSE which shows a smooth crossover between multi-site and single-site solitons, the VA predicts a discontinuous transition. In the lower panel of Fig. 7, we plot the axial length, , of the ground-state bright soliton as a function of strength , for and . The jump predicted by the variational calculation at , is a consequence of the inadequacy of of ansatz (21), which assumes the simple Gaussian waveform for the axial wave function, to describe multi-peaked states, as was also recently shown in the study of the quasi-1D model with the self-attraction and axial optical lattice we ().

As it is well known, cold Bose atoms with attractive interactions collapse if the interaction strength exceeds a threshold value, . Obviously, as increases towards , the profile becomes narrower and narrower, and therefore close to collapse in the presence of a transverse modulation, its shape should correspond to a single-site soliton, where the VA is appropriate. It is thus interesting to predict the collapse threshold, , as a function of parameters and of the transverse modulation by using our gaussian variational ansatz. In Fig. 8, we display the dependence predicted by the VA at five fixed values of modulation depth . For given , the critical value has its maximum at , which is natural, as Eq. (2) yields, in this case, the smallest constant value of the transverse-trapping frequency. Equation (2) also helps to understand another feature observed in Fig. 8., viz., that slowly diverges at as approaches (for instance, at ). Note also that there exists a modulation wavenumber, , at which attains its minimum, i.e., the collapse has the lowest threshold. Figure 8 shows that this minimum decreases with the increase of , which may be understood too: as mentioned above, the strong potential tends to squeeze the entire condensate into a single cell of the modulation structure, which facilitates the onset of the collapse. On the other hand, at large values of , becomes asymptotically constant, as the interaction of the condensate with the short-period modulation becomes exponentially weak, hence it produces little effect on the collapse threshold.

In Table 1, we display the comparison of the VA-predicted critical value, , to results following from the numerical solution of NPSE (10) (for ). The variational and numerically found critical values are seen to be in qualitative agreement. Note that, at , the critical value for the NPSE is , which is an obvious extension of found before at sala1 (); sala2 ().

0.00 | 2.12 | 2.76 | 2.16 | 0.86 |

0.10 | 1.99 | 2.52 | 1.20 | 0.72 |

0.25 | 1.74 | 2.17 | 0.91 | 0.65 |

0.50 | 1.55 | 1.85 | 0.75 | 0.56 |

1.00 | 1.41 | 1.57 | 0.76 | 0.57 |

1.50 | 1.34 | 1.47 | 1.19 | 0.79 |

2.00 | 1.32 | 1.48 | 1.25 | 0.83 |

2.50 | 1.45 | 1.54 | 0.99 | 0.75 |

3.00 | 1.44 | 1.55 | 0.93 | 0.73 |

3.50 | 1.42 | 1.55 | 0.92 | 0.73 |

Table 1. Properties of the bright soliton in the model with attraction in a proximity to the collapse, as found from the numerical solution of Eq. (10) at and different values of modulation wavenumber : is the critical value of the nonlinearity coefficient at the collapse point, is the axial width of the soliton, and its transverse width (at ) at the same critical point. For comparison, also included are values of predicted by the variational approximation, .

Table 1 demonstrates that, with the increase of , the critical value, , drops from the largest value, corresponding to , to a minimum at , and then slightly increases with the further increase of (the minimum at seems quite shallow from the side of ). This feature and the related ones can be explained if one notices that corresponds to a constant value of the modulation factor in Eqs. (10) and (11),

(28) |

and, on the other hand, for large (formally, for ), the modulation factor should be replaced by its average,

(29) |

where is the complete elliptic integral of the second kind. Then, Eqs. (10) and (11) with the constant modulation factor (and without the extra potential, ), take the form of

(30) |

(recall ). Further, the coefficient can be eliminated from Eq. (30) by means of an obvious rescaling [which also takes into regard the condition that the normalization of the solution must keep the form of Eq. (7)]:

As a consequence, the critical values of , together with the respective length scale, which correspond to the different constant values of , are related as follows:

(31) |

For (the value for which numerical data are collected in Table 1), Eqs. (28) and (29) yield and , hence the ratios predicted by Eq. (31) take the value . On the other hand, the numerical data from Table 1 yield , and . Thus, approximating at by (recall the change of at is insignificant), one can explain (at least, in a crude approximation) effects caused by the transition from long-scale to short-scale spatial modulation.

### iv.3 Solitons in the first finite bandgap

The above analysis was dealing with solitons (in the attractive model) whose chemical potential belongs to the semi-infinite bandgap in the linear spectrum of Eq. (10). On the other hand, it is known that the cubic self-attractive nonlinearity may also give rise to solitons located in higher-order (finite) bandgaps. Here, we report soliton solutions of the latter type, found from the NPSE by means of the self-consistent numerical method with periodic boundary conditions, which was outlined in the previous section.

In Fig. 9 we plot the first eigenvalues , as found from the numerical solution of the stationary nonlinear NPSE, versus the nonlinearity strength . The first eigenstates form the first band, and the other , which are well separated by a gap, cluster into the second band. Figure 9 shows that the lowest eigenvalue and the 33rd one split off from the first and second bands and move down (till the onset of the collapse) with the increase of , thus giving rise to localized states in the semi-infinite and first finite gaps, respectively. It is worthy to note that the second eigenvalue, which originally belonged to the first band, also splits off from it at larger values of . We have verified that the corresponding nonlinear eigenstate become localized, as one may expect. Qualitatively similar findings were reported in Ref. efremidis (), which was dealing with a numerical solution of the ordinary cubic GPE in one dimension, and, more recently, in Ref. we () which was dealing with the self-attractive BEC in the quasi-1D trap with an axial OL potential.

The density profiles, , of the 33th state which develops into a soliton belonging to the first finite bandgap, are displayed, for and six different values of interaction strength , in Fig. 10. The figure shows that this state is still fully delocalized (being similar to a Bloch wave) at , while at it becomes localized, with the width much smaller than the length of the periodic box, featuring many local maxima and minima (zeros). With further increase of , the gap soliton keeps compressing itself. Comparing Figs. 4 and 10, we conclude that the ground-state solitons, which reside in the semi-infinite gap, are drastically different from the gap solitons, i.e., ones found in the first finite bandgap. First, for the same parameters, the ground-state solitons are localized much stronger, and their local density minima are not zeros, unlike those of the gap solitons. In addition, it is noteworthy that local density maxima of the ground-state solitons correspond to minima of the effective periodic potential, while the local maxima of the gap soliton correspond to maxima of the periodic potential.

Although the gap solitons cannot play the role of ground states, they, as well as the ground-state solitons, are stable in direct simulations of the time-dependent variant of the NPSE equation. Therefore, the gap-soliton states are relevant to the experiment.

## V Conclusions

In this work, we have derived the effective one-dimensional NPSE (nonpolynomial Schrödinger equation) for a cigar-shaped trap whose transverse confining frequency is periodically modulated along the axial direction, thus inducing an effective periodic axial potential. Besides the usual quasi-1D linear geometry, the model may also be relevant as a means of creating an effective periodic potential in toroidal traps. In both cases of the repulsive and attractive nonlinearity, delocalized states and solitons were found, by means of numerical methods (which were applied to both the effective NPSE and the underlying 3D Gross-Pitaevskii equation) and VA (variational approximation; this method was applied to ground-state solitons in the model with attraction, and to the prediction of the collapse threshold in this case). It was found that the numerical solution to the NPSE is always extremely close to the full 3D solutions. The VA yields quite reasonable results too, except for the description of the crossover from single-site to multi-site solitons: numerical results reveal that the crossover is smooth and does not include a jump, contrary to the prediction of the VA. This shortcoming of the VA is explained by the fact that the simple Gaussian ansatz, on which the approximation is based, cannot adequately grasp the transition that alters the shape of the soliton, giving it the multi-peaked structure. The transition from delocalized states to gap solitons was studied in detail (by means of numerical methods) in the first finite bandgap, for both cases of the repulsive and attractive nonlinearities.

The above results, presented in terms of the dimensionless equations, can be easily translated into physical units. For instance, by considering an attractive Bose-Einstein condensate made of Li atoms, with scattering length nm, and choosing the transverse confining frequency as Hz, we have a typical value of the modulation period, m, for (recall it was a typical value for which the results were reported above). In this case the critical number of atoms at the collapse threshold is .

The results may be quite useful to design new experiments in toroidal traps, where an azimuthal periodic potential can be induced, as explained above, by the spatially-modulated transverse confinement.

The work of B.A.M. was partially supported by the Israel Science Foundation through Excellence-Center Grant No. 8006/03, and by German-Israel Foundation through Grant No. 149/2006. L.S. acknowledges GNFM-INdAM for partial support. L.S. and L.R. thank Prof. Alberto Parola for numerous discussions and suggestions.

## References

- (1) K. E. Strecker, G. B. Partridge, A. G. Truscott and R. G. Hulet, Nature 417, 150 (2002); K. E. Strecker, G. B. Partridge, A. G. Truscott, and R. G. Hulet, New J. Phys. 5, 73 (2003).
- (2) L. Khaykovich, F. Schreck, G. Ferrari, T. Bourdel, J. Cubizolles, L. D. Carr, Y. Castin, and C. Salomon, Science 256, 1290 (2002).
- (3) S. L. Cornish, S. T. Thompson and C. E. Wieman, Phys. Rev. Lett. 96, 170401 (2006).
- (4) V. M. Pérez-García, H. Michinel, and H. Herrero, Phys. Rev. A 57, 3837 (1998).
- (5) L. Salasnich, Laser Phys. 12, 198 (2002); L. Salasnich, A. Parola, and L. Reatto, Phys. Rev. A 65, 043614 (2002).
- (6) L. Salasnich, A. Parola, and L. Reatto, Phys. Rev. A 66, 043603 (2002).
- (7) L. Salasnich, A. Parola and L. Reatto, Phys. Rev. A 69, 045601 (2004).
- (8) L. Salasnich and B. A. Malomed, Phys. Rev. A 74, 053610 (2006).
- (9) A. E. Muryshev, G. V. Shlyapnikov, W. Ertmer, K. Sengstock, and M. Lewenstein, Phys. Rev. Lett. 89, 110401 (2002).
- (10) S. De Nicola, B. A. Malomed, and R. Fedele, Phys. Lett. A 360, 164 (2006).
- (11) D. E. Pelinovsky, Y. S. Kivshar, and V. V. Afanasjev, Physica D 116, 121 (1998); L. Khaykovich and B. A. Malomed, Phys. Rev. A 74, 023607 (2006).
- (12) O. Morsch and M. Oberthaler, Rev. Mod. Phys. 78, 179 (2006).
- (13) L. Salasnich, A. Cetoli, B. A. Malomed, F. Toigo, Phys. Rev. A 75, 033622 (2007).
- (14) S. Gupta, K. W. Murch, K. L. Moore, T. P. Purdy, and D. M. Stamper-Kurn, Phys. Rev. Lett. 95, 143201 (2005).
- (15) L. Salasnich, Laser Phys. 14, 291 (2004).
- (16) D. Mihalache, D. Mazilu, B. A. Malomed, F. Lederer, Phys. Rev. A 73, 043615 (2006).
- (17) E. Cerboneschi, R. Mannella, E. Arimondo, and L. Salasnich, Phys. Lett. A 249, 495 (1998); L. Salasnich, A. Parola, and L. Reatto, Phys. Rev. A 64, 023601 (2001).
- (18) G. L. Alfimov, V. V. Konotop, and M. Salerno, Europhys. Lett. 58, 7 (2002); B. B. Baizakov, V. V. Konotop, and M. Salerno, J. Phys. B 35, 51015 (2002); P. J. Y. Louis, E. A. Ostrovskaya, C. M. Savage, and Y. S. Kivshar, Phys. Rev. A 67, 013602 (2003).
- (19) A. Parola, L. Salasnich, R. Rota, and L. Reatto, Phys. Rev. A 72, 063612 (2005); L. Salasnich, A. Parola, and L. Reatto Phys. Rev. A 74, 031603(R) (2006).
- (20) B. A. Malomed, in Progress in Optics, edited by E. Wolf (North Holland, Amsterdam, 2002), Vol. 43, p. 71.
- (21) Equations (24) and (25) are solved by choosing as an independent parameter. Then, defining and , one finds, after some algebraic manipulations, , and .
- (22) M. G. Vakhitov and A. A. Kolokolov, Radiophys. Quantum Electron. 16, 783 (1973).
- (23) A. Gammal, L. Tomio, and T. Frederico, Phys. Rev. A 66, 043619 (2002).
- (24) B. A. Malomed, Z. H. Wang, P. L. Chu, and G. D. Peng, J. Opt. Soc. Am. B 16, 1197 (1999); G. Alfimov and V. V. Konotop, Physica D 146, 307 (2000).
- (25) B. B. Baizakov, B. A. Malomed, and M. Salerno, Europhys. Lett. 63, 642 (2003); Phys. Rev. A 70, 053613 (2004).
- (26) A. Gubeskys, B. A. Malomed, and I. M. Merhasin, Stud. Appl. Math. 115, 255 (2005).
- (27) N. K. Efremidis and D. N. Christodoulides, Phys. Rev. A 67, 063608 (2003).