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

We analyse the dynamical properties of three-dimensional solitary waves in cylindrically trapped Bose-Einstein condensates. Families of solitary waves bifurcate from the planar dark soliton and include the solitonic vortex, the vortex ring and more complex structures of intersecting vortex-line known collectively as Chladni solitons. The particle-like dynamics of these guided solitary waves provides potentially profitable features for their implementation in atomtronic circuits, and play a key role in the generation of metastable loop currents. Based on the time-dependent Gross-Pitaevskii equation we calculate the dispersion relations of moving solitary waves and their modes of dynamical instability. The dispersion relations reveal a complex crossing and bifurcation scenario. For stationary structures we find that for $\mu/\hbar\omega_\perp>2.65$ the solitonic vortex is the only stable solitary wave. More complex Chladni solitons still have weaker instabilities than planar dark solitons and may be seen as transient structures in experiments. Fully time-dependent simulations illustrate typical decay scenarios, which may result in the generation of multiple separated solitonic vortices.


Introduction
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 [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 Bose-Einstein condensates 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, 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 vortex rings [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 Ref. [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 subsection 2.1, and the solitonic vortex and vortex ring in subsection 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 localized 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 realization of these solitary waves.

Energy and inertial and physical masses
Solitary waves often exhibit particle-like dynamics. One manifestation of such particlelike 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 Bose-Einstein condensates: The solitonic vortex is stable when µ/( ω ⊥ ) > 2.65 while the dark soliton is stable below this value, where µ is the chemical potential and ω ⊥ 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 [23,6]. Unstable solitary waves like vortex rings 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 (µ, 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 =Ż. Requiring the energy to be a constant of motion then leads to where 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 µ(Z) = µ 0 − 1 2 mω 2 z Z 2 , we arrive at Newton's equation of motion in the form where the physical mass defined by is to be interpreted characteristic parameter of the solitary wave that gives rise to the buoyancy-like force on the right hand side of Eq. (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 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 Bose-Einstein condensates with tight transverse confinement [27,22]. 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) ∝ sin(Ωt) with [26,28] Such harmonic oscillations have already been observed in experiments and the frequency ratio measured for dark solitons [29] and for solitonic vortices [30] in Bose-Einstein condensates and for solitonic vortices in the superfluid Fermi gas [23,6]. In the remainder of this section we present numerical data of dispersion relations and mass parameters evaluated for the dark soliton, solitonic vortex and vortex rings, and compare with approximate analytical expressions.

Planar dark solitons
We will model solitary waves within the mean field theory given by the Gross-Pitaevskii equation for the condensate order parameter Ψ(r, t) where g = 4π 2 a/m is the interaction strengh determined by the positive s-wave scattering length a and the bosonic mass m, V (r) = mω 2 ⊥ r 2 ⊥ /2 is an external, cylindrically symmetric, harmonic potential in the transverse coordinate r 2 ⊥ = x 2 +y 2 , and the condensate particle number N follows from normalization N = dr 3 |Ψ| 2 . Our starting point is the search for stationary solutions to Eq. (6), Ψ(r, t) = e −iµt/ ψ(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 Ref. [16] ψ(r) = χ T F (r ⊥ ) tanh(z/ξ(r ⊥ )), which is valid in the Thomas-Fermi regime, where χ T F = µ loc (r ⊥ )/g is the transverse ground state, µ loc (r ⊥ ) = µ − V (r ⊥ ) is the local chemical potential, and a local healing length is defined by ξ(r ⊥ ) = / mµ loc (r ⊥ ). Employing Eq. (7) as initial ansatz, the numerical solution of Eq. (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 ψ 0 , by means of with the energy defined by the functional Substituting the ansatz (7) into Eq. (8) and neglecting the derivatives of χ T F (r ⊥ ) in Eq. (9), according to the Thomas-Fermi approximation, we get the energy of a planar dark soliton confined by a transverse harmonic trapping as where the quantities with tilde are measured in the characteristic units of the trap, µ = µ/ ω ⊥ ,ã = a/a ⊥ , and a ⊥ = /mω ⊥ . Normalization of Eq. (7) gives the missing number of particles in the soliton N s = N − N 0 : which can also be obtained from the relation Figure 1 shows a comparison of the analytical expressions (10) and (11) derived from the ansatz (7), and the exact numerical solution to Eq. (6). As expected, the agreement improves with increasing values of the chemical potential, and forμ > 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ãẼ s andãN s are fixed, as reflected in equations (10) and (11). A further improvement can easily be introduced in previous expressions for average quantities. Following Ref. [31], by properly incorporating the zero point energy of the harmonic oscillator in the calculations we obtain the improved expressions and which have the correct limits both in the Thomas-Fermi and in the quasionedimensional regimes. They also interpolate in between with very good accuracy as can be seen in Fig. 1.

Solitonic vortex and vortex ring
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 3 reports the resulting mass ratios for dark solitons (red curve), solitonic vortices (green curve with open squares), and vortex rings (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 vortex rings. As will be explained in later sections, at these points the kink solutions exist only for the static dark soliton, and moving vortex rings emerge with the same value of the chemical potential. Approximate formulas for the solitonic vortex properties from hydrodynamic theory in logarithmic accuracy appeared in Ref. [6]. The energy expression for the stationary solitonic vortex in a BEC reads and the Thomas Fermi approximation has been used for the one-and two-dimensional densities, in particular 4an 1 =μ 2 . The missing particle number is obtained by differentiation where γ = µ n ∂n ∂µ is a polytropic index characterising the equation of state, which evaluates to γ = 1 for the case of a BEC. For the inertial mass the following expression was obtained 3. Snaking instability of the dark soliton

Bogoliubov stability analysis
After finding stationary kink solutions to Gross-Pitaevskii equation (6), we look for elementary excitations {u(r), v(r)} with angular frequency ω around every equilibrium state ψ with chemical potential µ. The perturbed state can be written as Ψ(r, t) = e −iµt/ ψ + ω (u e −iωt + v * e iωt ) , and the excitation modes are the solutions to the linear Bogoliubov equations, where . Dynamical instabilities are related to solutions to Eq. (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 Eq. (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.

Approximate separation of variables
For a real-valued stationary state, as it is the case of the dark soliton or kink of Eq. (7), the Bogoliubov equations (21) can be transformed using f ± (r) = u(r) ± v(r) and For ω = 0 the two equations decouple. This was the problem solved in Ref. [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 Eq. (22) with finite imaginary or complex values of ω. Aiming at an approximate separation of transverse and longitudinal degrees of freedom we introduce the rescaled variablez = z/ξ(r ⊥ ).
For the dark soliton state of Eq. (7) we can write gψ 2 = µ loc (r ⊥ ) tanh 2z , where µ loc = µ − V (r ⊥ ) and may thus rewrite the Bogoliubov equation (22) as with The operators A ± 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 ± , since we expect the unstable modes to be localised near the dark soliton plane. The ground state eigenfunction of A + is ϕ 0 2 sech 2z 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 ϕ 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 Ignoring any x, y derivatives of the functions ϕ ± (z) and projecting onto the respective ground states by multiplying from the left with (ϕ 0 + , 0) and (0, ϕ −1/2 − ), and integrating over z, we obtain the matrix equation where 962 is a numerical constant close to one.

Homogeneous background
The transverse Bogoliubov equation (25) is easily solved in the absence of a transverse trapping potential, where µ loc = µ = 2 mξ 2 = const, with plane wave solutions, e.g., For the unstable eigenvalues we find For small q the growth rate Im(ω) grows linearly with the slope 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 Refs. [34,13,33] , where the exact slope for the Gross-Pitaevskii equation is The snaking instability is suppressed and eigenvalues become real-valued for wave numbers larger than q crit = 1/ξ, which is the exact value. For intermediate values 0 < q < 1/ξ the growth rate has previously only been obtained numerically, and Eq. (26) reproduces the results of Refs. [34,13] very closely.

Harmonically trapped kink state
We now proceed with solving the transverse Bogoliubov equation (25) for an isotropic transverse trapping potential with µ loc (r ⊥ ) = µ − 1 2 mω 2 ⊥ r 2 ⊥ . Assuming an azimuthal dependence ∝ cos(lθ) with the quantum number l = 0, 1, 2, . . . reduces Eq. (25) to a set of ordinary differential equations in the radial coordinate r ⊥ . It is now convenient to move to harmonic oscillator units. Introducing the rescaled radial coordinatẽ r ⊥ = r ⊥ /a ⊥ and dividing the equation by ω ⊥ , we obtain where where H l 1 represents a two-dimensional Laplacian and H l 2 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 Even though this represents a non-hermitian eigenvalue problem, we have only found real eigenvaluesω 2 ζ 2 in numerical investigations.

Chladni soliton bifurcation points
Solutions of Eq. (27) withω = 0 have special significance as they indicate the transition of a specific mode from representing an instabilityω 2 < 0 to a stable small amplitude oscillationω 2 > 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 branchoff points for Chladni solitons. The solutions are found in terms of the scaled Fock-Darwin radial eigenfunctions χ l p (r ⊥ ) = 2 − 1 4 R l p (2 − 1 4r ⊥ ) of the 2D harmonic oscillator Hamiltonian H l 2 with eigenvalues l p = (2p + l + 1)/ √ 2 and where L l p (x) is the generalised Laguerre polynomial. It is easily seen that χ l p solves the Bogliubov equation (30) withω = 0 when l p =μ/2, which translates into the condition for the bifurcation points of Chladni soliton solutions from the dark soliton, as found previously in Ref. [16].

Finite growth rates
For finite instability rates Im(ω) the eigenvalue equation (30) can be expanded in a basis of the normalised eigenfunctions |p, l) ≡ χ l p of H l 2 , which transforms it into a tridiagonal matrix eigenvalue equation with B l p,p = (p, l|H l For the matrix elements we find and B l p,p = 0 for |p − p | > 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ω l n ∼μ for largeμ and values of the azimuthal quantum number l > 1. All these eigenvalues have a zero crossing for a finite value ofμ corresponding to Eq. (32), where n = p can be identified. These results were obtained by diagonalising a truncated matrix B l p,p with p, p < p c . Changing the cutoff value p c affects the large-μ regime but leaves the zero crossings for n p c unaffected. The asymptotic behaviour reported above represents the limit of p c → ∞.
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 Eq. (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 Eq. (33) in a truncated basis with 0 < p < p c produces the asymptotic behaviour where still n = p can be identified at the zero crossings ofω. The results of the truncated eigenvalue problem (33) with the asymptotic results (37) and (38) closely resemble the full numerical results.

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 Fig. 3. The insets represent axial views of phase-colored 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 localized around the plane of the kink. Their emergence follow in a very good approximation the analytical prediction for bifurcations given by Eq. (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 comes into existence at µ 0,1 = 2.65 ω ⊥ , very close to the 2.8 ω ⊥ 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 pl = 2p + l + 1, with the characteristic degeneracy of the two-dimensional harmonic oscillator.

Stationary Chladni solitons
Every unstable mode of the kink is associated with a stationary Chladni soliton. The zero crossings of the unstable mode frequencies in Fig. 3 indicate bifurcation points of the Chladni solitons from the dark soliton. The sign patterns of the unstable mode functions at the bifurcation points χ l p (r ⊥ ) cos(θ) 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 Fig. 3. Increasing the nonlinearity parameterμ = µ/( ω ⊥ ) 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 Fig. 4. For instance, the first excited mode | p = 0, l = 1 generates the family of solitonic vortex states, while the next two unstable modes, | 1, 0 and | 0, 2 , that  branch off near µ 1,0 = 4 ω ⊥ in Fig. 3, generate the families of vortex rings and two crossed vortices states, respectively. The linear solutions of a general, isotropic or anisotropic, twodimensional 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 can be constructed as superposition of Laguerre modes |p, l with degenerate energies nxny = pl . For example,|n x = 0, n y = 2 ∝ |p = 0, l = 2 − |p = 1, l = 0 , and |n x = 1, n y = 2 ∝ |p = 0, l = 3 − |p = 1, l = 1 . Because of the energy splitting of the bundle of linear degenerate states presented in Fig. 3, the linear superposition mechanism is not acting in the nonlinear case. However, we have found that Hermitelike 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 solitonic vortex, 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.

Moving Chladni solitons
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 vortex rings) [19,20] and for the solitonic vortex [18]. In order to construct the full dispersion relations for moving solitary waves, we numerically search for states Ψ(r, t) = e −i(µ t+mvzz)/ ψ(r), moving along the z-axis with a constant velocity v z = (0, 0, v z ), which are solutions of the stationary Gross-Pitaevskii equation for a co-moving reference frame: where µ = µ + mv 2 z /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 Ref. [35]: whereμ = µ/ ω ⊥ . 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 5b shows the excitation energy E s as a function of the axial velocity v z for moving solitary waves with chemical potential µ = 5 ω ⊥ . 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 | < 0.24 c. Below them, at lower excitation energies, only solitary waves bifurcating at energies less or equal thanμ = 5 can be found in the graph, as per Eq. (32). As can also be deduced from the unstable frequencies of Fig. 3, indeed, only vortex rings (blue solid line), double solitonic vortex (dashed black) and the single solitonic vortex (solid green) are available atμ = 5 and appear in Fig. 5. At small velocity, and very close in energy to the vortex ring and the static cross soliton, a new type of solution emerges (see the inset of Fig. 5b). It is composed of a couple of almost parallel vortices, indeed 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 6a shows the density configuration of these states around zero velocity, which can be cross-checked with their features extracted from the three panels of Fig. 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 Fig. 6b we show density isocontours for some of such modes withμ = 5, and 10.
It is also interesting to look at the phase jump along the axial z-direction ∆φ = φ(z → ∞) − φ(z → −∞) created by the different solitary waves (Fig.  5a). 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π as the velocity approaches +c, or alternatively reduces to zero as v z → −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 vortex rings (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 Fig. 5b, on the curve corresponding to the vortex-antivortex pair. Between turning points, labeled as B and C in Fig. 5c, there  . Dispersion curves (non-exhaustive diagram) for Chladni solitons with µ = 10 ω ⊥ , some of which correspond to the isocontour diagrams in Fig. 4 for static solitons, and in Fig. 8 for moving solitons.
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 Along with the axial physical momentum p z = −i dr(ψ * ∂ z ψ − ψ∂ z ψ * ), carried by the particles traversing the plane of the moving soliton, the canonical momentum includes the contribution coming from the phase jump ∆φ between the axial boundaries of the condensate where n 1 is the axial density of the ground state of the system [20,26]. Figure 5c shows the dispersion relation, excitation energy versus axial canonical momentum, for solitonic states with µ = 5 ω ⊥ 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: vortex rings 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 vortex ring family, here an axisymmetric solitary wave very much like a grey soliton with a vortex ring phase singularity outside the Thomas-Fermi density of the trapped BEC, as is apparent in Figs. 5a,b. In this regime, vortex rings 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, Fig. 7 displays some of the curves making the dispersion diagram for µ = 10 ω. For this value of the chemical potential, the family of double vortex rings [(p = 2, l = 0), violet lines] is available, and produces a new couple of turning points compared to the case for µ = 5 ω. Following the different curves away from their maximum (corresponding to zero velocity v z = ∂E s /∂P 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 µ = 10 ω. It is apparent how their symmetry changes when compared to their static counterparts in 4.

Stability of Chladni solitons
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 solitonic solutions (as those represented in Fig. 3). 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).

Linear analysis
As mentioned before, and aside from the stability of vortex rings moving near to the speed of sound, our results indicate that the solitonic vortex (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 solitonic vortex 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 vortex ring (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 Fig. 9a. The red curve with open circles corresponds to the unstable frequencies for vortex rings 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 (gray dashed lines). It is worth to mention that, as can be seen in Fig. 9a, for intermediate values of the chemical potential, between 4 and 5 ω ⊥ , the cross soliton (2SV, dashed black curve) has smaller values of unstable frequencies than vortex rings. This fact suggests that cross solitons are good candidates for being experimentally realised in elongated BECs. Vortex rings have already been observed in experiments [3,37]. In this regard, we have noticed that in Ref. [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.
Additional remarks about vortex ring states are in order. As shown in Fig. 9a for the static case, vortex rings are unstable against decay modes with quantum numbers (p = 0, l = 1). Our numerical results (open circles in the inset of Fig. 9a) show that this instability decreases at slow rate with increasing chemical potential, in agreement with the analytical prediction of Ref. [36] (solid red line in the inset): where κ is the curvature of the vortex, ξ is the healing length, and R ⊥ is the Thomas-Fermi radius. This expression is valid for vortex rings 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 VR 2SV DS Figure 10.
Density isocontours (at 5% of maximum density) during real time evolution showing the decay dynamics of a kink (DS), a vortex ring (VR) and a cross soliton (2SV), with µ/ ω ⊥ = 5 and vz = 0, obtained by solving numerically the time-dependent Gross-Pitaevskii equation starting from the stationary configuration seeded with a small amount of numerical noise. The typical cross-section radius is 3.5a ⊥ .
for the life time of vortex rings. For instance, for µ = 21 ω ⊥ both the numerical and analytical methods predict an unstable frequency ω 0.16 ω ⊥ 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 9b, generated for µ = 5 ω ⊥ , 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 vortex rings (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 c = 0.2 π n 1 , (and P c = 1.8 π 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 vortex rings. This may seem surprising considering the cylindrical symmetry of the system, and gives support to their possible detection in experiments.

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 Ψ(r, t = 0), we have added a random noise perturbation δΨ(r), which typically amounts to 2% of the wave function Sling shot events in the real time evolution of a 2SV state (up) with chemical potential µ/ ω ⊥ = 8, and a 2VR state (down) with µ/ ω ⊥ = 15. Density isocontours at 5% of maximum density are shown in both cases, with a transparent external surface in the axial view of the 2VR.
amplitude. Afterwards, we have allowed these wave functions to evolve in time, without dissipation, at constant chemical potential according to Eq. (6). For example, we have followed this procedure for the static Chladni solitons with µ = 5 ω ⊥ , namely dark soliton, vortex ring, cross soliton and solitonic vortex, which have been previously characterized in Figs. 5 and Figs. 9b. As expected, we have observed the decay of all solitonic states except the soltionic 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.
For higher values of the interaction parameterμ, 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 Fig. 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 Fig. 11, a static vortex-antivortex pair, after a fairly long time scale (∼ 50ms), releases one of the vortices, whose acceleration can be noticed by its increasing curvature at later times. As another example, the lower panels of Fig. 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.

Conclusions
We have analysed the dynamical properties of static and moving Chladni solitons in cylindrically symmetric Bose-Einstein condensates, 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 object of a detailed analysis, and a formula for the prediction of the unstable frequencies has been proposed.
The linear excitations and the real time evolution of Chladni solitons have also been addressed. Our results suggest that these states, even other than the recently realised solitonic vortex, could be observed in current experiments with ultracold gases, given that their lifetimes are expected to be of tens of milliseconds. In this regard, several procedures could be followed. In particular, in Ref. [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 Ref. [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].