Bound solitonic states in trapped multidimensional Bose-Einstein condensates

We report on the existence and stability of multidimensional bound solitonic states in harmonically-trapped scalar Bose-Einstein condensates. Their equilibrium separation, as a measure of the strength of the soliton-soliton or the solitonic vortex-vortex interaction, is provided for varying chemical potential $\mu$. Static bound dark solitons are shown to be dynamically stable in elongated condensates within a range of intermediate (repulsive) interparticle-interaction strength. Beyond this range the snaking instability manifests during the time evolution of the planar solitons and produces the decay into non-stationary vortex states. A subsequent dynamical recurrence of solitons and vortices can be observed at low $\mu$. At equilibrium, the bifurcations of bound dark solitons are bound solitonic vortices. Among them, both two-open and two-ring vortex lines are demonstrated to exist with both counter- and co-rotating steady velocity fields. The latter flow configurations evolve, for high chemical potential, into a stationary 3D-chain-shaped vortex and a three vortex-antivortex-vortex ring sequence that arrest the otherwise increasing angular or linear momentum respectively. As a common feature to the bifurcated vortex states, their excitation spectra present unstable modes with associated oscillatory dynamics.


I. INTRODUCTION
Solitary waves are localized, nondispersive excitations that can transport energy and momentum in a medium. They usually show a particle-like dynamics that features a characteristic interaction with other solitary waves. Both the type of solitary waves and their interactions are determined by the underlying medium. Among the latter, superfluid systems made of ultracold quantum gases stand up as excellent playgrounds for the study of solitary waves. In particular, in bosonic systems with repulsive interparticle interactions, Bose-Einstein condensates supporting dark solitons are readily generated in current experiments. They have been explored by means of various methods, for instance by phase imprinting techniques [1], by engineering the hyperfine states of the same atomic species [2], by producing the interference of two condensates [3], or by the Kibble-Zurek mechanism after crossing the phase transition to Bose-Einstein condensation [4]. The availability of plenty of experimental data has allowed, on the one hand, to confront the early theoretical predictions arisen from one-dimensional integrable systems [5,6], and, on the other hand, to go beyond in order to explore multidimensional solitary waves [7][8][9][10][11][12][13].
Dark solitons are unstable structures in systems with dimensions higher than one [14,15]. But these instabilities can be arrested by the presence of an external trapping. In this case, there generally exists a range of parameters at low interparticle interaction where the dynamical stability of multidimensional dark solitons can be ensured [16], and therefore where they can be experimentally realized [3]. Inside the trap, the spatial extension of a dark soliton increases with the interparticle interac-tion. Beyond an interaction threshold the excitation of long wavelength modes on the soliton surface is energetically favorable, by means of which a planar soliton decays into localized, long-life vortex lines [9][10][11]17]. Although this fact is a drawback for the study of multidimensional dark solitons, it has provided a useful mechanism for the generation of different solitary waves [18,19]. In this way, vortex rings have been observed after the decay of dark solitons in three-dimensional systems with isotropic trapping [7]. In elongated systems instead, the unstable planar solitons lead eventually to a single vortex line, a solitonic vortex [20].
In order to understand the interactions between solitary waves, the realization of states containing multiple solitary waves is required. To this end, states containing two solitons in scalar BECs have been previously considered in one dimensional [21,22] and quasi-onedimensional [23] settings. In these works it is assumed that the soliton decay has been suppressed by selecting a tight enough transverse trapping [15,16,19]. Under these conditions, it has been experimentally demonstrated [3] that several solitons can survive moving and colliding in elongated harmonic traps during long times (up to seconds).
In scalar condensates, the force experienced by a dark soliton due to the presence of another nearby soliton is repulsive when the inter-atomic interactions are local [24][25][26]. As a consequence, in the absence of external trap, there is no static solution of two bound solitons. Still, analytical solutions for two moving solitons have been found in one-dimensional (1D) settings [27,28]. However, when a harmonic confinement is acting, two static dark solitons conveniently situated at symmetric positions around the arXiv:1807.10843v2 [cond-mat.quant-gas] 23 Aug 2018 center of the trap can form a bound state. This equilibrium is due to the fact that the bouyancy-like force experienced by the solitons in the inhomogeneous density background (produced by the confinement) balances the repulsive interaction between them [23]. Interestingly, it has been recently reported that in systems with two overlapped BECs the inter-soliton forces can show either repulsive or attractive character depending on the sign of the inter-atomic interactions between condensates [29].
Bound two-soliton states in 1D trapped systems show dynamical instabilities in a regime of inter-atomic interactions that is adjacent to the non-interacting limit (see Fig. 1) [23]. However, if the inter-atomic interaction is high enough the instability is suppressed and the bound solitons become stable against small perturbations. This situation is expected to change in multidimensional systems, since the increasing interaction leads to larger extensions in the spatial dimensions of the soliton that eventually allow for the excitation of unstable modes.
Inter-vortex interactions has comparatively received much more attention, mainly in two-dimensional (2D) settings (see for instance [30][31][32][33][34][35]). Contrary to the soliton case, vortex interactions are long range due to the velocity field that they create. In the referred 2D systems, and in the limit of negligible vortex cores, it has been shown that vortices behave as charged particles, with the circulation of their velocity fields playing the role of charges (see e.g. [35]). For this reason, two vortices experience an attractive force when their circulation are opposite and a repulsive force when they have equal sense of circulation. Bound two-straight-vortex states have been observed in 2D systems with isotropic trapping [36,37]. Both corotating and counter-rotating pairs of vortices have been realized. While the former state does not stay static inside a circular trap, the counter-rotating pair (or vortex dipole) does. The latter configuration, have been reported to support dynamically stable states in traps with small anysotropic aspect ratios [38,39], and the typical vortex separation has been numerically computed for the pan-cake-shaped BEC [40]. A more complex dynamics is expected in three-dimensional (3D) systems, where vortex lines can bend in response to external forces [41][42][43][44][45][46].
In the present work we address the study of steady bound states made of two solitary waves in 2D and 3D elongated condensates inside a harmonic trap. In particular, we look for double nonlinear excitations along the weak axis of a system with isotropic transverse trapping. This is a common arrangement in current experiments, where moving solitary waves can be tracked [9,10,12,45,46]. Our starting point is the simplest excitation of the mentioned type that corresponds to a two-soliton state symmetrically situated around the center of the trap. As we will see, two vortex lines, with both straight and ring configurations and both co-and counter-rotating velocity fields, can be found as bifurcations for increasing values of the inter-atomic interaction. As special cases at high interaction, it is worth mentioning the 3D chain-shaped-vortex or the three vortex-ring configuration presented by the states in the families of two co-rotating straight vortices and two vortex rings, respectively. Regarding the stability properties, we report on small windows of chemical potential where stable twodark-soliton states can be found. Two-straight-parallel counter-rotating 3D vortices inherit this stability close to the bifurcation, and, beyond a region of oscillatory instabilities, again at higher interaction.
The rest of the paper is organized as follows. In section II we introduce the mean-field theoretical framework: the Gross-Pitaevskii equation for the condensate wave function, and the Bogoliubov equations for the linear excitations of the stationary bound solitons. Sect. III describes 1D bound solitons in harmonic traps. In Sect. III A we investigate the 2D systems that provide the main features of multidimensional bound solitonic states; we discuss the forces acting on the solitary waves and show characteristic time evolutions of two-bound-dark solitons. In Sect. IV we study the 3D bound dark solitons and bound vortex lines in elongated condensates; we elaborate on the excitation frequencies, soliton bifurcations, and decay dynamics. To sum up, we present our conclusions and perspectives for future work in section V.

