text:
# Solitons in Bose-Einstein condensates trapped in a double-well potential

## 1 Introduction

Experimental observation of Bose-Einstein condensates (BECs) in trapped dilute gases  had opened new exciting possibilities for manifestations of nonlinear phenomena in various geometries. Indeed, in the mean-field approximation (which usually applies to a great accuracy), the order parameter, which can be identified with the single-atom wave function, obeys the nonlinear Schrödinger (NLS) equation with an external potential, also known as the Gross-Pitaevskii (GP) equation . By varying the trap potential, the shape of the condensate can be tailored to sundry geometries, in particular, the shape of the condensate may be approximately circular, i.e., the condensate itself may be effectively two-dimensional (2D) or one-dimensional (1D, alias cigar-shaped) .

It was recently shown that not only the geometry of the condensate, but also the magnitude and the sign of the scattering length, which determines interactions between atoms in the condensate and, thus, the nonlinear term in the corresponding GP equation, can be manipulated by varying the external magnetic field near the Feshbach resonance . This opens additional possibilities to control the quantum macroscopic dynamics of BECs.

The NLS equation is a fundamental model for many physical media. A well-known application of the 1D NLS equation is the description of the pulse propagation in nonlinear optical fibers . Similar to optics, where bright and dark solitons are supported by the focusing and defocusing nonlinearity, respectively, in BECs the $`s`$-wave scattering interaction between atoms is a similar determining factor. Therefore, condensates constrained to the one-dimensional shape can form dark or bright solitons. Dark solitons were found in condensates with repulsive interactions , while condensates with attractive interaction were recently experimentally shown to form stable bright solitons . The appearance of solitons is one of the most interesting manifestations of nonlinear dynamics, and in the case of BECs it is of paramount interest, as in this case the solitons represent self-localized “waves of matter”.

The similarity of the NLS and GP equations suggests that many nonlinear phenomena in BECs may have their counterparts in nonlinear optics, in particular, in fibers. However, due to unique manageability of the properties of BECs, new setups can be studied, which were not realized in nonlinear optics. One of such setups is related to the above-mentioned possibility of the effective control of the scattering length in BECs, i.e., the sign of the nonlinearity in the governing equations.

Anticipating such experiments, in the present paper we find and study, analytically and numerically, a new family of solitons in two weakly coupled effectively 1D condensates, trapped in a double-well magnetic potential, when the $`s`$-scattering lengths of the two condensates have different magnitudes, with particular emphasis on the case when the scattering lengths have opposite signs.

The spatially separated condensates can be created by focusing a far-off-resonant intense laser beam, which generates a repulsive optical dipole force, into the center of a magnetic trap. We assume that the condensates have the cigar-like shapes, with the transversal dimensions being strongly constrained by the trap, see Fig. 1. The chemical-potential difference between the traps, $`\mu _0`$, can be managed, moving the position of the barrier-generating laser beam by means of electro-optic or acousto-optic modulators (see also Ref. ). Such a setup was realized in experiments (see, for instance, Ref. ). We neglect a variation of the trap potential along the longitudinal (say, $`x`$) direction assuming it to be weak, so that for localized solutions well inside the trap the external potential can be considered as flat. Finally, as we are interested in the soliton solutions, we take into account the kinetic energy contribution.

The corresponding coupled-mode equations, after obvious scaling transformations, may be cast in the following dimensionless form

$`i_tu+_x^2u+|u|^2u+v=0,`$ (1)
$`i_tv+_x^2v(\mu _0+a|v|^2)v+u=0,`$ (2)
where $`\mu _0`$ accounts for a difference of the chemical potentials between the two condensates, described by the order parameters $`u`$ and $`v`$. It is assumed that the nonlinear interaction is always attractive in the $`u`$-condensate, while in the $`v`$-condensate the strength and sign of the nonlinearity are controlled by the coefficient $`a`$ (e.g., repulsive nonlinearity if $`a>0`$). The derivation of the couple-mode equations similar to Eqs. (1)-(2) from the corresponding GP equation was discussed in many works (consult, for instance, Refs. ), the outline is placed in appendix A. Thus, the system (1)-(2) realizes an interesting interplay between the dispersion (i.e., the kinetic energy contribution), self-focusing and self-defocusing nonlinearities, linear coupling and the potential shift, which deserves the study.

Dynamics of two BECs in a magnetic trap was studied before. For instance, existence of a macroscopic quantum-phase difference, experimentally demonstrated in Ref. , was used to study the coherent atomic tunnelling between two weakly coupled BECs confined in a double-well potential . In Ref.  it was shown that the coherent oscillations due to tunnelling are suppressed when the number of atoms exceeds a critical value. The effect of the trap oscillations on the atomic tunnelling between two BECs was analyzed in Ref. . In Ref.  an analogy with the nonlinear-optical directional fiber couplers was used to account for the kinetic terms in governing coupled NLS/GP equations (for a review of nonlinear dynamics in optical fiber couplers see section 6 of Ref. ; this analogy will play an important in the present work too). Further use of the analogy with guided-wave optics led to the development of a nonlinear collective-mode theory for BECs trapped in an external potential . In Ref. , nonlinear modes in the form of chains of bright and dark solitons, which have no linear counterparts, were shown to be stationary solutions to the NLS equation with a multi-well external potential. A related model, based on a multi-component (i.e., vector) NLS equation, that describes repulsively interacting two-component BECs, possesses a coupled dark-bright soliton solution .

For $`a<0`$ (i.e., when the nonlinearities are focusing in both subsystems) equations similar to the system (1)-(2) have been studied in connection with the nonlinear dynamics in the dual-core optical fibers (see the above-mentioned review  and original papers ). However, the system (1)-(2) with $`a0`$, i.e., with opposite signs of the nonlinearity in the two condensates, was not considered before. Besides its relevance to the description of the coupled BECs, as argued above, the study of the model with positive $`a`$ is of general interest: as we demonstrate in this work, it gives rise to interesting soliton solutions which exhibit some unusual properties (see figures 4 and 6 below, for instance). Whereas letting the parameter $`a`$ to have *a*rbitrary negative values we also uncover novel features of the solitons and bifurcations studied before for the case of $`a=1`$.

The paper is organized as follows. In the next section, we aim to identify a domain of the soliton existence (mainly for the case of $`a0`$) by deriving equations and inequalities which must be satisfied by the soliton solutions. In this section, we also find particular (sech-type) soliton solutions for $`a<0`$. Generic numerically found soliton solutions are presented in section 3. In particular, a noteworthy result reported in this section is the occurrence of an unusual bifurcation (which was found for $`a=1`$): a soliton may have its amplitude vanishing and width simultaneously diverging, while the number of atoms in this configuration remains finite. In section 4, we derive a criterion of the Vakhitov-Kolokolov (VK) type for the soliton stability, and study the fate of unstable solitons subject to small perturbations by means of direct numerical simulations. As the result, it is found that, depending on the values of the system parameters and on the form of the initial small perturbation, the unstable soliton either gives rise to a breather (with persistent intrinsic vibrations whose amplitude is dependent on the slope of the number of trapped atoms vs. the chemical potential) or completely decays into radiation. The concluding section summarizes results obtained in the work.

## 2 The region of soliton existence and analytical solutions

To search for the stationary solitary-pulse solutions to Eqs. (1) and (2) we set

$$u(t,x)=e^{i\mu t}U(x),v(t,x)=e^{i\mu t}V(x),$$
(3)
where $`\mu `$ is the normalized (dimensionless) chemical potential, and $`U(x)`$ and $`V(x)`$ are real functions vanishing as $`x\pm \mathrm{}`$. Thus we arrive at a system of two real ordinary differential equations for the pulse profiles:

$`{\displaystyle \frac{\mathrm{d}^2U}{\mathrm{d}x^2}}+(\mu +U^2)U+V=0,`$ (4)
$`{\displaystyle \frac{\mathrm{d}^2V}{\mathrm{d}x^2}}+(\mu \mu _0aV^2)V+U=0.`$ (5)
Equations (4)-(5) were numerically solved, to look for solitons in a wide domain in the parameter space ($`a`$, $`\mu _0`$, $`\mu `$). However, before proceeding to the numerical solution, some results can be obtained in an analytical form, which will provide for some insight into the existence of solitons. The analytical results will make it possible to narrow a domain in the parameter space where it makes sense to look for solitons numerically. Moreover, the analytical results will be used to check the numerical solutions.

First of all, considered as a dynamical system, equations (4)-(5) possess a Hamiltonian,

$$=\left(\frac{\mathrm{d}U}{\mathrm{d}x}\right)^2+\left(\frac{\mathrm{d}V}{\mathrm{d}x}\right)^2+\left(\mu +\frac{U^2}{2}\right)U^2+\left(\mu \mu _0\frac{aV^2}{2}\right)V^2+2UV.$$
(6)
Evidently, solutions vanishing as $`|x|\mathrm{}`$ correspond to $`=0`$. We specify the class of the soliton solutions we look for as even solutions, $`U(x)=U(x)`$, $`V(x)=V(x)`$, with a single maximum at $`x=0`$ (we aim to consider single-humped, i.e., fundamental solitons). For the definiteness’ sake, we set $`U(x=0)U_0>0`$. As the first derivatives of the fields vanish at the central point, we can derive the following quartic equation for the soliton amplitudes $`U_0`$ and $`V_0V(0)`$ from Eqs. (4) and (5):

$$\left(\mu +\frac{U_0^2}{2}\right)U_0^2+\left(\mu \mu _0\frac{aV_0^2}{2}\right)V_0^2+2U_0V_0=0.$$
(7)
Further, multiplication of Eq. (4) by $`\mathrm{d}U/\mathrm{d}x`$, and (5) by $`\mathrm{d}V/\mathrm{d}x`$, and the integration from $`x=0`$ to $`x=\mathrm{}`$ leads to the following identities:

$$\left(\mu +\frac{U_0^2}{2}\right)U_0^2=2\underset{0}{\overset{\mathrm{}}{}}dxV\frac{\mathrm{d}U}{\mathrm{d}x},\left(\mu \mu _0\frac{aV_0^2}{2}\right)V_0^2=2\underset{0}{\overset{\mathrm{}}{}}dxU\frac{\mathrm{d}V}{\mathrm{d}x}.$$
(8)
In principle, both in-phase ($`V_0>0`$) as out-of-phase ($`V_0<0`$) fundamental solitons are possible (recall we have set $`U_0>0`$). Evidently, the right-hand sides of Eqs. (8) are negative for the in-phase, and positive for the out-of-phase solitons. Thus, we conclude that the fundamental solitons must obey two inequalities:

$`U_0^2`$ $`<`$ $`2\mu ,\mathrm{if}V_0>0,`$ (9)
$`aV_0^2`$ $`<`$ $`2(\mu \mu _0),\mathrm{if}V_0<0.`$ (10)
For $`a>0`$ (the opposite signs of the nonlinearities in the two subsystems), from these inequalities it follows that, for the in-phase fundamental solitons, the chemical potential of the $`U`$-condensate is negative ($`\mu <0`$), while, for the out-of-phase solitons, the chemical potential of the $`V`$-condensate is positive ($`\mu \mu _0>0`$).

