## Abstract

We investigated the stability properties of breather soliton trains in a three-dimensional Bose–Einstein condensate (BEC) with Feshbach-resonance management of the scattering length. This is done so as to generate both attractive and repulsive interaction. The condensate is confined only by a one-dimensional optical lattice and we consider strong, moderate and weak confinement. By strong confinement we mean a situation in which a quasi two-dimensional soliton is created. Moderate confinement admits a fully three-dimensional soliton. Weak confinement allows individual solitons to interact. Stability properties are investigated by several theoretical methods such as a variational analysis, treatment of motion in effective potential wells, and collapse dynamics. Armed with all the information forthcoming from these methods, we then undertake a numerical calculation. Our theoretical predictions are fully confirmed, perhaps to a higher degree than expected. We compare regions of stability in parameter space obtained from a fully three-dimensional analysis with those from a quasi two-dimensional treatment, when the dynamics in one direction are frozen. We find that in the three-dimensional case the stability region splits into two parts. However, as we tighten the confinement, one of the islands of stability moves toward higher frequencies and the lower frequency region becomes more and more like that for the quasi two-dimensional case. We demonstrate these solutions in direct numerical simulations and, importantly, suggest a way of creating robust three-dimensional solitons in experiments in a BEC in a one-dimensional lattice.

## 1. Introduction

The creation of Bose–Einstein condensates (BEC) in vapours of alkali metals offers a splendid opportunity to investigate the wide range of nonlinear effects involving atomic matter waves. A particular challenge is to develop methods for creating and controlling matter-wave solitons. Both dark and bright one-dimensional solitons in harmonic traps have been observed (Burger *et al*. 1999; Denschlag *et al*. 2000; Khaykovich *et al*. 2002; Strecker *et al*. 2002; Eiermann *et al*. 2004). A promising approach to obtaining multi-dimensional solitons consists in varying the scattering length (SL) of interatomic collisions. This can be achieved by means of sweeping an external magnetic field through the zero-SL point close to the Feshbach resonance (Fedichev *et al*. 1996; Inouye *et al*. 1998; Donley *et al*. 2001; Saito & Ueda 2002; Theis *et al*. 2004). The application of an ac magnetic field may induce a periodic modulation of the SL, opening the way to so called ‘Feshbach-resonance management’ (FRM) (Kevrekidis *et al*. 2003). A noteworthy FRM-induced effect is the possibility of creating self-trapped oscillating BEC solitons (breathers) without an external trap in the two-dimensional case. The underlying idea is that fast modulations create an effective potential on a slower timescale. This potential can stabilize the soliton.

The mathematical model for a BEC based on the Gross–Pitaevskii equation (GPE) with harmonic modulation of the SL was investigated by Saito & Ueda (2002, 2004), Abdullaev *et al*. (2003), and Montesinos *et al*. (2004). The conclusion was that FRM enables us to stabilize two-dimensional breather solitons even *without* the use of an external trap. According to these references, three-dimensional breathers conversely require at least a tight, one-dimensional harmonic trap, practically reducing the problem to two-dimensional (Görlitz *et al*. 2001). Later on we will call this approach a quasi two-dimensional treatment (Q2D). It will be defined more precisely below. Stabilization of three-dimensional droplets is also possible by including dissipative effects (Saito & Ueda 2004).

Trippenbach *et al*. (2005) demonstrated that three-dimensional solitons can be stabilized by a combination of FRM and a one-dimensional optical lattice (1D OL), instead of a one-dimensional harmonic trap. This issue has practical relevance, as a one-dimensional OL can easily be created by illuminating the BEC by a pair of counterpropagating laser beams so that they form a periodic interference pattern (Stecher *et al*. 1997). Incidentally, it is easier to realize a tight confinement configuration in an optical lattice than in a harmonic trap. Hence this environment may well be more friendly for creating quasi two-dimensional solitons. The lattice will be weak or strong depending on the intensity of the laser light. The combined OL–FRM stabilization of three-dimensional solitons is possible even in a weak lattice, when atoms confined in different cells do interact.

