Stability of single and multiple matter-wave dark solitons in collisionally inhomogeneous Bose-Einstein condensates

TL;DR

This paper investigates the spectral stability of single and multiple matter-wave dark solitons in Bose-Einstein condensates with periodically modulated scattering length, revealing new instability phenomena not present in homogeneous cases.

Contribution

It provides a combined analytical and numerical study of dark soliton stability in inhomogeneous BECs, highlighting novel instability mechanisms caused by scattering length modulation.

Findings

Oscillatory instability of single dark solitons due to modulation

Absence of instabilities near the linear limit for two dark solitons

Exponential instability in three dark solitons with in-phase motion

Abstract

We examine the spectral properties of single and multiple matter-wave dark solitons in Bose-Einstein condensates confined in parabolic traps, where the scattering length is periodically modulated. In addition to the large-density limit picture previously established for homogeneous nonlinearities, we explore a perturbative analysis in the vicinity of the linear limit, which provides good agreement with the observed spectral modes. Between these two analytically tractable limits, we use numerical computations to fill in the relevant intermediate regime. We find that the scattering length modulation can cause a variety of features absent for homogeneous nonlinearities. Among them, we note the potential oscillatory instability even of the single dark soliton, the potential absence of instabilities in the immediate vicinity of the linear limit for two dark solitons, and the existence of an exponential instability associated with the in-phase motion of three dark solitons.

Equations

$$i{U}_t=-\frac{1}{2}{U}_{xx}+V(x)U+g(x)|U{|}^{2}U,$$

$$V(x)=\frac{1}{2}{\Omega }^2{x}^2,$$

$$g(x)=1+g_0{\sin }^2(kx).$$

$$U\propto H_n(\sqrt{\Omega }x)\exp (-\Omega x^2/2)\exp (-i{E}_{n}t),$$

$$u=\sqrt{ϵ}{u}_0+{ϵ}^{3/2}{u}_1+\dots ,\mu ={\mu }_0+ϵ{\mu }_1+\dots .$$

$${\mu }_1=\int |u_0{|}^{4}dx.$$

$$U(x,t)=e^{-i\mu t}\left[u(x)+a(x)e^{\lambda t}+b^{\star }(x)e^{{\lambda }^{\star }t}\right],$$

$$M=\left(\begin{array}{cc}{M}_{11}& M_{12}\\ M_{21}& M_{22}\end{array}\right),$$

$$\begin{array}[]{rcl}M_{11}&=&\left(-\frac{1}{2}\frac{d^{2}}{dx^{2}}+V+2g(x)|u|^{2}-\mu\right),\\[8.61108pt] M_{12}&=&g(x)u^{2},\\[4.30554pt] M_{21}&=&-g(x)u^{\star 2},\\[4.30554pt] M_{22}&=&-\left(-\frac{1}{2}\frac{d^{2}}{dx^{2}}+V+2g(x)|u|^{2}-\mu\right).\end{array}$$

$$Mv=({\mathcal{H}}_0+ϵ{\mathcal{H}}_1)v=\omega v,$$

$$v$$

$${\mathcal{H}}_0$$

$$\displaystyle{\cal H}_{1}=\left(\begin{array}[]{cc}2g(x)|u_{0}|^{2}-\mu_{1}&g(x)u_{0}^{2}\\[4.30554pt] -g(x)(u_{0}^{2})^{\star}&\mu_{1}-2g(x)|u_{0}|^{2}\end{array}\right).$$

$${\mu }_0{u}_1+{\mu }_1{u}_0=\mathcal{L}{u}_1+g(x)u_0^3$$

$$u_1=\sum _{n\ne m}\frac{⟨v_n,F(x)⟩}{{\omega }_n-{\mu }_0}{v}_n,$$

$${\mathcal{M}}_{ij}=⟨W_i|{\mathcal{H}}_1|W_j⟩,$$

$$\omega =(n-m)\Omega +ϵ{\omega }^{\mathrm{cor}}.$$

$$\omega =\Omega -\frac{1}{6}\left(\mu -\frac{3}{2}\Omega \right).$$

$$\omega =2\Omega -\frac{1}{12}\left(\mu -\frac{3}{2}\Omega \right),$$

$$\omega =3\Omega -\frac{7}{32}\left(\mu -\frac{3}{2}\Omega \right).$$

$${\omega }_n=\sqrt{\frac{n(n+1)}{2}}\Omega .$$

$$\omega =\Omega -\frac{5}{41}\left(\mu -\frac{5}{2}\Omega \right).$$

