Dark-in-bright solitons in Bose-Einstein condensates with attractive interactions

We demonstrate a possibility to generate localized states in effectively one-dimensional Bose-Einstein condensates with a negative scattering length in the form of a dark soliton in the presence of an optical lattice (OL) and/or a parabolic magnetic trap. We connect such structures with twisted localized modes (TLMs) that were previously found in the discrete nonlinear Schrödinger equation. Families of these structures are found as functions of the OL strength, tightness of the magnetic trap and chemical potential, and their stability regions are identified. Stable bound states of two TLMs are also found. In the case when the TLMs are unstable, their evolution is investigated by means of direct simulations, demonstrating that they transform into large-amplitude fundamental solitons. An analytical approach is also developed, showing that two or several fundamental solitons, with the phase shift π between adjacent ones, may form stable bound states, with parameters quite close to those of the TLMs revealed by simulations. TLM structures are also found numerically and explained analytically in the case when the OL is absent, the condensate being confined only by the magnetic trap.


Introduction
The current experimental and theoretical studies of Bose-Einstein condensates (BECs) [1] have attracted a great deal of interest to nonlinear patterns which can exist in them, including dark [2] and bright [3] solitons. Two-dimensional (2D) excitations in the form of vortices were also realized experimentally [4]. Other nonlinear excitations, such as Faraday waves [5], ring dark solitons and vortex necklaces [6], were predicted to occur in BECs.
The behaviour of the condensate crucially depends on the sign of the atomic interactions: dark (bright) solitons can be created in BECs with repulsive (attractive) interactions, resulting from the positive (negative) scattering length. Recent experiments have demonstrated that 'tuning' of the effective scattering length, including a possibility to change its sign, can be achieved by means of a Feshbach resonance [7].
In this work we demonstrate that another type of localized nonlinear excitation, which may be regarded as dark solitons embedded in bright ones, can be created in attractive BECs. They may also be considered as bound states of two bright solitons with a phase shift π between them. In some respects, we follow the path of very recent results of [8]. These, in turn, were based on earlier works on the so-called twisted localized modes (TLMs) in discrete systems [9]; that is why we apply the same term, TLM, to these solitons of the combined dark-bright type. TLMs in a discrete system are bright solitons in which the field vanishes, changing its sign, at the central point. An alternative way to view a TLM is as a concatenation of an 'up-pulse' (i.e., sech(x −x 1 )) with a 'down pulse' (i.e., −sech(x − x 2 )), where x 1 and x 2 are appropriately separated, to form an up-down, two-pulse configuration. Notice, however, that our work is essentially different from [8], since we focus on the case of BECs with negative scattering length (such as 7 Li [10] and 85 Rb [11]), and demonstrate robustness and stability of TLMs in the latter context. These structures, if regarded as effectively dark solitons in the attractive BEC, are the inverse of recently predicted [12] bright solitons in repulsive BECs with an optical lattice (OL). It is expected that a characteristic longitudinal length of the TLM soliton in physical units will be on the order of the wavelength of light which induces the OL, i.e., ∼1 µm, and the number of atoms in the soliton may typically be ∼10 3 -10 4 .

64.3
We consider a quasi-one-dimensional (1D) (cigar-shaped) BEC with negative scattering length, the corresponding normalized Gross-Pitaevskii (GP) equation being [1,13,14] where u(x, t) is the macroscopic wavefunction, measures the strength of the parabolic magnetic trap and the last term accounts for the OL potential created by interference of two counterpropagating coherent light beams [15]- [20]. As we consider only (effectively) 1D configurations, the attractive interactions between the atoms (characterized by the negative scattering length) cannot produce collapse in the framework of the present model, hence various solitons have a chance to be stable.
It is obvious that the scaling invariance of equation (1) makes it possible to fix the optical wavelength, setting k ≡ 2π/λ = 1. This essentially means that all the remaining length scales of the problem are measured by comparison to the period of the OL which has been fixed to π. We follow this convention below. The number of atoms in the condensate, given by is a dynamical invariant of equation (1).
In the next section we report detailed numerical results for TLM solutions in attractive BECs, as modelled by equation (1). In section 3, we will support some of the numerical findings by analytical considerations. Section 4 concludes the paper.

