Stability and dispersion relations of three-dimensional solitary waves in trapped Bose–Einstein condensates

Nonlinear waves in superfluids are the subject of intense theoretical and experimental research. The exquisite control achieved in manipulating ultracold atomic gases has enabled the creation, manipulation, and detection of dark solitons [1, 2], vortex rings (VRs) [3, 4] and solitonic vortices in Bose–Einstein condensates (BECs) [5] and Fermi gases along the BEC–BCS crossover [6]. All these structures are strongly influenced by the confinement of the particle cloud and represent solitary waves in the sense that they are characterised by an excitation energy density above the condensate ground state that is localised with respect to the long axis of the confining geometry. These multi-dimensional solitary waves distinguish themselves by their non-trivial topology associated with the constituting superfluid currents.

The combination of confinement and superfluid currents is also the main constituent in the development of atomtronic devices, and, to this end, an in-depth understanding of the nonlinear phenomena involved in such dynamics is required. In particular, the role played by nonlinear waves deserves special attention. On the one hand, the shape-preserving evolution of solitary waves, in both repulsively and attractively interacting systems, could be a useful feature to be implemented in future applications, in a similar way as optical solitons are being currently used in optical fibers [7]. On the other hand, the necessity to better understand the role played by solitary waves in the generation of superfluid currents has manifested itself in a series of experiments with superfluid rings at NIST [8, 9]. Therein, vortices of solitonic nature, due to the transverse trapping along the radius of the rings, have been found after driving the superfluid into motion. Additionally, the generation of such metastable loop currents has been demonstrated to be mediated by the existence of solitary waves that produce an energy barrier preventing phase slips [10].

Dark solitons, or kinks, are density dips with an associated jump in the phase of the order parameter, and represent nonlinear excitations in BECs with repulsive interparticle interactions [11]. In one-dimensional rings, kink excitations represent intermediate stages connecting states with different winding numbers [12]. A one-dimensional dark soliton can be understood as a vortex which is crossing the ring, and hence providing a characteristic density depletion and phase slip that depends on the position of the vortex. In higher dimensions, the structure of a vortex line crossing the ring, a solitonic vortex (SV), can be more easily identified on the cross section of the system. Alternatively, other transverse states containing vortex lines can be excited in order to produce a given phase jump along the ring circumference. In general, multidimensional stationary kinks are dynamically unstable [13], unless a tight trap could keep the system in the quasi-one-dimensional regime so that higher energy transverse excitations were excluded [14, 15].

In trapped superfluids, Chladni solitons [16] emerge from the decay of three-dimensional kinks, as a result of the excitation of standing waves on the nodal plane of the kink. Such waves produce patterns of vorticity along the nodal lines of the transverse modes in analogy to the Chladni figures visualising the nodal lines of plate vibration modes [17]. In traps with cylindrical symmetry, the different families of solitary waves can be described by the radial p and azimuthal l quantum numbers, indicating the number of transverse nodal lines along their respective directions. Solitonic vortices, belonging to the family $(p=0,l=1)$, are the lowest energy states in elongated condensates [15, 18], while VRs [19, 20], presenting higher excitation energies, belong to the family $(p=1,l=0)$. Along with these previously known states, there exist a sequence of more complex stationary solitary waves with all possible combinations of p and l quantum numbers, as was pointed out by the authors [16]. Very recently evidence for the observation of the Φ-shaped Chadni soliton with quantum numbers $(p=1,l=1)$ in a superfluid Fermi gas at unitary was reported in [21]. The relative strength of the decay modes that can produce Chladni solitons from the decay of the kink, as well as the robustness of the Chladni solitons are key points that remain to be clarified.

Motivated by the previous considerations, in the present work we study the dynamics of solitary wave excitations within the framework of the time-dependent Gross–Pitaevskii equation. Section 2 is devoted to characteristic mass parameters relevant for the Landau quasiparticle dynamics of solitary waves. Numerical data is presented and compared to analytical approximations for energy, inertial, and physical masses of the dark soliton in section 2.1, and the SV and VR in section 2.2. The snaking instability of the kink state is analysed in detail in section 3 by solving the Bogoliubov equations of linearised excitations for the trapped kink state numerically and by developing a semi-analytical theory of the unstable modes. Numerical results for stationary Chladni solitons are reported in section 4, while the dispersion relations and phase step of moving Chladni solitons are considered in section 5. A stability analysis of Chladni solitons—stationary and moving—is performed in section 6, where also results from real-time evolution beyond the linear response regime are reported. We identify two characteristic scenarios in the fate of Chladni solitons: either a chain of decay episodes into single vortex lines which are localised around a transverse plane of the system, or the generation of secondary travelling waves. Finally, the dynamical features pointed out in our study are used to propose feasible protocols for the experimental realisation of these solitary waves.

Solitary waves often exhibit particle-like dynamics. One manifestation of such particle-like dynamics occurs when solitary waves move across a slowly varying background where energy radiation is being suppressed. As a consequence of energy conservation, the solitary wave will then adapt adiabatically to the changing environmental conditions, adjusting its internal parameters as to maintain its local energy constant, and acting as a Landau quasiparticle [22]. The nonlinear wave solutions considered in this article are all solitary waves in this sense, because they are localised with respect to the long axis of a trapped geometry and thus may perform guided quasiparticle motion, even though their properties are significantly influenced by the presence of a transverse confining potential. In sections 3 and 6 we present evidence to show that only two types of solitary waves are dynamically stable in cylindrically trapped BECs: the SV is stable when $\mu /({\hslash }{\omega }_{\perp })\gt 2.65$ while the dark soliton is stable below this value, where μ is the chemical potential and ${\omega }_{\perp }$ the frequency of the transverse harmonic trapping potential. Dynamically stable solitary wave can be expected to perform near-hamiltonian quasiparticle dynamics for a long time and have been observed in experiments with trapped superfluid Fermi gases for several seconds [6, 23]. Unstable solitary waves like VRs may still exhibit quasiparticle like dynamics if the competing decay dynamics is slow enough or suppressed by symmetry constraints [24, 25].

In the framework of Landau quasiparticle dynamics, the equations of motion of a solitary wave in a trapped quantum gas can be derived from knowing the excitation energy ${E}_{s}(\mu ,{v}_{s})$ of the solitary wave as a function of the chemical potential μ and its velocity v_s [22]. In a trapped gas, the chemical potential is then treated as a (slowly varying) function of the position Z of the solitary wave, while the velocity is the time derivative of position ${v}_{s}=\dot{Z}$. Requiring the energy to be a constant of motion then leads to

$$\begin{eqnarray}&&0=\displaystyle \frac{{\rm{d}}{E}_{s}}{{\rm{d}}t}=\left({\left.\displaystyle \frac{\partial {E}_{s}}{\partial \mu }\right|}_{{v}_{s}}\displaystyle \frac{{\rm{d}}\mu }{{\rm{d}}Z}+{M}_{\mathrm{in}}\ddot{Z}\right)\dot{Z},\end{eqnarray} \tag{ 1 }$$

where

$$\begin{eqnarray}&&{M}_{\mathrm{in}}=\displaystyle \frac{1}{{v}_{s}}{\left.\displaystyle \frac{\partial {E}_{s}}{\partial {v}_{s}}\right|}_{\mu },\end{eqnarray} \tag{ 2 }$$