The domain where the soliton solutions are permitted is determined by the dispersion law. Assuming the exponential decay $`Ue^{kx}`$, $`Ve^{kx}`$ for $`\left|x\right|\mathrm{}`$, and linearizing equations (4)-(5), we find two branches of the dispersion relation for the solitons:

$$k_1^2=\mu +\frac{\mu _0}{2}\sqrt{\left(\frac{\mu _0}{2}\right)^2+1},k_2^2=\mu +\frac{\mu _0}{2}+\sqrt{\left(\frac{\mu _0}{2}\right)^2+1}.$$
(11)
It is seen that the condition $`k_1^2>0`$ implies $`\mu <\mu _0`$, which excludes the out-of-phase solitons for this branch for positive $`a`$. The appealing conclusion that the $`k_1`$-branch corresponds to the in-phase soliton solutions is confirmed by our numerical results (for both positive and negative $`a`$), see the next section. This fact is also used in the proof of the Vakhitov-Kolokolov stability criterion for the in-phase solitons in section 4.

Finally, analysis of solutions of Eq. (7) by means of the inequalities (9) and (10) shows that, for fixed $`a>0`$, $`\mu _0`$, $`\mu `$, and $`U_0`$, there is one negative \[satisfying the condition (10)\] and, at most, two positive \[satisfying the condition (9)\] real solutions for $`V_0`$. Thus, there may be at most two branches of the in-phase solitons (however, we were able to find, in a numerical form, only one in-phase soliton solution corresponding to a given $`U_0`$ for fixed $`a`$, $`\mu _0`$, and $`\mu `$, see the next section).

For $`a<0`$, when the system (4)-(5) describes two coupled condensates with attractive interactions, it also provides for a straightforward generalization of the model describing the (mismatched) nonlinear dual-core optical fiber. In this context, the soliton solutions were studied in detail for the case $`a=1`$, i.e., equal Kerr coefficients in both cores . In this case, the parameter $`\mu _0`$ accounts for the phase-velocity difference between the cores .

For negative $`a`$, Eqs. (4)-(5) have special exact soliton solutions, both in- and out-of-phase ones, which generalize, respectively, the well-known symmetric and antisymmetric solitons in the standard model of the dual-core optical fiber . The exact solutions for the in-phase solitons are

$$\left(\begin{array}{c}U\\ V\end{array}\right)=\left(\begin{array}{c}A_1\\ A_1/\sqrt{a}\end{array}\right)\mathrm{sech}\left(\frac{A_1x}{\sqrt{2}}\right),A_1=\sqrt{2}\left(\mu \frac{1}{\sqrt{a}}\right)^{1/2},$$
(12)
which exist in the special case, when $`\mu _0`$ and $`a`$ are not independent parameters, but are related as follows:

$$\mu _0=\sqrt{a}\frac{1}{\sqrt{a}}\mu _0\left(a\right),\mu <\mu _{\mathrm{max}}^{(1)}\frac{1}{\sqrt{a}}.$$
(13)
The out-of-phase solitons are

$$\left(\begin{array}{c}U\\ V\end{array}\right)=\left(\begin{array}{c}A_2\\ A_2/\sqrt{a}\end{array}\right)\mathrm{sech}\left(\frac{A_2x}{\sqrt{2}}\right),A_2=\sqrt{2}\left(\mu +\frac{1}{\sqrt{a}}\right)^{1/2},$$
(14)
and they exist in the case

$$\mu _0=\mu _0\left(a\right),\mu <\mu _{\mathrm{max}}^{(2)}\frac{1}{\sqrt{a}}$$
(15)
\[recall $`\mu _0\left(a\right)`$ is defined in Eq. (13)\]. Note that the in-phase sech-type solitons correspond to the first branch of the dispersion law (11), $`A_1^2/2=k_1^2`$, while the out-of-phase solitons correspond to the second branch, $`A_2^2/2=k_2^2`$. Unfortunately, the solutions (12) and (14) cannot be continued to positive values of $`a`$.

Some additional exact results can be obtained concerning families of soliton solutions to Eqs. (4)-(5) with $`a<0`$. For instance, by using the standard bifurcation analysis it is easy to show that there is another family of the in-phase solitons, bifurcating from the solutions (12) at the point

$$\mu =\mu _{\mathrm{bif}}\frac{4a}{3\sqrt{a}},$$
(16)
if $`\mu _0`$ and $`a`$ are related as in equation (13) (this bifurcation is considered in detail numerically in section 3.2 and analytically in section 4). Setting $`a=1`$, one recovers the bifurcation which was found, in terms of the dual-core fiber, in Ref.  (in this case $`\mu _0=0`$). Note that due to the relation $`\mu _{\mathrm{bif}}<\mu _{\mathrm{max}}^{(1)}`$ (which can be easily checked to hold) this bifurcation is present for any $`a<0`$. Actually, the bifurcating in-phase solitons generalize the so-called $`A`$-type solitons (defined in Ref.  for $`a=1`$, see also Refs. ) for arbitrary (negative) values of $`a`$.

Another family of soliton solutions, the so-called $`B`$-type solitons, bifurcates from the out-of-phase solitons in a non-perturbative way . However, the $`B`$-type solitons were shown to be unstable in the whole domain of their existence, while the out-of-phase solitons could be stable only in a very narrow interval . For this reason, we discard the out-of-phase solitons from further consideration, so that all the solitons considered below are implied to be of the in-phase type.

## 3 Soliton solutions (numerical analysis)

We will discuss only the normalized quantities corresponding to the system (1)-(2). The numbers of the trapped atoms are also normalized accordingly, so that our discussion below is general instead of being tied to a particular shape and size of the external trap. The actual numbers of trapped atoms depend heavily on the particulars of the trap according to formula (53) given in appendix A. The transformation coefficient between the physical and “normalized” numbers of atoms is estimated there to be on the order of $`10^3`$ for the current experimental setups. Thus the soliton solutions we discuss below correspond to the interval from thousands to a hundred thousands of the trapped atoms. Those are the characteristic numbers of atoms in the current experiments. On the other hand, it is not surprising to come up with precisely this interval in our model, since it is derived under the assumption of the weak nonlinearity, in which case the upper bound on the number of trapped atoms is determined by the external trap. To avoid confusion, we note here that below by the “number of trapped atoms” we mean only the number of particles in the model (1)-(2).

To solve Eqs. (4)-(5) numerically the following iterative scheme was adopted

$`{\displaystyle \frac{\mathrm{d}^2U^{(n)}}{\mathrm{d}x^2}}+\left(U^{(n1)}\right)^2U^{(n)}+V^{(n)}=\mu ^{(n)}U^{(n)},`$ (17)
$`{\displaystyle \frac{\mathrm{d}^2V^{(n)}}{\mathrm{d}x^2}}\left[\mu _0+a\left(V^{(n1)}\right)^2\right]V^{(n)}+U^{(n)}=\mu ^{(n)}V^{(n)},`$ (18)
i.e. we solve the eigenvalue problem for $`\mu ^{(n)}`$ at the $`n`$-th iterative step. Selecting the lowest eigenvalue at each iterative step leads to convergence to a soliton solution111In some cases for $`a<0`$, when there are several soliton solutions corresponding to the same $`a`$, $`\mu _0`$, and $`\mu `$, we used the shooting method to obtain the branches of solitons, to which the iterative scheme (17)-(18) has poor convergence. (the use of this scheme, and the selection of the eigenvalue, are prompted by the well-known facts from the one-dimensional eigenvalue theory). We used the Fourier spectral (collocation) method with up to 256 grid points (see Refs.  for introduction into the spectral methods). Geometric convergence of the numerical solution to a soliton was noted for quite arbitrary initial profiles. As our method required only one of the two amplitudes to be specified (by appropriate normalization of the eigenfunction at each step), while the other was the result of the computation, we used the exact relations (7) to check the correctness of the numerical solutions. We also tested our approach, using the explicit soliton solutions (12).

### 3.1 Solitons for $`a0`$

In the model with the opposite signs of the nonlinearity in the two condensates, $`a>0`$, we have found a family of soliton solutions existing for *all* values of $`a`$ and $`\mu _0`$, and *all* the values of the chemical potential $`\mu `$ permissible by the dispersion law, i.e., which satisfy the condition

$$\mu <\mu _{\mathrm{max}}\frac{\mu _0}{2}\sqrt{\left(\frac{\mu _0}{2}\right)^2+1}.$$
(19)
It should be stressed that, in comparison with the previous works done in the context of nonlinear optics, these are essentially novel solitons, as all the previously studied cases  assumed the same sign of the nonlinearity in both cores of the system (although a possibility of the existence of bright gap solitons was shown in the case when the signs were opposite in front of the second derivatives ; in optics, this case is quite possible, corresponding to opposite signs of the group-velocity dispersion in a dual-core fiber with asymmetric cores, while in the case of BEC this case makes no sense).

Case I. $`a=0`$

We start with the special case $`a=0`$, when the second condensate is characterized by the zero scattering length (no nonlinearity). In figure 2 we plot the total number of atoms, and the numbers of atoms in each condensate as functions of the chemical potential $`\mu `$ for several values of the chemical-potential difference $`\mu _0`$. These dependencies, besides being important to quantity the BECs, determine the stability of the soliton solution: for fixed $`a`$ and $`\mu _0`$ ($`a0`$), the solitons are stable if the slope of the curve $`N=N(\mu )`$ ($`N`$ is the total number of atoms) is negative, $`\mathrm{d}N/\mathrm{d}\mu <0`$. This is a stability condition of the Vakhitov-Kolokolov (VK) type  (proof for the case under consideration is given in section 4). At $`\mu =\mu _{\mathrm{max}}`$, where the curves in Fig. 2 start from the zero number of atoms, the solitons bifurcate from the trivial solution $`U=V=0`$. As $`\mu `$ approaches $`\mu _{\mathrm{max}}`$, the width of the solitons tends to infinity, while their amplitudes decrease to zero.

Decrease in the chemical potential difference $`\mu _0`$ causes the curve $`N=N(\mu )`$ to sag in and develop a local maximum close to the upper limit value $`\mu =\mu _{\mathrm{max}}`$, see Figs. 2(b)-(d) \[in Fig. 2(d) we show only the sagging part of the curve, the shape of the whole curve being similar to that in Fig. 2(c)\]. It is interesting to note that for $`\mu `$ close to the value that corresponds to the maximum of $`N=N(\mu )`$ almost all the atoms are trapped in the linear $`v`$-condensate. The share of the number of atoms trapped in the condensate with no atomic interaction grows even further (for $`\mu `$ around the local maximum of $`N=N(\mu )`$) with further decrease of $`\mu _0`$ towards minus infinity. For instance, for $`a=0`$ and $`\mu _0=10`$, the number of atoms in the linear condensate can be up to $`99\%`$ (the respective function $`N=N(\mu )`$ is similar to that in Fig. 2(d), but with the local maximum at $`N2500`$). On the contrary, for any given $`\mu _0`$, for large negative values of the chemical potential $`\mu `$ almost all the atoms are trapped in the $`u`$-condensate.