Stationary TLM solutions and their linear stability
We seek stationary solutions to equation (1) in the form where µ is the chemical potential, and a real function v(x) obeys the equation supplemented with free boundary conditions. We have checked that the boundary conditions do not significantly affect the results. The corresponding boundary-value problem was solved by means of a finite-difference discretization. A second-order difference scheme with spacing x = 0.2 was used to approximate v xx . The resulting set of nonlinear algebraic equations was solved by means of a Newton iteration. More details on finite-difference schemes and Newton type methods can be found in [21]. The solutions of interest were obtained by choosing, as the initial condition for the Newton iterations, the following functional form: The ansatz (4) approximates a superposition of nonlinear Schrödinger (NLS) solitons. The presence of the OL suggests choosing positions of the solitons' centres, ξ j , as ξ j = j π/2, with some integer j . The sign factor (−1) j implies the phase difference π between adjacent solitons. This choice of the phase pattern is suggested by the known result that, in the discrete NLS equation, only such a configuration may be stable, see [22] and references therein. Similar

64.4
results for a BEC model with elliptic-function potentials were obtained in [23]. The number n of solitons in the initial ansatz (4) will be taken as n = 2, in order to find the TLM soliton proper, or n = 3, with the objective to find a bound state of two TLMs, see below. However, it will become clear below that, concatenating more such 'building block' structures (such as the ones with n = 2 or 3), one can, in principle, construct multi-pulse configurations with a larger number (n) of elementary pulses.
Once TLM states were found as solutions to equation (3) (starting with the initial approximation (4)), their stability was subsequently analysed by examining the corresponding linearized problem for small perturbations. To this end, the perturbed solution was taken as where δ is an infinitesimal amplitude of the perturbation, a(x) and b(x) are eigenfunctions of a perturbation mode and ω is the corresponding eigenfrequency, * standing for the complex conjugation. Using the ansatz of equation (5) to O(δ), one obtains the linearization problem as an eigenvalue problem where ω is the corresponding eigenfrequency, while (a, b ) T (where the T denotes transpose) is the corresponding eigenfunction. Using the same finite-difference scheme as the one used for computing v(x), we convert this to a matrix eigenvalue problem and find all the eigenvalues of the corresponding matrix. This is done by standard matrix eigenvalue solvers built into Matlab [24]. ω can, in general, be complex, i.e., ω = ω r + iω i , where the subscripts denote the real and imaginary parts of the eigenfrequency. If eigenfrequencies with non-zero imaginary parts exist, they lead to exponential instabilities, while their absence implies linear stability of the solution under consideration. This is monitored in the spectral plane pictures of the imaginary (ω i ) versus the real (ω r ) part of the eigenfrequency (see, e.g., the bottom right panels of figures 1-3 below). The presence of even a single ω with ω i = 0 implies instability. In the case when stable stationary solitons were identified according to this criterion (see below), the stability was also checked and confirmed by direct simulations of the full equation (1). For more details on the set-up and solution of the eigenvalue problem, the interested reader is referred to [25,26].
To systematically trace the evolution of numerically exact TLM solutions of equation (1), we performed one-parameter continuations in the three-dimensional parameter space of equation (3), (µ, , V 0 ). Recall that we have set λ ≡ 2π, hence λ is not a free parameter. We start by varying the OL potential's strength V 0 for fixed = 0 and µ = −1.5. We note that negative values of the chemical potential are typical of cases when the corresponding wavefunction pattern is localized, although µ can sometimes be positive (see below), which is explained by the presence of the OL potential in equation (3).
Numerical results corresponding to this case (varying V 0 ) are presented in figure 1. In the top subplot of the upper part of the figure, the number of atoms in the soliton, defined as per equation (2), is shown versus V 0 , and the largest instability growth rate, found from the associated linearized problem, is shown in the bottom panel of the upper part of figure 1 (zero instability growth rate implies that the soliton is stable). Three characteristic examples of the TLM solution are displayed, together with their linear-stability eigenfrequencies, in the lower part of figure 1.
This solution branch exists only for V 0 > 0.075. The presence of this cut-off point is expected, as the OL is crucial for the existence of TLM solutions, which are not found in the usual NLS equation (at least, in the absence of the parabolic-trap potential). Nevertheless, below we will present a case in which a TLM/dark soliton solution is possible even in the absence of OL, being supported solely by the magnetic trap. In the interval 0.275 < V 0 < 0.7, this soliton 64.5 family is subject to oscillatory instability which is accounted for by the so-called Hamiltonian Hopf bifurcation [27]. In the latter case, the collision of two pairs of eigenvalues of opposite Krein signature [27] results in the generation of a quartet of genuinely complex eigenvalues. The strongest instability, with Im ω ≈ 0.107 occurs at V 0 ≈ 0.5.
Similarly, we examined the evolution of the solutions, setting V 0 = 1 and µ = −1.5 in equation (3) and performing continuation in , starting from = 0. In this case too, oscillatory instability was found in a finite interval, 0.05 < < 1.61. The maximum instability growth rate (which is very large in this case) is 0.48, and occurs at ≈ 0.6. The results for this case are shown in figure 2.
In these cases, the oscillatory instability arises quite naturally. Indeed, similar instabilities had been observed in the discrete counterpart of the system and, in fact, in the next section we will give a theoretical justification for their occurrence. Another interesting numerical finding is that the TLM solutions appear to be much more robust than one would initially suspect. In particular, there seem to exist TLM solutions even when the OL has been almost completely smeared out by the harmonic trap (see e.g., the last panel of figure 2). We will return to this point in the next section.
We also performed the continuation of the TLM solutions with respect to the chemical potential µ, which yielded the results displayed in figure 3. In this case, = 0 and V 0 = 1 were chosen. The solutions exist in the region µ < µ max ≈ 0.125. The fact that the cut-off value µ max is positive is explained by a contribution of the mean value of the OL potential V 0 cos 2 x in equation (3). In this case, oscillatory instability is found in two finite intervals, −0.375 < µ < −0.1 and −0.075 < µ < 0.