By analysing the stability charts in configuration space we find two distinct regions where stable solutions exist. The first of these regions has its counterpart in the Q2D treatment. The other region appears when the frequency of modulation exceeds the lowest excitation frequency of the confining potential. It is not present in the Q2D treatment; and so can only correspond to fully three-dimensional solitons. In the limit of tight confinement, the latter region moves up to extremely high frequencies and the Q2D stability chart is recovered.

## 2. Theoretical approach

We describe our system by the GPE in reduced units, including a time-dependent (FRM-controlled) nonlinear coefficient *g*(*t*) and an external potential *U*(* r*,

*t*):(2.1)Initially the BEC is in the ground state of a radial (two-dimensional) parabolic trap with frequency

*ω*

_{⊥}, supplemented, in the longitudinal direction, by ‘end caps’, induced by transverse light sheets. The configuration is very similar to that used to create soliton trains in a Li

^{7}condensate (Khaykovich

*et al*. 2002). A one-dimensional lattice potential in the axial direction is adiabatically turned on from

*ϵ*=0 to

*ϵ*=

*ϵ*

_{f}, see figure 1. Thus, the full potential is, with period normalized to

*π*,(2.2)where

*ϱ*is the radial variable in the plane transverse to

*z*, and the axial ‘end-cap’ potential,

*U*

_{0}(

*z*) is approximated by a sufficiently deep one-dimensional rectangular potential well. The width of the well determines the number of peaks in the final structure. The

*f*(

*t*) is a switch function (see figure 1).

The nonlinear interaction coupling is described by(2.3)Initially *g*_{1}(0)=0 and *g*(0)=*g*_{0}(0)>0. At some moment *t*_{1}, we begin decreasing *g*_{0}(*t*) linearly. It vanishes at time *t*_{2}, and remains zero up to *t*_{3}, when we start to gradually switch on the rapid FRM modulation of the SL. In the interval [*t*_{3},*t*_{4}], *g*_{0}(*t*) decreases linearly from zero to a negative value *g*_{0f} and the amplitude of the modulation *g*_{1}(*t*) increases from zero to *g*_{1f}, see figure 1. Simultaneously, both the radial confinement and end-caps are gradually switched off (see the behaviour of the *f*(*t*) function in figure 1). At times *t*>*t*_{4}, *g*(*t*) oscillates with a constant amplitude *g*_{1f} around a negative average value *g*_{0f}. Consequently, a soliton so created, if any, is supported by the combination of the one-dimensional lattice and FRM. Note that we choose a set of specific ramp functions *g*_{0} and *g*_{1} and they grow rather rapidly in time. In general the stability diagram may depend on the shape and rapidity of the ramp functions. We did not investigate this issue in detail, but believe that our results are to some extend universal. The reason we start from a positive value of *g*_{0} is so that the atoms can be uniformly distributed in the cells. This will allow stabilization of solitons in all filled cells simultaneously.

Numerical experiments following the path outlined in figure 1 indicate that it is possible to create stable three-dimensional solitons (see inset to figure 2). Before displaying the results, we first resort to the variational approximation (VA) in order to predict conditions on the modulation frequency and the size of the negative average nonlinear coefficient *g*_{0f}, necessary to support three-dimensional solitons.

## 3. Variational analysis

The VA can be applied to the description of BEC dynamics under diverse circumstances (e.g. Baizakov *et al*. 2003; Saito & Ueda 2003; Montesinos *et al*. 2004). Equation (2.1) is derived from the Lagrangian density(3.1)We use the VA for *t*>*t*_{3} and choose a hybrid ansatz composed of an hyperbolic secant function and a Gaussian for the solution for one lattice cell (calculations with two Gaussians give slightly inferior results, in terms of agreement with numerics). The amplitude is *A*(*t*), the overall phase is *ϕ*, the radial and axial widths are *W*(*t*) and *V*(*t*), respectively, and *b*(*t*) and *β*(*t*) are the corresponding chirps(3.2)

