Accurate one-dimensional effective description of realistic matter-wave gap solitons

Bose–Einstein condensates (BECs) of dilute atomic gases are an invaluable source of insight into the quantum behavior of matter at macroscopic scales [1, 2]. Topics like vortex structures or nonlinear matter waves are being elucidated through BEC investigations [3]. In particular, considerable effort has been invested in recent years in the study of solitary matter waves (solitons) [4–6] due, in part, to their similarity with well-known equivalent structures in the field of nonlinear optics [7–9]. A great variety of solitons have been experimentally observed on the background of condensed atomic clouds. Dark solitons have been found [10–14] as excitations in systems with repulsive interatomic interactions (the analogue of self-defocusing nonlinear media), while bright solitons have been generated [15–17] as ground states when the interaction is attractive (self-focusing nonlinear media). A particularly important special case of bright solitons is that of gap solitons [6]. These matter-wave nonlinear structures are localized states that can occur in BECs with repulsive interatomic interactions if a periodic potential (such as that produced by an optical lattice) is present. Although gap solitons have been observed and studied in nonlinear optics since some time ago, only more recently have they been found in BEC experiments [18]. This has triggered a renewed theoretical interest in the study of these nonlinear structures.

Strictly speaking, solitons in a one-dimensional (1D) setting are exact solutions of the one-dimensional nonlinear Schrödinger equation (1D NLSE). Accordingly, in the framework of the mean-field approximation, matter-wave solitonic structures in 1D optical lattices have usually been investigated theoretically by means of the 1D Gross–Pitaevskii equation (GPE) [19–24], which is nothing but a nonlinear Schrödinger equation describing the evolution of the condensate wavefunction in the presence of the confining potentials. However, soliton-like states realized in BECs always have an intrinsic three-dimensional (3D) character, even when the trapping potential imposes a tight radial confinement. In fact, typical experimental parameters for the generation of gap solitons in BECs lie in a region of parameter space where the quasi-1D approximation is not strictly valid [6, 18] and, as a consequence, very often the 1D GPE cannot quantitatively reproduce the dynamics of realistic matter-wave gap solitons.

Various approaches have been put forward with the aim of incorporating in the description of condensed gases in elongated geometries the consequences of their intrinsic 3D character. Effective models for studying the stationary properties of highly anisotropic BECs beyond the analytically solvable perturbative and Thomas–Fermi regimes were proposed in [25–27]. An analytic expression for the local chemical potential μ_⊥(n₁) of elongated condensates as a function of the axial linear density was independently proposed in [27] and [28]. Such expression remains valid in the 3D–1D crossover regime and, in the case of [27], is also valid in the presence of an axially symmetric vortex. Effective equations have also been introduced to study topics such as solitons [29], the influence of the transverse confining on the dynamics of BECs in 1D optical lattices [30], or the emergence of Faraday waves [31]. Of particular relevance because of their simplicity and usefulness are the effective 1D equations derived in [32] and [27, 33], which govern the condensate dynamics accounting properly for the contributions from excited transversal modes through a nonpolynomial nonlinear term. Several variants of the above equations aimed at improving their accuracy to the price of incorporating adjustable [34] or a larger number of variational parameters [35, 36] have also been proposed. However, the intended increase in accuracy hardly compensates for the complexity of the resulting equations or the arbitrariness in the choice of free parameters. On the other hand, a number of generalizations of the above equations have been derived to study issues such as the expansion of elongated BECs [37], spin-1 atomic condensates [38], matter-wave vortices [39], BECs in anisotropic transverse traps [40], binary condensates in mixed dimensions [41], dark–bright solitons in two-component BECs [42], dipolar BECs [43], or localized modes in condensates with spin-orbit and Rabi couplings [44]. The equations of [32] and [33] have also been used to study BECs in nonlinear lattices [45], bright solitons in condensates with inhomogeneous defocusing nonlinearities [46], or the Bogoliubov spectrum of elongated condensates [47]. Since these effective equations are extensions of the standard 1D GPE (to which they reduce in the proper limit) and can be numerically solved with similar computational effort, they may offer a clear advantage over the latter. In this work, we will focus on the effective equation proposed in [27, 33]. This equation has been used in the description of matter-wave dark solitons in elongated BECs and has proven to give results in good agreement with the experiments [48–51]. As was demonstrated in [52], it can also be applied to the study of stationary gap solitons in the different physically relevant mean-field regimes. This is a nontrivial issue because this kind of systems are characterized by a small spatial scale (the lattice period d) that can represent a rather restrictive condition. However, to establish the utility of the model in realistic situations and determine whether it may be a real alternative to a fully 3D treatment one has to analyze its ability to reproduce a number of physical properties of fundamental interest in the theoretical description of stationary matter-wave gap solitons that were not studied in [52]. One of these properties is the dependence of the soliton chemical potential μ on the number N of constituent atoms, a quantity that can be used to define trajectories in the μ − N plane which are essential for the classification of gap solitons into different families. Another one is the frequency spectrum of the elementary excitations determining the dynamical stability of the gap soliton family. To provide an answer to the previous question, in this work we analyze to what extent the 1D GPE and the effective equation of [27, 33] can reproduce the above-mentioned fundamental physical properties in realistic situations. To this end we consider condensed atomic clouds confined by transverse harmonic potentials in a range of parameters typical of current experiments and compare with exact numerical results obtained from the full 3D GPE.