The deformation of the solitons with the variation of the chemical potential for fixed $`a`$ and $`\mu _0`$ is illustrated by Fig. 3. The soliton solutions are plotted for three values of the chemical potential, which are marked by stars in Fig. 2(c), that correspond to the same total number of atoms. Note that, while at point 1, where $`\mu =3.9894`$, the $`u`$-component of the soliton is significantly higher, the $`v`$-component takes over as one moves to the right along the curve $`N=N(\mu )`$ in Fig. 2(c): for $`\mu =2.0236`$ and $`\mu =1.6576`$ (the points 2 and 3, respectively) the $`v`$-component of the soliton solution has a larger amplitude.

II. Small positive $`a`$

For small positive values of $`a`$ (roughly, for $`a0.1`$) and positive or small negative $`\mu _0`$, the numbers of atoms in the two condensates vs. the chemical potential have the forms similar to those of Figs. 2(a) through (c). For instance, for $`a=0.1`$ and for $`\mu _0=1`$, the curves are found to be similar to those in Fig. 2(c). Thus, for such values of $`a`$ and $`\mu _0`$, the solitons bifurcate from the trivial solution at $`\mu =\mu _{\mathrm{max}}`$ (19), as in the case of zero $`a`$. Moreover, for very small values of $`a`$ and some negative $`\mu _0`$, almost all the atoms can be trapped in the condensate with the (weak) repulsive inter-atomic interaction. As an example, a zoomed-in part of the functions $`N=N(\mu )`$, $`N_u=N_u(\mu )`$ and $`N_v=N_v(\mu )`$ around the local maximum, together with the soliton solutions, are displayed in Fig. 4 for $`a=0.001`$ and $`\mu _0=5`$. At point 2 in this figure, about $`96\%`$ of all the atoms are trapped in the $`v`$-condensate with the weak repulsive interaction \[a zoomed-in part of Fig. 2(d) would be similar to left panel of Fig. 4 with a slightly lower local maximum\].

However, the positive scattering length (self-repulsion) in the $`v`$-condensate brings new features too, as compared to the linear $`v`$-condensate. Most importantly, the grows of the local maximum of $`N=N(\mu )`$ with decrease of $`\mu _0`$ eventually saturates and changes for decrease, i.e., the maximum reaches its peak value at some negative $`\mu _0`$. Further decrease of $`\mu _0`$ causes the local maximum to disappear. Instead, a divergence develops at the boundary: $`\mathrm{d}N/\mathrm{d}\mu \mathrm{}`$ as $`\mu \mu _{\mathrm{max}}`$. For instance, for $`a=0.1`$ and $`\mu _0=2`$ there is no local maxima at all, while the share of atoms in the $`v`$-condensate is large for $`\mu `$ close to $`\mu _{\mathrm{max}}`$ due to the above mentioned divergence. Recall that the VK stability criterion demands the negative slope of $`\mathrm{d}N/\mathrm{d}\mu `$, hence the solitons with a large share of atoms in the repulsive condensate are unstable for such values of $`a`$ and $`\mu _0`$. This case is illustrated by Fig. 5, where we give the numbers of atoms vs. chemical potential and typical $`u`$\- and $`v`$-shapes of the corresponding unstable soliton with larger number of atoms in the repulsive condensate.

The following property of the soliton solutions is observed as $`\mu \mu _{\mathrm{max}}`$ (for curves similar to those in the left panel of figure 5): the soliton’s amplitude remains bounded, while its width is not. Therefore, our conjecture is that the solitons develop a “pedestal” (long shelf) of an increasing width, as the chemical potential approaches its limit value $`\mu _{\mathrm{max}}`$ (this limit is very difficult to study numerically, precisely for the same reason). In any case, all such solitons are unstable, since the slope of the corresponding curve $`N=N(\mu )`$ is positive, see Fig. 5.

We have checked that, for various small positive $`a`$ and sufficiently large negative values of $`\mu _0`$, the dependence of number of atoms on the chemical potential is quite similar to what is shown in Fig. 5. Vice versa, for arbitrary negative $`\mu _0`$ there is a threshold value of $`a`$ such that, for $`a`$ larger than this value, the derivative $`\mathrm{d}N/\mathrm{d}\mu `$ grows to infinity as $`\mu \mu _{\mathrm{max}}`$. For example, for $`\mu _0=1`$ it was found that, for $`a=0.25`$, the dependence of the numbers of atoms on the chemical potential has essentially the same form as in Fig. 5.

The results on the solitons for small (positive) values of $`a`$ and various $`\mu _0`$ can be summarized as follows: the soliton solutions with $`\mu `$ close to $`\mu _{\mathrm{max}}`$ and positive or small negative $`\mu _0`$ are similar to those displayed in Fig. 3, and for large negative $`\mu _0`$ are similar to the soliton displayed in the right panel of Fig. 5. On the other hand, for $`\mu \mu _{\mathrm{max}}`$ the solitons are similar to the solution with $`\mu =3.9894`$ displayed in Fig. 3 (tagged by 1).

III. Large positive $`a`$

For $`a1`$ and negative values of $`\mu _0`$, the curves $`N=N(\mu )`$, $`N_u=N_u(\mu )`$, and $`N_v=N_v(\mu )`$ are still similar to those in Fig. 5, while for positive $`\mu _0`$ they are similar to the curves in Fig. 2(a). Increasing $`a`$ further (i.e., making the repulsion between atoms in the $`v`$-condensate still stronger), we have found that the numbers of atoms in the condensates vs. the chemical potential take the form of what is displayed in the left panel of Fig. 5 also for positive $`\mu _0`$.

The value separating the two different types of behavior of $`N=N(\mu )`$ as $`\mu \mu _{\mathrm{max}}`$, i.e., given by Fig. 2(a) and the left panel of Fig. 5, was found to be $`a=1`$. At $`a=1`$ and $`\mu _0=0`$, the soliton has its amplitude gradually vanishing, as $`\mu \mu _{\mathrm{max}}`$, but, nevertheless, the number of atoms approaches a *non-zero* value, which is evident from the upper panel of Fig. 6. This seemingly strange bifurcation implies that in the limit $`\mu \mu _{\mathrm{max}}`$ the soliton’s width is coupled to its amplitude, thus the integral that gives the number of atoms remains constant (see the bottom panel in Fig. 6). Such a bifurcation may be naturally called a discontinuous one, as the soliton bifurcates from the zero solution without a visible jump in the shape, but with a jump in the number of atoms.

It may be interesting to check if the numerically found solitons for $`a0`$ are approximated by a sech-based ansatz. Let us assume, for example, that the soliton solutions can be approximated by the $`\mathrm{sech}^\alpha `$-functions, i.e., the solitons have the shape of $`f=A\mathrm{sech}^\alpha (x/d)`$, where $`A`$ (amplitude) and $`d`$ (width), in general, different for the two components, are the dynamic parameters (determined through the variational analysis) and the power $`\alpha `$ is fixed. This is an improvement of the usual sech-type ansatz. Given a numerical soliton solution we need to determine $`\alpha `$ such that the corresponding sechα-function approximates the soliton in some sense. This can be done in many different ways. Here we adopt as the proximity measure the following functional

$$P_f=\frac{\underset{\mathrm{}}{\overset{\mathrm{}}{}}dxx^2f\underset{\mathrm{}}{\overset{\mathrm{}}{}}dxf^2}{\underset{\mathrm{}}{\overset{\mathrm{}}{}}dxf\underset{\mathrm{}}{\overset{\mathrm{}}{}}dxx^2f^2}.$$
(20)
Inverting the function $`P=P_{\mathrm{sech}^\alpha }P(\alpha )`$ and using the numerically computed $`P_S`$ (for a numerically found soliton solution $`S=S(x)`$) we get the corresponding $`\alpha `$: $`\alpha =\alpha (P_S)`$. \[The proximity measure was chosen in the ad hoc way. However, there are two conditions to be satisfied: ($`i`$) the functional must be independent of both $`A`$ and $`d`$; ($`ii`$) the function $`P=P(\alpha )`$ must be monotonous (it is decreasing for $`P_f`$).\] If the solitons could be approximated by the above ansatz, the expected result would be a sharply peaked distribution of the power parameter values $`\alpha `$ computed for the numerically found solitons. Then one could assume the value of $`\alpha `$ corresponding to the peak as the approximation value.

We have computed the power parameter $`\alpha `$ in the above approximation for the numerically found soliton solutions for various $`a`$, $`\mu _0`$ and $`\mu `$. For instance, for the solitons shown in Fig. 3 the corresponding values of $`\alpha `$ are as follows (the subscripts refer to the $`u`$\- and $`v`$-components of the soliton)

$$\begin{array}{cc}\alpha _u^{(1)}=0.825,& \alpha _v^{(1)}=1.61,\\ \alpha _u^{(2)}=0.218,& \alpha _v^{(2)}=1.01,\\ \alpha _u^{(3)}=0.711,& \alpha _v^{(3)}=0.922.\end{array}$$

Evidently the soliton solutions cannot be approximated by the sechα-type functions with the same common power $`\alpha `$ due to a broad distribution of the power parameters corresponding to the solitons. Here we note that, in contrast, the variational approach based on the sech-type ansatz with the amplitudes and widths as the dynamic variables was successfully applied before to the solitons for $`a<0`$ (see Refs.  where the case of $`a=1`$ was considered).

### 3.2  Solitons for $`a<0`$

In the case of $`a<0`$, when the interactions are attractive in both condensates, it is enough to consider the interval $`1a<0`$, as Eqs. (4)-(5) are invariant against the following substitution:

$$U\stackrel{~}{U}=\sqrt{a}V,V\stackrel{~}{V}=\sqrt{a}U,\mu _0\stackrel{~}{\mu }_0=\mu _0,a\stackrel{~}{a}=\frac{1}{a},$$
(21)
and $`\mu `$ $`\stackrel{~}{\mu }=\mu \mu _0`$.

The solitons for $`a<0`$ were studied before in the context of the nonlinear dual-core optical fibers , chiefly in the particular case $`a=1`$. Though our new results, presented below, pertain to $`a1`$, we also consider in detail the previously studied case of $`a=1`$ and obtain some results, which were not known before.

We will frequently refer to the special case when $`a`$ and $`\mu _0`$ are related as follows $`\mu _0=\sqrt{a}1/\sqrt{a}\mu _0\left(a\right)`$. Note that only in this case Eqs. (4)-(5) admit (real) sech-type soliton solutions.

In section 2 it was shown that there are two branches of the in-phase solitons, one given by Eq. (12) for arbitrary (negative) $`a`$ and $`\mu _0=\mu _0(a)`$, and the other one bifurcating from these sech-type solitons at $`\mu =\mu _{\mathrm{bif}}`$, see Eq. (16). We have found numerically that these branches  *coexist* for arbitrary values of $`a`$ and $`\mu _0`$, although they do not always intersect, hence do not always undergo a collision bifurcation. Thus, the sech-type soliton (12) is a special case of a broader family of soliton solutions. Moreover, we have verified numerically that this broad family of (in-phase) solitons corresponds to the $`k_1`$-branch of the dispersion relation, see Eq. (11), i.e. the solitons exist for $`\mu `$ satisfying the condition (19).