is the inertial mass of the solitary wave, often also called the effective mass. Equation (1) already looks like Newton's law. For weak harmonic trapping potential along the z axis where the Thomas–Fermi approximation demands that $\mu (Z)={\mu }_{0}-\frac{1}{2}m{\omega }_{z}^{2}{Z}^{2}$, we arrive at Newton's equation of motion in the form

$$\begin{eqnarray}&&{M}_{\mathrm{in}}\ddot{Z}=-{M}_{\mathrm{ph}}{\omega }_{z}^{2}Z,\end{eqnarray} \tag{ 3 }$$

where the physical mass defined by

$$\begin{eqnarray}&&{M}_{\mathrm{ph}}=m{\left.\displaystyle \frac{\partial {E}_{s}}{\partial \mu }\right|}_{{v}_{s}},\end{eqnarray} \tag{ 4 }$$

is to be interpreted as the characteristic parameter of the solitary wave that gives rise to the buoyancy-like force on the right-hand side of equation (3). The interpretation of the physical mass to be related to a buoyancy phenomenon is further supported by it being closely related to number of missing particles N_s, i.e. particles depleted from the background density due to the presence of the solitary wave, where ${M}_{\mathrm{ph}}={{mN}}_{s}$ holds in many cases [26].

Solitary waves in repulsively interacting quantum gases typically have negative inertial and physical masses, which leads to oscillatory motion of the solitary waves in a trapped gas. This, e.g. is the case for the one-dimensional Gross–Pitaevskii equation describing BECs with tight transverse confinement [22, 27]. If the physical and inertial masses are independent of position and velocity, or in the limit of small-amplitude motion, we obtain simple harmonic oscillations $Z(t)\propto \mathrm{sin}({\rm{\Omega }}t)$ with [26, 28]

$$\begin{eqnarray}&&\displaystyle \frac{{\omega }_{z}^{2}}{{{\rm{\Omega }}}^{2}}=\displaystyle \frac{{M}_{\mathrm{in}}}{{M}_{\mathrm{ph}}}.\end{eqnarray} \tag{ 5 }$$

Such harmonic oscillations have already been observed in experiments and the frequency ratio measured for dark solitons [29] and for solitonic vortices [30] in BECs and for solitonic vortices in the superfluid Fermi gas [6, 23]. In the remainder of this section we present numerical data of dispersion relations and mass parameters evaluated for the dark soliton, SV and VRs, and compare with approximate analytical expressions.

2.1. Planar dark solitons

We will model solitary waves within the mean field theory given by the Gross–Pitaevskii equation for the condensate order parameter ${\rm{\Psi }}({\bf{r}},t)$

$$\begin{eqnarray}&&{\rm{i}}{\hslash }\displaystyle \frac{\partial {\rm{\Psi }}}{\partial t}=\left(-\displaystyle \frac{{{\hslash }}^{2}}{2m}{{\rm{\nabla }}}^{2}+V({\bf{r}})+g{\left|{\rm{\Psi }}\right|}^{2}\right){\rm{\Psi }},\end{eqnarray} \tag{ 6 }$$

where $g=4\pi {{\hslash }}^{2}a/m$ is the interaction strengh determined by the positive s-wave scattering length a and the bosonic mass m, $V({\bf{r}})=m{\omega }_{\perp }^{2}{r}_{\perp }^{2}/2$ is an external, cylindrically symmetric, harmonic potential in the transverse coordinate ${r}_{\perp }^{2}={x}^{2}+{y}^{2}$, and the condensate particle number N follows from normalisation $N=\int {\rm{d}}{{\bf{r}}}^{3}{\left|{\rm{\Psi }}\right|}^{2}$.

Our starting point is the search for stationary solutions to equation (6), ${\rm{\Psi }}({\bf{r}},t)={{\rm{e}}}^{-{\rm{i}}\mu t/{\hslash }}\psi ({\bf{r}})$, with chemical potential μ, having the form of planar kinks across the axial coordinate z. This task has been carried out numerically because no analytical solution is known for the 3D Gross–Pitaevskii equation (6). Nevertheless, we have been guided by the asymptotic analytical solution proposed in [16]

$$\begin{eqnarray}&&\psi ({\bf{r}})={\chi }_{\mathrm{TF}}({r}_{\perp })\mathrm{tanh}(\frac{z}{\xi ({r}_{\perp })}),\end{eqnarray} \tag{ 7 }$$

which is valid in the Thomas–Fermi regime, where ${\chi }_{\mathrm{TF}}=\sqrt{{\mu }_{\mathrm{loc}}({r}_{\perp })/g}$ is the transverse ground state, ${\mu }_{\mathrm{loc}}({r}_{\perp })=\mu -V({r}_{\perp })$ is the local chemical potential, and a local healing length is defined by $\xi ({{\bf{r}}}_{\perp })={\hslash }/\sqrt{m{\mu }_{\mathrm{loc}}({r}_{\perp })}$. Employing equation (7) as initial ansatz, the numerical solution of equation (6) is obtained without difficulty either using a Newton method or by imaginary time evolution.

The ansatz (7) also provides an excellent description of relevant properties of the kink, such as excitation energy or missing number of particles. We define the excitation energy of a soliton ψ, relative to the ground state ${\psi }_{0}$, by means of

$$\begin{eqnarray}&&{E}_{s}=E[\psi ]-\mu N-(E[{\psi }_{0}]-\mu {N}_{0}),\end{eqnarray} \tag{ 8 }$$

with the energy defined by the functional

$$\begin{eqnarray}&&E[\psi ]=\int {\rm{d}}{\bf{r}}\left(\displaystyle \frac{{{\hslash }}^{2}}{2m}{\left|{\rm{\nabla }}\psi \right|}^{2}+V({\bf{r}}){\left|\psi \right|}^{2}+\displaystyle \frac{g}{2}{\left|\psi \right|}^{4}\right).\end{eqnarray} \tag{ 9 }$$

Substituting the ansatz (7) into equation (8) and neglecting the derivatives of ${\chi }_{\mathrm{TF}}({r}_{\perp })$ in equation (9), according to the Thomas–Fermi approximation, we get the energy of a planar dark soliton confined by a transverse harmonic trapping as

$$\begin{eqnarray}&&\displaystyle \frac{{E}_{s}}{{\hslash }{\omega }_{\perp }}=\displaystyle \frac{4}{15}\displaystyle \frac{{\tilde{\mu }}^{\frac{5}{2}}}{\tilde{a}},\end{eqnarray} \tag{ 10 }$$

where the quantities with tilde are measured in the characteristic units of the trap, $\tilde{\mu }=\mu /{\hslash }{\omega }_{\perp }$, $\tilde{a}=a/{a}_{\perp }$, and ${a}_{\perp }=\sqrt{{\hslash }/m{\omega }_{\perp }}$ . Normalisation of equation (7) gives the missing number of particles in the soliton ${N}_{s}=N-{N}_{0}$ :

$$\begin{eqnarray}&&{N}_{s}=-\displaystyle \frac{2}{3}\displaystyle \frac{{\tilde{\mu }}^{\frac{3}{2}}}{\tilde{a}},\end{eqnarray} \tag{ 11 }$$

which can also be obtained from the relation

$$\begin{eqnarray}&&\partial {E}_{s}/\partial \mu ={-N}_{s}.\end{eqnarray} \tag{ 12 }$$

Figure 1 shows a comparison of the analytical expressions (10) and (11) derived from the ansatz (7), and the exact numerical solution to equation (6). As expected, the agreement improves with increasing values of the chemical potential, and for $\tilde{\mu }\gt 10$ errors are below 1%. The preceding analysis shows that for a harmonically trapped kink, given a chemical potential in the trap units, the parameters $\tilde{a}\tilde{{E}_{s}}$ and $\tilde{a}{N}_{s}$ are fixed, as reflected in equations (10) and (11).