Since in this work we are interested in the mean-field regime, we shall restrict our analysis to physical configurations with a large number of atoms per lattice well, which are amenable to a treatment in terms of a mean-field macroscopic wavefunction satisfying the GPE.

The paper is organized as follows. Section 2 briefly reviews the theoretical models considered in this work. In sections 3 and 4 we analyze the numerical results: in section 3 we focus on the wavefunctions and the distinctive trajectories of the different gap soliton families, and in section 4 we study the frequency spectrum of the elementary excitations and the stability properties of these nonlinear structures. Finally, in section 5 we summarize our results and present our conclusions.

When a periodic potential is imposed on a BEC, the chemical potential of the system develops a band-gap structure. In these circumstances, it is possible to find either extended states with chemical potential inside the bands (nonlinear Bloch waves) or localized gap solitons between bands. At zero temperature and in the mean-field regime, these coherent states are accurately described by the 3D GPE [1, 2]:

$$\begin{equation} {\rm i}\hbar \frac{\partial \psi }{\partial t}=\left( -\frac{\hbar ^{2}}{2m}\nabla ^{2}+V(\mathbf {r})+gN|\psi |^{2}\right) \psi , \end{equation} \tag{ 1 }$$

where g = 4πℏ²a/m is the interaction parameter, with a being the s-wave scattering length, N is the number of atoms, and V(r) is the potential of the trap. In this work we will consider an external potential of the form V(r) = V_⊥(r_⊥) + V_z(z), which besides the axial term V_z(z) = V₀sin ²(πz/d) accounting for a 1D optical lattice with period d, includes a transverse harmonic trapping $V_{\bot }(r_{\bot })=m\omega _{\bot }^{2}r_{\bot }^{2}/2$, where $r_{\bot } ^{2}=x^{2}+y^{2}$.

A description in terms of the GPE corresponds to a classical mean-field description that neglects quantum correlations. It represents the classical limit of the underlying quantum field theory and becomes an excellent approximation in the limit N → ∞ with the nonlinear interaction parameter gN held constant. In practice, quantum fluctuations decay rapidly with N and one expects the GPE to be a good approximation for configurations with a few hundreds of atoms per lattice well.

Although the 3D GPE can be solved numerically with no approximations [53, 54], it demands a considerable computational effort when either the 3D character of the system or the nonlinearity increases. In order to simplify the computational problem, simpler models with reduced dimensionality are commonly used in the study of gap solitons in 1D optical lattices. Whenever it is possible to make the assumption that the condensate remains in the ground state of the transverse harmonic trap (so that higher transverse modes are not excited), a good approximation to equation (1) is given by the 1D GPE with rescaled nonlinearity:

$$\begin{equation} {\rm i}\hbar \frac{\partial \phi }{\partial t}=-\frac{\hbar ^{2}}{2m}\frac{\partial ^{2}\phi }{\partial z^{2}}+V_{z}(z)\phi +g_{1D}N|\phi |^{2}\phi , \end{equation} \tag{ 2 }$$