Bound states of two TLM solitons
Complexes which may be considered as a bound state of two TLMs, or simply as a concatenation of three fundamental solitons with phase difference π between adjacent ones, were also studied. In this case, for symmetry purposes, cos x in equations (1) and (3) was replaced by sin x (recall we have set λ = 2π), and three pulses in the initial ansatz (4) were set at ξ j = −π, 0, π. We performed the continuation in , which demonstrated that the bound state persists to very large values of , see figure 4. The branch becomes unstable, through a quartet of complex eigenfrequencies, at > 0.015. A second instability, accounted for by another eigenvalue quartet, sets in at > 0.075. This second oscillatory instability can be easily explained: as we argue below, for each pair of anti-phase pulses there exists a pair of eigenvalues with negative Krein signature, which is known to give rise to a Hamiltonian Hopf bifurcation [27]. The strongest instability is found at = 0.23, when the two instability growth rates are 0.39 and 0.35. The first quartet returns to stability at > 0.29, while the second quartet is stabilized for > 0.33. The bound states are completely stable for > 0.33.

Evolution of unstable TLM solitons
In the cases for which the TLM solitons or their bound states were found to be unstable (due to the Hamiltonian Hopf bifurcations), we simulated the nonlinear development of the instability in the framework of the full GP equation (1). To this end, a finite-difference scheme with dx = 0.2 (the same step-size as in dealing with the stationary equation (3)) was used, and time integration was performed by means of the fourth-order Runge-Kutta scheme with a time step     The two top ones correspond, respectively, to the strong and weak oscillatory instability of the TLM soliton. In both of these cases, the TLM state evolves into a single fundamental soliton. A difference between the two cases is that, in the former one, with the relatively weak OL, the newly formed soliton keeps moving between two minima of the potential, while in the latter case, which corresponds to a stronger OL, the fundamental soliton becomes pinned. Finally, the bottom subplot of figure 5 demonstrates that the oscillatory instability of a bound state of two TLM solitons transforms it into a single fundamental soliton with a large amplitude.