Zoom In Zoom Out Reset image size

A further improvement can easily be introduced in previous expressions for average quantities. Following [31], by properly incorporating the zero point energy of the harmonic oscillator in the calculations we obtain the improved expressions

$$\begin{eqnarray}&&\displaystyle \frac{{E}_{s}}{{\hslash }{\omega }_{\perp }}=\displaystyle \frac{4}{15}\displaystyle \frac{{(\tilde{\mu }-1)}^{\frac{5}{2}}}{\tilde{a}}+\displaystyle \frac{2}{3}\displaystyle \frac{{(\tilde{\mu }-1)}^{\frac{3}{2}}}{\tilde{a}},\end{eqnarray} \tag{ 13 }$$

and

$$\begin{eqnarray}&&{N}_{s}=-\displaystyle \frac{2}{3}\displaystyle \frac{{(\tilde{\mu }-1)}^{\frac{3}{2}}}{\tilde{a}}-\displaystyle \frac{{(\tilde{\mu }-1)}^{\frac{1}{2}}}{\tilde{a}},\end{eqnarray} \tag{ 14 }$$

which have the correct limits both in the Thomas–Fermi and in the quasi-onedimensional regimes. They also interpolate in between with very good accuracy as can be seen in figure 1.

2.2. SV and VR

We have determined the inertial and physical masses for solitary waves by computing the energy of the fully numerical solution of the 3D Gross–Pitaevskii equation (6) and taking numerical derivatives with respect to chemical potential and velocity near the stationary point. Figure 2 reports the resulting mass ratios for dark solitons (red curve), solitonic vortices (green curve with open squares), and VRs (blue curve with open circles) in 3D condensates confined by isotropic radial harmonic traps. For the kink, one can observe discontinuities arising at the bifurcation points of $(p,0)$ modes, corresponding to VRs. As will be explained in later sections, at these points the kink solutions exist only for the static dark soliton, and moving VRs emerge with the same value of the chemical potential.

Zoom In Zoom Out Reset image size

Approximate formulas for the SV properties, obtained from hydrodynamic theory in logarithmic accuracy, appeared in [6]. The energy expression for the stationary SV in a BEC reads

$$\begin{eqnarray}{E}_{\mathrm{SV}} & = & \displaystyle \frac{\pi {{\hslash }}^{2}m}{{m}_{B}^{2}}{n}_{2}\mathrm{ln}\displaystyle \frac{{R}_{\perp }}{\xi }\\ & = & \displaystyle \frac{\sqrt{2}}{3}{\tilde{\mu }}^{\frac{3}{2}}\mathrm{ln}(2\tilde{\mu })\displaystyle \frac{{a}_{\perp }}{a}{\hslash }{\omega }_{\perp },\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&=\displaystyle \frac{4\sqrt{2}}{3\sqrt{\tilde{\mu }}}\mathrm{ln}(2\tilde{\mu }){a}_{\perp }{n}_{1}{\hslash }{\omega }_{\perp },\end{eqnarray} \tag{ 16 }$$

and the Thomas–Fermi approximation has been used for the one- and two-dimensional densities, in particular $4{{an}}_{1}={\tilde{\mu }}^{2}$. The missing particle number is obtained by differentiation

$$\begin{eqnarray}&&{N}_{\mathrm{SV}}=-\displaystyle \frac{m}{{m}_{B}}\pi {n}_{2}{\xi }^{2}\left(\displaystyle \frac{2\gamma +1}{\gamma }\mathrm{ln}\displaystyle \frac{{R}_{\perp }}{\xi }+\displaystyle \frac{2}{\gamma }\right)\end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}&&=-\displaystyle \frac{\sqrt{2}}{2}\displaystyle \frac{{a}_{\perp }}{a}\sqrt{\tilde{\mu }}\left(\mathrm{ln}(2\tilde{\mu })+\displaystyle \frac{2}{3}\right),\end{eqnarray} \tag{ 18 }$$

where $\gamma =\frac{\mu }{n}\frac{\partial n}{\partial \mu }$ is a polytropic index characterising the equation of state, which evaluates to $\gamma =1$ for the case of a BEC. For the inertial mass the following expression was obtained

$$\begin{eqnarray}&&{M}_{\mathrm{in}}^{{SV}}=-\displaystyle \frac{4\pi }{2\gamma +1}\displaystyle \frac{{n}_{2}{R}_{\perp }^{2}}{\mathrm{ln}\displaystyle \frac{{R}_{\perp }}{\xi }}m,\end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray}&&=-\displaystyle \frac{8\sqrt{2}}{9}\displaystyle \frac{{\tilde{\mu }}^{\frac{5}{2}}}{\mathrm{ln}(2\tilde{\mu })}\displaystyle \frac{{a}_{\perp }}{a}m.\end{eqnarray} \tag{ 20 }$$

3.1. Bogoliubov stability analysis

After finding stationary kink solutions to Gross–Pitaevskii equation (6), we look for elementary excitations $\{u({\bf{r}}),v({\bf{r}})\}$ with angular frequency ω around every equilibrium state ψ with chemical potential μ. The perturbed state can be written as ${\rm{\Psi }}({\bf{r}},t)={{\rm{e}}}^{-{\rm{i}}\mu t/{\hslash }}\left[\psi +{\displaystyle \sum }_{\omega }(u\;{{\rm{e}}}^{-{\rm{i}}\omega t}+{v}^{*}{{\rm{e}}}^{{\rm{i}}\omega t})\right]$, and the excitation modes are the solutions to the linear Bogoliubov equations

$$\begin{eqnarray}&&\left({H}_{0}-\mu +2\;g| \psi {| }^{2}\right)u+g{\psi }^{2}v={\hslash }\omega u,\end{eqnarray} \tag{ 21a }$$

$$\begin{eqnarray}&&-g{{\psi }^{*}}^{2}u-\left({H}_{0}-\mu +2\;g| \psi {| }^{2}\right)v={\hslash }\omega v,\end{eqnarray} \tag{ 21b }$$

where ${H}_{0}={-{\hslash }}^{2}{{\rm{\nabla }}}^{2}/2m+V({{\bf{r}}}_{\perp })$. Dynamical instabilities are related to solutions to equation (21) with complex frequencies ω, where the imaginary part of ω is a rate of exponential growth of the corresponding unstable mode. The inverse of the imaginary part of the unstable frequency can thus be interpreted as the lifetime of the particular mode. Below we will report fully numerical solutions of equation (21) for the dark soliton and Chladni solitons. Here we will proceed with finding analytical solutions to the Bogoliubov equations for the dark soliton based on the Thomas–Fermi approximation.

3.2. Approximate separation of variables

For a real-valued stationary state, as it is the case of the dark soliton or kink of equation (7), the Bogoliubov equations (21) can be transformed using ${f}_{\pm }({\bf{r}})=u({\bf{r}})\pm v({\bf{r}})$ and ${g}_{\pm }=(2\pm 1)g$ into

$$\begin{eqnarray}&&\left({H}_{0}-\mu +{g}_{\pm }{\psi }^{2}\right){f}_{\pm }={\hslash }\omega {f}_{\mp }.\end{eqnarray} \tag{ 22 }$$