where g_1D = 2aℏω_⊥ is the 1D interaction parameter, obtained after averaging over the 'frozen' transverse degrees of freedom, and ϕ(z) is the axial wavefunction of the condensate. It is important to bear in mind that this equation is only valid when the axial energies are small enough in comparison with the radial quantum, ℏω_⊥, a condition that, in general, is not satisfied in realistic situations.

A simple dimensional analysis shows that the physical problem is characterized by three dimensionless parameters which can be conveniently chosen as

$$\begin{equation} \frac{E_{R}}{\hbar \omega _{\bot }},\;\frac{V_{0}}{E_{R}},\;\frac{Na}{d} \end{equation} \tag{ 3 }$$

where E_R ≡ ℏ²k²/2m is the recoil energy of the lattice, with k = π/d being the lattice wavevector, s = V₀/E_R is the lattice depth, and Na/d is the nonlinear coupling constant determining the mean-field interaction energy. The derivation of equation (2) relies on the assumption that the 3D wavefunction ψ can be factorized in independent radial and axial components, with the radial component corresponding to the Gaussian wavefunction of the ground state. Such assumption entails an implicit adiabatic approximation according to which the transverse degrees of freedom can instantaneously follow the axial dynamics. In order for this to be true, the characteristic time scale of the axial motion must be much longer than that of the radial motion (${\sim }\omega _{\bot }^{-1}$). Moreover, since the validity of the 1D GPE (2) requires the radial configuration to remain in the ground state, the axial energies involved must be sufficiently small in comparison with the characteristic energy of the transverse harmonic oscillator, ℏω_⊥. Taking into account that the axial energy of the underlying linear problem is of the order of the recoil energy E_R and that the nonlinear interaction energy is of the order of an₁ℏω_⊥ (with n₁ = N/d being the linear density of the gap soliton along the axial direction), the above requirements can be stated as

$$\begin{equation} E_{R}\ll \hbar \omega _{\bot },\;\;an_{1}\ll 1. \end{equation} \tag{ 4 }$$

These conditions, which ensure the permanence of the system in the transverse ground state, also guarantee automatically the validity of the adiabatic approximation. However, in order for the mean-field approximation implicit in equations (1) and (2) to be applicable, it is also necessary that N ≫ 1. Incorporating this latter requirement, the conditions for the validity of the 1D GPE can be rewritten as

$$\begin{equation} d\gg a_{\bot },\;\;1\ll N\ll d/a, \end{equation} \tag{ 5 }$$

where $a_{\bot }=\sqrt{\hbar /m\omega _{\bot }}$ is the radial harmonic-oscillator length. These conditions are very restrictive and usually are not met in real condensates. For instance, for ⁸⁷Rb atoms a = 5.3 nm and $a_{\bot }=10.78/\sqrt{\omega _{\bot }/2\pi }$ μm (with ω_⊥/2π given in Hz) while, typically, the lattice parameter d varies between 0.4 and 1.6 μm [55], the lattice depth s between 0 and 30 [56] and, for elongated condensates, ω_⊥/2π varies between 85 and 1000 Hz. Thus, typically, d ≲ 1.6 μm, a_⊥ ≳ 0.34 μm, and d/a ≲ 300 which makes the conditions (5) hard to satisfy. In particular, in the experiment of [18], where gap solitons were observed in a ⁸⁷Rb condensate, d = 0.39 μm, a_⊥ = 1.17 μm, s = 0.7, and N ∼ 900 atoms. In that of [56], d = 0.42 μm, a_⊥ = 1.14 μm, s = 0 − 30, and N ∼ 3 × 10^{5 87}Rb atoms. It is important to note, however, that even though most experiments have been carried out in the above parameter regimes, larger values for d and ω_⊥ (which can realize the 1D mean-field regime) are well within experimental reach.

Next, we briefly recall the 1D effective model proposed in [27, 33], which is an extension of the standard 1D GPE. As before, the starting point is the adiabatic approximation, which allows the factorization of the 3D wavefunction in the form

$$\begin{equation} \psi (\mathbf {r},t)=\varphi (\mathbf {r}_{\bot };n_{1})\phi (z,t), \end{equation} \tag{ 6 }$$

where the radial wavefunction φ is assumed to be normalized to unity and n₁ is the local linear density