To present the results in more detail, we first consider small (negative) values of $`a`$. We note that the bifurcation involving the two branches of solitons, with one branch corresponding to the sech-type solitons (12), belongs to the tangential type, i.e., it is a one-sided cusp bifurcation (consult, for instance, Ref. ). This is clearly seen from Fig. 7, where we take the values $`a=0.1`$ and $`\mu _0=\mu _0(a)`$ (this statement is also proven analytically in section 4). For small deviations $`\mu _0\mu _0(a)`$, and small $`a`$ (roughly for $`0.5a<0`$), the bifurcation is similar. For larger deviations $`\mu _0\mu _0(a)`$, the two branches develop a trend to collide, which gives rise to a picture resembling a collision bifurcation (see Fig. 8, where this type of behavior is shown, though for a larger negative value of $`a`$).

Finally, for small negative $`a`$ and for small but finite deviations $`\mu _0\mu _0(a)`$ we have checked that the solitons belonging to the branch which corresponds to the sech-type solitons at $`\mu _0=\mu _0(a)`$ can be uniformly approximated by the sech ansatz for *all* values of the chemical potential $`\mu `$. This is, however, not so for larger negative values of $`a`$ (see the discussion below for $`a=1`$).

We have found that the collision bifurcation is replaced by the coexisting branches when the deviation $`\mu _0\mu _0(a)`$ has large enough modulus. For instance, when $`a=0.8`$, the two branches are clearly separated for $`\mu _0=0.2`$, see Fig. 8, which is $`10\%`$ away from the corresponding value $`\mu _0(a=0.8)0.22`$ (the case of the small deviation $`\mu _0\mu _0(a)`$ is analyzed in detail for $`a=1`$ below).

In Fig. 9 we display the bifurcation diagrams for $`a=1`$. Here it is worth noting that, due to the symmetry against the substitution (21) the solution components $`U=U(x)`$ and $`V=V(x)`$ are interchangeable in this case. For instance, the sech-type solitons have identical components, $`U=V`$ (the symmetric solitons), while the solitons belonging to the other branch (the asymmetric ones) admit two possibilities: $`U>V`$, or $`V>U`$. Only the diagrams with $`U>V`$ are displayed in Figs. 9(b) and (d).

First of all, we notice that Figs. 9(a) and 9(b) reproduce the known bifurcation diagram for $`\mu _0=0`$  (note that in this case $`\mu _0=\mu _0(a=1)`$). Although it is demonstrated in section 4 that the asymmetric solitons bifurcate tangentially from the symmetric branch, we were not able to show this feature in the picture, because a region where it takes place for $`a=1`$ and $`\mu _0=0`$ is extremely narrow (for the same reason, this fact was left unnoticed in Ref. ).

As soon as the chemical-potential difference $`\mu _0`$ deviates from zero (in terms of the dual-core-fiber model, this corresponds to a phase mismatch between the cores), the collision bifurcation undergoes a significant transformation, see figure 9(c). It is interesting to note that, as soon as $`\mu _0`$ goes away from zero (we take $`\mu _0=0.01`$), the former sech-type branch divides itself into two different branches of soliton solutions, while the former asymmetric branch splits into two close curves, which pertain to different branches as well.

Consider the top part of the branch with the arrow in Fig. 9(c). In this region, the power parameters $`\alpha _u`$ and $`\alpha _v`$, numerically computed using Eq. (20), differ from $`1`$ by less than $`2\%`$ . Although this curve corresponds to the solitons which are well approximated by the sech ansatz and the corresponding soliton solutions have nearly identical components, $`UV`$, it nevertheless can be identified with the deformed asymmetric branch of solitons. This statement is supported by comparison of Fig. 9(d) with Fig. 9(b): in both cases, the bifurcation diagram has a pitchfork-type form, but in the case of Fig. 9(d) the stem and the middle prong of the pitchfork, both lying very close to the sech-type branch for $`\mu _0=0`$, belong to different branches of solutions. We were not able to trace the bifurcation point $`\mu _{\mathrm{bif}}`$ by continuing Fig. 9(c) in the direction of negative $`\mu `$, i.e., for $`\mu \mathrm{}`$, and checking stability of the solitons (the bifurcation point is the threshold of the soliton instability, see section 4). Thus, when $`a=1`$, the collision bifurcation is replaced by two coexisting branches of solitons already for $`\mu _0=0.01`$.

Finally, we observe that the branch of the soliton solutions denoted by the thin curves in Figs. 7 through 9, which corresponds to the asymmetric solitons for $`a=1`$ and $`\mu _0=0`$, undergoes a turning-point bifurcation. The turning point is clearly seen in these figures.

As it was mentioned above, we would not consider bifurcations of the out-of-phase (e.g., antisymmetric) solutions, as they are always unstable (see also Ref. ). For the same reason, we do not consider of higher-order multi-humped solitons, which are also found in the numerical solution, but turn out to be unstable in all the cases.

## 4 Soliton stability and evolution under perturbations

After having found the soliton solutions, the next necessary step is to analyze the soliton stability under action of perturbations, as well as the fate of unstable solitons. Accordingly, we split this section into two parts. In the first part, we prove the above-mentioned VK stability criterion, $`\mathrm{d}N/\mathrm{d}\mu <0`$, for the in-phase soliton solutions of Eqs. (4)-(5) with $`a0`$ and show where the VK criterion applies to the solitons for $`a<0`$. For instance, we demonstrate that the bifurcation point (16), discussed in the previous section, is (quite naturally) an instability threshold for the sech-type solitons (12). Our proof of the stability criterion for the soliton solutions of the system (1)-(2) generalizes a similar proof for the solitons of a single (scalar) equation, which can be found in Refs. . In the second part of this section, we numerically study the evolution of unstable solitons subject to initial perturbations.

Here we point that for $`a=1`$ (with $`\mu _0=0`$) the sech-type soliton solutions (12) discussed in the present paper, and the in-phase solitons in general, reduce to the solitons which were studied in Ref. , in the context of the dual-core optical fibers. The stability properties of the in-phase solitons for $`a=1`$ were investigated in Ref. , and evolution of unstable solitons under perturbations was simulated in Ref. . Below, we (in particular) generalize the stability results of Ref.  for arbitrary $`a<0`$ and $`\mu _0`$, which is relevant to the model (1)-(2) for trapped BECs.

### 4.1 Stability analysis

I. Stability of the in-phase solitons for $`a0`$

The soliton stability with respect to small perturbations is determined by the system (1)-(2) linearized about the soliton solution. Thus we need to consider a perturbed soliton solution:

$$u=e^{i\mu t}\left[U(x)+u_1(t,x)\right],v=e^{i\mu t}\left[V(x)+v_1(t,x)\right],$$

where $`U(x)`$ and $`V(x)`$ are the stationary soliton solution components, while $`u_1(t,x)𝒰_R(t,x)+i𝒰_I(t,x)`$ and $`v_1(t,x)𝒱_R(t,x)+i𝒱_I(t,x)`$ represent the perturbation. Linearizing equations (1)-(2) with respect to the perturbation we arrive at the system

$$_t\left(\begin{array}{c}𝒰_R\\ 𝒰_I\\ 𝒱_R\\ 𝒱_{}\end{array}\right)=\left(\begin{array}{cccc}0& L_0^{(u)}& 0& 1\\ L_1^{(u)}& 0& 1& 0\\ 0& 1& 0& L_0^{(v)}\\ 1& 0& L_1^{(v)}& 0\end{array}\right)\left(\begin{array}{c}𝒰_R\\ 𝒰_I\\ 𝒱_R\\ 𝒱_{}\end{array}\right),$$
(22)
where the linear operators are defined as follows:

$$L_0^{(u)}\left(\frac{\mathrm{d}^2}{\mathrm{d}x^2}+\mu +U^2(x)\right),L_0^{(v)}=\left(\frac{\mathrm{d}^2}{\mathrm{d}x^2}+\mu \mu _0aV^2(x)\right),$$

$$L_1^{(u)}\left(\frac{\mathrm{d}^2}{\mathrm{d}x^2}+\mu +3U^2(x)\right),L_1^{(v)}=\left(\frac{\mathrm{d}^2}{\mathrm{d}x^2}+\mu \mu _03aV^2(x)\right).$$
(23)
The solitons are unstable if the matrix operator on the right-hand side of equation (22) has an eigenvalue $`\lambda `$ with the positive real part \[in fact, the eigenvalues are real, since $`\mathrm{\Lambda }_0`$ is non-negative, see Eqs. (24)-(25) and the discussion below\]. The system (22) can be written in an equivalent second-order form, hence the corresponding eigenvalue problem can be reformulated for a fourth-order differential operator acting on a two-component vector $`(X_u,X_v)`$, with $`X_u`$ and $`X_v`$ defined as follows (consult also appendix B)

$$𝒰_R(t,x)=e^{i\mathrm{\Omega }t}X_u(x)+\mathrm{c}.\mathrm{c}.,𝒱_R(t,x)=e^{i\mathrm{\Omega }t}X_v(x)+\mathrm{c}.\mathrm{c}..$$

Here, $`\mathrm{\Omega }=i\lambda `$ has the meaning of a frequency. The linear eigenvalue problem then takes the form

$$\mathrm{\Lambda }_0\mathrm{\Lambda }_1\left(\begin{array}{c}X_u\\ X_v\end{array}\right)=\mathrm{\Omega }^2\left(\begin{array}{c}X_u\\ X_v\end{array}\right),$$
(24)
where we have introduced two symmetric matrix operators:

$$\mathrm{\Lambda }_0\left(\begin{array}{cc}L_0^{(u)}& 1\\ 1& L_0^{(v)}\end{array}\right),\mathrm{\Lambda }_1\left(\begin{array}{cc}L_1^{(u)}& 1\\ 1& L_1^{(v)}\end{array}\right).$$
(25)
The proof of the stability criterion relies on properties of the factorization operators $`\mathrm{\Lambda }_0`$ and $`\mathrm{\Lambda }_1`$. We first aim to show that $`\mathrm{\Lambda }_0`$ is non-negative, i.e., the scalar product $`\mathrm{\Psi }|\mathrm{\Lambda }_0|\mathrm{\Psi }`$, defined as $`_{\mathrm{}}^+\mathrm{}dx\mathrm{\Psi }^{}(x)\mathrm{\Lambda }_0\mathrm{\Psi }(x)`$ (here $``$ stands for the Hermitian conjugation), is non-negative for any real two-component vector $`\mathrm{\Psi }(x)=(\psi _1(x),\psi _2(x))^T`$. Indeed, the following inequalities follow from the definitions the operators $`L_0^{(u)}`$ and $`L_0^{(v)}`$:

$$dx\psi ^{}(x)L_0^{(u)}\psi (x)dx\frac{V(x)}{U(x)}|\psi (x)|^2,dx\psi ^{}(x)L_0^{(v)}\psi (x)dx\frac{U(x)}{V(x)}|\psi (x)|^2.$$
(26)
To arrive at Eq. (26) it is enough to note that, due to Eqs. (4)-(5), these operators can be rewritten as

$$L_0^{(u)}=\frac{1}{U(x)}\frac{\mathrm{d}}{\mathrm{d}x}U^2(x)\frac{\mathrm{d}}{\mathrm{d}x}\frac{1}{U(x)}+\frac{V(x)}{U(x)},L_0^{(v)}=\frac{1}{V(x)}\frac{\mathrm{d}}{\mathrm{d}x}V^2(x)\frac{\mathrm{d}}{\mathrm{d}x}\frac{1}{V(x)}+\frac{U(x)}{V(x)}.$$

Taking advantage of the inequalities (26) and of the fact that, for the in-phase soliton solutions, the product $`U(x)V(x)`$ is positive, we get

$$\mathrm{\Psi }|\mathrm{\Lambda }_0|\mathrm{\Psi }=dx\left\{\psi _1^{}(x)L_0^{(u)}\psi _1(x)+\psi _2^{}(x)L_0^{(v)}\psi _2(x)\psi _1^{}(x)\psi _2(x)\psi _1(x)\psi _2^{}(x)\right\}$$

$$dx\left\{\frac{V(x)}{U(x)}|\psi _1(x)|^2+\frac{U(x)}{V(x)}|\psi _2(x)|^2\psi _1^{}(x)\psi _2(x)\psi _1(x)\psi _2^{}(x)\right\}$$

$$=dx\left|\sqrt{\frac{V(x)}{U(x)}}\psi _1(x)\sqrt{\frac{U(x)}{V(x)}}\psi _2(x)\right|^20.$$

The operator $`\mathrm{\Lambda }_0`$ has, in fact, a zero eigenvalue, with the soliton solution proper $`(U,V)^T`$ being the corresponding eigenfunction:

$$\mathrm{\Lambda }_0\left(\begin{array}{c}U(x)\\ V(x)\end{array}\right)=0.$$
(27)
In its turn, the other factorization operator $`\mathrm{\Lambda }_1`$ also has a zero mode,

$$\mathrm{\Lambda }_1\left(\begin{array}{c}\mathrm{d}U(x)/\mathrm{d}x\\ \mathrm{d}V(x)/\mathrm{d}x\end{array}\right)=0,$$
(28)
which, however, has a node (zero) at $`x=0`$. Thus, according to the Sturm oscillation theorem, the operator $`\mathrm{\Lambda }_1`$ may have negative eigenvalues (it is found numerically that $`\mathrm{\Lambda }_1`$ has either one or two negative eigenvalues). In general, it may have two nodeless eigenfunctions corresponding to two negative eigenvalues:

$$\mathrm{\Lambda }_1\left(\begin{array}{c}f^{(1)}(x)\\ g^{(1)}(x)\end{array}\right)=\lambda _1\left(\begin{array}{c}f^{(1)}(x)\\ g^{(1)}(x)\end{array}\right),\mathrm{\Lambda }_1\left(\begin{array}{c}f^{(2)}(x)\\ g^{(2)}(x)\end{array}\right)=\lambda _2\left(\begin{array}{c}f^{(2)}(x)\\ g^{(2)}(x)\end{array}\right),$$

where $`f^{(1)}g^{(1)}>0`$ and $`f^{(2)}g^{(2)}<0`$. However, it is shown in appendix C that $`\mathrm{\Lambda }_1`$ cannot have an eigenfunction of the second type for $`a0`$ (corresponding to a negative eigenvalue).

Return now to the eigenvalue problem (24). We have

$$\mathrm{\Omega }^2\mathrm{\Lambda }_0^1\left(\begin{array}{c}X_u\\ X_v\end{array}\right)=\mathrm{\Lambda }_1\left(\begin{array}{c}X_u\\ X_v\end{array}\right)+C_0\left(\begin{array}{c}U\\ V\end{array}\right),$$
(29)
where $`C_0`$ is a scalar constant. In deriving (29), we have taken into account that for $`\mathrm{\Psi }=(X_u,X_v)^T`$

$$(U,V)|\mathrm{\Psi }=0,$$
(30)
since $`\mathrm{\Lambda }_0`$ is a symmetric operator. Thus, the relation (29) is justified. The minimum eigenfrequency $`\mathrm{\Omega }_0`$ can be also found from the following minimization problem:

$$\mathrm{\Omega }_0^2=\mathrm{min}\frac{\mathrm{\Psi }|\mathrm{\Lambda }_1|\mathrm{\Psi }}{\mathrm{\Psi }|\mathrm{\Lambda }_0^1|\mathrm{\Psi }},$$
(31)
with $`\mathrm{\Psi }=(\psi _1,\psi _2)^T`$ subject to the constraint (30). To evaluate the sign of $`\mathrm{\Omega }_0^2`$, one has to evaluate the sign of the numerator in the expression (31). Due to the constraint (30), the latter problem is equivalent to evaluation of the sign of the lowest eigenvalue of the following generalized eigenvalue problem,

$$\mathrm{\Lambda }_1\left(\begin{array}{c}X\\ Y\end{array}\right)=\lambda \left(\begin{array}{c}X\\ Y\end{array}\right)+\beta \left(\begin{array}{c}U\\ V\end{array}\right),$$
(32)
where $`\beta `$ is a constant and $`\mathrm{\Psi }=(X,Y)^T`$ obeys Eq. (30). From the above consideration, we know that $`\mathrm{\Lambda }_1`$ has only one negative eigenvalue, and the eigenfunction cannot be orthogonal to $`(U,V)`$ for $`UV>0`$, since $`f^{(1)}g^{(1)}`$ is positive too. Therefore, the expansion of $`(X,Y)^T`$ over the eigenfunctions of $`\mathrm{\Lambda }_1`$,

$$\mathrm{\Lambda }_1\left(\begin{array}{c}\varphi _1^{(n)}\\ \varphi _2^{(n)}\end{array}\right)=\lambda _n\left(\begin{array}{c}\varphi _1^{(n)}\\ \varphi _2^{(n)}\end{array}\right),$$

has the following form (for simplicity of the presentation, the summation below also denotes the integration over the continuous spectrum)

$$\left(\begin{array}{c}X\\ Y\end{array}\right)=\beta \underset{n,n2}{}\frac{(\varphi _1^{(n)},\varphi _2^{(n)})|\left(\begin{array}{c}U\\ V\end{array}\right)}{\lambda _n\lambda }\left(\begin{array}{c}\varphi _1^{(n)}\\ \varphi _2^{(n)}\end{array}\right),$$
(33)
where the zero eigenvalue $`\lambda _2=0`$ does not enter the sum. In deriving Eq. (33), we have made use of Eq. (32) and the orthogonality of the eigenfunctions of $`\mathrm{\Lambda }_1`$. The constraint (30) requires that the generalized eigenvalue $`\lambda `$ from Eq. (32) satisfies

$$\beta \underset{n,n2}{}\frac{(\varphi _1^{(n)},\varphi _2^{(n)})|\left(\begin{array}{c}U\\ V\end{array}\right)^2}{\lambda _n\lambda }\beta F(\lambda )=0.$$
(34)
Consider the function $`F=F(\xi )`$ defined in (34). It is discontinuous at the points $`\xi =\lambda _n`$ and is a monotonically growing function elsewhere. Note that the summation involves one negative eigenvalue $`\lambda _1`$, the next (smallest) eigenvalue $`\lambda _3`$ being positive. Then, since the generalized eigenvalue $`\lambda `$ defined in Eq. (32) is a zero of the function $`F`$, we have $`\mathrm{sgn}\lambda =\mathrm{sgn}\{F(0)\}`$. Therefore, we must demand $`F(0)<0`$ for the soliton stability. Note that

$$F(0)=(U,V)|\mathrm{\Lambda }_1^1|\left(\begin{array}{c}U\\ V\end{array}\right).$$
(35)
On the other hand, the differentiation of Eqs. (4)-(5) with respect to $`\mu `$ yields (note that the operator $`\mathrm{\Lambda }_1`$ is even with respect to the substitution $`xx`$ )

$$\mathrm{\Lambda }_1\left(\begin{array}{c}U/\mu \\ V/\mu \end{array}\right)=\left(\begin{array}{c}U\\ V\end{array}\right),$$
(36)
or

$$\mathrm{\Lambda }_1^1\left(\begin{array}{c}U\\ V\end{array}\right)=\left(\begin{array}{c}U/\mu \\ V/\mu \end{array}\right).$$

Finally, substitution of the latter result into Eq. (35) leads to the VK criterion for the soliton stability:

$$F(0)=(U,V)|\mathrm{\Lambda }_1^1|\left(\begin{array}{c}U\\ V\end{array}\right)=(U,V)|\left(\begin{array}{c}U/\mu \\ V/\mu \end{array}\right)=\frac{1}{2}\frac{N}{\mu }<0.$$
(37)
It should be pointed out that, in the course of the derivation of the stability criterion (37) for the in-phase soliton solutions of Eqs. (4)-(5), we used the following two facts about the operator $`\mathrm{\Lambda }_1`$: ($`i`$) there is only one (non-degenerate) negative eigenvalue, and ($`ii`$) the corresponding eigenfunction is not orthogonal to the soliton solution $`(U,V)`$ (otherwise the lowest generalized eigenvalue $`\lambda `$ coincides with the negative eigenvalue of $`\mathrm{\Lambda }_1`$). With these two conditions satisfied, one can use the VK stability criterion for other in-phase soliton solutions of Eqs. (1)-(2). Below, this fact is used for the stability analysis of the sech-type solitons (12), and the solutions bifurcating from the sech-type solitons at $`\mu =\mu _{\mathrm{bif}}`$.

II. Stability analysis of the in-phase solitons for $`a<0`$

First of all, let us consider the special case when $`\mu _0`$ is not an independent parameter, but is the function of $`a`$ considered above, i.e. $`\mu _0=\sqrt{a}1/\sqrt{a}=\mu _0(a)`$, see Eq. (13). To derive the expression for the bifurcation point (16), we note that, in general, at such a point two branches of the soliton solutions $`(U^{(1)},V^{(1)})`$ and $`(U^{(2)},V^{(2)})`$ coincide, hence due to Eq. (36) we have

$$\mathrm{\Lambda }_1\left[\left(\begin{array}{c}U^{(1)}/\mu \\ V^{(1)}/\mu \end{array}\right)\left(\begin{array}{c}U^{(2)}/\mu \\ V^{(2)}/\mu \end{array}\right)\right]=0.$$
(38)
Therefore, the bifurcation point is characterized by the appearance of the second (nodeless) zero mode of the operator $`\mathrm{\Lambda }_1`$, in addition to the zero mode given by Eq. (28). Such a zero mode $`(X_1,X_2)^T`$ can be easily found analytically for the sech-type solitons (12); in this case, it is a solution to the system