The reduced Lagrangian can be found upon substituting (3.2) into (3.1) and integrating over space. It is(3.3)where *I*_{1}=2*π* ln 2, *I*_{2}=(9/4)*πζ*(3), and *I*_{3}=(*π*/3)(4 ln 2−1). By varying this reduced Lagrangian with respect to *ϕ* we obtain the constant , where the integral extends over one cell of the lattice, and *n* is the number of occupied lattice cells. Notice that the total number of atoms is included in the definition of the nonlinear coupling *g*(*t*), and the total wavefunction is normalized to unity. When the other four variational equations are derived, we can deduce two dynamical equations for the widths by substituting for the chirps(3.4)(3.5)to obtain(3.6)(3.7)where *J*_{1}=(*I*_{1}−*I*_{3})/*I*_{2}, , and . All *J*s are positive. These equations describe the dynamics of a single peak, and so with small corrections can be applied to the problem of a BEC confined in a one-dimensional harmonic trap (Abdullaev *et al*. 2003; Saito & Ueda 2003; Montesinos *et al*. 2004). Note that for *f*(*t*)=0 one can derive a simple condition: . This equation depends on time only implicitly, and is used when describing collapse scenarios (§6).

In the corresponding Q2D treatment we drop the *z* dimension in equation (2.1). We assume that in this direction the profile of the wavefunction is fixed and reproduces the ground state *ψ*_{0} of the single lattice cell or harmonic potential, as in Saito & Ueda (2003) and Montesinos *et al*. (2004). The reduced potential in two-dimensional GP will take the form . In the VA we take *V* constant, *V*_{0} as found from the equation , and only solve equation (3.6). In numerical simulations we rescale the nonlinear coupling coefficient .

## 4. Numerical results

We simulated both the full GPE, equation (2.1), using an axisymmetric code (for the three-dimensional case), a Cartesian code (for Q2D), and the variational equations for comparison. Numerical simulations followed the path outlined in figure 1.

An example of the numerical results in the moderate confinement regime is shown in figure 2. The parameters used in the simulations would correspond, for ^{87}Rb atoms, to an OL period of *λ*=1.5 μm, an initial radial-confinement frequency of *ω*_{⊥}=2*π*×160 Hz, an FRM frequency of *Ω*=2*π*×19 kHz, a lattice depth of *ϵ*_{f}=25 recoil energies, an effective nonlinear coefficient of *N**a*=±2×10^{−5} m, where *N* is the number of atoms per lattice cell, with a total number of atoms in the range of 10^{4}–10^{6}. The corresponding values of the normalized parameters are given in the figure caption. This figure shows the evolution of the central-peak's amplitude versus time in dimensionless units, defined by equation (2.1). After an initial transient, a stable breathing structure is established. The difference between the dynamics in the three-dimensional treatment (lower curve) and the corresponding Q2D treatment (upper curve) is obvious. This is an example of a moderate confinement, when only the three-dimensional treatment describes the dynamics of the system properly. The inset shows the three-dimensional soliton structure for a fixed moment of time.

The dynamics in the weak confinement regime are displayed in figure 3. The figure reveals the influence of neighbouring solitons on the amplitude of the central peak. The thick curve corresponds to the evolution of the full multi-peak structure, similar to that shown in the inset of figure 2. The interaction between neighbouring solitons can be seen in the form of a beat. This phenomenon can be explained by the fact that the oscillation periods of the adjacent solitons are slightly different, as there is a small difference in the number of atoms in neighbouring cells. To obtain the thin curve we repeated the above calculations up to the time *t*=7500 and then removed all but the central peak. The structure so formed can be seen in the inset of figure 3. The beat is no longer visible. The difference between these two cases proves that the dynamics of stable breathers as described here is a collective multi-peak phenomenon.

In figure 4 we collected results of a systematic scan of parameter space based on GPE simulations and compared them with the predictions of the VA (a similar analysis can be performed if we replace one-dimensional OL with a one-dimensional harmonic trap—the conclusion does not depend on the form of confinement in the *z* direction). The agreement between VA and direct simulations is very good. The borders of the VA stability regions were found analytically, e.g. upon considering effective potentials. A detailed derivation is presented in §5. As seen from figure 4*a*, in a fully three-dimensional treatment we have found *two islands of stability*. Note the similarity of the lower regions of figure 4*a*,*c* to that of figure 4*b*, which portrays the results of the Q2D treatment. This region corresponds to Q2D solitons. On the other hand, the upper region contains fully three-dimensional solitons. They appear when the frequency of modulation exceeds the lowest excitation frequency of the confining potential. As we see in figure 4, when the strength of the lattice *ϵ*_{f} is increased, this region moves towards higher frequencies (and has risen beyond the scope of figure 4*c*), the Q2D region expands upwards, and the whole picture becomes more and more like figure 4*b*.

## 5. Stability analysis

We base our calculations on the variational equations (3.6) and (3.7) with *f*(*t*)=0. In the context of these equations, we will derive approximate formulas for borders of stability regions in parameter space. Several types of instability will be considered. Some of them can only be described by the full three-dimensional treatment. All stability conditions derived below lead to limiting curves in figure 4, with surprising accuracy. They are in good agreement with the results of numerical simulations of both variational equations (3.6) and (3.7) and the GPE (2.1).

*Q*2*D treatment.* In the Q2D treatment we neglect the dynamics of *V*, assuming *V*=*V*_{0} and only consider the variational equation for *W*, (3.6),(5.1)Here we introduced the constants *A*=*J*_{1}(*g*_{0}/*g*_{c}−1) and *B*=−*J*_{1}*g*_{1}/*g*_{c}. The critical nonlinearity *g*_{c}=−*V*_{0}*J*_{1}/*J*_{2} is equal to the threshold for collapse in the absence of rapid modulations.

In the adiabatic approximation we separate *W* into slowly and rapidly oscillating parts *W*=*w*+*δ*, and assume that the latter is relatively small, *δ*≪*w*. Equation (5.1) takes the form(5.2)In the adiabatic limit of small we obtain approximately(5.3)

The second formula is obtained using the first and can be interpreted as an equation of motion of a particle in an effective potential(5.4)where *C*=*J*_{1}(*g*_{0}/*g*_{c}−1)/2 and *D*=[*J*_{1}*g*_{1}/(2*g*_{c}*Ω*)]^{2}.

We denote the initial width *w*_{0} and choose . From (5.4) we see that the particle will be trapped by the potential provided that *U*_{ef}(*w*_{0})<0, as *U*_{ef}(∞)=0. This gives a stability condition(5.5)which leads to the lower limiting curve in figure 4*b*. Note that the above requires *C*>0, or equivalently *g*_{0}/*g*_{c}>1. This was the *only* relevant stability condition known so far (Saito & Ueda 2003; Montesinos *et al*. 2004).

We now consider if the adiabatic assumption used when deriving equations (5.3) is consistent with the result, i.e. whether *ξ*, defined above, remains small. Its maximum value can easily be estimated. If the oscillations in the potential (5.4) are small, then is small and so is *ξ*. If on the other hand these oscillations are large, one cycle can be split into two stages, (a) motion in the potential *U*_{ef}≈−*C*/*w*^{2} for large *w*, and (b) motion in the potential *U*_{ef}≈*D*/*w*^{6} for small *w*. The *ξ* coefficient reaches its maximum value during the (b) stage. All possible trajectories in the (b) stage potential are approximately described by a function, which scales with *D* and the initial conditions in the following way:(5.6)where *w*_{min} is the value of *w* at the smaller turning point and can be calculated from *w*_{0}, using equation (5.4). It turns out that if the right-hand side of (5.6) exceeds a critical value, an instability occurs in the numerical simulations. We observe an increase in the particle energy (equation (5.3) is no longer valid) around the smaller turning points, and consequently particle escape from the potential well (5.4). We stress that, for *ξ*_{max} less than critical, the system is stable for a very long time. Thus we obtain the upper limit in figure 4*b*.

*Three-dimensional treatment.* Now we extend our considerations to the case when the axial width of the peak, *V* is no longer constant and this influences stabilization. We keep the assumption that *V* is close to the width of the localized Wannier function of our lattice potential. Thus we write *V*=*V*_{0}+*η*, where *η*≪*V*_{0}.

The value of *V*_{0} is taken such that the first two terms on the right hand side of equation (3.7) cancel (there are in fact two such values and we take the lower one, known to be stable). From the third, nonlinear term we take only the oscillating part (the first two non-oscillating linear terms are assumed to dominate). In first order in *η*/*V*_{0} we end up with a harmonic oscillator type formula with a small driving term,(5.7)where and . The solution is(5.8)Now we rewrite equation (3.6), neglecting terms of higher order in *η*/*V*_{0},(5.9)and substitute *η* from (5.8). In the adiabatic limit, we repeat the reasoning that led to equation (5.4), to obtain an equation of motion for a particle in an effective potential, including a new term(5.10)where . If the last term in (5.10) can be neglected as compared to the middle one (this is true if *Ω* is close to *Ω*_{0}), another stability condition can be found(5.11)This gives the lower limit of the upper region in figure 4*a*. Strictly speaking, all terms in (5.10) should be taken into account when calculating these stability conditions. However, approximate formulas (5.5) and (5.11) prove to be surprisingly accurate for the range of parameters studied here. A condition analogous to (5.6) will follow:(5.12)However, another, lower border of the stability region can be found when considering resonance with double frequency of the forcing term. In this case we have to include higher, anharmonic terms in equation (5.7). In the following simple reasoning we will just take the first, *η*^{2} term,(5.13)where . We define . The linear solution of (5.13) is(5.14)and the next order correction *η*^{(1)}, obtained when the quadratic term is included, satisfies(5.15)For a solution in the form we get(5.16)Thus the threshold for appearance of a strong resonance when the driving frequency *Ω* comes close to 2*Ω*_{0} is given by(5.17)where we take and the minimal value of the width during the evolution as . In fact *w*_{min} cannot be any smaller, as it would imply another instability (compare to equation (5.5)). A quadratic in *Ω* is obtained. Simple manipulations show that *ϵ* tends to zero as *C* tends to zero when *g*_{0f}=*g*_{0c}≈20 and *Ω*=2*Ω*_{0} on figure 4*a*.

A fuller analysis would include an *η*^{3} term in (5.13) and the result would depend on *b*, the amplitude of *η*^{(1)}. This is why our *ϵ* (5.17) is just a threshold value. For a full analysis, see Landau & Lifshitz (1960). Thus we obtain the upper limiting curve figure 4*a*. The *η*^{3} term indicates that the resonance in question appears *above* this limiting curve in figure 4*a*.

Interestingly, the instability discussed in the Q2D analysis, associated with a breakdown of the adiabatic approximation (5.6), is strongly suppressed in the three-dimensional treatment in the range *Ω*_{0}<*Ω*<2*Ω*_{0} (see equation 5.12). This is due to flattening of the effective potential (5.10) by the additional term, and a decrease of the value of the *ξ* coefficient during the evolution.

## 6. Collapse and spreadout

When considering possible collapse scenarios we look at the leading terms in equation (3.6), when *τ*=*t*_{col}−*t* tends to zero. The sin(*Ωt*) term either tends to a constant or else to ±*Ωτ*, in which latter case it can be ignored in our considerations. We consider the leading small *τ* terms of equation (3.6) to be(6.1)and of equation (3.7)(6.2)

Now *Γ* and *Γ*′ can take either common sign, but we will see that only *Γ*, *Γ*′>0 can lead to collapse. As *τ*→0 either all three terms in equation (6.1) cancel, or else two predominate and cancel in lowest order *τ*. When all three cancel, *W*∼*τ*^{1/2}, *V*→*V*_{col} (1). When two terms cancel, either the first and third are of lowest order and *W*∼*V*∼*τ*^{2/5} (2), or else the second and third cancel because *V*_{col}=*Γ*/*J*_{1} and we find that *W*∼*V*−*V*_{col}∼*τ*^{2/3} (3). We now look at all three possibilities in some detail.