$$\omega =2\Omega +\left(\mu -\frac{5}{2}\Omega \right)\frac{-55\pm 3\sqrt{23}i}{656}.$$

$$\omega =3\Omega -\frac{87}{656}\left(\mu -\frac{5}{2}\Omega \right).$$

$$\omega =\Omega -\frac{19}{156}\left(\mu -\frac{7}{2}\Omega \right)$$

$$\omega =2\Omega +\left(\mu -\frac{7}{2}\Omega \right)\frac{-163\pm 7\sqrt{71}i}{2352}.$$

$$\omega =3\Omega +\left(\mu -\frac{7}{2}\Omega \right)\frac{-769\pm \sqrt{124031}i}{16384}.$$

$${\mu }_1=\frac{3(2+g_0){\Omega }^2-e^{-\frac{k^2}{2\Omega }}{g}_0(k^4+6k^2\Omega +{\Omega }^2)}{8\sqrt{2\pi }{\Omega }^{3/2}}.$$

$$\omega =\Omega +ϵ{\omega }_1^{\mathrm{cor}}$$

$${\omega }_1^{\mathrm{cor}}=\frac{{\Omega }^{1/2}}{32\sqrt{2\pi }}\left(A\pm \sqrt{B}\right),$$

$$A$$

Full text

I Introduction

Over the past 20 years, the physics of atomic Bose-Einstein condensates (BECs) has enabled the examination of numerous physical concepts book1 ; book2 (see also Ref. dumitr1 for a recent review on the progress in this setting). One of the key directions that have been explored lies at the interface between nonlinear wave dynamics and such atomic (as well as optical) systems, concerning the study of so-called matter-wave solitons emergent ; book_new . These coherent structures have been not only theoretically predicted but also in numerous cases experimentally verified. Some of the most notable relevant examples are bright expb1 ; expb2 ; expb3 , dark djf and gap gap matter-wave solitons. Higher dimensional analogues of these structures have also been studied including vortices fetter1 ; fetter2 , solitonic vortices and vortex rings komineas_rev .

In the one-dimensional (1D), single-component repulsive BEC setting, arguably, the most prototypical nonlinear excitations are dark solitons. These states have received considerable attention, perhaps in large measure due to their accessibility in a diverse array of experiments featuring a wide variety of methods. These have led to their (progressively more) well-controlled creation chap01:denschl ; chap01:dark1 ; ourmarkus2 ; ourmarkus3 ; seng1 ; seng2 . Additionally, numerous works have examined dark solitons in higher-dimensional settings observing experimentally their instability in the latter settings. This, in turn, leads to the formation of vortices, vortex rings and vortex lines, as illustrated, e.g., in Refs. chap01:dark ; engels ; jeff ; rcg:90 ; pitasv . One of the important associated recent developments has been the realization of such nonlinear excitations in fermionic superfluids martinz ; fprl . A relatively comprehensive summary of the relevant activity, encompassing both experiment and theory —although not including some of the most recent developments— can be found in Ref. djf and also in Chap. 2 of Ref. book_new .

A context that has recently received considerable interest as regards the dynamics of solitary waves is that of spatially inhomogeneous nonlinearities. A review exploring diverse aspects of this topic can be found in Ref. borisreview . In BECs, the nonlinearity stems from interatomic interactions; in the mean-field approach, where BECs are described by a macroscopic BEC wavefunction obeying the Gross-Pitaevskii (GP) equation, the nonlinearity coefficient is proportional to the $s$-wave scattering length book1 ; book2 . The fact that the sign and magnitude of the interatomic interactions can be controlled using Feshbach resonances Koehler ; feshbachNa ; ofr , has led to a wide array of theoretical theor and experimental possibilities, including the realization of matter-wave bright solitons expb1 ; expb2 ; expb3 or the revelation of the BEC-BCS crossover exp . For instance, Feshbach resonances have been experimentally used to induce spatial inhomogeneities in the scattering length of Yb BECs Tak . Such collisional inhomogeneities, lead to an effective nonlinear potential, in addition to the customary external (e.g., parabolic) potential. As a result, this may lead to features absent in spatially uniform condensates NLPots ; Chang ; borisreview . Some such examples include, but are not limited to, adiabatic compression of matter waves our1 , enhancement of the transmission of matter waves through barriers our2 , dynamical trapping of solitary waves our2 , delocalization transitions of matter waves LocDeloc , emission of atomic solitons in4 ; fot_scat , generation of solitons 21 and vortex rings 22 , control of Faraday waves alex , and others. These types of nonlinear potentials have also led to interesting insights in the context of photonic structures kominis .

