Gap solitons in elongated geometries: the one-dimensional Gross-Pitaevskii equation and beyond

Matter-wave gap solitons (GSs) are localized modes that can be created in elongated Bose-Einstein condensates (BECs) loaded into one-dimensional (1D) optical-lattice (OL) potentials, with the intrinsic nonlinearity induced by repulsive interactions between atoms. GSs have been the topic of a large number of original works. Results produced by these studies were summarized in several reviews Kon1 ; Mor1 ; Carre1 ; Carre2 ; Barcelona . Within the framework of the mean-field approximation, the description of matter-wave patterns is based on the Gross-Pitaevskii equation (GPE) for the macroscopic wave function $\psi$ GPE :

which has proven to be very successful in reproducing experimental results for the zero-temperature BEC (for instance, for multiple vortices, as shown in detail in Refs. Vort ; Huht ). In this equation, $M$ is the atomic mass, the norm of the wave function $\psi$ is unity, $N$ is the number of atoms, and $g=4\pi\hbar^{2}a/M$ is the interaction strength with $a$ being the s-wave scattering length. In this work, we consider only repulsive interactions (i.e., $a>0$). Furthermore, $V_{\bot}(\mathbf{r}_{\bot})=(1/2)M\omega_{\bot}^{2}r_{\bot}^{2}$ is the radial-confinement potential and $V_{z}(z)$ is the axial potential, which may include the axial harmonic trap and the 1D OL, with depth $V_{0}$ and period $d$:

The energy scale in the underlying (no-interaction) linear problem is set by the OL’s recoil energy, $E_{R}=\hbar^{2}\pi^{2}/(2Md^{2})$.

Most commonly, theoretical studies of GSs in elongated geometries are carried out in terms of the 1D GPE Kon1 ; Mor1 ; Kiv1 ; Abd1 ; Mal1 ; Wu1 ; Wu2 :

where $g_{\mathrm{1D}}=2a\hbar\omega_{\bot}$. This is an effective evolution equation for the axial dynamics —described by the 1D wave function $\phi(z,t)$— which can be derived from the full three-dimensional (3D) GPE after averaging out the radial degrees of freedom under the assumption that the radial confinement is so tight that the transverse dynamics are frozen to zero-point oscillations. These conditions, however, are not easy to realize using typical experimental parameters Mor1 , and when such conditions are not met the transverse excitations can no longer be neglected, making an essentially 3D analysis necessary. In this work, we aim to investigate different physically relevant regimes and capture 3D effects in the generation and stability of matter-wave GSs in elongated BECs. This analysis should make it possible to determine the range of validity of the 1D GPE for the description of the GSs under typical experimental conditions, as well as characteristic features of the fundamental and higher-order GSs in realistic situations.

In principle, when the transverse excitations are relevant, Eq. (3) fails and one has to resort to the full 3D GPE (1), as recently done in Ref. PRA82 , or, alternatively, use extended 1D models that, within the framework of certain assumptions, take into account effects of higher-order radial modes on the axial dynamics of the condensate, leading to effective 1D equations with the cubic-quintic CQ or nonpolynomial Reatto1 ; PRA77 nonlinearities. In this work, we will consider both the full 3D GPE and the effective 1D model with a nonpolynomial nonlinearity, which was derived, for the case of the self-repulsive nonlinearity, in Ref. PRA77 . The latter one, which represents a simple generalization of the usual 1D GPE, that reduces to Eq. (3) in the appropriate limit, has demonstrated an excellent quantitative agreement with experimental observations Kev1 ; Wel1 ; Kev2 ; Carre3 (chiefly, for delocalized dark solitons, which are natural patterns in the case of the self-repulsion; however, the comparison was not reported before for localized GS modes). We will demonstrate that, while the range of applicability of the 1D GPE (3) is severely limited in realistic situations, the above-mentioned generalization gives a good description of stationary matter-wave GSs in virtually all cases of practical interest.

The paper is organized as follows. The model is formulated in Section II, where we also recapitulate the derivation of the effective nonpolynomial 1D equation following the lines of Ref. PRA77 , as this derivation is essential for the presentation of the results for the GSs. The main findings are collected in Sections III, where families of GS solutions are reported in several physically relevant regimes for strong, intermediate, and weak transverse confinement and shallow or deep OL potential. Both fundamental GSs and their multipeak bound states are considered. Only in the case of tight confinement combined with a shallow OL does the ordinary 1D GPE provide for a sufficiently accurate description of the GSs in comparison with results of full 3D computations. On the other hand, the effective nonpolynomial 1D equation provides for good accuracy in all cases; for the fundamental solitons and their bound states alike. The stability of the GSs is studied in Section IV by means of systematic simulations of the evolution of perturbed solitons. The conclusion is that the fundamental GSs are stable (except in a narrow region close to the upper edge of the band gaps), even in the case of strong deviation from the usual 1D description. Multisoliton bound states are stable if the OL potential is deep enough. Conclusions following from results of this work, including applicability limits for the mean-field approximation and 1D approximations, are formulated in Section V.