II. MEAN FIELD MODEL
In the mean field regime, the dynamics of a scalar Bose-Einstein condensate (BEC) can be accurately described by the Gross-Pitaevskii (GP) equation for the wave function Ψ(r, t): where g = 4π 2 a/m is the interparticle interaction strength characterized by the scattering length a and the particle mass m. We consider a BEC confined by a cylindrically symmetric harmonic trap; V trap = m( ω 2 ⊥ (x 2 + y 2 ) + ω 2 z z 2 )/2 in 3D, and V trap = m( ω 2 ⊥ y 2 + ω 2 z z 2 )/2 in 2D, with aspect ratio λ = ω ⊥ /ω z . The stationary states of Eq. (1) are Ψ(r, t) = exp(−iµt/ ) ψ(r), where µ is the chemical potential. The number of particles N is fixed by normalization Ψ * (r, t)Ψ(r, t) dr = N .
The dynamical stability of the stationary states is checked by introducing linear modes {u(r), v(r)} with energy µ ± ω around the equilibrium state, i.e. Ψ(r, t) = e −iµt/ ψ + ω (u e −iωt + v * e iωt ) . After substitution in the GP Eq. (1), and keeping terms up to first order in the modes, one obtains the Bogoliubov equations for the linear excitations of the condensate where H L = − 2 ∇ 2 /2m+V trap −µ, is the linear Hamiltonian. The existence of frequencies ω with non-vanishing imaginary parts indicates the presence of dynamical instabilities that can produce the decay of the stationary state.
For the sake of comparison with current experiments, in what follows we report on 2D and 3D condensates made of 87 Rb atoms in the hyperfine state |F = 1, m F = 0 . The corresponding scattering length is a = 5.14 × 10 −9 m. The numerical results are obtained for a transverse harmonic trap of frequency ω ⊥ = 2π × 200 Hz and different aspect ratios.