Our aim in the present work is to explore states involving one or more dark solitons in both collisionally homogeneous and collisionally inhomogeneous BECs. In the former setting, such states have already been explored, e.g., in Ref. ourmarkus3 and their spectral analysis in the limit where the solitary waves can be considered as particles was theoretically studied in Ref. coles . However, even in that setting we offer a twist by examining the opposite (near-linear) limit by virtue of the perturbative analysis developed in Ref. mprizola . The two limits, when combined, constitute – in our view – the full set of parameter ranges that can be analytically studied, and the intermediate parameter range between the two is supplemented by means of numerical computations. We also extend both our analytical methodology and our numerical investigations to the case of collisionally inhomogeneous condensates. There we find a rich variety of features that are distinct from the homogeneous case. For instance, the single dark soliton, which is generically stable in 1D BECs, can become subject to a Hopf bifurcation and an oscillatory instability in the inhomogeneous setting. The two-soliton state, on the other hand, is immediately subject to this type of instability beyond the linear limit for constant nonlinearity. The associated degeneracy can be broken for inhomogeneous nonlinearities, and this instability may be absent in the immediate vicinity of the limit. Finally, for three dark solitons, an instability that did not exist in the homogeneous limit arises featuring exponential in-phase divergence of the three solitons from their equilibrium solution. In this case too, the theoretical analysis is supplemented by numerical bifurcation results. We also use direct numerical simulations in order to showcase some of the above instabilities.

Our presentation is structured as follows. First, in Sec. II, we introduce the model and present our analytical approach. Next, in the first part of Sec. III, we provide the numerical results for the homogeneous case, while in the second part, we consider the inhomogeneous case. Finally, in Sec. IV, we summarize our findings and present a number of directions for future work.

II Model and Computational Setup

Our model of choice will be the quasi-1D Gross-Pitaevskii (GP) equation of the form book1 ; book2 ; book_new :

[TABLE]

where $U(x,t)$ is the macroscopic BEC wavefunction, $V(x)$ is the external potential, and $g(x)$ is the, potentially, spatially-dependent nonlinearity coefficient. While improved models have been devised to more adequately capture the effect of BEC dimensionality – and, particularly, the interplay between longitudinal and transverse directions of the BEC (see, e.g., Refs. salasnich ; delgado ) – the methods we propose here would be equally applicable to the latter variants. Nevertheless, for simplicity and clarity of the exposition (and of the analytical results provided), as well as to connect to earlier literature on the subject, we opt to present the results in the simpler GP setting.

The linear external potential will be assumed to be parabolic, in the functional form

[TABLE]

with $\Omega$ being the normalized trap strength, while the collisional inhomogeneity, when present, will be assumed to have a typical periodic modulation around a constant value (see, e.g., Ref. borisreview ), namely:

[TABLE]

In principle, for $g_{0}$ sufficiently negative, $g$ may even change sign, yet here we will restrict our considerations to cases where this does not occur. The linear limit where the wavefunction $U\rightarrow 0$ amounts to the quantum harmonic oscillator with eigenfunctions

[TABLE]

and corresponding eigenenergies $E_{n}=(n+1/2)\Omega$, where $H_{n}$ is the $n$-th Hermite polynomial.

We now focus on stationary states of the nonlinear problem in the form: $U(x,t)=e^{-i\mu t}u(x)$. It is well known that such states bifurcate from the corresponding linear states kat ; zezy . Following the methodology of Ref. mprizola , we can use a series expansion of the solution in the vicinity of this linear limit, in the form:

[TABLE]

where $(\mu_{0},u_{0})$ correspond, respectively, to the eigenvalue and eigenfunction of a state at the linear limit. As a result, we find from Eq. (1), at O$(\epsilon)$, the solvability condition:

[TABLE]

From this formula, we can specify $\epsilon=(\mu-\mu_{0})/\mu_{1}$ (to leading order).

The next step is to consider the spectral stability i.e., the so-called Bogolyubov-de Gennes (BdG) analysis associated with the linearization around a stationary state. We use the ansatz

[TABLE]

where $\lambda$ denotes the relevant eigenvalue, and $(a,b)^{T}$ is its eigenfunction (star denotes complex conjugate). This leads to the eigenfrequency matrix:

[TABLE]

with

[TABLE]

Here, the eigenfrequency $\omega$ is connected to the eigenvalue $\lambda$ through $\lambda=i\omega$.

Now, using the expansion in powers of $\epsilon$ within the stability matrix, we obtain

[TABLE]

with

[TABLE]