The minimum size of the TLM state
The existence of the TLM state in the GP equation with OL can be explained, by means of perturbation theory, as a bound state of two fundamental solitons with phase difference π between them. Here, we consider the case of = 0 (no parabolic trap), but the analysis can be easily generalized to include the trap. If the present model is considered as a perturbed NLS equation, we can use an ansatz in the form of expression (4) for two solitons. The system's Hamiltonian reads Recall that we have set = 0; here, for generality's sake, we do not fix λ = 2π . Then, known methods (see a review in [28]) yield an effective net potential of the interaction between two solitons and interaction of both solitons with the OL. This potential is realized as a part of the Hamiltonian that depends on the coordinates ξ 1 and ξ 2 of the solitons and phase difference φ between them (cf a similar potential for solitons in the discrete-NLS model found in [22]) and is of the form The second term in this expression is the Peierls-Nabarro potential (see, e.g., [25] and references therein) induced by the OL. Equilibrium positions of the two-soliton system are determined by equations Straightforward consideration of equations (7) and (6) shows that a bound state of two solitons with φ = π , which corresponds to a TLM-like pattern, may exist if the separation L ≡ |ξ 1 −ξ 2 | between them exceeds the minimum value, This prediction for the minimum separation between fundamental-soliton components of the TLM state was tested numerically. As a typical example, in figure 6  In the latter panel, |u(x, t)| 2 in final and initial states and the potential are shown, respectively, by solid, dashed and dashed-dotted curves. The first example is displayed for = 0.6 and V 0 = 1; recall that the instability growth rate is ≈0.48 in this case, i.e., it is a case of strong instability. A random-noise perturbation of amplitude 10 −5 was added to the initial configuration in order to excite the instability. The second example corresponds to = 0 and V 0 = 1 for µ = −0.25. As the instability growth rate has a very small value in this case, 0.032, a randomperturbation seed with a larger amplitude, ∼10 −3 , was used. The last example pertains to an unstable bound state of two TLM solitons in the case of = 0.23 and V 0 = 1.0. In this case, a uniformly distributed random perturbation of amplitude 10 −5 was used to excite the instability.  (8) for the minimum separation between two π-out-of-phase fundamental solitons in the bound state (stars), and the separation between two density maxima in the TLM found from numerical data (circles), versus V 0 . Notice that the numerical results should be considered with an error bar of ±0.2, due to the discretization used.
of L min with V 0 for λ = 2π and η = √ 3. The distance between the centres of the solitons was extracted from numerical data as the distance between two local maxima of |u(x)|, an inherent discretization error of ±0.2 being imposed by the finite-difference scheme. As is seen, the agreement between the theoretical prediction and the numerically approximated values of L min is quite good.
A more detailed analytical consideration of dynamical properties of a perturbed bound state (following the lines of [22]) shows that, within the framework of the present approximation, the bound state is stable. Moreover, the existence of stable bound states of three fundamental solitons, which correspond to the bound state of two TLMs considered above, can also be demonstrated by means of this approach.

Oscillatory instabilities of the single-TLM and double-TLM states
It is well known (see, e.g., [25] and references therein) that the linearized equations for the perturbations defined in equation (5) can be rewritten, in terms of new variables U = a + b * and W = a − b * , as

64.13
The Krein signature of the eigenvalue can be written as [26] Analytical consideration is possible for ω that are real; then U and W are also real, hence K amounts to the expression As was done above, we approximate a TLM solution as a superposition, v(x) = P (x)+N(x), of two stationary fundamental-soliton waveforms, positive P (x) and negative N(x), the underlying assumption being a weak overlap between P (x) and N(x). Since P and N are, essentially, individual solitons, in the lowest-order approximation they obey equations [V (x) − µ]P = P xx + P 3 and [V (x) − µ]N = N xx + N 3 . Summing them up, and taking into account that the product |P N| is everywhere much smaller than P 2 + N 2 (due to the assumed weak overlap between P and N), we conclude that, up to the same lowest-order accuracy, the function P (x) + N(x) satisfies the equation L − (P + N) = 0 (recall that the operator L − is defined in equation (9)). Thus, in the lowest approximation, P (x) + N(x) is a zero mode of the operator L − . In a similar way, the following approximate result can be obtained for the function P (x) − N(x): The linearization around each of the TLM-constituent pulses (if considered in isolation) carries a pair of zero eigenvalues due to the phase invariance. The corresponding eigenfunction, W , has the shape of the pulse itself (P and N, respectively). There is also a generalized eigenfunction [26] for each pair. We expect to have two pairs of TLM eigenvalues originating from these pairs of the individual soliton ones, the corresponding eigenfunctions being (in the lowest approximation) W = P + N and W = P − N . The former one is related to the phase invariance of the system (in other words, to the conservation of the total number of atoms), which means that the corresponding eigenvalue pair remains exactly equal to zero, while the other one may become different from zero. Then, it follows from equations (12) and (9) that, in the lowest approximation, U = −P N(P − N). In turn, equation (11) yields Since, by definition, the signs of P (x) and N(x) are opposite, equation (13) implies that the Krein signature of this phase eigenmode is negative. A known consequence of this [27] is that an oscillatory instability will be generated by collision of this nonzero eigenvalue with other discrete ones (and/or those belonging to the continuous spectrum) that carry a positive Krein signature.
The above consideration indicates that, for every pair of the interacting π -out-of-phase pulses, there exists an eigenvalue with negative Krein signature and, consequently, a potentially ensuing Hamiltonian Hopf bifurcation. This conclusion fully agrees with our numerical results presented above, that have identified a single oscillatory instability for the single TLM soliton, and two such instabilities in the case of bound states of two TLM solitons.