All three terms in equation (6.1) are of equal strength. We find from both equations (6.1) and (6.2),(6.3)(6.4)where , ,

*Γ*>*V*_{col}*J*_{1}.When the first and the third term in equation (6.1) cancel, we find(6.5)(6.6)where

*α*=(25/6)^{1/5}(*Γ*^{3}/*Γ*′)^{1/10},*β*=[25*Γ*′^{2}/(6*Γ*)]^{1/5},*Γ*,*Γ*′>0, and(6.7)(6.8)

The second and third term in equation (6.1) cancel when

*V*=*Γ*/*J*_{1}+*βτ*^{2/3},*W*=*ατ*^{2/3}. However, we find that*α*_{4}=−9*J*_{1}*β*/(2*Γ*) and . As*J*_{1}>0, and*Γ*,*Γ*′ take the same sign, no real*α*,*β*pair satisfies these equations. This possibility is ruled out on the level of the coefficients.

On the other hand, spreadout, when it occurs, is such that *W* and *V* are proportional to *t* at large times. Alternatively, just *W*∼*t* and *V*_{0}→const.

## 7. Conclusions

Only for strong confinement is the Q2D analysis an adequate approach. For moderate confinement the three-dimensional region of stability is richer from that of Q2D (see figure 4). A large, new region of stability appears in the three-dimensional treatment as compared to the Q2D treatment. This is proof of the fact that our solitons are truly three-dimensional. We are more used to the effect of adding a new dimension simply shrinking or abolishing the basin of stability of solitons or waves. For example, water waves are unstable with respect to perturbations along their direction of propagation only when the depth exceeds a critical value (Benjamin & Feir 1967). When, however, two-dimensional perturbations are allowed, there will always be an unstable angle regardless of the depth (Hayes 1973). Another example is that of one-dimensional solitons of the Nonlinear Schrödinger equation (NLS) with constant coefficients. They are stable in the one-dimensional case, but unstable in the two-dimensional or three-dimensional case. This is also true for some NLS waves (Anderson *et al*. 1979; Infeld & Rowlands 1980, 2000). In the problem treated here this is not the case. For situations corresponding to much of the quasi two-dimensional stability chart adding a degree of freedom *stabilizes* the soliton solution. The key to this dichotomy is the presence of a periodic modulation, absent in the above mentioned classical examples. This can be illustrated by a simple case involving oscillators. Take as the one-dimensional version a forced oscillator problem:(7.1)If *y* is fixed, the solution has a secular component *yt*/(2*ω*_{0})sin(*ω*_{0}*t*), and so the amplitude will grow as *t*. If, however, we allow a second degree of freedom, such that *y* also oscillates (for instance ) the solution stabilizes, unless . In general this can be the case when there are periodic modulations. This fact, obvious in oscillator theory, is perhaps less well known in the soliton context.

The main result of this paper is pointing out the possibility of creating fully three-dimensional breather solitons in a BEC confined by a one-dimensional optical lattice potential, corresponding to the upper region in figure 4. The stable patterns may feature a multi-cell structure. Depending on the strength of the confinement, we identified three different cases: (a) strong confinement; practically no evolution in the *z* direction, thus Q2D can be applied, (b) moderate confinement; fully three-dimensional soliton-train, but no interaction between cells, (c) weak confinement; solution in form a set of weakly interacting fundamental solitons. The scheme proposed here is based on a combination of FRM and a one-dimensional optical lattice, and could be implemented in an experiment relatively easily. This would open the way to the creation of robust three-dimensional solitons (breathers) in BECs.

## Acknowledgments

M.M. acknowledges support of the KBN grant 2P03 B4325, M.T. was supported by the Polish Ministry of Scientific Research and Information Technology under grant PBZ MIN-008/P03/2003, and E.I. acknowledges support of grant 2P03B09722. The authors are deeply indebted to Professor Boris Malomed for valuable discussions.

## Footnotes

- Received February 2, 2005.
- Accepted June 16, 2005.

- © 2005 The Royal Society