where ${\cal L}=-\frac{1}{2}\frac{d^{2}}{dx^{2}}+V(r)$, while

[TABLE]

A nonzero imaginary part of $\omega$ (or, equivalently, a nonzero real part of $\lambda$) in this Hamiltonian system signals the presence of a dynamical instability.

We should note in passing that for the existence problem the leading-order correction $u_{1}$ satisfies

[TABLE]

Decomposing $u_{1}$ into modes of the quantum harmonic oscillator $u_{1}=\sum_{n\neq m}a_{n}v_{n}$ with frequency (energy) $\omega_{n}$, and considering the state of interest $u_{0}$ to be a multi-dark soliton state composed by $m$ dark solitons, we obtain

[TABLE]

where $\langle f,g\rangle=\int_{-\infty}^{\infty}{f(x)g(x)\,dx}$. While the denominator acquires a particularly simple form given that $\omega_{n}=(n+1/2)\Omega$ and $\mu_{0}=(m+1/2)\Omega$, hence their difference is $(n-m)\Omega$, we will not pursue the existence problem of Eq. (20) at higher order further. We will instead focus on the stability problem of Eqs. (12)–(18) where this difference also appears, characterizing the eigenvalues $(\pm(n-m)\Omega)$ of the stability problem associated with the operator ${\cal H}_{0}$. In the case where $n>m$, these are referred to as positive energy eigenvalues, while if $n<m$, they are referred to as negative energy (or anomalous) eigenvalues/eigenmodes.

To identify the dependence of the eigenvalues on the chemical potential parameter $\mu$ that we consider in our continuations, one uses the degenerate perturbation theory of Ref. mprizola , constructing the matrix ${\cal M}$ with elements

[TABLE]

for all the pairs of $i,j$ (and eigenmodes $W_{i},W_{j}$) which correspond to degenerate eigenvalues in the linear limit. The eigenfrequencies of the matrix ${\cal M}$ will yield the corrections $\omega^{{\rm cor}}$ to the linear limit of the stability problem according to:

[TABLE]

Recall that here $n$ indexes the $n$-th eigenmode, while $m$ is the index of the state of interest, i.e., $m=1$ for a single dark soliton, $m=2$ for a double dark soliton, $m=3$ for a triple dark soliton, and so on.

Having laid out this theoretical formulation, we are now ready to obtain specific results for the two cases at hand, namely the collisionally homogeneous case of $g(x)=1$ and the collisionally inhomogeneous, sinusoidally varying $g(x)$.

III Numerical Results and Comparison to Theory

III.1 Collisionally Homogeneous Case

We start with the case of the single dark soliton, which it is well known to be spectrally stable at the level of the GP equation book_new ; djf . Since $m=1$, there will be two modes of the linearization that will be degenerate at $\omega=\Omega$, a positive energy one, with $n=2$, and a negative energy one, with $n=0$. However, one of these modes (the $n=2$ one) is associated with the so-called dipolar or Kohn oscillation book1 ; book2 and hence remains invariant under variations of $\mu$. The other one (the anomalous mode) becomes associated with the oscillation frequency of the dark soliton in the parabolic trap, predicted in the asymptotic limit of large chemical potential to be $\omega=\Omega/\sqrt{2}$ (this prediction was first reported in Ref. busch and was subsequently confirmed by numerous additional works – cf. discussion in Refs. book_new ; djf and references therein). Our theoretical analysis near the linear limit leads to $\mu_{1}=3\Omega^{1/2}/(4\sqrt{2\pi})$ (and from this, we express $\epsilon=(\mu-3\Omega/2)/\mu_{1}$ for this branch). Then, using the degenerate perturbation theory, we obtain that the frequency associated with the anomalous mode near the linear limit is:

[TABLE]

Importantly, all higher-order modes in this case are non-degenerate and, hence, the corresponding eigenfrequency corrections can be computed as scalars with $i=j=1$ in Eq. (21). Following this path, we obtain

[TABLE]

for the case of the mode with $n=2$, while finally for the mode with $n=3$, we have:

[TABLE]

The large chemical potential case also enables analytical consideration. In this setting, the dark solitons become narrow as their characteristic scale, the healing length, gets smaller, and can be assumed as having no internal dynamics over that scale. However, they still possess their own internal modes (the anomalous or negative energy modes) that pertain to the in-trap soliton motion. The corresponding eigenvectors encompass evolution of the coherent structures, but not that of the background. In this large chemical potential, so-called Thomas-Fermi (TF) limit, the spectrum consists of two separate ingredients. The above mentioned internal modes of the solitons and those of the background, i.e., the ground state on top of which the solitary waves are “supported”. The latter modes have been characterized since the early work of Stringari sandro2 (see Ref. pelirecent for a recent discussion in different dimensionalities) as pertaining, through a suitable transformation, to Legendre differential equation and being eventually given by