The model that was developed in Ref. PRA77 for a BEC in the absence of OL potentials resorted to the adiabatic approximation, to neglect correlations between the transverse and axial motions. This approximation assumes that the axial density varies slowly enough to allow the transverse wave function to follow these slow variations. Because the OL imposes a new spatial scale, which may be more restrictive, it is necessary to find out if the adiabatic approximation remains valid in the presence of the OL. To this end, we will now briefly recapitulate the derivation of the effective 1D equation based on this approximation.

The starting point is the ansatz based on the factorized 3D wave function

with $n_{1}(z,t)$ being the local condensate density per unit length characterizing the axial configuration:

To derive Eq. (5), we have assumed that the transverse wave function $\varphi$ is normalized to unity. Next, the substitution of Eq. (4) into the 3D equation (1) leads to

Because of the very different time scales of the axial and radial motions, it is natural to assume that the slow axial dynamics may be accurately described by averaging Eq. (6) over the fast (radial) degrees of freedom. Doing this, one obtains

Since $n_{1}(z,t)$ enters the last term of Eq. (8) merely as an external parameter, it is clear that, whenever the axial kinetic energy associated with the transverse wave function may be neglected; namely,

the radial dynamics decouple and $\mu_{\bot}(n_{1})\!$ can be determined without the knowledge of the axial wave function. Actually, the right-hand side of Eq. (9) is the chemical potential of an axially homogeneous condensate characterized by the density per unit length $n_{1}$ and wave function $\varphi$. For a sufficiently small value of the linear density ($an_{1}\ll 1$), the chemical potential of the lowest-energy state of this homogeneous condensate can be readily obtained perturbatively. In this case, to the lowest order, $\varphi(\mathbf{r}_{\bot})$ is given by the Gaussian wave function of the ground state of the radial harmonic oscillator, with the corresponding chemical potential being $\hbar\omega_{\bot}(1+2an_{1})$. As the linear density increases, more radial modes become excited and, in general, the ground state of the corresponding homogeneous condensate involves many harmonic-oscillator modes.

In previous works, it was shown using different approaches that, also in this case, an analytical solution can be constructed PRA77 . In particular, by using a variational approach based on the direct minimization of the chemical-potential functional, it was shown that, for any (dimensionless) linear density $an_{1}$, a very accurate estimate for the chemical potential of the ground state is given by $\hbar\omega_{\bot}\sqrt{1+4an_{1}}$ Ger1 . This expression can be easily extended to incorporate the case in which the condensate contains an axisymmetric vortex PRA77 (see also Ref. Luca , where the intrinsic vorticity was included into the derivation of the 1D nonpolynomial nonlinear Schrödinger equation in the case of self-attraction). However, in this work we restrict the consideration to zero vorticity. Using Eq. (9), we see that a sufficient condition for the second derivative in $z$ appearing in Eq. (8) to be negligible is

Taking into account an estimate, $K_{z}[\varphi]\sim\hbar^{2}/(2M\Delta_{z}^{2})$, where $\Delta_{z}$ is the characteristic length scale in the axial direction, Eq. (10) can be rewritten as

When this condition holds, $\mu_{\bot}(n_{1})$ coincides, to a good approximation, with the transverse local chemical potential of the stationary radial GPE:

and is given by

Finally, substituting Eqs. (13) and (5) into Eq. (7) and taking into account that $\varphi$ is a real normalized wave function (which implies the vanishing of the integral in the last term), we arrive at the following effective 1D equation to govern the slow axial dynamics of the condensate:

We note that the contribution from interatomic interactions enters the above equation through term $\hbar\omega_{\bot}(\sqrt{1+4aN|\phi|^{2}}-1)\phi$, which vanishes for $N\rightarrow 0$. Thus, Eq. (14) incorporates an energy shift $\hbar\omega_{\bot}$, which is irrelevant for the dynamics but simplifies the form of the equation and makes the global chemical potential that follows from this effective 1D equation exactly coinciding with the corresponding 3D result.

The above derivation demonstrates that the validity of Eq. (14) relies on two conditions:

(i) The typical time scale $\Delta_{t}$ of the axial motion must be much larger than the typical time scale of the radial motion ($\sim\omega_{\bot}^{-1}$).

(ii) The axial kinetic energy $K_{z}[\varphi]$ associated with the transverse wave function must be negligible.

The former requirement is necessary to allow the radial wave function to adiabatically follow the axial dynamics. While the specific temporal scale $\Delta_{t}$ depends on the particular initial conditions, typically, in the presence of an OL, the fulfillment of this condition gets more difficult as the lattice period $d$ decreases, or the linear density $an_{1}$ increases. In particular, a sufficient condition for the validity of the adiabatic approximation is: $d\gg a_{\bot}$, $an_{1}\ll 1$, with $a_{\bot}\equiv\sqrt{\hbar/(M\omega_{\bot})}$ being the radial harmonic-oscillator length. Nevertheless, such a constraint, which is hard to satisfy in realistic situations, is not a necessary condition. Actually, for stationary states $\Delta_{t}\rightarrow\infty$ and condition (i) always holds. In the present work we are interested in this case, which is the most relevant 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 an OL, the characteristic length scale $\Delta_{z}$ is typically on the order of the lattice period $d$, hence Eq. (11) becomes

where $E_{R}$ is the corresponding recoil energy. We will demonstrate that this condition is well satisfied for stationary GSs in most cases of interest.

It is clear that, when the linear density is small enough ($an_{1}\ll 1$), Eq. (14), which exhibits the nonpolynomial nonlinearity, reduces to the usual 1D GPE (3) with the cubic nonlinearity. This is the quasi-1D mean-field regime, which corresponds to condensates whose radial state, as given by a solution to Eq. (12), may be well approximated by the Gaussian ground state of the corresponding harmonic oscillator PRA77 ; PRA74 . Thus, Eq. (14) represents an extension of the 1D GPE for condensates with larger linear densities or, equivalently, larger mean-field interaction energies. In such condensates, the radial ground state satisfying Eq. (12) is, in general, a linear combination of many harmonic-oscillator modes (a similar situation takes place in the derivation of the effective 1D equation for fermionic gases by means of the density-functional approach SKA ).

Inspection of Eqs. (1) or (14) demonstrates that the dynamical problem is fully controlled by four dimensionless parameters, which may be defined as

As said above, the recoil energy $E_{R}$ sets the energy scale of the underlying linear problem, while the nonlinear coupling constant $Na/a_{\bot}$ determines the order of magnitude of the mean-field interaction energy in units of the radial quantum $\hbar\omega_{\bot}$. Note that the dimensionless linear density $an_{1}$ is proportional to $Na/a_{\bot}$.

In this work, we assume the axial confinement to be so weak that the corresponding harmonic-oscillator length, $a_{z}\equiv\sqrt{\hbar/(M\omega_{z})}$, is much larger than the lattice period $d$. Under these conditions, it is safe to neglect the modulation induced in the condensate density by the axial harmonic trap and set $\omega_{z}=0$, so that we are left with a uniform 1D lattice potential acting along the axial direction.

The stationary solutions of the 3D GPE are wave functions of the form $\psi(\mathbf{r},t)=\psi_{0}(\mathbf{r})\exp(-i\mu t)$, with $\psi_{0}$ obeying the time-independent GPE:

where $\mu$ is the chemical potential. When Eq. (14) is applicable, one can instead generate the stationary axial wave function $\phi_{0}(z)$ by solving the effective 1D equation

which, in the limit of $an_{1}\ll 1$, reduces to the time-independent version of the usual 1D GPE (3):

Matter-wave GSs are characterized by a chemical potential lying in a band gap of the energy spectrum of the underlying linear system. To construct solutions for realistic 3D matter-wave GSs in the effectively 1D geometry, we looked for numerical solutions to Eqs. (17) and (18) in different physically relevant regimes, and compared the obtained results with those produced by the usual 1D GPE (19) to determine to what extent 3D contributions are relevant. In particular, the comparison allows us to estimate the accuracy and range of applicability of the effective 1D equations (18) and (19).