$`{\displaystyle \frac{\mathrm{d}^2X_1}{\mathrm{d}x^2}}+\left[\mu +3A_1^2\mathrm{sech}^2\left({\displaystyle \frac{A_1x}{\sqrt{2}}}\right)\right]X_1+X_2=0,`$
$`{\displaystyle \frac{\mathrm{d}^2X_2}{\mathrm{d}x^2}}+\left[\mu \mu _0+3A_1^2\mathrm{sech}^2\left({\displaystyle \frac{A_1x}{\sqrt{2}}}\right)\right]X_2+X_1=0.`$
Solving these equations, we obtain the expression (16) for the bifurcation point $`\mu _{\mathrm{bif}}`$, and the corresponding zero mode

$$\left(\begin{array}{c}X_1\\ X_2\end{array}\right)=\left(\begin{array}{c}1\\ \sqrt{a}\end{array}\right)\mathrm{sech}^2\left(\frac{A_1x}{\sqrt{2}}\right),A_1=\left(\frac{2(1a)}{3\sqrt{a}}\right)^{1/2}.$$
(39)
Now, we aim to show that the sech-type solitons (12) are stable in the interval

$$\mu _{\mathrm{bif}}<\mu <\mu _{\mathrm{max}}^{(1)},$$
(40)
where $`\mu _{\mathrm{max}}^{(1)}=1/\sqrt{a}`$ \[see Eq. (13)\], and unstable otherwise \[we remind that here $`\mu _0=\mu _0(a)`$\]. We need to determine the direction of the shift of the zero eigenvalue for $`\mu \mu _{\mathrm{bif}}`$ as we move along a given branch of the soliton solutions. To this end, we can use the perturbation theory for eigenvalues of linear operators. A shift of the chemical potential from the bifurcation point by a small amount $`ϵ`$, $`\mu =\mu _{\mathrm{bif}}+ϵ`$, results in a perturbation of the eigenfunction (39) and the corresponding eigenvalue $`\stackrel{~}{\lambda }=0+ϵl_1+𝒪(ϵ^2)`$. Substitution of the perturbed eigenfunction,

$$\left(\begin{array}{c}\stackrel{~}{X}_1\\ \stackrel{~}{X}_2\end{array}\right)=\left(\begin{array}{c}X_1\\ X_2\end{array}\right)+ϵ\left(\begin{array}{c}\chi _1\\ \chi _2\end{array}\right)+𝒪(ϵ^2),$$

into the eigenvalue problem

$$\mathrm{\Lambda }_1\left(\begin{array}{c}\stackrel{~}{X}_1\\ \stackrel{~}{X}_2\end{array}\right)=\stackrel{~}{\lambda }\left(\begin{array}{c}\stackrel{~}{X}_1\\ \stackrel{~}{X}_2\end{array}\right),$$
(41)
and keeping only the first-order terms in $`ϵ`$, we derive the following equation for $`l_1`$,

$$\mathrm{\Lambda }_1\left(\begin{array}{c}\chi _1\\ \chi _2\end{array}\right)=l_1\left(\begin{array}{c}X_1\\ X_2\end{array}\right)+\left(\begin{array}{cc}1+6U(U/\mu )& 0\\ 0& 16aV(V/\mu )\end{array}\right)\left(\begin{array}{c}X_1\\ X_2\end{array}\right),$$
(42)
where the soliton solutions $`U`$ and $`V`$ are given by Eq. (12), and the operator $`\mathrm{\Lambda }_1`$, together with $`U`$, $`V`$ and their derivatives, are taken at $`\mu =\mu _{\mathrm{bif}}`$. The left multiplication of Eq. (42) by $`(X_1,X_2)`$ and integration over $`x`$ leads to the following expression for the sign of $`l_1`$:

$$\mathrm{sgn}l_1=\mathrm{sgn}dx\left\{\left(1+6U\frac{U}{\mu }\right)X_1^2+\left(16aV\frac{V}{\mu }\right)X_2^2\right\}.$$
(43)
The expressions (42) and (43) are valid for the two branches of the soliton solutions which collide at the bifurcation point $`\mu _{\mathrm{bif}}`$ (16).

We now apply Eq. (43) to the sech-type solitons of Eq. (12). We have

$$U=A_1\mathrm{sech}y,\frac{U}{\mu }=2\mathrm{s}\mathrm{e}\mathrm{c}\mathrm{h}^2y(1y\mathrm{tanh}y),V=\frac{U}{\sqrt{a}},y\frac{A_1x}{\sqrt{2}},$$

thus

$$\mathrm{sgn}l_1=\mathrm{sgn}dy(16\mathrm{s}\mathrm{e}\mathrm{c}\mathrm{h}^2y[1y\mathrm{tanh}y])\mathrm{sech}^4y>0.$$

Therefore, noticing that $`\mathrm{sgn}\stackrel{~}{\lambda }=\mathrm{sgn}(ϵl_1)`$, we arrive at the following formula for the sign of the eigenvalue as we move along the sech-type branch of the soliton solutions:

$$\mathrm{sgn}\stackrel{~}{\lambda }|_{\mathrm{sech}}=\mathrm{sgn}(\mu \mu _{\mathrm{bif}}).$$
(44)
Thus, we have $`\stackrel{~}{\lambda }>0`$ for $`\mu >\mu _{\mathrm{bif}}`$. As the result, the nodeless eigenfunction from Eq. (41) with the property $`\stackrel{~}{X}_1\stackrel{~}{X}_2<0`$ corresponds to a positive eigenvalue, and $`\mathrm{\Lambda }_1`$ has only one non-degenerate negative eigenvalue for $`\mu _{\mathrm{bif}}<\mu <\mu _{\mathrm{max}}^{(1)}`$.

To complete the proof of the stability of the sech-type solitons in the interval (40), we note that the eigenfunction corresponding to the (only) negative eigenvalue $`\lambda _1`$ of $`\mathrm{\Lambda }_1`$ cannot be orthogonal to the soliton solution $`(U,V)`$ (12). Indeed, such an eigenfunction $`(Z_1,Z_2)^T`$ has no nodes and satisfies the condition $`Z_1Z_2>0`$ at $`\mu =\mu _{\mathrm{bif}}`$, as it is orthogonal to the zero mode (39). Thus, to become orthogonal to the soliton solution $`(U,V)`$, one of the components of this eigenfunction should first pass through zero at some value of $`\mu `$ which is impossible.

For $`\mu _0=\mu _0(a)`$, the stability properties of the in-phase solitons bifurcating from the sech-type soliton solutions at $`\mu =\mu _{\mathrm{bif}}`$ are affected by the turning-point bifurcation mentioned in section 3 (see, for instance, the top panel of figure 7). To find where the VK criterion applies to these solitons, we note that, at the collision bifurcation point $`\mu _{\mathrm{bif}}`$, these soliton solutions share the operator $`\mathrm{\Lambda }_1`$ with the sech-type solitons and at the turning point one of the (two) negative eigenvalues of $`\mathrm{\Lambda }_1`$ passes through zero and becomes positive, as we go from the upper part of the branch to the lower one \[this is due to the fact that the function $`F(\xi )`$ makes a jump from $`\mathrm{}`$ to $`+\mathrm{}`$ at $`\xi =0`$ as we pass the turning point downwards, due to $`N/\mu =2F(0)`$\]. Therefore, as we move along this (i.e., bifurcating) branch from the bifurcation point $`\mu _{\mathrm{bif}}`$, the zero eigenvalue corresponding to the nodeless zero mode (39) of $`\mathrm{\Lambda }_1`$ becomes negative, but at the turning point it passes through zero and becomes positive. Therefore, the VK criterion applies to the lower part of the branch, while the solitons of the upper part are always unstable (there is one negative generalized eigenvalue for them).

To complete the consideration of the special case with $`\mu _0=\mu _0(a)`$, we note that, at the bifurcation point $`\mu =\mu _{\mathrm{bif}}`$, the numbers of atoms for the two branches of the in-phase solitons vs. chemical potential, say $`N_1=N_1(\mu )`$ and $`N_2=N_2(\mu )`$, have equal slopes: $`N_1/\mu =N_2/\mu `$, which is a straightforward corollary of Eq. (36). Indeed, the zero mode (39) of the operator $`\mathrm{\Lambda }_1`$, also defined as

$$\left(\begin{array}{c}X_1\\ X_2\end{array}\right)\left(\begin{array}{c}U^{(1)}/\mu \\ V^{(1)}/\mu \end{array}\right)\left(\begin{array}{c}U^{(2)}/\mu \\ V^{(2)}/\mu \end{array}\right),$$

is orthogonal to the soliton solution $`(U,V)^T`$ taken at the bifurcation point $`\mu _{\mathrm{bif}}`$, as it is the right-hand side of Eq. (36). Thus,

$$\frac{N_1}{\mu }\frac{N_2}{\mu }=dx\left(UX_1+VX_2\right)=0.$$

Geometrically, it means that the in-phase solitons undergo a  tangential (one-sided cusp-like) bifurcation. This fact was already illustrated numerically in the previous section.

The stability properties of the soliton solutions for $`\mu _0\mu _0(a)`$ can be summarized as follows. According to their presentation in Figs. 7, 8 and 9(c)-(d), we will refer to the two branches of the soliton solutions as the thin and the thick branch, respectively, where the latter one corresponds to the sech-type solitons (12) for $`\mu _0=\mu _0(a)`$.

First of all, the VK criterion applies to the solitons belonging to the thick branch for all values of $`\mu `$, thus they are stable when the number of atoms decreases with increase of the chemical potential. Second, to understand the stability of the solitons belonging to the thin branch, we note that, at the turning point, one of the negative eigenvalues of $`\mathrm{\Lambda }_1`$ passes through zero as one goes from the upper part of the branch to the lower part, similar as in the special case $`\mu _0=\mu _0(a)`$. Thus, the solitons belonging to the upper branch are unstable, as the operator $`\mathrm{\Lambda }_1`$ has two negative eigenvalues. The VK criterion applies to the lower branch if the (only) negative eigenvalue of $`\mathrm{\Lambda }_1`$ there has the eigenfunction which is non-orthogonal to the soliton solution proper. It is enough to verify this condition numerically just at one point. The outcome is that the VK criterion is indeed applicable to the lower part of the thin branch.

The analytical results on the soliton stability were checked by direct numerical solution of the linear eigenvalue problem (24), and by counting the negative eigenvalues of the operator $`\mathrm{\Lambda }_1`$ (we used the Fourier spectral discretization method with up to 256 grid points in the LAPACK routines of Matlab). As the outcome, it has been verified that the analytical results completely agree with numerical ones in all the cases.

### 4.2 Evolution of perturbed unstable solitons