[TABLE]

Hence, for a single dark soliton all the TF limit modes consist of those of Eq. (26), for non-negative integers $n$, and the $\Omega/\sqrt{2}$ mode pertaining to the single soliton in-trap oscillation. The result of Eq. (26) is obtained directly by using the (approximate) TF solution profile $u_{0}=[\max(\mu-V(x),0)]^{1/2}$ in the linearization equations that get considerably simplified in this case and, in the limit, approach the Legendre differential equation.

Both sets of predictions, namely those of Eqs. (23)–(25), as well as those of Eq. (26) (and the soliton internal mode) are shown in Fig. 1. The former modes can be seen to provide a good description for small chemical potentials near the limit of $3\Omega/2$, while the latter provide the proper asymptotic limit for large $\mu$. The numerical computations interpolate between these two asymptotic, analytically tractable limits, providing the full spectral picture of the BdG analysis. It is interesting to mention here that in the limit of large $\mu$, we systematically observe that the convergence of the numerical BdG results (solid blue line) to the analytical ones of Eq. (26) (dashed orange line) occurs “sooner”, i.e., for lower values of $\mu$, for lower modes. This trait will also be discernible in the multiple soliton cases that follow.

In the case of two dark solitons, the picture is fairly similar, however bearing the following differences. Since now we are dealing with the $m=2$ mode, bifurcating out of $\mu_{0}=5\Omega/2$, it is the modes with $n=3$ (positive energy) and $n=1$ (negative energy) that will be resonant at $\omega=\Omega$. Again, here the $n=3$ mode corresponds to the dipolar oscillation leading to an invariant eigenfrequency in the stability analysis. However, the $n=1$ mode represents the lowest vibrational mode that can be tracked via our degenerate perturbation theory, which in this case predicts:

[TABLE]

A similar degeneracy, for $m=2$, can be diagnosed for the frequency $2\Omega$. This is due to the fact that the modes with $n=4$ (positive energy) and $n=0$ (negative energy) are degenerate. In this case, the degenerate perturbation theory is necessary for these two modes, resulting in a complex eigenfrequency of the form:

[TABLE]

Finally, in the case of the mode starting from frequency $3\Omega$, there is no degeneracy and the relevant mode with $n=5$ can be found to have the frequency:

[TABLE]

In the large chemical potential limit, the two-soliton case naturally still bears all the modes associated with the ground state, as per Eq. (26). However, now there are two solitonic modes, the lowest one still given by $\Omega/\sqrt{2}$ and associated with the in-phase vibration of the two dark solitons (tantamount to the oscillation of a single one inside the trap, hence bearing the same frequency). The higher one among the two modes is associated with the out-of-phase motion of the two solitons, and is provided by Eq. (21) in Ref. coles , illustrating both the single and the double logarithmic dependence on the chemical potential $\mu$.

The numerical computations associated with this two-soliton state are illustrated in Fig. 2. While the large density limit has been explored before (e.g., in Ref. coles ), and its accuracy is perhaps expected, it is relevant to highlight the success of the degenerate perturbation theory near the linear limit. The latter captures the decreasing dependence of all of the first few modes as the chemical potential increases. Moreover, it accurately tracks down not only the motion of the first and third mode towards decreasing frequencies, but also the instability caused by the degenerate second mode at $2\Omega$. In fact, it provides a very good quantitative handle of the growth rate of the associated instability, as shown in the right panel of Fig. 2. Once again, the numerical results interpolate between the analytically tractable limits in this two-soliton case.

Finally, in Fig. 3, we consider the case example of a three dark soliton state. Now all three of the first frequencies at $\Omega$, $2\Omega$ and $3\Omega$ are resonant. The first pair however, associated with $n=4$ (positive energy) and $n=2$ (negative energy) for the case of $m=3$ does not lead to instability, as again one of these modes is the dipolar one ($n=4$) and remains invariant. The anomalous mode, on the other hand, moves according to:

[TABLE]

and decreases in frequency approximating in the TF limit, once again, $\Omega/\sqrt{2}$ being associated with the in-phase motion of all three dark solitons. The other two pairs indeed do lead to instabilities near the linear limit. Indeed, the modes with $n=5$ (positive energy) and $n=1$ (negative energy) lead to a resonance, resolved by degenerate perturbation theory, according to the result:

[TABLE]

Finally, the modes with $n=6$ (positive energy) and $n=0$ (negative energy) will be degenerate in this case of $m=3$ at $3\Omega$, and can be captured via the degenerate perturbation methodology as:

[TABLE]

We note that for higher chemical potential, these complex eigenfrequency quartets split into two pairs. In each of these, the positive energy one follows the asymptotics prescribed by Eq. (26), while the negative energy modes become associated with vibrations of the dark solitons. The second anomalous mode relates to the out-of-phase vibration of the outer dark solitons, while the middle one remains quiescent (cf. also a relevant experimental result in Ref. ourmarkus3 ). The third anomalous mode corresponds to a more complex motion where the two outer solitons are in-phase, while the middle one is out-of-phase with respect to them. These are both captured accurately by Eq. (31) of Ref. coles .

Figure 3 confirms that the above theoretical predictions are again in good agreement with numerical observations. The decreasing tendency of the real part of the first three eigenfrequencies (or eigenfrequency pairs) is accurately captured by Eqs. (30)–(32). Equally importantly, the growth rates relating to the two instabilities due to the resonances at $2\Omega$ and $3\Omega$ that are also accurately represented near the linear limit.

III.2 Collisionally Inhomogeneous Case

We now extend our considerations to the collisionally inhomogeneous case, to examine the variations on the spectral picture that are induced by the presence of $g(x)$. As indicated above, we will primarily focus on a periodic variation, cf. Eq. (3), as a principal building block towards more complex such variations. Interestingly, the analytical considerations of the previous section can still be carried out in the present setting. However, the explicit formulae resulting are far more tedious. For instance, in the case of the single dark soliton the chemical potential correction $\mu_{1}$ is given by:

[TABLE]

It is even more interesting to examine in this case the corrections to the spectral frequencies. Recall that, as before, at $\omega=\Omega$ two frequencies collide, the positive energy one, with $n=2$, and the negative energy one, with $n=0$, for this case of $m=1$. However, importantly, the presence of the periodic variation destroys the invariance associated with the dipolar mode. As a result, the resonance may arise in this setting and is associated with an eigenvalue given by:

[TABLE]

where $\epsilon=(\mu-3\Omega/2)/\mu_{1}$ and $\omega_{1}^{{\rm cor}}$ is given by the following lengthy expression:

[TABLE]

where

[TABLE]

Furthermore, it is straightforward to observe that in the homogeneous limit of $g_{0}\rightarrow 0$, this result retrieves Eq. (23), as well as the dipolar mode with $\omega=\Omega$. It is relevant to also point out that the above expression, due to the presence of the radical, provides implicitly conditions (through the zero crossing of $B$) under which the two real pairs colliding at this frequency lead to a complex eigenfrequency quartet. However, it is also evident that while these expressions are available in explicit form, they are particularly tedious and, hence, we will only provide them for the single soliton case. The other two modes associated with $n=3$ and $n=4$, at $2\Omega$ and $3\Omega$ respectively, are non-resonant and yield eigenfrequency corrections, as follows:

[TABLE]

Let us now compare these predictions with the numerical results, as shown in Fig. 4. The bottom left panel of the figure shows an example of the single soliton profile, at a relatively large value of $\mu=3.962$, in order to clearly illustrate the effect of the periodic variation towards the solitonic structure inside the parabolic trap (also shown). The rest of the parameters are chosen as $\Omega=0.1$, $g_{0}=0.1$ and $k=0.5$. These three parameters will indeed be fixed hereafter to those values, although it is clear from Eqs. (33)–(III.2) that our methodology can tackle general parametric combinations. In this case, it turns out that the quantity $B$ in Eq. (35) is negative, hence the relevant pair of modes predicted in Eq. (35) is complex for values of $\mu$ near the linear limit. This, as well as the overall trend of the eigenfrequency variations in Eq. (35) is accurately predicted by the theory (cf. top left panel of Fig. 4). Moreover, the imaginary part of the relevant eigenfrequency (i.e., the growth rate of the associated instability) is also accurately captured, as shown in the top right panel of the figure.