$$\begin{equation} n_{1}(z,t)=N|\phi (z,t)|^{2}=N\int {\rm d}^{2}\mathbf {r}_{\bot }|\psi (\mathbf {r} _{\bot },z,t)|^{2}. \end{equation} \tag{ 7 }$$

After substituting equation (6) in equation (1) and averaging over the radial degrees of freedom one obtains the following 1D effective equation

$$\begin{equation} {\rm i}\hbar \frac{\partial \phi }{\partial t}\!=\!-\frac{\hbar ^{2}}{2m}\frac{\partial ^{2}\phi }{\partial z^{2}}+V_{z}\phi +\mu _{\bot }\phi , \end{equation} \tag{ 8 }$$

where μ_⊥(n₁) is the transverse local chemical potential which follows from the stationary 2D GPE

$$\begin{equation} \left( \!-\frac{\hbar ^{2}}{2m}\nabla _{\bot }^{2}+V_{\bot }(\mathbf {r}_{\bot })+gn_{1}|\varphi |^{2}\!\right) \varphi =\mu _{\bot } (n_{1})\varphi . \end{equation} \tag{ 9 }$$

When the dimensionless linear density is small enough (an₁ ≪ 1), this latter equation can be solved perturbatively and, to the lowest order (i.e., for φ(r_⊥) given by the Gaussian ground state of the corresponding harmonic oscillator), one obtains μ_⊥ = ℏω_⊥(1 + 2an₁). Substituting back in equation (8) and taking into account equation (7) it is thus clear that one arrives at the 1D GPE (2). However, a more general solution to equation (9) can be found. Using two independent approaches it was demonstrated in [27, 33] that, for condensates with no vorticity (the case we are interested in here), a very accurate solution, valid for arbitrary an₁, is given by

$$\begin{equation} \mu _{\bot }=\hbar \omega _{\bot }\sqrt{1+4an_{1}}. \end{equation} \tag{ 10 }$$

Substituting this expression in equation (8), one finally finds

$$\begin{equation} {\rm i}\hbar \frac{\partial \phi }{\partial t}=-\frac{\hbar ^{2}}{2m}\frac{\partial ^{2}\phi }{\partial z^{2}}+V_{z}(z)\phi +\hbar \omega _{\bot }\sqrt{1+4aN|\phi |^{2}}\phi . \end{equation} \tag{ 11 }$$

This effective 1D equation governs the axial dynamics of the condensate, having a much wider range of applicability than the standard 1D GPE (2). In particular, now the linear density an₁ can take arbitrary values and in the limit an₁ ≪ 1 one recovers equation (2). Yet, the applicability of the above equation still requires the fulfilment of the adiabatic approximation and, as far as the dynamical problem is concerned, this condition can be hard to satisfy in the presence of an optical lattice. Things are different, however, for the stationary problem, which is, in fact, the most relevant in the characterization of matter-wave gap solitons. In this case, the time scale of the axial dynamics tends to infinity and the adiabatic approximation always holds. As a result, unlike the usual 1D GPE (2), the effective model (11) will be able to correctly predict the trajectories in μ − N plane as well as the wavefunctions and stability properties of the different gap soliton families.

Next, we will apply both equations (2) and (11) to the study of stationary gap solitons in a regime corresponding to typical experimental parameters. To this end, by using a Newton continuation method we look for numerical solutions of the form ϕ(z, t) = ϕ₀(z) e^−iμt/ℏ, with μ being the condensate chemical potential. These solutions define distinctive trajectories in μ − N plane which, as already said, are essential for the classification of gap solitons into different families. To investigate how accurately the above 1D equations can predict such trajectories and, thus, the physical properties of the various families, we compare the predictions from both equations with those obtained from the numerical solution of the full 3D GPE.

In what follows we consider a ⁸⁷Rb condensate subject to a rather tight transverse potential of frequency ω_⊥/2π = 800 Hz and in the presence of a 1D optical lattice with period d = 1.5μm and depth s = 20. For these parameters, the recoil energy of the lattice takes the value E_R = 0.32ℏω_⊥.