III. BOUND SOLITONIC STATES IN LOW DIMENSIONAL SYSTEMS
In the non-interacting regime, the second excited axial eigenmode of the 1D harmonic oscillator is the second Hermite polynomial π, which has two axial nodes symmetrically situated around the trap center at distance q ≡ z node = a z / √ 2. The nonlinear continuation (for g > 0) of this state keeps the topology (the same number of modes) and builds the family of two-bound-dark solitons in the harmonic trap. The equilibrium distance between nodes is maximum at g = 0. Two representative examples are depicted in the top panel of Fig. 1, containing data, density (lines) and phase (symbols), obtained from the numerical solution of the time-independent GP equation for chemical potentials µ =10 and 50 ω z .
In a system with repulsive inter-atomic interactions, dark solitons experience also repulsive forces between themselves. In the absence of external trap, it has been shown that this force decays exponentially with the soliton separation 2q as F int ∝ exp (−4q/ξ) [23,24], where ξ = / √ mg n is the healing length, and prevents the existence of bound states. In harmonically-trapped 1D systems, however, two equal dark solitons symmetrically situated at a distance z = q from the trap center find an equilibrium configuration. From a particle-like approach for the soliton dynamics [47,48], the soliton separation is determined at equilibrium by the balance between the inter-soliton force F int (z) = −16N s 2 exp (−4z/ξ)/mξ 3 and the buoyancy-like force F b (z) = N s mω 2 z z due to the inhomogeneous density background induced by the trap, where ξ is evaluated at maximum density, and N s (<0) is the number of particles depleted by the soliton. The equilibrium leads to a transcendental equation forq = q/ξ which, as can be seen in the bottom panel of Fig. 1, provides a very good estimate of the soliton separation as obtained from the numerical solution of the 1D GP equation. As the chemical potential (hence the inter-atomic interaction) increases, the healing length decreases and so does the soliton separation in absolute (harmonic oscillator, a z ) units. In healing length units, however, the inter-soliton distance q increases with the chemical potential, which reflects a lower underlying force of solitonsoliton interaction.