For $\omega =0$ the two equations decouple. This was the problem solved in [16] in order to find the bifurcation points of Chladni solitons from the dark soliton. Here, we are interested in the more general problem of finding solutions to equation (22) with finite imaginary or complex values of ω. Aiming at an approximate separation of transverse and longitudinal degrees of freedom we introduce the rescaled variable $\bar{z}=z/\xi ({r}_{\perp })$. For the dark soliton state of equation (7) we can write $g{\psi }^{2}={\mu }_{\mathrm{loc}}({r}_{\perp }){\mathrm{tanh}}^{2}\bar{z}$, where ${\mu }_{\mathrm{loc}}=\mu -V({r}_{\perp })$ and may thus rewrite the Bogoliubov equation (22) as

$$\begin{eqnarray}&&\left(-\displaystyle \frac{{{\hslash }}^{2}}{2m}{{\rm{\nabla }}}_{\perp }^{2}+{\mu }_{\mathrm{loc}}{A}_{\pm }\right){f}_{\pm }={\hslash }\omega {f}_{\mp },\end{eqnarray} \tag{ 23 }$$

with

$$\begin{eqnarray}&&{A}_{-}=-\displaystyle \frac{1}{2}\displaystyle \frac{{\partial }^{2}}{\partial {\bar{z}}^{2}}-1+{\mathrm{tanh}}^{2}\bar{z},\end{eqnarray} \tag{ 24a }$$

$$\begin{eqnarray}&&{A}_{+}=-\displaystyle \frac{1}{2}\displaystyle \frac{{\partial }^{2}}{\partial {\bar{z}}^{2}}-1+3{\mathrm{tanh}}^{2}\bar{z}.\end{eqnarray} \tag{ 24b }$$

The operators ${A}_{\pm }$ are both Hamilton operators of one-dimensional Schrödinger equations with a shifted Rosen-Morse potential, whose eigenfunctions are known analytically [32]. Here we use the localised ground states of the respective operators as a restricted basis for the z-dependence of ${f}_{\pm }$, since we expect the unstable modes to be localised near the dark soliton plane. The ground state eigenfunction of ${A}_{+}$ is ${\varphi }_{+}^{0}(\bar{z})=\frac{\sqrt{3}}{2}{\mathrm{sech}}^{2}\bar{z}$ with eigenvalue 0. It corresponds to the well-known Goldstone mode of translation of the dark soliton in z direction. This does not constitute an instability in itself but the mode will be relevant for constructing the decaying modes with imaginary omega. The operator ${A}_{-}$ has the ground state wave function ${\varphi }_{-}^{-1/2}(\bar{z})=\frac{1}{\sqrt{2}}\mathrm{sech}\;\bar{z}$ with eigenvalue $-\frac{1}{2}$. It is this mode that is responsible for the existence of unstable Bogoliubov modes.

The z dependence can now be removed from the Bogoliubov equation (22) by starting from the ansatz ${({f}_{+},{f}_{-})}^{t}={\chi }_{+}(x,y){({\varphi }_{+}^{0},0)}^{t}+{\chi }_{-}(x,y){(0,{\varphi }_{-}^{-1/2})}^{t}$. Ignoring any $x,y$ derivatives of the functions ${\varphi }_{\pm }(\bar{z})$ and projecting onto the respective ground states by multiplying from the left with $({\varphi }_{+}^{0},0)$ and $(0,{\varphi }_{-}^{-1/2})$, and integrating over $\bar{z}$, we obtain the matrix equation

$$\begin{eqnarray}\left(\begin{array}{cc}-\displaystyle \frac{{{\hslash }}^{2}}{2m}{{\rm{\nabla }}}_{\perp }^{2} & -{\hslash }\omega \zeta \\ -{\hslash }\omega \zeta & -\displaystyle \frac{{{\hslash }}^{2}}{2m}{{\rm{\nabla }}}_{\perp }^{2}-\displaystyle \frac{1}{2}{\mu }_{\mathrm{loc}}\end{array}\right)\left(\begin{array}{c}{\chi }_{+}\\ {\chi }_{-}\end{array}\right)=0,\end{eqnarray} \tag{ 25 }$$

where $\zeta =\int {\varphi }_{-}^{-1/2}{\varphi }_{+}^{0}\;d\bar{z}=\frac{\pi }{4}\sqrt{\frac{2}{3}}\approx 0.962$ is a numerical constant close to one.

3.3. Homogeneous background

The transverse Bogoliubov equation (25) is easily solved in the absence of a transverse trapping potential, where ${\mu }_{\mathrm{loc}}=\mu =\frac{{{\hslash }}^{2}}{m{\xi }^{2}}=\mathrm{const}$, with plane wave solutions, e.g., ${\chi }_{\pm }\propto \mathrm{exp}({\rm{i}}{qx})$. For the unstable eigenvalues we find

$$\begin{eqnarray}&&{\omega }^{\mathrm{hom}}=\displaystyle \frac{{\rm{i}}q}{2\zeta }\sqrt{\displaystyle \frac{\mu }{m}}\sqrt{1-{q}^{2}{\xi }^{2}}.\end{eqnarray} \tag{ 26 }$$

For small q the growth rate $\mathrm{Im}(\omega )$ grows linearly with the slope $\sqrt{\frac{\mu }{4{\zeta }^{2}m}}\approx 0.520\sqrt{\frac{\mu }{m}}$ as is demanded by general hydrodynamic arguments [33]. Although being approximate due to the restricted basis expansion of the z dependence, this result compares very well with the previously obtained ones in [13, 33, 34], where the exact slope for the Gross–Pitaevskii equation is $\sqrt{\frac{\mu }{3m}}\approx 0.577\sqrt{\frac{\mu }{m}}$. The snaking instability is suppressed and eigenvalues become real-valued for wave numbers larger than ${q}_{\mathrm{crit}}=1/\xi$, which is the exact value. For intermediate values $0\lt q\lt 1/\xi$ the growth rate has previously only been obtained numerically, and equation (26) reproduces the results of [13, 34] very closely.

3.4. Harmonically trapped kink state

We now proceed with solving the transverse Bogoliubov equation (25) for an isotropic transverse trapping potential with ${\mu }_{\mathrm{loc}}({r}_{\perp })=\mu -\frac{1}{2}m{\omega }_{\perp }^{2}{r}_{\perp }^{2}$. Assuming an azimuthal dependence $\propto \mathrm{cos}(l\theta )$ with the quantum number $l=0,1,2,\ldots$ reduces equation (25) to a set of ordinary differential equations in the radial coordinate ${r}_{\perp }$. It is now convenient to move to harmonic oscillator units. Introducing the rescaled radial coordinate ${\tilde{r}}_{\perp }={r}_{\perp }/{a}_{\perp }$ and dividing the equation by ${\hslash }{\omega }_{\perp }$, we obtain

$$\begin{eqnarray}\left(\begin{array}{cc}{H}_{1}^{l} & -\tilde{\omega }\zeta \\ -\tilde{\omega }\zeta & {H}_{2}^{l}-\displaystyle \frac{\tilde{\mu }}{2}\end{array}\right)\left(\begin{array}{c}{\chi }_{+}^{l}\\ {\chi }_{-}^{l}\end{array}\right)=0,\end{eqnarray} \tag{ 27 }$$

where $\tilde{\omega }=\omega /{\omega }_{\perp }$ and