Figure 1 shows the different theoretical predictions the above equations make for the family of fundamental gap solitons. This figure depicts the dimensionless chemical potential $\tilde{\mu }\equiv (\mu -\hbar \omega _{\bot })/E_{R}$ of the solitons versus the number N of constituent atoms. The 3D band-gap structure of the underlying linear problem is shown in the background, with shaded stripes corresponding to the allowed energy bands. The two narrow stripes located at $\tilde{\mu }=10.45$ and 16.7 are purely 3D Bloch bands associated with excited transversal states and, thus, cannot be accounted for by any 1D model. They are replicas of the lowest energy band (located at $\tilde{\mu }=4.2$) shifted up in energy, respectively, by two and four quanta of the radial harmonic oscillator, ℏω_⊥/E_R [57, 58]. Open symbols in this figure are exact results obtained by solving numerically the full 3D GPE (1), the solid red curve is the prediction from the 1D effective equation (11), and the dashed blue curve is that from the standard 1D GPE (2). As is evident, this latter equation clearly fails. In fact, even in the first (1D) gap, it makes errors larger than 100% and thus can hardly be considered valid. On the contrary, the effective model of equation (11) is in good quantitative agreement with the 3D results over the entire curve. Gap solitons in this family are nonlinear stationary states exhibiting a single peak well localized at a lattice site. A representative example is displayed in the inset, which shows the 3D atom density of a gap soliton of this family containing 325 particles (point A marked by an open triangle in the figure). Its chemical potential, μ = 4.7 ℏω_⊥ > ℏω_⊥, reflects the fact that the corresponding transversal state is composed of multiple harmonic-oscillator modes. The prediction the effective model (11) makes for this soliton has an error in the number of particles of only 2%. Figure 1 also shows that while the results from the 1D GPE get worse as the number of atoms increase, those from the effective equation (11) remain in good quantitative agreement with the 3D results.

Zoom In Zoom Out Reset image size

This behavior is corroborated by figure 2. Panels A and B in this figure display the dimensionless axial densities an₁ (left panels) and wavefunctions $\tilde{\phi }(z)=\sqrt{a N}\phi _{0}(z)$ (right panels) of the fundamental gap solitons corresponding, respectively, to points A and B in figure 1. As before, solid red curves have been obtained from the effective equation (11), dashed blue curves from the standard 1D GPE, and open symbols are numerical results obtained from the full 3D GPE. For reference, the lattice potential is also shown in arbitrary units (dotted lines). In order to deduce the axial wavefunction from the corresponding 3D results we have used that $\tilde{\phi }(z)=\sqrt{an_{1}(z)}\,{\rm e}^{{\rm i}\varphi (z)}$, where the linear density can be readily obtained by performing the integral in equation (7). As for the axial phase φ(z), we have defined it as the average local phase at every z-plane. Again the agreement between the predictions of equation (11) and the 3D results is very good, while the 1D GPE clearly fails.

Zoom In Zoom Out Reset image size