Here we report results of direct numerical simulations of Eqs. (1)-(2) with perturbed soliton solutions as the initial conditions. We have used the Fourier spectral discretization method in the $`x`$ coordinate combined with the leap-frog time-stepping scheme. The stability of the scheme itself was guaranteed by selecting a small enough time step (we had $`\mathrm{\Delta }t=0.001`$), such that the stability domain of the leap-frog time-stepping method contains the eigenvalues of the spacial discretization operator multiplied by $`\mathrm{\Delta }t`$. The radiation was absorbed by introducing a smooth distributed damping (smooth damping is necessary for stability of a numerical scheme based on the spectral methods), which had a negligible effect on the localized solution.

Previously, the results of numerical simulations for the case of negative $`a`$ in the context of the nonlinear dual-fiber optical model were reported in Ref. , where $`a=1`$ and $`\mu _0=0`$ were used. Here, our main interest is the evolution of unstable solitons in the most interesting case, when $`a0`$, corresponding to nonlinearities of the opposite signs in the two condensates.

First of all, we note that the in-phase solitons may have only one unstable mode, i.e., the linear eigenvalue problem (24) can have only one negative eigenvalue $`\mathrm{\Omega }^2`$ (or equivalently, one imaginary eigenfrequency $`\mathrm{\Omega }`$). For $`a=1`$ and $`\mu _0=0`$, similar result was reported in Ref. .

In the context of the stability analysis presented above, it is easy to understand why there is just one unstable mode. Indeed, where the VK criterion applies, and the solitons become unstable due to the positive slope, $`\mathrm{d}N/\mathrm{d}\mu >0`$, this means that one generalized eigenvalue of the problem (32) has passed zero and become negative with the change of sign of $`F(0)`$, see Eq. (37). For the soliton solutions corresponding to the top part of the thin branches in Figs. 7-9 the VK criterion does not apply, but note, however, that the slope is negative. Thus, the only possible unstable mode is due to the appearance of the second negative eigenvalue in the spectrum of $`\mathrm{\Lambda }_1`$ and, hence, a single negative generalized eigenvalue of the linear problem (32). The fate of the corresponding (second) negative eigenvalue of $`\mathrm{\Lambda }_1`$ in this case depends on whether there is a collision bifurcation or not. In the former case, the negative eigenvalue must go to zero as we approach the bifurcation point $`\mu _{\mathrm{bif}}`$, see Fig. 7, since at $`\mu =\mu _{\mathrm{bif}}`$ this eigenvalue is zero (the stability criterion applies to the branch of the sech-type solitons on one side from the bifurcation point). In the latter case, when there is no collision bifurcation, the negative eigenvalue tends to $`\mathrm{}`$ as we move along the top part of the thin branch away from the turning point, see for instance, the top part of figure 8. In both cases, this negative eigenvalue also goes to zero as we approach the turning point (i.e., as we move in the opposite direction).

The fact that there is just one unstable mode, i.e., one unstable direction for the growth of a small perturbation, simplifies the task of understanding the fate of the unstable soliton solutions subject to perturbations: it is sufficient to solve the initial-value problem for Eqs. (1)-(2), taking as the initial condition the soliton plus a perturbation proportional to its unstable mode.

Here we note that though the perturbation based on the unstable mode is complex (consult appendix B for more details), it is sufficient to use just its real part. This is due to the fact that the soliton solutions are given by the real functions and, hence, only the real part of the perturbation, e.g., proportional to $`(X_u,X_v)^T`$ from equation (24), affects the number of atoms in the first order approximation. Indeed, for a perturbation $`uU=ϵ(X_u+iY_u)`$ and $`vV=ϵ(X_v+iY_v)`$, where $`U`$ and $`V`$ are the soliton components, the corresponding variation of the number of atoms reads

$$\mathrm{\Delta }N=2ϵ\underset{\mathrm{}}{\overset{\mathrm{}}{}}dx\left(UX_u+VX_v\right)+𝒪(ϵ^2).$$
(45)
Below by the small perturbation (or the unstable mode) we mean the real part.

We now focus on the case of $`a0`$. We found that the evolution of an unstable soliton is determined by the shape of the stability curve $`N=N(\mu )`$ and the sign of the small perturbation in the form of the unstable mode. There are two distinct cases, which correspond to the shapes shown in Figs. 2(c) and the left panel of Fig. 5, respectively.

We first consider the evolution of the unstable solitons in the former case. Figure 10 shows the (unstable) soliton and its unstable mode for this case. We chose $`a=0`$ (no nonlinearity in the $`v`$-subsystem) and $`\mu _0=1`$, but the evolution is similar for other values of these parameters, provided that they give rise to a similar shape of the curve $`N=N(\mu )`$. For example, it is so for $`a`$ and $`\mu _0`$ of Fig. 4 (and in any other case of $`a>0`$, when the $`v`$-condensate has atoms with repulsion, but the shape of $`N(\mu )`$ is similar to the one given by Fig. 2(c)).

Figure 11 illustrates the evolution of the unstable soliton solution, where the top picture indicates the direction of the evolution if one adds (the case denoted by 1) or subtracts (2) the unstable mode, as it was specified above. For instance, in scenario 1 there were, initially, more atoms trapped in the $`u`$-condensate, but this distribution of atoms in the condensates is unstable. For such initial conditions, the corresponding attracting configuration has more atoms in the $`v`$-condensate, as the arrow 1 indicates in the top part of Fig. 11.

The picture in the bottom part of Fig. 11 shows the evolution of the two numbers of atoms, initiated by the instability of the original soliton. As it follows from this picture, in both cases the instability gives rise to a breather-like state with persistent internal vibrations. However, the amplitude of the vibrations is dependent on the slope of the stability curve: for the steeper slope the amplitude is smaller.

Now let us consider what happens to the unstable soliton when the stability curve $`N=N(\mu )`$ has the shape of that shown in the left panel of Fig. 5. In this case, a scenario of the type 2 in terms of Fig. 11 (with formation of a breather) is also observed. However, a scenario of the type 1 does not take place in this case. It is replaced by complete decay of the soliton into radiation, as is shown in Fig. 12, which pertains to the case $`a=0.25`$ and $`\mu _0=1`$, the curve $`N(\mu )`$ being in this case similar to the one in the left panel of Fig. 5.

## 5 Conclusion

We have found novel soliton states for the system of two linearly coupled NLS equations, which describe two Bose-Einstein condensates trapped in an asymmetric double-well potential, in particular, when the atomic interactions have opposite signs in the two condensates. This system realizes an interplay between dispersion (corresponding to the kinetic energy in the condensates), attractive and repulsive nonlinearities, linear coupling, and the chemical-potential difference between the two traps.

We have found stable soliton solutions with almost all atoms ($`96\%`$) being trapped in the condensate with weak repulsive interaction, which is a novel self-trapping phenomenon in a system of two weakly coupled condensates having the scattering lengths of the opposite signs. An unusual bifurcation was found too, when a soliton merges with the zero background, having its amplitude vanishing and width diverging, while the number of atoms trapped in the soliton remains finite.

The stability of solitons was studied in detail, and evolution of the unstable ones was investigated by means of direct numerical simulations. The outcome is that the unstable solitons either give rise to breathers or completely decay into incoherent waves (radiation). Our results for the two-component solitons are also relevant to the study of coherent atomic tunnelling (see, for instance, Refs.), which has attracted a lot of attention in the context of the Bose-Einstein condensation.

New results concerning the two-component solitons and their stability were also obtained for the asymmetric model in which the self-interaction is attractive in both components. This model has direct applications to the dual-core nonlinear optical fibers. Although it was studied in many works, new results for this model have been obtained here.

The investigation of the solitons in linearly coupled quasi-one-dimensional condensates in a double-well potential can be continued in several directions. These include the study of breathers, and collisions between moving solitons. In particular, boundary effects will have to be taken into account for moving solitons, i.e., the form of the longitudinal potential will enter the stage. Moreover, the antisymmetric (out-of-phase) and higher-order (multi-humped) solitons, although unstable, may play an important role in the quantum tunnelling phenomena, since the solution may oscillate about these solitons during a finite time. These aspects of the soliton dynamics will be addressed elsewhere.

## 6 Acknowledgements

The research by V.S. was supported by the Fundação de Amparo à Pesquisa do Estado de São Paulo (FAPESP), Brazil. B.A.M. appreciates hospitality of Instituto de Fisica Teórica at Universidade Estadual Paulista (São Paulo, Brazil).

## Appendix A Derivation of the coupled-mode equations

Let us outline the derivation of the coupled-mode equations (1)-(2). We will use the tilde for the physical variables in the GP equation to distinguish them from the corresponding dimensionless variables of the system (1)-(2). The crucial observation is that the modulus of the order parameter (single-atom wave function) is exponentially small in the barrier region, see figure 1 of section 1. This allows us to approximate the solution to the corresponding GP equation with the double-well potential $`V_{\mathrm{ext}}(\stackrel{}{\rho })`$ strongly confining in the transverse directions $`\stackrel{}{\rho }(y,z)`$ (we neglect the longitudinal variation of $`V_{\mathrm{ext}}`$),

$$i\mathrm{}\frac{\mathrm{\Psi }}{\stackrel{~}{t}}=\frac{\mathrm{}^2}{2m}\left(\frac{^2}{\stackrel{~}{x}^2}+_\stackrel{}{\rho }^2\right)\mathrm{\Psi }+\left(V_{\mathrm{ext}}+g_0|\mathrm{\Psi }|^2\right)\mathrm{\Psi },$$
(46)
as a sum of the factorized wave functions for the two wells:

$$\mathrm{\Psi }(\stackrel{~}{t},\stackrel{~}{x},\stackrel{}{\rho })=\psi _u(\stackrel{~}{t},\stackrel{~}{x})\mathrm{\Phi }_u(\stackrel{}{\rho })+\psi _v(\stackrel{~}{t},\stackrel{~}{x})\mathrm{\Phi }_v(\stackrel{}{\rho }),$$
(47)
where $`\mathrm{\Phi }_u(\stackrel{}{\rho })`$ and $`\mathrm{\Phi }_v(\stackrel{}{\rho })`$ would be the ground-states (for the transversal degrees of freedom) if the two wells were isolated. Here we point that $`g_0`$ is different in the two wells: $`g_0^{(v)}/g_0^{(u)}`$ has arbitrary value and sign. This quotient can be easily managed by technique based on the Feshbach resonance  with the magnetic field applied to only one of the wells. Inserting the ansatz (47) in the GP equation (46) and using the conditions

$$\mathrm{\Phi }_u\mathrm{\Phi }_v\mathrm{d}^2\stackrel{}{\rho }0,\mathrm{\Phi }_u^2\mathrm{d}^2\stackrel{}{\rho }=\mathrm{\Phi }_v^2\mathrm{d}^2\stackrel{}{\rho }=1,$$
(48)
one arrives at the system of two equations for the wave functions describing the longitudinal evolution of the condensates in the wells:

$$i\mathrm{}\frac{\psi _u}{\stackrel{~}{t}}=\frac{\mathrm{}^2}{2m}\frac{^2\psi _u}{\stackrel{~}{x}^2}+(_u+g_u|\psi _u|^2)\psi _uK\psi _v,$$
(49)
$$i\mathrm{}\frac{\psi _v}{\stackrel{~}{t}}=\frac{\mathrm{}^2}{2m}\frac{^2\psi _v}{\stackrel{~}{x}^2}+(_v+g_v|\psi _v|^2)\psi _vK\psi _u.$$
(50)
Here the parameters are defined as follows (by changing the sign of either $`\mathrm{\Phi }_u`$ or $`\mathrm{\Phi }_v`$, if necessary, one can set $`K>0`$)

$$_{u,v}=\left(\frac{\mathrm{}^2}{2m}(_\stackrel{}{\rho }\mathrm{\Phi }_{u,v})^2+V_{\mathrm{ext}}\mathrm{\Phi }_{u,v}^2\right)\mathrm{d}^2\stackrel{}{\rho },g_{u,v}=g_0\mathrm{\Phi }_{u,v}^4\mathrm{d}^2\stackrel{}{\rho },$$

$$K=\left(\frac{\mathrm{}^2}{2m}(_\stackrel{}{\rho }\mathrm{\Phi }_u)(_\stackrel{}{\rho }\mathrm{\Phi }_v)+\mathrm{\Phi }_uV_{\mathrm{ext}}\mathrm{\Phi }_v\right)\mathrm{d}^2\stackrel{}{\rho }.$$
(51)
(Note that though the overlap of the ground states $`\mathrm{\Phi }_u`$ and $`\mathrm{\Phi }_v`$ as in Eq. (48) is negligible, the coupling coefficient $`K`$ is *not*, due to the local maximum of the potential $`V_{\mathrm{ext}}(\stackrel{}{\rho })`$ in the overlap region, see Fig. 1 of section 1). Finally, assuming that $`g_u<0`$ (i.e., attractive atomic interaction in the $`u`$-condensate) and setting

$$\psi _u=\sqrt{\frac{K}{g_u}}\mathrm{exp}\left(\frac{i_ut}{\mathrm{}}\right)u,\psi _v=\sqrt{\frac{K}{g_u}}\mathrm{exp}\left(\frac{i_ut}{\mathrm{}}\right)v,$$
(52)
in Eqs. (49)-(50) one arrives at the dimensionless system (1)-(2) with $`a=g_v/g_u`$ and $`\mu _0=(_v_u)/K`$, where the space and time variables of the system (1)-(2) are given as follows $`x=(\sqrt{2mK}/\mathrm{})\stackrel{~}{x}`$ and $`t=(K/\mathrm{})\stackrel{~}{t}`$.

Recalling the well-known expression for the coupling coefficient $`g_0^{(u)}=4\pi \mathrm{}^2a_s^{(u)}/m`$, where $`a_s^{(u)}`$ is the scattering length in the $`u`$-condensate, and using the definition of $`g_u`$ (51) in the relation (52) we conclude that the actual (i.e., physical) numbers of atoms are related to the numbers of particles in the model (1)-(2) as follows

$$\stackrel{~}{N}_u=\frac{\left[\mathrm{\Phi }_u^4\mathrm{d}^2\stackrel{}{\rho }\right]^1}{8\pi d_0|a_s^{(u)}|}dx|u|^2,\stackrel{~}{N}_v=\frac{\left[\mathrm{\Phi }_u^4\mathrm{d}^2\stackrel{}{\rho }\right]^1}{8\pi d_0|a_s^{(u)}|}dx|v|^2,$$
(53)
where we have introduced the characteristic overlap distance by $`d_0\mathrm{}/\sqrt{2mK}`$.

In the above derivation it was assumed that the transverse profile of the order parameter is determined by the linear part of the r.h.s. in equation (46). Such approximation is valid for the weak atomic interaction:

$$\frac{\mathrm{max}(|g_0^{(u)}|,|g_0^{(v)}|)\stackrel{~}{N}}{d_{}^2d_{}}\frac{\mathrm{}\omega _{}}{2},$$
(54)
where $`d_{}=(\mathrm{}/m\omega _{})^{1/2}`$ and $`d_{}`$ are the transverse and longitudinal sizes of the condensate, with the transverse size being given by the oscillator length, $`\stackrel{~}{N}`$ is the (physical) total number of atoms, and $`\omega _{}`$ is the external trap frequency (in the transverse dimensions) without the separation barrier. The r.h.s. of formula (54) is an estimate on the energy of the linear transverse terms in the GP equation (46). Substitution of the expression for the coupling coefficients $`g_0^{(u,v)}=4\pi \mathrm{}^2a_s^{(u,v)}/m`$ in the condition (54) leads to a more convenient equivalent form:

$$8\pi \mathrm{max}(|a_s^{(u)}|,|a_s^{(v)}|)\frac{\stackrel{~}{N}}{d_{}}1.$$
(55)
To make an estimate we use the usual values: $`a_s10^9`$m and $`d_{}10^4`$m (the latter corresponds to the asymmetry parameter $`\gamma d_{}/d_{}0.01`$ for $`d_{}10^6`$m, usual for experiments on the cigar-shaped condensates). Thus $`\stackrel{~}{N}10^6`$ thereby allowing for up to a hundred thousands of atoms in the condensates.

Due to the scaling $`d_{}N`$, where $`N`$ is the total number of particles in the model (1)-(2) (compare the widths of the solitons from figs. 3 and 4-5), the condition (55) is satisfied for all values of the parameters $`\mu _0`$ and $`a`$ in the system (1)-(2) and all values of the chemical potential $`\mu `$ if it is satisfied just at one point $`(a,\mu _0,\mu )`$. Therefore, though the total numbers of atoms do vary by an order of magnitude, nevertheless, the model system (1)-(2), in fact, applies uniformly in the whole parameter space.

A very rough estimate on the numbers of trapped atoms can be easily carried out. Note that the coupling coefficient $`g_u`$ depends on the geometry of the trap through the multiplier

$$\mathrm{\Phi }_u^4\mathrm{d}^2\stackrel{}{\rho }\frac{1}{(d_{}^{(u)})^2}.$$

Here $`d_{}^{(u)}`$ is the characteristic transverse size of the $`u`$-condensate. To make an estimate on the numbers of trapped atoms we assume: $`|a_s^{(u)}|10^9`$m, $`d_{}^{(u)}d_{}=10^6`$m, and $`d_00.1d_{}^{(u)}=10^7`$m, what leads to the transformation coefficients in formula (53) of the order $`10^3`$. Thus this estimate shows that all the stable solitons we discuss in the present paper correspond to numbers of trapped atoms lying below the estimated bound of $`10^6`$.

## Appendix B Unstable mode and the related perturbation of the soliton solution

Here we provide some details on the unstable eigenfunction of the linear stability problem and the related initial perturbation. Looking for eigenmodes of the linear system (22) we write

$$\left(\begin{array}{c}𝒰_R\\ 𝒰_I\\ 𝒱_R\\ 𝒱_{}\end{array}\right)=e^{i\mathrm{\Omega }t}\left(\begin{array}{c}X_u\\ Y_u\\ X_v\\ Y_v\end{array}\right)+c.c.,$$
(56)
since we need a real solution. The linear problem (22) then reduces to an eigenvalue problem for $`\lambda =i\mathrm{\Omega }`$, which can be written as a system of two matrix equations:

$$i\mathrm{\Omega }\left(\begin{array}{c}X_u\\ X_v\end{array}\right)=\mathrm{\Lambda }_0\left(\begin{array}{c}Y_u\\ Y_v\end{array}\right),i\mathrm{\Omega }\left(\begin{array}{c}Y_u\\ Y_v\end{array}\right)=\mathrm{\Lambda }_1\left(\begin{array}{c}X_u\\ X_v\end{array}\right),$$
(57)
from which one derives equation (24). Note that $`(X_u,X_v)^T`$ is real (as the eigenfunction of the eigenvalue problem (24)). Then from the second equation in (57) we conclude that the eigenfunction $`(X_u,Y_u,X_v,Y_v)^T`$ corresponding to the real eigenvalue (i.e. imaginary eigenfrequency $`\mathrm{\Omega }`$) is real. Taking this into account and using equations (56)-(57) it is straightforward to deduce that the initial perturbation $`(u_1,v_1)^T`$ of the soliton solution which is based on the unstable eigenfunction $`(X_u,Y_u,X_v,Y_v)^T`$ is given as

$$\left(\begin{array}{c}u_1\\ v_1\end{array}\right)=\sigma \left(1\frac{i\mathrm{\Lambda }_1}{\mathrm{Im}\mathrm{\Omega }}\right)\left(\begin{array}{c}X_u\\ X_v\end{array}\right),$$
(58)
where $`\sigma `$ is arbitrary real constant.

## Appendix C Negative eigenvalues of the operator $`\mathrm{\Lambda }_1`$ for $`a0`$

To prove that $`\mathrm{\Lambda }_1`$ has only one nodeless eigenfunction (and, hence, only one negative eigenvalue) for $`a0`$, we use the fact that the in-phase soliton solutions correspond to the $`k_1`$-branch of the dispersion relation (11). The eigenfunction $`(X_1,X_2)^T`$ of the operator $`\mathrm{\Lambda }_1`$ corresponding to an eigenvalue $`\lambda `$ satisfies the following system of equations,

$`{\displaystyle \frac{\mathrm{d}^2X_1}{\mathrm{d}x^2}}+\left(\mu +\lambda +3U^2\right)X_1+X_2=0,`$ (59)
$`{\displaystyle \frac{\mathrm{d}^2X_2}{\mathrm{d}x^2}}+\left(\mu +\lambda \mu _03aV^2\right)X_2+X_1=0,`$ (60)
which follow from the definition of $`\mathrm{\Lambda }_1`$. We will prove that, among two possible combinations of nodeless eigenfunctions, $`X_1X_2>0`$ and $`X_1X_2<0`$, the latter one is not possible for the in-phase solitons when $`a0`$.

The multiplication of Eq. (60) by $`\mathrm{d}X_2/\mathrm{d}x`$ and integration from $`x=0`$ to $`x=\mathrm{}`$ give

$$(\mu \mu _0+\lambda )X_2^2(0)=\underset{0}{\overset{\mathrm{}}{}}dx\left\{6aV^2X_2\frac{\mathrm{d}X_2}{\mathrm{d}x}+2X_1\frac{\mathrm{d}X_2}{\mathrm{d}x}\right\}.$$
(61)
Here we have used that the eigenfunction has no nodes, hence $`\mathrm{d}X_2/\mathrm{d}x=0`$ at $`x=0`$. Taking into account that the condition $`k_1^2>0`$ (see equation (11)) implies $`\mu \mu _0<0`$, and that $`X_2(\mathrm{d}X_2/\mathrm{d}x)<0`$ for $`x>0`$, we conclude that Eq. (61) cannot be satisfied for a negative eigenvalue $`\lambda `$ if $`X_1(\mathrm{d}X_2/\mathrm{d}x)>0`$ for $`x>0`$, i.e., the operator $`\mathrm{\Lambda }_1`$ cannot have an eigenfunction corresponding to a negative eigenvalue and, simultaneously, having the property $`X_1X_2<0`$.