$$\begin{eqnarray}&&{H}_{1}^{l}=-\displaystyle \frac{1}{2}\left(\displaystyle \frac{{\partial }^{2}}{\partial {\tilde{r}}_{\perp }^{2}}+\displaystyle \frac{1}{{\tilde{r}}_{\perp }}\displaystyle \frac{\partial }{\partial {\tilde{r}}_{\perp }}-\displaystyle \frac{{l}^{2}}{{\tilde{r}}_{\perp }^{2}}\right),\end{eqnarray} \tag{ 28 }$$

$$\begin{eqnarray}&&{H}_{2}^{l}={H}_{1}^{l}+\displaystyle \frac{1}{4}{\tilde{r}}_{\perp }^{2},\end{eqnarray} \tag{ 29 }$$

where H₁^l represents a two-dimensional Laplacian and H₂^l is the Hamiltonian of a two-dimensional harmonic oscillator with weakened trap potential compared to the one experienced by the atoms. Equation (27) can also be rewritten in the form

$$\begin{eqnarray}&&{H}_{1}^{l}\left({H}_{2}^{l}-\displaystyle \frac{\tilde{\mu }}{2}\right){\chi }_{-}^{l}={\tilde{\omega }}^{2}{\zeta }^{2}{\chi }_{-}^{l}.\end{eqnarray} \tag{ 30 }$$

Even though this represents a non-hermitian eigenvalue problem, we have only found real eigenvalues ${\tilde{\omega }}^{2}{\zeta }^{2}$ in numerical investigations.

3.5. Chladni soliton bifurcation points

Solutions of equation (27) with $\tilde{\omega }=0$ have special significance as they indicate the transition of a specific mode from representing an instability ${\tilde{\omega }}^{2}\lt 0$ to a stable small amplitude oscillation ${\tilde{\omega }}^{2}\gt 0$. At the same time they indicate a bifurcation of the stationary nonlinear solutions of the Gross–Pitaevskii equation (6) and here relate to the branch-off points for Chladni solitons. The solutions are found in terms of the scaled Fock–Darwin radial eigenfunctions ${\chi }_{p}^{l}({\tilde{r}}_{\perp })={2}^{-\frac{1}{4}}{R}_{p}^{l}({2}^{-\frac{1}{4}}{\tilde{r}}_{\perp })$ of the 2D harmonic oscillator Hamiltonian H₂^l with eigenvalues ${\epsilon }_{p}^{l}=(2p+l+1)/\sqrt{2}$ and

$$\begin{eqnarray}&&{R}_{p}^{l}(r)=\sqrt{\displaystyle \frac{2p!}{(p+l)!}}{r}^{l}{L}_{p}^{l}\left({r}^{2}\right){{\rm{e}}}^{-\frac{{r}^{2}}{2}},\end{eqnarray} \tag{ 31 }$$

where ${L}_{p}^{l}(x)$ is the generalised Laguerre polynomial. It is easily seen that ${\chi }_{p}^{l}$ solves the Bogliubov equation (30) with $\tilde{\omega }=0$ when ${\epsilon }_{p}^{l}=\tilde{\mu }/2$, which translates into the condition

$$\begin{eqnarray}&&\displaystyle \frac{\mu }{{\hslash }{\omega }_{\perp }}=\sqrt{2}(2p+l+1),\end{eqnarray} \tag{ 32 }$$

for the bifurcation points of Chladni soliton solutions from the dark soliton, as found previously in [16].

3.6. Finite growth rates

For finite instability rates $\mathrm{Im}(\tilde{\omega })$ the eigenvalue equation (30) can be expanded in a basis of the normalised eigenfunctions $| p,l)\equiv {\chi }_{p}^{l}$ of H₂^l, which transforms it into a tridiagonal matrix eigenvalue equation

$$\begin{eqnarray}&&\displaystyle \sum _{{p}^{\prime }}{B}_{p,{p}^{\prime }}^{l}{v}_{{p}^{\prime }}={\tilde{\omega }}^{2}{\zeta }^{2}{v}_{p}\end{eqnarray} \tag{ 33 }$$

with ${B}_{p,{p}^{\prime }}^{l}=\left(p,l| {H}_{1}^{l}\left({H}_{2}^{l}-\displaystyle \frac{\tilde{\mu }}{2}\right)| {p}^{\prime },l\right)={\displaystyle \int }_{0}^{\infty }{\tilde{r}}_{\perp }\;{\rm{d}}{\tilde{r}}_{\perp }{\chi }_{p}^{l}{H}_{1}^{l}({H}_{2}^{l}-\displaystyle \frac{\tilde{\mu }}{2}){\chi }_{{p}^{\prime }}^{l}$. For the matrix elements we find

$$\begin{eqnarray}&&{B}_{p,p}^{l}=\displaystyle \frac{{(2p+l+1)}^{2}}{4}-\tilde{\mu }\displaystyle \frac{2p+l+1}{4\sqrt{2}},\end{eqnarray} \tag{ 34 }$$

$$\begin{eqnarray}&&{B}_{p,p-1}^{l}=\displaystyle \frac{1}{4}[2(p-1)+l+1]\sqrt{{p}^{2}+{lp}}-\tilde{\mu }\displaystyle \frac{\sqrt{{p}^{2}+{lp}}}{4\sqrt{2}},\end{eqnarray} \tag{ 35 }$$

$$\begin{eqnarray}&&{B}_{p-1,p}^{l}=\displaystyle \frac{1}{4}(2p+l+1)\sqrt{{p}^{2}+{lp}}-\tilde{\mu }\displaystyle \frac{\sqrt{{p}^{2}+{lp}}}{4\sqrt{2}},\end{eqnarray} \tag{ 36 }$$

and ${B}_{p,{p}^{\prime }}^{l}=0$ for $| p-{p}^{\prime }| \gt 1$. The matrix is block-diagonal in the azimuthal quantum number l due to the azimuthal symmetry of the problem. We have solved the corresponding eigenvalue equations numerically and have found the approximate asymptotic behaviour

$$\begin{eqnarray}&&{\tilde{\omega }}_{n}^{l}\sim \displaystyle \frac{\tilde{\mu }-\sqrt{2}(2n+1)}{4\zeta },\end{eqnarray} \tag{ 37 }$$

for large $\tilde{\mu }$ and values of the azimuthal quantum number $l\gt 1$. All these eigenvalues have a zero crossing for a finite value of $\tilde{\mu }$ corresponding to equation (32), where n  =  p can be identified. These results were obtained by diagonalising a truncated matrix ${B}_{p,{p}^{\prime }}^{l}$ with $p,{p}^{\prime }\lt {p}_{{\rm{c}}}$. Changing the cutoff value p_c affects the large-$\tilde{\mu }$ regime but leaves the zero crossings for $n\ll {p}_{{\rm{c}}}$ unaffected. The asymptotic behaviour reported above represents the limit of ${p}_{{\rm{c}}}\to \infty$.

In the l  =  0 sector a special case occurs, where the $| 0,0)$ state needs to be excluded from the basis in order to eliminate an unphysical unstable eigenvalue with p  =  0 that otherwise occurs in solving equation (33). The mode with $p=0,l=0$ corresponds to zero point motion of the kink, which is not captured correctly by the underlying Thomas–Fermi approximation. Indeed, no unstable mode with these quantum numbers is found in the full numerical solution of the Bogoliubov equations. Diagonalising equation (33) in a truncated basis with $0\lt p\lt {p}_{{\rm{c}}}$ produces the asymptotic behaviour

$$\begin{eqnarray}&&{\tilde{\omega }}_{n}^{l=0}\sim \displaystyle \frac{\tilde{\mu }-\sqrt{2}(2n-1)}{4\zeta },\end{eqnarray} \tag{ 38 }$$