However, in the same panel, it can be seen that the existence of a resonance and an instability of even the single dark soliton (entirely contrary to what is the case in the collisionally homogeneous environment) is not the only feature of the spectrum. In addition to that instability (occurring, in this case, for $0.15<\mu<1.16$), there is another instability that arises for higher values of the chemical potential i.e., for $3.16<\mu<4.76$. This instability can be straightforwardly inferred from the bottom right panel, showing the spectral plane $(\lambda_{r},\lambda_{i})$ of the eigenvalues $\lambda=\lambda_{r}+i\lambda_{i}$. There, it can be seen that the anomalous mode as it increases collides with the third pair of eigenfrequencies in the vicinity of $3\Omega$ producing this instability. It is important to highlight once again that, for $g_{0}=0$, the single dark soliton is generically spectrally stable, hence all these instability features do not arise.

In the case of the two soliton state, for the above mentioned choice of the parameters, the stability and existence results are shown in Fig. 5. Here, the correction to the chemical potential is given by:

[TABLE]

yet the rest of the analytical expressions is, arguably, too complex to be presented here. Instead, for reasons of completeness, we provide the theoretical predictions for the numerically selected values of $\Omega=0.1$, $g_{0}=0.1$ and $k=0.5$. We find that:

[TABLE]

with $\epsilon=(\mu-0.25)/0.0834$. The first two pairs in Eq. (36) are meant to indicate that the two eigenfrequency pairs for these parameter values are not resonant immediately after the linear limit. In fact, in one of the cases their sign of variation is opposite (for those bifurcating at $\omega=\Omega$). The other pair at $2\Omega$ leads to an instability shortly after the linear limit [although in agreement with Eq. (36) this does not happen immediately]. It can be seen in the top left panel of Fig. 5 that, once again, the near linear stability predictions are fairly accurate here. Naturally, for larger values of $\mu$, more complex phenomena may arise as, e.g., the collision of the first positive and the first negative energy modes observed in the interval $1.312<\mu<2.313$. This is also showcased in the bottom panel examples of the solution profile and its BdG spectrum, shown for $\mu=1.813$.

However, it is important to highlight again here the differences that the collisional inhomogeneity may bring to the collisionally homogeneous picture. Namely, in the collisionally homogeneous setting, the dark soliton pair would be unstable in the immediate vicinity of the linear limit, while here the presence of a spatial modulation of the nonlinearity may (parametrically) “delay” this instability, while it may induce other ones in different parametric ranges (similarly to what we found also in the single dark soliton case).

Finally, we also consider the continuation over the chemical potential of the branch of three solitons in a collisionally inhomogeneous setting. In this case, the results are presented in Fig. 6. The corresponding theoretical predictions near the linear limit are

[TABLE]

with $\epsilon=(\mu-0.35)/0.0746$. Despite the fairly complicated nature of the relevant structure, once again we can see in Fig. 6 that the analytical predictions do well in capturing the numerical results. This is true both for the first two sets of modes, where despite the potential for resonance, it is instead predicted that for the considered parameter values, no such resonance exists. It is also true for the third set of pairs that still leads to a resonance and to a complex eigenvalue quartet, accurately captured both in the real part (top left panel of the figure) and in the imaginary part (top right panel showing the instability growth rate) of the corresponding eigenfrequency.

However, in this case too, there are features that are arguably somewhat unexpected and, in any case, fundamentally distinct from the collisionally homogeneous case. In particular, we observe that the lowest anomalous eigenfrequency decreases in value until it collides with the origin, resulting in the relevant eigenfrequency pair exiting as purely imaginary for the parametric interval $1.36<\mu<2.04$. This type of instability was never observed for the configurations considered in the collisionally homogeneous case. The instability is also illustrated in the bottom panels of the figure representing the solution profile and the BdG stability analysis for $\mu=1.688$.

The above results clearly show that the collisionally inhomogeneous case presents a number of features (illustrated in Figs. 4-6) that are absent in the context of the collisionally homogeneous setting. Some of the prototypical ones of these features are worth dynamically exploring, as is depicted in Fig. 7. Our earlier work in the collisionally homogeneous setting ourmarkus3 demonstrated that the oscillatory instabilities of the two- and three-dark solitons there corresponded to modes of energy exchange between the solitons and the background – cf. Figs 2 and 3. In the top panel of Fig. 7, we encounter an oscillatory instability even for a single dark soliton in line with the results of Fig. 4 (an instability absent for homogeneous settings). This instability is seen to lead to an oscillation of the solitary wave and an energy exchange with a corresponding mode of the background (as is implied by the associated resonance). Nevertheless, the dynamics of the solitary wave remains oscillatory for a long time, without any other major features arising in the dynamics for the soliton shown in Fig. 4, for $\mu=3.962$.