Experimentally relevant values for the OL period $d$ range from $0.4$ to $1.6$ $\operatorname{\mu m}$ Mor1 , which implies that, for the condensate of ⁸⁷Rb, $E_{R}/(2\pi\hbar)$ ranges from $3.6$ $\operatorname{kHz}$ to $220$ $\operatorname{Hz}$. Taking into account that, typically, $\omega_{\bot}/2\pi\lesssim 1$ $\operatorname{kHz}$, we conclude that $E_{R}/\hbar\omega_{\bot}\gtrsim 1/4$. On the other hand, for the condensate of ²³Na one has $E_{R}/\hbar\omega_{\bot}\gtrsim 1$ Na ; therefore, in this case the GS energy is sufficiently large to excite higher modes of the radial confinement, which makes 3D contributions always relevant. Since the ⁸⁷Rb condensate is most relevant for the experimental realization of GSs in the quasi-1D setting Ober1 ; Isa1 , in what follows below we use particular parameters of this condensate to illustrate our results. Nevertheless, the results of the present work, which are always expressed in terms of relevant dimensionless parameters, are valid for any BEC. Actually, this is a direct consequence of the scaling properties of Eqs. (1), (3), and (14).

We have investigated the stability of the GS solutions by monitoring their long-time behavior after the application of a random perturbation Wu2 ; JYang . Specifically, we have perturbed the corresponding wave functions with a small-amplitude ($\sim 2\%$) additive Gaussian white noise and have monitored the subsequent nonlinear evolution for $1$ $\operatorname{s}$. To this end we have numerically solved the full 3D GPE (1), as well as the effective 1D equation (14), using a Laguerre-Fourier pseudospectral method with a third-order Adams-Bashforth time-marching scheme. Numerical integration of the 3D GPE for such a long time is a very demanding computational task. While it is possible for a certain GS to be metastable and decay on a time scale still longer than $1$ $\operatorname{s}$, in practice this time is long enough in comparison with the lifetime of the condensate per se, hence the soliton surviving for $1$ $\operatorname{s}$ may be categorized as stable.

We have found that both the full 3D GPE (1) and the effective nonpolynomial 1D equation (14) lead to the same conclusions regarding the stability of the GSs. However, once an unstable soliton begins to decay, one cannot expect, in general, the latter equation to reproduce the dynamical behavior correctly. The reason for this problem is that, as already mentioned in Sec. II, in the presence of the OL potential the above-mentioned adiabatic condition (i) is hard to fulfill in time-dependent settings. Of course, the same is true for the 1D equation (3), which is a limit form of Eq. (14) and, consequently, has a much narrower range of applicability. The results presented below have been obtained using the full 3D GPE (1). For each GS we have used at least two different basis sets and time steps to check that the results do not depend on details of the numerical procedure.

Our simulations demonstrate that, in general, $(1,0,0)$ fundamental GSs are stable, except in a narrow region close to the top edge of the band gaps. In particular, the solitons corresponding to points A, D, K, and P in Figs. 1, 6, and 9 remain stable up to $t=1$ $\operatorname{s}$. On the contrary, GSs in a shallow OL ($V_{0}/E_{R}=2$), located close to the top edge of a band gap, such as those corresponding to points C, H, J, and Q in Figs. 1, 6, and 9, turn out to be unstable. This instability manifests itself as a steady decay of the norm of the GS through emission of radiation, on a time scale that increases as one moves away from the top edge of the band gap. In this regard, one finds that the GS corresponding to point H decays much slower than the other ones. As the lattice depth $V_{0}/E_{R}$ increases, in general, GSs become more stable and the instability region shrinks. In particular, our simulations indicate that GSs such as those corresponding to points F, G, L, M and R in Figs. 1, 6, and 9 (which are fundamental GSs in the deep OL, with $V_{0}/E_{R}=20$, located close to the top edge of the band gap) also remain stable up to $t=1$ $\operatorname{s}$. This is in good agreement with previous analyses carried out in the context of deep optical lattices in terms of the ordinary 1D GPE (3) Kiv1 ; Wu2 .

An important result is that $(1,0,0)$ fundamental GSs remain stable even in the weak-radial-confinement regime ($E_{R}/\hbar\omega_{\bot}\simeq 1$). In this regime, GSs always have sufficient energy to excite higher radial modes, and, as a consequence, the 3D effects are always relevant and the usual 1D GPE (3) fails. Figure 11(a) displays the long-time behavior, after the application of the perturbation, of the fundamental GS shown in Fig. 10 [it corresponds to point P in Fig. 9(b)]. This figure shows (by means of a color map) the evolution of the axial density $an_{1}(z,t)$, which has been obtained from the 3D wave function $\psi(\mathbf{r}_{\bot},z,t)$ by integrating out the radial dependence, as per Eq. (5). The left panel in this figure shows the 3D condensate density at $t=0$ (i.e., just after the application of the perturbation) as an isosurface taken at $5\%$ of the maximum density, while the right panel represents the density at $t=1$ $\operatorname{s}$. This GS has $\widetilde{\mu}=1.59$, which implies $\mu=2.59\hbar\omega_{\bot}$. Note that, despite the fact that this chemical potential is greater than the quantum $\hbar\omega_{\bot}$, the GS remains stable up to $1$ $\operatorname{s}$. The same is true for the fundamental GS corresponding to point R in Fig. 9(c), which corresponds to a disk-shaped BEC containing more than $2300$ ⁸⁷Rb atoms with chemical potential $\mu=12.5\hbar\omega_{\bot}\gg\hbar\omega_{\bot}$. As Fig. 11(b) shows, this GS also remains stable.