where still n  =  p can be identified at the zero crossings of $\tilde{\omega }$.

The results of the truncated eigenvalue problem (33) with the asymptotic results (37) and (38) closely resemble the full numerical results.

3.7. Fully numerical results

We have also solved the Bogoliubov equation (21) numerically in full three-dimensions without making use of the approximations discussed in the previous paragraphs. The results for the unstable modes and associated frequencies are collected as a function of the chemical potential of the kink in figure 3. The insets represent axial views of phase-coloured density isocontours of the excitation modes (p, l) just after the appearance of bifurcation points. These modes present a structure of nodal lines, derived from the radial p and azimuthal l quantum numbers, characteristic of the linear excitations of the transverse trap. As can be seen in the lower left inset, the only one displaying a longitudinal view, they are strongly localised around the plane of the kink. Their emergence follow in a very good approximation the analytical prediction for bifurcations given by equation (32), which are indicated by red arrows below the horizontal axis. As the interaction energy increases from the quasi-onedimensional configuration, where kink states are stable structures and generate only real frequency excitations, the first bifurcation point denoting the appearance of a complex frequency for the mode $| \;p=0,l=1\rangle$ comes into existence at ${\mu }_{\mathrm{0,1}}=2.65{\hslash }{\omega }_{\perp }$, very close to the $2.8{\hslash }{\omega }_{\perp }$ value predicted by (32). From this point on, kinks are unstable states, and further increase of the chemical potential is accompanied by the emergence of new bifurcation points grouped around the integer energy values ${\epsilon }_{{pl}}=2p+l+1$, with the characteristic degeneracy of the two-dimensional harmonic oscillator.

Zoom In Zoom Out Reset image size

Every unstable mode of the kink is associated with a stationary Chladni soliton. The zero crossings of the unstable mode frequencies in figure 3 indicate bifurcation points of the Chladni solitons from the dark soliton. The sign patterns of the unstable mode functions at the bifurcation points ${\chi }_{p}^{l}({r}_{\perp })\mathrm{cos}(l\theta )$ determine the direction of flow along or counter the z direction, and the nodal lines translate into vortex lines. Numerically obtained mode functions are also shown in figure 3. Increasing the nonlinearity parameter $\tilde{\mu }=\mu /({\hslash }{\omega }_{\perp })$ above the bifurcation point, we obtain a family of stationary Chladni solitons with the same structure and symmetries, as endowed by the unstable mode at the bifurcation point. A number of Chladni soliton solutions is shown in figure 4. For instance, the first excited mode $| \;p=0,l=1\rangle$ generates the family of SV states, while the next two unstable modes, $| \;1,0\rangle$ and $| \;0,2\rangle$, that branch off near ${\mu }_{\mathrm{1,0}}=4{\hslash }{\omega }_{\perp }$ in figure 3, generate the families of VRs and two crossed vortices states, respectively.

Zoom In Zoom Out Reset image size

The linear solutions of a general, isotropic or anisotropic, two-dimensional harmonic oscillator can be written in terms of Hermite modes with cartesian symmetries instead of the Laguerre modes of cylindrical symmetry. The Hermite modes are characterised by the pair $({n}_{x},{n}_{y})$ of quantum numbers indicating the nodal lines along x and y directions. Corresponding Chladni solitons would be made of vortex lines in such rectangular patterns. Although we have found the stationary Chladni solitons bifurcating from the kink to all conform to the Laguerre type, one may still expect to find Hermite-type structures to be relevant in the decay process of the kink. In the linear regime, the Hermite modes $| {n}_{x},{n}_{y}\rangle$ can be constructed as superposition of Laguerre modes $| p,l\rangle$ with degenerate energies ${\epsilon }_{{n}_{x}{n}_{y}}={\epsilon }_{{pl}}$. For example,$| {n}_{x}=0,{n}_{y}=2\rangle \propto | p=0,l=2\rangle -| p=1,l=0\rangle$, and $| {n}_{x}=1,{n}_{y}=2\rangle \propto | p=0,l=3\rangle -| p=1,l=1\rangle$. Because of the energy splitting of the bundle of linear degenerate states presented in figure 3, the linear superposition mechanism is not acting in the nonlinear case. However, we have found that Hermite-like modes do emerge as nonlinear bifurcations at a second stage. In the first stage, families of stationary solitary waves bifurcate from the kink in close relation to its unstable Laguerre-like modes. Except for the SV, all these families are made of unstable states, the unstable excitation modes of which can generate new solitary waves at a second stage. It is then when the new solitonic families turn out to be composed of Hermite-like modes $({n}_{x},{n}_{y})$.The instabilities of Chladni solitons will be discussed in more detail in section 6.

Up to now we have been looking at static solitons. In order to understand how the different families of solitonic states appear and connect, it is convenient to consider a more general picture of moving solitary waves. Here, we extend upon the work of Komineas and Papanicolao, who already computed energy dispersion relations and phase steps for axisymmetric solitary waves (kinks and nested VRs) [19, 20] and for the SV [18].

In order to construct the full dispersion relations for moving solitary waves, we numerically search for states ${\rm{\Psi }}({\bf{r}},t)={{\rm{e}}}^{-{\rm{i}}({\mu }^{\prime }t+{{mv}}_{z}z)/{\hslash }}\psi ({\bf{r}})$, moving along the z-axis with a constant velocity ${{\bf{v}}}_{z}=(0,0,{v}_{z})$, which are solutions of the stationary Gross–Pitaevskii equation for a co-moving reference frame:

$$\begin{eqnarray}&&\displaystyle \frac{1}{2m}{(-{\rm{i}}{\hslash }{\boldsymbol{\nabla }}-m{{\bf{v}}}_{z})}^{2}\psi +V({\bf{r}})\psi +g| \psi {| }^{2}\psi ={\mu }^{\prime }\psi ,\end{eqnarray} \tag{ 39 }$$

where ${\mu }^{\prime }=\mu +{{mv}}_{z}^{2}/2$ is the shifted chemical potential. Moving solitons have nonzero density dips associated to reduced (smaller than π) phase jumps, and their velocities are limited by the speed of sound c, at which solitonic states become linear sound excitations. In order to calculate the speed of sound, which will be used as velocity unit in what follows, we will make use of the analytical expression given in [35]:

$$\begin{eqnarray}&&\displaystyle \frac{c}{\sqrt{{\hslash }{\omega }_{\perp }/m}}={\left(\displaystyle \frac{{\tilde{\mu }}^{2}-1}{2\tilde{\mu }}\right)}^{\frac{1}{4}},\end{eqnarray} \tag{ 40 }$$

where $\tilde{\mu }=\mu /{\hslash }{\omega }_{\perp }$. This expression provides the speed of sound for elongated, harmonically trapped condensates with arbitrary values of the interaction, and gives the exact limits both in the quasi-onedimensional and Thomas–Fermi regimes.