In the middle panel of Fig. 7, we see an instability of the two soliton state arising from the collision of the anomalous mode with $n=1$ with the “former dipolar” mode of $n=3$. Note that this instability for $\mu=1.813$ (the soliton is shown in the bottom panel of Fig. 5) would be absent in the collisionally homogeneous limit. Namely, there is no instability in the latter limit for the soliton in-phase motion. Yet, the presence of collisional inhomogeneity destroys the invariance associated with the dipolar motion and allows the resonant collision of the in-phase mode with the dipolar one and the corresponding emergence of an oscillatory growth of the in-phase motion, as depicted in the middle panel of Fig. 7. Eventually, this mode undergoes a dramatic modification in its character, as a result of the instability resulting in one of the dark solitons executing smaller amplitude oscillations near the center, while the other one executes large amplitude oscillations, reaching the rims of the condensate. This behavior is reminiscent, albeit involving a smaller number of solitons, of the quantum’s Newton cradle behavior observed in larger dark soliton chains, as described in Ref. NewtonCradle . It is interesting to note that some of dark soliton collisions in our setting render evident the repulsive nature of the interaction, while others illustrate the finite height of the repulsive barrier, given that the solitons effectively can go through each other (if they possess sufficient kinetic energy, as, e.g., in the collisions around $t=1340$ or $t=1500$).

Lastly, yet another example of an instability absent in the collisionally homogeneous limit is showcased for $\mu=1.688$ in the bottom panel of Fig. 7. This is the genuinely exponential instability arising in the case of the three-soliton solution, shown in the bottom panel of Fig. 6. Here, we can see that the mode is no longer resonant with another mode, but merely reflects the exponential growth eventually associated with the in-phase motion of the three solitons, as expected based on the spectrum of Fig. 6. As an aside, it is interesting to point out here the more curved trajectory of the outer solitons as they reach the outer rims of the condensate, presumably because they feel more intensely there the periodic variation of the nonlinearity coefficient.

IV Conclusions and Future Work

In the present work we explored the stability – and partly the associated dynamics – of single and multiple matter-wave dark solitons in BECs confined in parabolic traps. We focused on the contradistinction between the collisionally homogeneous case, of a constant nonlinearity coefficient, and the collisionally inhomogeneous one, of a spatially-varying nonlinearity coefficient. As our approach of choice, we utilized the vicinity of the linear limit (corresponding to the quantum harmonic oscillator problem) enabling us to employ a variant of the degenerate perturbation theory for this Hamiltonian system. Our analytical approach provided a systematic handle on the nature of the (positive energy, or negative energy –anomalous–) modes. This also allowed us to quantitatively characterize these modes, at least in the vicinity of the linear limit. In all the cases considered, this was found to be a useful tool for understanding both the motion of the relevant modes, as well as for assessing their potential for yielding oscillatory instabilities associated with complex bifurcations near this limit. Supplemented with an understanding (at least in the collisionally homogeneous case) of the spectral picture for large chemical potentials, our analytical approximations completed the full range of parametric regimes that are analytically tractable. On the other hand, the numerical computations provided a smooth interpolation between these asymptotic cases.

Our investigation revealed a significant wealth of differences between the collisionally homogeneous and the collisionally inhomogeneous cases. For instance, in the latter, even the single dark soliton structure might become unstable, while it is always stable in one-dimensional collisionally homogeneous setting. Also, two-soliton states may feature an in-phase resonant instability, which is again excluded by symmetry in the homogeneous case. Finally, for three solitons even an exponential instability of their in-phase mode is possible, again a trait that is not encountered in the homogeneous case. For these instabilities, given their previously unexplored nature, we examined the ensuing dynamics through direct numerical simulations.

There are numerous questions that are of interest to explore in future work. On the one hand, tackling the large chemical potential limit of the inhomogeneous case, and understanding the relevant Bogolyubov-de Gennes spectrum in that limit in as much generality as possible, would certainly be of interest. On the other hand, in the case of higher-dimensional structures (such as vortex states) in a single component, or in that of structures arising in multicomponent systems (even in one dimension, such as the dark-bright soliton), developing a perturbative count of the relevant eigenvalues and of their dependence on the chemical potential (or potentials, as in multiple components, there are multiple such parameters) would be a direction also very worthwhile to pursue. Work along some of these directions is currently in progress and will be reported in future publications.

Acknowledgements.

P.G.K. gratefully acknowledges the support of NSF-DMS-1312856, and NSF-PHY-1602994, the Alexander von Humboldt Foundation, as well as from the the ERC under FP7, Marie Curie Actions, People, International Research Staff Exchange Scheme (IRSES-605096). R.C.G. gratefully acknowledges the support of NSF-DMS-1309035.