TLM solitons in the absence of the OL
Finally, it is possible to present analytical arguments justifying the existence of TLM solitons (both single and multiple ones) even when the harmonic trap is overwhelmingly stronger than the OL potential (see the numerical results shown in the last panels in figures 2 and 4), or when the OL potential is simply absent. Within the framework of the same perturbation theory which was employed above to derive equation (6), each of the pulses constituting the TLM state is then subject to the action of two forces. One of them is exerted by the magnetic trap, while the other comes from the interaction between the pulses. Accordingly, the effective potential is Predictions following from this approximation were tested against direct numerical results obtained in the absence of the OL, i.e., with V 0 = 0 in equation (1). In figure 7, an example is displayed for = 0.2 and µ = −1.5. The top panel shows the effective potential (14), which, in this case, has a minimum at ξ 1,2 = ±1.35. The middle and bottom panels present a numerical solution, with maxima of |u(x)| at ξ 1,2 ≈ ±1.4. Thus, the perturbative approximation not only explains the existence of the TLM state in attractive BECs even in the absence of OL, but also accurately predicts the location of centres of the pulses whose concatenation generates the TLM structure. This, in turn, supports the numerical findings presented in section 2, which suggested that TLM solitons in the attractive BEC would indeed persist even for very tight magnetic traps. Thus, we conclude that the TLM solitons are a very robust feature of the system.

Conclusions
In this work we have demonstrated that counterparts of TLMs, which were earlier found in the discrete 1D NLS equation, can exist as robust objects in attractive, effectively 1D BECs, provided that the condensate is confined by the OL and/or parabolic magnetic trap. The TLMs in BECs may be considered as bound states of two fundamental solitons with a phase difference π between them, or, alternatively, as a dark soliton embedded in a broader bright one. The family of the TLM solitons and their stability were investigated in detail, by means of numerical continuation in the strengths of the OL or magnetic-trap parameters, as well as in the chemical potential. Stable bound states of two TLMs, alias a bound state of three fundamental solitons with a π phase shift between adjacent ones, were also found. In the cases when the TLMs are unstable, the development of the instability was monitored by means of direct simulations, with a conclusion that the TLM rearranges itself into a fundamental soliton.
We have also developed an analytical approach, which treats each fundamental soliton as a quasi-particle subject to the action of effective forces exerted by the OL, magnetic trap and interaction with another soliton. The analysis predicts the existence of the stable bound state of two fundamental solitons with a phase shift of π between them. The smallest possible distance between the bound solitons predicted by this approach is in good agreement with the numerically found size of the corresponding TLM. The analysis can be readily extended to predict that multipulse bound states with the phase shift π between adjacent solitons also exist and may be stable.
The origin of the instability of the TLMs, in the case when they are unstable, was also examined by means of an analytical approximation, and was found to stem from Hamiltonian Hopf bifurcations, which occur when eigenvalues with negative Krein signature collide with 64.15  Figure 7. The top panel shows the effective potential V eff , as a function of the coordinate ξ of the right-hand pulse, which is induced by the repulsion between the pulses and the restoring force from the magnetic trap, following the analytical expression (14). The potential minimum (equilibrium position) is located at ξ ≈ 1.35. The middle and bottom panels show the TLM configuration and its stability, as found in the numerical form. The equilibrium position corresponding to the numerical solution is inferred from the middle plot to be at ξ ≈ 1.4. The parameters are V 0 = 0, = 0.2 and µ = −1.5.
positive-signature ones. This result was based on the existence of a negative-Krein-sign eigenvalue set, that was constructed by means of the perturbative analysis.
It will be very interesting to generalize these considerations to higher dimensions-in particular to 2D settings. In that case, a relevant issue is the possible existence of a 2D counterpart of the TLM in the form of a bright vortex (which was found in the 2D discrete NLS equation, where it was shown to be stable under certain conditions [29]). A very recent preliminary result is that bright vortices, stabilized by a 2D OL, can indeed be found [30].