Similar conclusions can be drawn in the case of multiple gap solitons. Figure 3 displays the trajectories in μ − N plane corresponding to the composite family consisting of two in-phase peaks (see also figure 4). Solitons in this family are bound states formed by a symmetric superposition of two fundamental gap solitons of the family (1, 0, 0), i.e. those displayed in figure 1. Here we are using the notation introduced in [57], which enables us to classify the different fundamental gap solitons into families characterized by the quantum numbers (n, m, n_r) corresponding to the 3D Bloch band from which they bifurcate, with n = 1, 2, 3, ... being the band index of the corresponding 1D axial problem and m = 0, ±1, ±2, ... and n_r = 0, 1, 2, ... being, respectively, the angular-momentum and radial quantum numbers characterizing the transversal state. In this work, as is usually the case, we restrict ourselves to gap solitons consisting of fundamental constituents of the type (n, 0, 0), which are amenable to an effective treatment in terms of the 1D equations (2) and (11). Gap solitons with (m, n_r) ≠ (0, 0) have been rarely considered in the literature [57, 58]. They feature a nontrivial radial topology and require a more elaborate treatment [58]. The inset in figure 3 shows the atomic density (as an isosurface taken at 5% of the maximum density) of the soliton corresponding to point D (marked by an open triangle). This soliton contains 1765 atoms and has a chemical potential $\tilde{\mu }=16.6$. As can be seen, the two fundamental constituents occupy adjacent lattice sites and have the same phase (shown in the inset in terms of a color map). This is also apparent from figure 4, which depicts the axial linear densities an₁ and wavefunctions $\tilde{\phi }(z)$ of the gap solitons marked by C and D in figure 3. A simple comparison also shows that the wavefunctions of the composite solitons are indeed symmetric superpositions of those in figure 2.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Gap solitons consisting of two out-of-phase peaks are analyzed in figure 5. The upper trajectories in this figure (those bifurcating from the first Bloch band) correspond to composite solitons formed by an antisymmetric superposition of two fundamental (1, 0, 0) gap solitons. An example of a soliton of this family is displayed in the upper inset, which shows a representation of the 3D wavefunction of the soliton corresponding to point E in the figure. This wavefunction is given in terms of an isosurface of the atomic density (taken at 5% of its maximum) and a color map indicating the local phase at each point. As in the previous case, the two fundamental constituents of this compound soliton occupy adjacent lattice sites, but now they exhibit a π jump in their relative phase. A detailed comparison of the axial density and wavefunction of this soliton (given in the upper panels of figure 6) with those of the corresponding fundamental solitons (given in the lower panels of figure 2) confirms that this soliton is a bound state consisting of the antisymmetric superposition of two of the fundamental gap solitons previously studied.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The lower trajectories in figure 5, which bifurcate from the second 1D Bloch band, constitute a family that has been referred to as the subfundamental family [23]. Gap solitons in this family are associated with the first excited (n = 2) energy band of the corresponding 1D axial problem, which reflects in the fact that they exhibit two major peaks in the axial direction (i.e., one axial node) localized in a single lattice site (see the lower inset in figure 5 as well as panels F in figure 6). Unlike the previous case, these solitons (which can be regarded as excited states featuring an embedded dark soliton) are not bound states of two fundamental (1, 0, 0) gap solitons. They are fundamental solitons of the (2, 0, 0) family which can be used as elementary constituents of more complex compound structures.

As is apparent from the above results, the predictions from the 1D effective model (11) are always in good quantitative agreement with the 3D numerical results.

Stability is of primary importance for the experimental relevance of gap solitons. In this section we are interested in determining whether the effective equation (11) is able to correctly predict the stability properties of the gap solitons previously considered. To this end we perform a stability analysis based on the above effective equation and compare with the respective 3D results. We begin by deriving the Bogoliubov–de Gennes (BdG) equations corresponding to equation (11) and then obtain numerically the frequency spectrum of the elementary excitations. To derive the BdG equations, we perturb a given gap soliton stationary wavefunction ϕ₀(z) in the form [3]

$$\begin{equation} \phi (z,t)=[ \phi _{0}(z)+u(z)\,{\rm e}^{-{\rm i}\omega t}+v^{*}(z)\,{\rm e}^{{\rm i}\omega t}] \,{\rm e}^{-{\rm i}\mu t/\hbar } \end{equation} \tag{ 12 }$$

and substitute in equation (11) retaining only linear terms in the amplitudes u and v of the normal modes of the system. The numerical solution of the resulting linear equations provides the frequencies of the elementary excitations determining the dynamical stability of the soliton. A similar procedure has been followed to derive the BdG equations corresponding to the 3D GPE (1).