Figure 5(b) shows the excitation energy E_s as a function of the axial velocity v_z for moving solitary waves with chemical potential $\mu =5{\hslash }{\omega }_{\perp }$. Kinks, represented by the solid red line at the top of the figure, have the highest excitation energy among the solitary waves, and exist in this case only for low velocities, $| {v}_{z}| \lt 0.24\;c$. Below them, at lower excitation energies, only solitary waves bifurcating at energies less or equal than $\tilde{\mu }=5$ can be found in the graph, as per equation (32). As can also be deduced from the unstable frequencies of figure 3, indeed, only VRs (blue solid line), double SV (dashed black) and the single SV (solid green) are available at $\tilde{\mu }=5$ and appear in figure 5. At small velocity, and very close in energy to the VR and the static cross soliton, a new type of solution emerges (see the inset of figure 5(b)). It is composed of a couple of almost parallel vortices, more specifically a vortex–antivortex pair or vortex dipole, which is not coming directly from a bifurcation of the kink, but from a secondary bifurcation of the cross soliton. Indeed, the solution can be connected to a decay instability of the cross soliton produced by the Hermite modes $({n}_{x}=2,{n}_{y}=0)$ or $({n}_{x}=0,{n}_{y}=2)$. Figure 6(a) shows the density configuration of these states around zero velocity, which can be cross-checked with their features extracted from the three panels of figure 5, where clear differences in the associated phase step arise for the three states A, B and C. At higher velocities, the picture is slighly different. The structure of the kink changes and so do its excitation modes. Since the cross soliton is a static state, it can not be found between the bifurcations of moving kinks, and is substituted by the mentioned, moving Hermite modes. For growing values of the chemical potential, new, higher energy Hermite modes emerge by equivalent mechanisms. In figure 6(b) we show density isocontours for some of such modes with $\tilde{\mu }=5,{\rm{and}}\;10$.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

It is also interesting to look at the phase jump along the axial z-direction ${\rm{\Delta }}\phi =\phi$ $(z\to \infty )$ $-\phi (z\to -\infty )$ created by the different solitary waves (figure 5(a)). The phase shift of dark solitons of the one-dimensional nonlinear Schrödinger equation (dotted line) is shown as a reference. Its phase shows a π jump in the static configuration, and grows up to $2\pi$ as the velocity approaches $+c$, or alternatively reduces to zero as ${v}_{z}\to -c$. Similar behaviour, but with different variation rates, is in general followed by the phases of 3D solitonic states. However, a particular feature can be noted as characteristic of the 3D case. It is the existence of turning points with vertical phase slopes. Looking at the curves for VRs (blue lines with triangles), one can notice that there are two separated branches ending at respective turning points, and connected by kink states (red curve). A more striking manifestation of this phenomenon can be observed in the inset of figure 5(b), on the curve corresponding to the vortex–antivortex pair. Between turning points, labeled as B and C in figure 5(c), there exist a set of almost static and degenerate states of this type which produce different phase jumps. The turning points here indicate a transition between the Hermite-like symmetry of the vortex–antivortex pair and the Laguerre-like vortex cross.

The precedent analysis, in terms of phase and velocity, can now be completed by constructing the dispersion relations of Chladni solitons. To this end, as usual, we define the axial canonical momentum P_c of a solitonic state as the conjugate variable of the axial coordinate z, that fulfils

$$\begin{eqnarray}&&{v}_{z}=\displaystyle \frac{\partial {E}_{s}}{\partial {P}_{{\rm{c}}}}.\end{eqnarray} \tag{ 41 }$$

Along with the axial physical momentum ${p}_{z}=-{\rm{i}}{\hslash }\int {\rm{d}}{\bf{r}}({\psi }^{*}{\partial }_{z}\psi -\psi {\partial }_{z}{\psi }^{*})/2$, carried by the particles traversing the plane of the moving soliton, the canonical momentum includes the contribution coming from the phase jump ${\rm{\Delta }}\phi$ between the axial boundaries of the condensate

$$\begin{eqnarray}&&{P}_{{\rm{c}}}={p}_{z}-{\hslash }{n}_{1}{\rm{\Delta }}\phi ,\end{eqnarray} \tag{ 42 }$$

where n₁ is the axial density of the ground state of the system [20, 26]. Figure 5(c) shows the dispersion relation, excitation energy versus axial canonical momentum, for solitonic states with $\mu =5{\hslash }{\omega }_{\perp }$ moving along the axial coordinate. The turning points that have been described on the phase graph, appear here as the vertexes of cusps connecting states with the twofold symmetry: VRs and vortex–antivortex pairs. It is also worth noticing that the lowest excitation energy levels for fixed momentum are occupied by solitonic vortices, wherever they exist. This last remark accounts for the fact that there exist a small regime of soliton speeds, approaching the speed of sound, where the only solitonic state is the continuation of the VR family, here an axisymmetric solitary wave very much like a grey soliton with a VR phase singularity outside the Thomas–Fermi density of the trapped BEC, as is apparent in figures 5(a) and (b). In this regime, VRs are dynamically stable states.

When the chemical potential increases, and the number of bifurcations grows, the dispersion diagram of Chladni solitons becomes more complex, because of the emergence of new connections between solitary waves sharing symmetries. As an instance of this complexity, figure 7 displays some of the curves making the dispersion diagram for $\mu =10\;{\hslash }\omega$. For this value of the chemical potential, the family of double VRs $[(p=2,l=0)$, violet lines] is available, and produces a new couple of turning points compared to the case for $\mu =5\;{\hslash }\omega$. Following the different curves away from their maximum (corresponding to zero velocity ${v}_{z}=\partial {E}_{s}/\partial {P}_{{\rm{c}}}=0$) the density patterns of Chladni solitons can change dramatically from their static configuration. To illustrate this point, figure 8 shows the density isocontours of some of the moving Chladni solitons that can be found with $\mu =10\;{\hslash }\omega$. It is apparent how their symmetry changes when compared to their static counterparts in figure 4.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

It remains to know how robust the Chladni solitons are. To this end, we have numerically solved the Bogoliubov equations (21) for the linear excitation modes of the stationary and moving solitonic solutions (as those represented in figures 1 and 5). In addition, we have checked the stability of these nonlinear systems against perturbations by monitoring their evolution in real time through the full time dependent Gross–Pitaevskii equation (6).

6.1. Linear analysis

As mentioned before, and aside from the stability of VRs moving near to the speed of sound, our results indicate that the SV branch is the only one containing dynamically stable states. Solitons corresponding to other families are unstable, and decay through the instability channels opened by lower energy branches. Specifically for the stationary solitary waves, the SV as the lowest energy solitary wave has no channel of instability, since there is no other, lower energy solution bifurcating from it (or from the kink). However, the second excited state, which is the single VR, does present one instability channel associated to the bifurcation of solitonic vortices from the kink with lower energy. The next family is that of cross solitons (2SV), and is unstable through two channels, and so on. This analysis for static states is displayed in figure 9(a). The red curve with open circles corresponds to the unstable frequencies for VRs as a function of the chemical potential, and the two blue curves with open triangles indicate the two instability channels for the cross soliton. The dotted vertical lines mark the bifurcation points for the Chladni solitons considered (VR and 2SV), by intersecting the instability curves of kinks (grey dashed lines). It is worth to mention that, as can be seen in figure 9(a), for intermediate values of the chemical potential, between 4 and $5{\hslash }{\omega }_{\perp }$, the cross soliton has smaller values of unstable frequencies than VRs. This fact suggests that cross solitons are good candidates for being experimentally realised in elongated BECs. VRs have already been observed in experiments [3, 37]. In this regard, we have noticed that in [38], where the decay of dark solitons in anisotropic cigar-shaped condensates was observed in experiments, travelling solitary waves composed of vortex–antivortex pairs (see below) were clearly identified, and a cross soliton structure was found at the turning points of motion.

Zoom In Zoom Out Reset image size