Our simulations also demonstrate that the compound GS from Fig. 2, corresponding to point B in Fig. 1(b), is unstable. This is seen in Fig. 12(a), which shows how the soliton decays on a time scale $\simeq 120$ $\operatorname{ms}$ after the addition of a $2\%$ white-noise perturbation. This compound GS, which has $\widetilde{\mu}=1.75$ and $N=35$, was built in the shallow OL ($V_{0}/E_{R}=2$) as the symmetric combination of three fundamental constituents. We have found that, in such shallow lattices, GSs of this type are always unstable. However, it is not difficult to find stable solitons of this type in somewhat deeper lattices. An example is shown in Fig. 12(b), which represents a stable three-peak soliton with $\widetilde{\mu}=3.13$ and $N=59$, trapped in the OL with $V_{0}/E_{R}=4$ and $E_{R}/\hbar\omega_{\bot}=1/10$. Similar dynamics is observed for the GS displayed in Fig. 7, which corresponds to point I in Fig. 6(b). This is a five-peak compound generated in the shallow OL ($V_{0}/E_{R}=2$) in the regime of the intermediate radial-confinement strength ($E_{R}/\hbar\omega_{\bot}\simeq 1/4$). Our simulations indicate that this compound soliton is unstable. In fact, no stable five-peak solitons were found in such a shallow lattice. On the other hand, it is not difficult to find stable compounds of this kind for deeper OLs. An example is the five-peak pattern with $V_{0}/E_{R}=4$, $E_{R}/\hbar\omega_{\bot}=1/4$, $\widetilde{\mu}=2.83$, and $N=212$. In general, GSs naturally become more stable as the lattice depth increases. For deep OLs ($V_{0}/E_{R}=20$), the stability of compound GSs is essentially identical to that of their fundamental constituents. The GS in Fig. 4, which corresponds to point E in Fig. 1(c), is an example of a stable three-peak soliton realized in a deep OL.

To appraise the physical relevance of the results reported in this work, it is pertinent to recall that the mean-field treatment of the GSs (gap solitons), based on the GPE, is valid if $N\gg 1$ and the condensate is sufficiently dilute so that $a^{3}\overline{n}\ll 1$, where $a$ is the scattering length of the inter-atomic collisions and $\overline{n}$ is the atomic density. In shallow OLs (optical lattices) such conditions can be easily met, though the former one imposes a serious limitation on the range of applicability of the usual 1D GPE (19), which fails as $N$ increases. In deep OLs, where the tunneling between adjacent cells is strongly suppressed, the above conditions must be fulfilled at each site of the lattice. Under these circumstances, GSs may have a rather high local density. A good estimate for the 3D peak density $n(\mathbf{0})=N|\psi(\mathbf{0})|^{2}$ in terms of the peak axial density, $n_{1}(0)$, can be obtained from the following relation PRA77 :

Substituting the values of $an_{1}(0)$ obtained in Section III for different GSs into Eq. (25), one can easily verify that condition $a^{3}n(\mathbf{0})\ll 1$ is satisfied in all cases, which justifies the use of the mean-field description. The most critical situation occurs for the GS corresponding to point G in Fig. 1(c), which has $a^{3}n(\mathbf{0})\simeq 10^{-4}$. For the GS corresponding to point R in Fig. 9(c), one finds $a^{3}n(\mathbf{0})\simeq 5\times 10^{-5}$.

Despite the fact that applicability conditions for the mean-field treatment can be readily met, it is much harder to justify the validity of the usual 1D GPE (3), which requires both conditions $an_{1}\ll 1$ and $N\gg 1$ to hold simultaneously. Only in this case may the 3D effects be safety neglected. In most experimentally relevant situations, these conditions are not met, hence the applicability range of Eq. (3) turns out to be very limited. On the contrary, it has been shown in this work that the effective 1D GPE (14) with the nonpolynomial nonlinearity provides for an accurate description of the stationary matter-wave GSs in most cases of practical interest.

A relevant extension of the present analysis may be to GSs with intrinsic vorticity, as well as to mobility of stable solitons. It is also interesting to consider two-component GSs in the realistic 3D setting.