Results for gap solitons B, D, E, and F are given in the panels labeled with the same letters in figure 7. Left panels show the corresponding excitation spectra as obtained from the 1D effective model (11), while right panels display the eigenfrequencies of the BdG equations associated to the 3D GPE. The horizontal axis in each figure represents the real part of the excitation frequencies ω shifted by the chemical potential $\tilde{\mu }_{0}$ of the corresponding stationary soliton, $\tilde{\mu }=\tilde{\mu }_{0}+\mathrm{{\rm Re}}(\omega /E_{R})$. Vertical axes show the respective imaginary parts, $\mathrm{{\rm Im}}(\omega /\omega _{\bot })$. As is well known, the absence of complex eigenfrequencies in the spectrum implies the linear stability of the system. Of course, the excitation spectrum of the full 3D problem necessarily contains frequencies that cannot exist in the effective problem with reduced dimensionality (for instance, those associated to radial normal modes). This is an expected result. The point, however, is that, except in those cases where the solitons exhibit a nontrivial (unstable) transverse configuration, one might also reasonably expect that the stability properties of the different gap solitons will be dictated by the axial degrees of freedom, which ultimately are the ones responsible for the existence of self-trapping. If this expectation proves true, the effective 1D model (11) could be a reliable tool for the investigation of the stability of these nonlinear structures. As is apparent from a comparison between left and right panels in figure 7, this is indeed the case: except for the gap soliton F (lower panels), which is a bright soliton containing a dark soliton that extends considerably in the radial direction, the effective 1D model predicts accurately the complex eigenfrequencies determining the stability of the various gap solitons. According to the excitation spectra shown in the figure, the fundamental soliton B and the symmetric composite soliton D are linearly stable while soliton E (the antisymmetric counterpart of the latter) is unstable. Since in all cases considered $\tilde{\mu }_{0}=16.6$, it is clear from the figure that the spectrum of this soliton exhibits a pair of purely imaginary eigenfrequencies, indicating an exponential instability. As for soliton F, even though the effective model (11) still predicts correctly the (oscillatory) unstable nature of this soliton, it cannot give a quantitative account of its 3D excitation spectrum which now contains a quartet of complex frequencies not appearing in the corresponding 1D spectrum (which exhibits another quartet of complex frequencies, although only a pair is visible in the field of view of the figure). One expects these complex frequencies not appearing in the 1D spectrum to be associated with unstable transverse excitations of the embedded dark soliton (snake-like instability). To verify this and the previous conclusions we have introduced an additive Gaussian white noise to produce a random small perturbation in the stationary configuration of the different gap solitons and have followed the ensuing time evolution by integrating numerically both the effective 1D equation (11) and the 3D GPE (1). The results obtained are collected in figure 8. The upper left panel shows the dynamical evolution of soliton D in terms of a density map where brighter regions correspond to higher densities. This image, which has been obtained from the 1D effective equation (11), is indistinguishable from the corresponding 3D result and confirms the stability of the symmetric composite soliton D as predicted by the linear stability analysis of figure 7. The same holds true for the fundamental soliton B (not shown in the figure). For later convenience, the upper right panel displays the time evolution of the real part of the wavefunction of soliton D (in arbitrary units) obtained from both the 1D (top) and the 3D models (bottom). Since this soliton is stable, its dynamical evolution, apart from small local fluctuations induced by the Gaussian perturbation, is given essentially by a periodic amplitude modulation $\mathrm{{\rm Re} }(\phi (t))\simeq \phi (0)\cos (\mu t/\hbar )$ with period T ≃ 0.2 ms, as corresponds to a (quasi) stationary state. Thus stability implies, in particular, that the relative phase between any pair of points in the global amplitude envelope must remain (quasi) constant over time. Accordingly, even though an analysis of the density evolution of the antisymmetric composite soliton E (middle left panel in figure 8) seems to indicate that this soliton is stable (the same conclusions follow from both 1D and 3D results), a more detailed study including also phase information reveals that this is not true. Indeed, as can be seen from the middle right panel, in this case the time evolution of the real part of the wavefunction no longer is a periodic amplitude modulation. While the two major peaks are initially out of phase, both the 1D effective model (top) and the 3D model (bottom) predict that they become essentially in phase at t ≈ 4.5–5 ms, so that, this soliton turns out to be unstable, in agreement with the linear stability analysis of figure 7. This example clearly reflects the importance of the spectrum of elementary excitations in the analysis of the dynamical stability of these nonlinear structures.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the two lower panels we compare the evolution in time of the density of soliton F obtained from equation (11) (left panel) with that obtained from the 3D GPE (right panel). As is apparent, in this case the effective 1D model predicts a longer soliton lifetime than the 3D model, in good agreement with the linear stability analysis of figure 7. The inset in the right panel shows that the soliton decay occurs when the nodal plane of the embedded dark soliton begins to undergo a considerable deformation, thus, indicating that the decay is a consequence of a transverse modulational instability which clearly cannot be accounted for with a 1D model.