Additional remarks about VR states are in order. As shown in figure 9(a) for the static case, VRs are unstable against decay modes with quantum numbers $(p=0,l=1)$. Our numerical results (open circles in the inset of figure 9(a)) show that this instability decreases at slow rate with increasing chemical potential, in agreement with the analytical prediction of [36] (solid red line in the inset):

$$\begin{eqnarray}&&\omega =-\displaystyle \frac{3{\hslash }}{\sqrt{2}}\displaystyle \frac{\mathrm{ln}[{\xi }^{2}({R}_{\perp }^{-2}+{\kappa }^{2}/8)]}{2{{mR}}_{\perp }^{2}},\end{eqnarray} \tag{ 43 }$$

where κ is the curvature of the vortex, ξ is the healing length, and ${R}_{\perp }$ is the Thomas–Fermi radius. This expression is valid for VRs in harmonically trapped elongated condensates within the Thomas–Fermi regime, and gives nonzero unstable frequencies ω in the limit of very high chemical potential, thus providing an estimation for the life time of VRs. For instance, for $\mu =21\;{\hslash }{\omega }_{\perp }$ both the numerical and analytical methods predict an unstable frequency $\omega \simeq 0.16\;{\omega }_{\perp }$ corresponding to a life time of about 20 ms for 50 Hz transverse trap.

In the case of moving solitons, the linear stability analysis follows essentially the preceding procedure for static states. Figure 9(b), generated for $\mu =5{\hslash }{\omega }_{\perp }$, shows our numerical result for the unstable frequencies of moving Chladni solitons as a function of the canonical momentum. The unstable frequencies decrease rapidly for VRs (blue lines with open triangles) of increasing speed (and thus their life time increases), and indeed they become stable past the bifurcation with solitonic vortices close to ${P}_{{\rm{c}}}=0.2\;\pi {\hslash }{n}_{1}$, (and ${P}_{{\rm{c}}}=1.8\;\pi {\hslash }{n}_{1}$) where the speed approaches the sound speed. As anticipated, there are also no unstable frequencies for solitonic vortices. As it is the case for the cross soliton, moving vortex dipoles (black dashed curves) present lower unstable frequencies than VRs. This may seem surprising considering the cylindrical symmetry of the system, and gives support to their possible detection in experiments.

6.2. Real time evolution

In order to check the predictions given by the linear stability analysis, we have also tested the nonlinear stability of Chladni solitons by the real time evolution of their stationary configurations. For, on the initial states ${\rm{\Psi }}({\bf{r}},t=0)$, we have added a random noise perturbation $\delta {\rm{\Psi }}({\bf{r}})$, which typically amounts to 2% of the wave function amplitude. Afterwards, we have allowed these wave functions to evolve in time, without dissipation, at constant chemical potential according to equation (6). For example, we have followed this procedure for the static Chladni solitons with $\mu =5\;{\hslash }{\omega }_{\perp }$, namely dark soliton, VR, cross soliton and SV, which have been previously characterised in figures 5 and 9(b). As expected, we have observed the decay of all solitonic states except the solitonic vortex, which, as a result, emerges at the final stage of the time evolution in all cases. Figure 10 summarises the decay processes, showing snapshots of the evolution at intermediate times. In particular it shows complex patterns localised in the plane of the initial stationary state at intermediate times and the emergence of a single solitionic vortex at late times, while some small amplitude radiation moves away from the solitary wave at the speed of sound.

Zoom In Zoom Out Reset image size

For higher values of the interaction parameter $\tilde{\mu }$, different scenarios can be found in the decay of Chladni solitons. In particular more than one solitary wave can appear and move away from the location of the initial unstable soliton. The interaction energy released from the parent state can transform into translational kinetic energy of a descendant solitary wave, leaving a simpler structure at the initial position. This is the case shown at figure 11, which displays two types of 'sling-shot' events, similar to the one recently observed in experiments and simulations with elongated BECs [38]. In the upper panels of figure 11, a static vortex–antivortex pair, after a fairly long time scale (∼50 ms), releases one of the vortices, whose acceleration can be noticed by its increasing curvature at later times. As another example, the lower panels of figure 11 follow the evolution of a static double ring that has been perturbed with a radial noise imprint preserving the cylindrical symmetry. One can observe how the external ring escapes from the original transverse plane, and, as in the previous case, increases its velocity.

Zoom In Zoom Out Reset image size

We have analysed the dynamical properties of static and moving Chladni solitons in cylindrically symmetric BECs, within the mean-field regime described by Gross–Pitaevskii equation. These states, strongly influenced by the geometry of the trap, emerge from the excitation of standing waves on planar kink states, and inherit particle-like features characterised by lower excitation energies and higher inertial masses than the kink. We have calculated numerically such quantities, and presented analytical expressions for their evaluation. The unstable standing waves producing the decay of the kink have been the object of a detailed analysis, and a formula for the prediction of the unstable frequencies has been proposed.

It is an interesting question to consider how a small deviation from the cylindrically symmetric trapping potential considered here will influence the nature and stability of Chladni solitons. Clearly, the continuous degeneracy of symmetry-breaking Chladni solitons (with $l\geqslant 1$) with respect to azimuthal rotation will be lifted. In addition, the near degeneracies of the different Laguerre-type unstable modes of the dark soliton (e.g. seen in figure 3) can be overcome with sufficiently strong azimuthal asymmetry, which will favour Hermite-type Chladni solitons with cartesian symmetries similar to the structures found at finite velocities and shown in figure 6. Since Laguerre and Hermite-type structures co-exist already in a symmetric trap where they have similar excitation energies and stability properties (see e.g. figure 9(b)), we do not expect the stability properties to change dramatically with small trap asymmetry. For more strongly asymmetric traps, we expect trap distortion to eventually destroy the vortex-ring-type structures and lead to chains of solitonic vortices with alternating sign in the limit of a near-two-dimensional trap. The further exploration of this interesting topic is left for future work.

The stability of Chladni solitons was studied by a linear stability analysis of the stationary states as well as by real-time evolution. Even though the recently observed SV and fast moving grey solitons are the only stable solitary waves in the strongly nonlinear regime ($\mu \gt 2.65{\hslash }{\omega }_{\perp }$), the more complex structures with crossing vortex lines are expected to be observable in current experiments with estimated lifetimes of tens of milliseconds, comparable to the previously observed VRs. In this regard, several procedures could be followed. In particular, in [16], we have proposed a feasible protocol for seeding a particular Chladni soliton on a planar kink. By means of a dark–bright soliton [39] in a two component condensate in the immiscible regime, a proper density and phase pattern could be imprinted on the bright soliton of one of the components occupying the kink depletion in the other component. The subsequent transfer of the selected pattern into the kink component, through a controlled Raman pulse, could serve the purpose of seeding the decay into the corresponding Chladni soliton. Other procedures with scalar condensates relying on an adequate trap geometry, have already been demonstrated. This is the case in [38], where vortex dipoles and the cross soliton has been identified after the decay of kinks in anisotropic harmonic traps. Very recently also the Φ soliton $(p=1,l=1)$ has been identified in what appears to be a seeded decay of a kink in a unitary Fermi gas [21]. This existing experimental evidence for Chladni solitons indicates that small amounts of dissipation due to trap losses and finite temperature are not fatal for the existence of these structures. It is an interesting question for further study whether Chladni solitons or analogous structures may exist as attractors in inherently dissipative superfluids, such as exciton–polariton condensates.

We acknowledge useful discussions with Martin Zwierlein, Shih-Chuan Guo, and Péter Jeszenszki.