A. 2D Bound solitonic states
The family of 1D two-bound-dark solitons contains dynamically unstable states for low values of the chemical potential (see the inset in the top panel of Fig. 1), where out of phase excitation modes break the static configuration [23]. In multidimensional systems new instabilities appear. The equilibrium states are unstable against the bending of the dark soliton stripe (in 2D) or the dark soliton plane (in 3D) when their extensions are long enough to support long wavelength (above a few healing lengths) transverse modes [16]. However dynamically stable bound solitons can still be found in multidimensional settings, just in between the end of the out-ofphase 1D instability and the beginning of the snaking instability. As an example, the top panel of Fig. 2 shows a robust 2D state with µ = 2.2 ω ⊥ in a trap with aspect ratio λ = 10. After seeding a perturbative small amount of white noise on the stationary state, the equilibrium configuration is preserved for a long time evolution. This situation changes for higher chemical potentials (or equivalently for higher inter-atomic interactions), where the snaking instability can operate. To illustrate this process, we have prepared an initial 2D state with µ = 3 ω ⊥ having two dark solitons at their equilibrium distance. As can be seen in the bottom panel of Fig. 2, (again after adding a small white noise on the initial state) the real time evolution shows the soliton decay into counter-rotating vortices (inset (b)) that os- cillate around their center of mass. During the motion along the weak trap z-axis, the vortices transit through regions of lower local chemical potential where there is not enough energy to support vortex structures. As a result, close to the turning points, the vortices smoothly transform back into dark solitons (inset (c)). The opposite behavior, soliton to vortex conversion, occurs in the proximity to the trap center. This phenomenon of dynamic inter-conversion between stationary states [38,49] is another instance of a nonlinear recurrence in the GP equation that reflects the dynamical instability of the involved stationary states [50]. However, at higher chemical potential, away of this oscillatory instability, robust states made of counter-rotating vortices can be found again [38].
In elongated condensates, it has been shown that static states made of solitonic vortex lines bifurcate from a single dark soliton excitation [19]. A priori, one could expect a similar scenario to happen by starting with two dark solitons, but this is only the case if the dark solitons are far away from each other. Otherwise, the interaction between the emerging solitary waves play a decisive role in selecting the possible stationary states in the trap. Due to the different nature of vortex-vortex interaction, the distance of equilibrium between two bound vortices is generally different from the distance between two bound dark solitons (which is at the origin of the oscillations observed in the lower panel of Fig. 2). In addition, the dimensionality of the system itself, modulating the underlying density profiles, makes also this distance to change even for the same type of solitary waves.
The effect of both factors, type of solitary wave and di-mensionality, on the solitary-wave separation in 2D systems is shown across the panels of Fig. 3. First of all, three new static bound-vortex states are shown to bifurcate for increasing chemical potential. As depicted in the right panel, they correspond to two counter-rotating vortices (or vortex dipole, VD), two co-rotating vortices (2V), and a double vortex dipole (2VD). The axial separations between the corresponding topological defects (plotted in the top left panel of Fig. 3) measure the associated force of interaction. Similarly to 1D solitons, co-rotating vortices present repulsive interactions in homogeneous-density backgrounds. The difference resides in the range of the interaction, which is long range for vortices (Coulomb type in 2D [35]). However, in elongated settings the solitonic character of the vortices makes their features localized [20] and the interaction turns essentially local [45]. As shown in Fig. 3, the separation distance decreases with the increasing chemical potential (initially) faster along the co-rotating (2V) family than along the dark soliton family. Therefore the vortex-vortex interaction, for given chemical potential in this regime, is weaker than the soliton-soliton interaction. For higher chemical potential, a striking difference arises in the former family, which reaches a minimum vortex separation. This is due to the entry of a pair of new co-rotating vortices with opposite circulation to the original vortices. As can be seen in the profile 2V(a) of the right panel of Fig. 3, the entry is symmetric on the perpendicular bisector of the line joining the original vortices. As a consequence the angular momentum, which kept increasing with decreasing vortex separation, reaches a maximum value near the minimum separation. Note that the angular momentum is not a conserved quantity along each family of stationary solitonic states. Only at the bifurcation point, where two solitonic families merge, all the dynamical properties are common. Beyond this point, the states belonging to a particular family change generally their properties for varying chemical potential (which is the nonlinear continuation parameter).
The scenario is quite different for two bound counterrotating vortices (or vortex dipole, VD). In this case, the vortex-antivortex pair produces zero angular momentum, and presents attractive interactions on homogenous density backgrounds. Again, the finite size of both the system (due to the trap) and the vortex cores induces an effective repulsive interaction that allows for a bound configuration. Contrary to the states considered before, the vortex-antivortex separation increases with the chemical potential (see Fig. 3), denoting higher interaction forces. Interestingly, this behavior contrasts with the roughly constant (but slightly decreasing) vortex separation found in pan-cake geometries [40].
The effect of dimensionality on the bound state configuration (hence on the inter solitary wave interaction) is illustrated for bound dark solitons (DS) in the bottom left panel of Fig. 3. It can be seen that for a given chemical potential the inter-soliton distance increases with the number of dimensions. The cause is the inhomogeneous density profile along the transverse directions, which has larger healing lengths at lower local chemical potential and contributes with higher buoyancy forces to the overall interaction. More importantly, as we show below, dimensionality has striking consequences in 3D configurations due to the bending of vortex lines.