One may wonder why the effective model (11) is able to properly predict the stability properties of the gap solitons even though, in the presence of the optical lattice (and because of the high frequency axial motion it might induce), in general, it is unable to accurately reproduce the detailed evolution in time. To address this question, we begin by noting that the adiabatic (Born–Oppenheimer) approximation essentially consists in assuming that the axial dynamics is so slow in comparison to the typical time scale of the radial degrees of freedom that, at every instant t, the latter can adjust to their equilibrium configuration compatible with the axial configuration n₁ occurring at that instant. This is exactly what the stationary radial equation (9) reflects. In the presence of an optical lattice, one thus reasonably expects this approximation to work better in those situations where the time variable does not play a prominent role, such as occurs for the stationary solutions of the time-dependent GPE, which are completely determined from a (time-independent) eigenvalue problem. Something similar can be expected for the linear stability properties of the system. Indeed, such properties are completely characterized by the frequency spectrum of the elementary excitations, which follow from the solution of an eigenvalue problem (BdG equations) where time does not play any role. Even though an instability manifests itself in the evolution in time, the very existence of a certain unstable mode is a characteristic feature of the system which depends only on the interplay among the various energies and interactions contributing to the Hamiltonian. This explains why the effective model (11) can predict the stability properties of the gap solitons even though it is unable to properly account for their dynamical evolution.

Most theoretical studies of gap solitons realized in condensates loaded in one-dimensional (1D) optical lattices have been carried out in terms of the 1D Gross–Pitaevskii equation (GPE). This is motivated, in part, by the analogy of the latter equation with the 1D nonlinear Schrödinger equation governing the evolution of similar structures in nonlinear optics and, in part, by the fact that, as occurs with dark solitons, one would reasonably expect these structures to be more stable in a quasi-1D regime, where their potential decay would be inhibited by the strong radial confinement. However, it has been demonstrated that robust, long-lived gap solitons exist in a three-dimensional (3D) regime where many higher-order radial modes are excited [57]. Moreover, matter-wave gap solitons are intrinsically 3D and in realistic situations can hardly satisfy conditions (5), which are necessary for the validity of the 1D GPE. In these circumstances, effective equations like equation (11), which can take the 3D character of stationary Bose–Einstein condensates into account by incorporating contributions from higher-order radial modes, offer a clear advantage over the latter. However, the question remains open whether the effective 1D equation (11) represents a reliable alternative to a fully 3D treatment. In this work we have addressed this question by analyzing the ability of the above effective equation to reproduce the fundamental physical properties of stationary matter-wave gap solitons in realistic conditions. Our results demonstrate that the predictions from the effective model (11) are in good quantitative agreement with the exact numerical results from the full 3D GPE. In particular, unlike the standard 1D GPE (which fails in most cases of practical interest), the 1D effective model (11) correctly predicts the distinctive trajectories characterizing the different gap soliton families as well as the corresponding axial wavefunctions along the entire band gaps. Our results also show that this effective model can predict correctly the stability properties of the different gap soliton families. This follows from both a linear stability analysis in terms of the corresponding Bogoliubov–de Gennes (BdG) equations and a representative set of numerical computations monitoring the long time response of the system to a sudden small random perturbation. In particular, by numerically solving the BdG equations we have analyzed the prediction the effective model makes for the spectrum of elementary excitations, which proves to be essential for unambiguously determining the stability properties of certain matter-wave gap solitons. Our results show that, except in those cases where the soliton features a nontrivial (unstable) transverse configuration, the model is able to correctly predict the spectrum of complex eigenfrequencies responsible for the dynamical stability of the system. These results thus indicate that the above effective 1D model can be a useful tool for the physical description of realistic matter-wave gap solitons in 1D optical lattices.

As a natural extension of the present work, one might consider to quantitatively investigate the effects of quantum fluctuations on the above stationary matter-wave gap solitons. In the tight-binding regime (s ≫ 1), the tunneling of atoms between adjacent lattice sites can be strongly inhibited. Under these conditions, intersite atom-number fluctuations become significantly suppressed while the corresponding phase fluctuations become enhanced, which can lead to a reduced relative phase-coherence. While for fundamental gap solitons (those located at a single lattice well) with a few hundreds atoms these quantum effects are not expected to be very significant, they can be more important in the case of composite in-phase or out-of-phase solitons. The problem of quantifying the contribution of these quantum fluctuations can be addressed in terms of a stochastic phase-space description based on the truncated Wigner approximation [59]. In this respect, we note that since equation (11) properly captures the Bogoliubov excitation spectrum of the system, it represents a particularly convenient starting point for such an analysis. This will be the subject of a future publication.