IV. 3D BOUND SOLITONIC STATES
It is convenient to start the search for families of bound solitary waves by considering the static state made of two bound dark solitons. In the non-interacting 3D case such state is given by where the axial part corresponds to the second Hermite function H 2 (z/a z ) of the 1D system. The continuation of this solution in the nonlinear regime gives the family of 3D bound solitons. Figure 4 shows the frequencies of unstable modes in their excitation spectrum, obtained from the numerical solution of the Bogoliubov equations (2). The modes are grouped by the azimuthal polarity l, a positive integer that indicates the number of nodal diameters on the transverse section. As can be seen, the 3D bound solitons inherit the 1D (l = 0) outof-phase instability [23], but it is strongly suppressed at high aspect ratios λ 1. More relevant for the bifur- cation of new solitonic states, there are transverse excitations leading to the so-called snaking instability that bends the soliton plane [15,16]. This instability appears beyond a chemical potential threshold (or equivalently an inter-atomic interaction threshold) that marks the excitation of the lowest energy transverse mode at the soliton planes. Such mode has l = 1 and presents a single transverse nodal line at each soliton plane [19]. It can be viewed as the linear superposition (with equal weight) of a transverse vortex and an antivortex with angular momentum L z = and L z = − respectively. The subsequent growth of this unstable mode produces a solitonic vortex [20] on the corresponding soliton planes. For this reason, the threshold for the excitation of the lowest transverse mode coincides with the bifurcation of a nonlinear wave made of solitonic vortices. The threshold increases with the condensate aspect ratio, tending to the limit case of no axial trapping, where it takes the value µ ≈ 2.65 ω ⊥ [19,51]. In between the two mentioned instabilities, for an intermediate range of the chemical potential, our results demonstrate that there exist dynamically stable states made of two bound dark solitons. The stability window enlarges with the trap aspect ratio due to the shift of the snaking instability threshold towards higher µ and the suppression of the 1D-out-of-phase instability. Therefore, the states made of two static bound dark solitons are susceptible of observation in current experiments.
Once the snaking instability starts to operate in a state with high enough chemical potential, the solitons decay into vortex lines. The vortices are generated from the excitation of the unstable transverse modes available at that value of µ, as depicted in Fig. 4. Apparently, this is the same scenario found for single dark solitons [19]. However, the interaction between the emerging solitary waves introduces additional features that lead to a different outcome. First, not all the unstable frequencies of a bound-two-soliton state are pure imaginary, which causes the decay into non-static, oscillatory states [52]. Some of the unstable modes with l = 0 (dashed lines) depicted in Fig. 4 belong to this set. Second, for given l the instabilities appear in pairs, corresponding to the even and odd axial parity of the associated modes, with the even parity modes emerging at lower chemical potential than the odd ones (that present an axial node). The prototypical example is the instability with l = 1, which gives rise to bifurcations of a solitonic vortex in each soliton plane (see Fig. 5). In an isolated dark soliton all the transverse directions with the two possible circulations of the solitonic vortex are degenerate. In two bound solitons the orientations of the two emerging vortices are parallel, and the axial parity breaks the degeneracy between the counter-rotating configuration (VD, with even parity and lower energy) and the co-rotating one (2V, odd parity). As can be seen in Fig. 4, the energy differences induced by the axial parity are reduced for larger aspect ratios. tary waves from bound two solitons in a system with λ = 5. The chemical potential of the solitonic states is shown versus interaction, which is parametrized by χ = N a/a ⊥ √ λ. The bifurcations coincide with the emergence of unstable modes of same azimuthal polarity l that possess pure imaginary frequencies (as shown in Fig. 4 Representative examples of bound solitary waves are shown in Figs. 6 and 7, obtained from the numerical solution of the 3D GP equation with interaction parameters χ = 8.94 and 44.7, respectively. The top panels represent the density isocontours (colored by phase) of the BEC at 5% of its maximum density, and the bottom panels depict the (non-dimensional) axial density profiles an 1 (z) = a dx dy |ψ(x, y, z)| 2 , after integration along the transverse coordinates. Similar features to those discussed for 2D systems can be observed in the twostraight-vortex states of Fig. 6, and in the two-vortexrings shown in Fig. 7, where the inter-vortex separation is clearly larger for the vortex dipole configurations (VD and VRD).
Further comments about the configurations of the states containing co-rotating vortices are in order. These states evolve, for increasing chemical potential, into complex 3D structures. Two examples are shown in Fig. 8  FIG. 9. Real time evolution of the two co-rotating-vortex state depicted in Fig. 6, after the addition of perturbative white noise. The side views (top row) and the inner axial views (bottom row) of the density isocontours (at 5% of maximum density) are shown, from left to right, at times t = 19, 21.5, 23.5 and 26 ms.
for the families of co-rotating open vortices (2V) and corotating vortex rings (2VR). As it was the case in 2D systems, these solitonic families present a minimum of the inter-vortex distance when measured on the axial density profile. However, the 3D states evolve with the chemical potential by increasing the bending of the initially straight vortex lines, in such a way that at the center of the condensate (thus at higher density) the inter-vortex separation is greater than close to the condensate surface (at lower density). Eventually, near the chemical potential value where the axial inter-vortex distance reaches the minimum, the end points of the two co-rotating vortices merge with two new perpendicular vortices (one at each end point of the initial vortices). This 3D configuration (visible in the perspective view (2V) of Fig. 8) reduces the otherwise increasing angular momentum for increasing inter-atomic interaction. From this point on, the axial density profile (represented with a solid red line in the bottom panel of Fig. 7) can not capture the features of the 3D structure, and a single thick density notch opens in the middle of the condensate.
In a similar way, the family of two co-rotating vortex rings incorporate new vortex lines (a vortex ring larger and with opposite circulation) at the trap center. In this case, the new configuration reduces the axial linear momentum generated by the approaching off center rings. The resulting vortex pattern (visible in the perspective view (2VR) of Fig. 8) shows a three vortex-antivortex-vortex sequence of rings.
A relevant difference with respect to the case of single solitary waves concerns the stability and dynamics of vortex states. We have only found stable states in the family of counter-rotating-strainght vortices (VD). Stability is somehow expected just after bifurcation, inheriting this property from the parent bound-soliton states [39]. However, the bound-vortex states also inherit an oscillatory instability at higher chemical potential. As a result, for instance at λ = 5, a state with µ = 2.45 ω ⊥ is stable while a state with µ = 3.25 ω ⊥ is not. After this instability region, stability is recovered at even higher chemical potential. The VD state with µ = 4 ω ⊥ , shown in Fig. 6, is a stable example. In general, the typical decay of bound states show a complex vortex dynamics of rebounds or reconnections [10,46]. We show an example of bound vortex decay in Fig. 9, where selected snapshots during the real time evolution illustrate the process. The initial stationary state corresponds to the two-co-rotating-vortex configuration shown in Fig.  6, upon which a peturbative white noise has been added. As can be seen, a small variations in the vortex positions are manifested at t ≈ 20 ms. This initial departure from equilibrium leads eventually to the breakdown of the stationary flow created by the vortices, that show strong bending (visible in the inner views at the bottom of Fig. 9). As a result, the particle density distribution becomes strongly distorted and the two solitonic vortices move apart.

V. CONCLUSIONS
In this work, we have analyzed static configurations of bound states made of solitary waves in harmonically trapped BECs. Motivated by a common setting in current experiments, we have focused on elongated, multidimensional condensates with realistic parameters. Therein solitonic states bifurcating from two bound dark solitons have been considered in 2D and 3D settings. We report on states containing either two-co-rotating or twocounter-rotating vortex lines having both a straight and a ring configuration. Among the families of solitonic states, only bound dark solitons and bound vortex dipoles have been found to be dynamically stable, and then feasible to experimental realization, in a range of small chemical potentials. In the unstable regime of two bound dark solitons, the emergence of unstable modes with pure imaginary excitation frequencies marks the bifurcation of the static vortex states. Contrary to the case of a single dark soliton, we have also found (genuine 3D) unstable modes with non-vanishing real frequencies that give rise to oscillatory dynamics.
The analyzed states shed light on the nature of the soliton-soliton and (solitonic) vortex-vortex interactions. Contrary to the untrapped systems, the inhomogeneous density profile induces repulsive inter-vortex interactions irrespective of the vortex circulations. In all cases, the associated repulsive forces balance the buoyancy forces dragging the solitons and vortices towards the trap center. In the range of chemical potential analyzed, counter rotating vortices show equilibrium distances that increase with the inter-atomic interaction. For co-rotating vortices this quantity shows a non-monotonic behavior that reflects a change in the vortex configuration: a 3D vortexchain is shaped in the family of open vortex lines in order to arrest the increasing angular momentum, whereas a new counter-rotating vortex ring is introduced at the trap center in the family of vortex rings to reduce the axial linear momentum.
Natural extensions of the present work are envisaged.
The evolution of the bound vortex families at higher chemical potential, or the search for additional analytical expressions for the equilibrium distances of solitons and vortices are issues that deserve further exploration.