Bound solitonic states in trapped multidimensional Bose-Einstein condensates
I. Morera M. Guilleumas R. Mayol A. Muñoz Mateo


 


Departament de Física Quàntica i Astrofísica, 
Facultat de Física, Universitat de Barcelona, E–08028
Barcelona, Spain
Institut de Ciències del Cosmos de la Universitat de 
Barcelona, ICCUB, E–08028 Barcelona, Spain

Departament de Física Quàntica i Astrofísica, 
Facultat de Física, Universitat de Barcelona, E–08028 Barcelona, Spain
Institut de Ciències del Cosmos de la Universitat de 
Barcelona, ICCUB, E–08028 Barcelona, Spain

Departament de Física Quàntica i Astrofísica, Facultat 
de Física, Universitat de Barcelona, E–08028 Barcelona, Spain
Institut de Ciències del Cosmos de la Universitat de 
Barcelona, ICCUB, E–08028 Barcelona, Spain

Departament de Física Quàntica i Astrofísica, 
Facultat de Física, Universitat de Barcelona, E–08028
Barcelona, Spain




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 μ. 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 
μ. 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. In spite of this,
the family of two-open counter-rotating vortices support dynamically stable 3D 
states.
  


03.75.Lm, 03.75.Mn, 67.85.-d

Bound states of dark solitons and vortices in trapped
multidimensional Bose-Einstein condensates
    A. Muñoz Mateo
    December 30, 2023
================================================================================================





§ 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 <cit.>, by 
engineering the hyperfine states of the same atomic species 
<cit.>, or by producing the interference of two 
condensates 
<cit.>. 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 <cit.>, and, on the 
other hand, to go beyond in order to explore  multidimensional solitary waves
<cit.>.

Dark solitons are unstable structures in systems with dimensions higher than 
one <cit.>. 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 
<cit.>, and therefore where they can be experimentally realized 
<cit.>. Inside the trap, the transverse extension of a 
dark soliton (measured on the surface of minimum density) increases 
with the interparticle interaction. 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 (with either open or ring 
configurations) <cit.>. 
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 <cit.>. In this way, vortex rings 
have been observed after the decay of dark solitons in three-dimensional 
systems with isotropic trapping <cit.>. In elongated systems 
instead, due to the presence of perturbative noise, the unstable 
planar solitons lead eventually to a single vortex line, 
a solitonic vortex that provides a perturbative amount of angular momentum 
<cit.>.  


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 
<cit.> and quasi-one-dimensional <cit.> 
settings. In these works it is assumed that the soliton decay has been 
suppressed by selecting a tight enough transverse trapping 
<cit.>. Under these conditions, it 
has been experimentally demonstrated <cit.> 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 <cit.>. 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 <cit.>. 
However, when a harmonic confinement is acting, two static dark solitons 
conveniently situated at symmetric positions around the 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 <cit.>. 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 <cit.>.

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. <ref>) <cit.>. 
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 
<cit.>). 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. <cit.>). 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 <cit.>. Both co-rotating 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) can. The latter configuration,  
have been reported to support dynamically stable states in traps with small 
anisotropic aspect ratios <cit.>, and the 
typical vortex separation has been numerically computed for the pancake-shaped 
BEC <cit.>. A more complex dynamics is expected in 
three-dimensional (3D) systems, where vortex lines can bend in response to 
external forces 
<cit.>.

In the present work we address the study of steady bound states 
made of dark solitons or vortices in 2D and 3D elongated condensates inside a 
harmonic trap. From now on we will use the terms bound solitonic 
states or bound solitary waves to denote such generic states. In 
particular, we look for double nonlinear excitations along the weak axis of a 
system with isotropic transverse trapping. The latter is a common arrangement 
in current experiments, where moving solitary waves can be tracked 
<cit.>. 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 two-dark-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, recover it at higher interaction. 

The rest of the paper is organized as follows. In section 
<ref> 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. <ref> describes 1D bound solitons 
in harmonic traps. In Sect. <ref> 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. <ref> 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 
<ref>.




§ 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 Ψ(𝐫,t): 


    iħ∂/∂ tΨ =(-ħ^2/2m ∇^2+V_ trap(𝐫)+g|Ψ|^2 )Ψ ,


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)+ω_z^2 z^2)/2 
in 3D, 
and V_ trap= m ( ω_⊥^2 y^2+ω_z^2 z^2)/2 in 2D,
with aspect ratio λ=ω_⊥/ω_z.
 The stationary states of Eq. (<ref>) are Ψ(𝐫,t)=exp(-i 
μ 
t/ħ) ψ(𝐫), where μ is the chemical potential. The 
number of particles N is fixed by normalization ∫Ψ^*(𝐫,t) 
Ψ(𝐫,t) d𝐫=N.


The dynamical stability of the stationary states is checked by introducing 
linear modes {u(𝐫),v(𝐫)} with energy 
μ±ħω around the equilibrium state, i.e. 
Ψ(𝐫,t)= e^-iμ t/ħ[ψ+∑_ω(u  e^-iω 
t+ v^*e^iω t)]. After substitution in the GP
Eq. (<ref>), and keeping terms up to first order in the modes, 
one obtains the Bogoliubov equations for the linear excitations of the 
condensate


    (  H_L + 2g|ψ|^2) u +
    gψ^2 v = ħω  u   ,
    
     
    -g(ψ^*)^2 u -(  H_L + 2g|ψ|^2) v
     = ħω  v  ,

 

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.

The stationary states have been numerically obtained by using a 
Newton continuation method based on a pseudo-spectral approach with an adapted 
basis of Fourier modes, along the axial z-coordinate, and Laguerre functions, 
in the transverse sections. The same pseudo-spectral approach has been used to 
solve the Bogoliubov equations from exact diagonalization. For the time 
evolution we have used a third-order Adams-Bashford scheme.
For the sake of comparison with current experiments, in what follows we
report on 2D and 3D 
condensates made of ^87Rb atoms in the hyperfine state 
|F=1,m_F=0⟩. The corresponding scattering length is a=5.14× 
10^-9m. The numerical results are obtained for a transverse harmonic trap 
of frequency ω_⊥=2π×200 Hz and different aspect ratios. 




§ BOUND SOLITONIC STATES IN LOW DIMENSIONAL SYSTEMS


It is instructive to start with the analysis of pure 1D systems.
In the non-interacting regime, the second excited axial eigenmode
of the 1D harmonic oscillator is the second Hermite polynomial
H_2=[2(z/a_z)^2-1]exp(-z^2/2a_z^2)/2√(a_z√(π)), which has
two axial nodes symmetrically situated around the trap center at distance q 
≡  z_node= a_z/√(2), where 
a_z=√(ħ/mω_z). The nonlinear continuation (for g>0) of this 
state keeps the topology (the same number of nodes) 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. <ref>, 
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(-4 q/ξ) <cit.>, where 
ξ=ħ/√(m g 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 <cit.>, the soliton separation is 
determined at equilibrium by the balance between the inter-soliton force 
F_int(z)=-16 N_s ħ^2 exp(-4 z/ξ)/m ξ^3 and
the buoyancy-like force F_b(z)=N_s m  ω_z^2 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 for q̃=q/ξ

    q̃ e^4q̃=(4μ/ħω_z)^2,

which, as can be seen in the bottom panel of Fig. <ref>, 
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 
soliton-soliton interaction.





 §.§ 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. <ref>), where out of phase excitation modes 
break the static configuration <cit.>. 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 <cit.>. However dynamically 
stable bound solitons can still be found in multidimensional settings, just in 
between the end of the out-of-phase 1D instability and the beginning of the 
snaking instability.
As an example, the top panel of Fig. <ref> 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. <ref>, 
(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 oscillate 
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 <cit.> is 
another instance of a nonlinear 
recurrence in the GP equation that reflects the dynamical instability 
of the involved stationary states <cit.>. However, at higher 
chemical potential, away of this oscillatory instability, robust states made of 
counter-rotating vortices can be found again <cit.>.


In elongated condensates, it has been shown that static states made 
of solitonic vortex lines bifurcate from a single dark soliton excitation 
<cit.>. 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. <ref>). 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 dimensionality, on the 
solitary-wave separation in 2D systems is shown across the panels of Fig. 
<ref>. 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. <ref>) 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 
<cit.>). 
However, in elongated settings the solitonic character of the vortices makes 
their features localized <cit.> and the interaction turns essentially 
local <cit.>.  




As shown in Fig. <ref>, 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. 
<ref>, 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 counter-rotating 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. 
<ref>), denoting higher interaction forces. 
Interestingly, this behavior contrasts with the roughly constant 
(but slightly decreasing) vortex separation found in pancake geometries 
<cit.>.

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. <ref>. 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.



§ 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 
Ψ_2(x,y,z,t)=exp(-ε_2 
t/ħ) exp[-(x^2+y^2)/2a_⊥^2]/√(2π a_⊥^2)H_2(z/a_z), with  
energy ε_2=ħ(ω_⊥ +ω_z/2), 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 <ref> shows the 
frequencies of unstable modes in their excitation spectrum, obtained from the 
numerical solution of the Bogoliubov equations (<ref>). 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) out-of-phase instability 
<cit.>, but it is strongly suppressed at high aspect ratios 
λ≫1.
More relevant for the bifurcation of new solitonic states, there are
transverse excitations leading to the so-called snaking instability that 
bends the soliton plane <cit.>. 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
<cit.>. 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 <cit.> 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ħω_⊥ <cit.>.

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. <ref>.
Apparently, this is the same scenario found for single dark solitons 
<cit.>. 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 <cit.>. Some of the unstable modes 
with l=0 (dashed lines) depicted in Fig. <ref> 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. <ref>). 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. <ref>, the 
energy differences induced by the axial parity are reduced for larger aspect 
ratios.


Figure <ref> shows the bifurcations of the first bound solitary 
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. <ref>), whereas 
the complex frequencies lead to oscillatory dynamics.
Apart from the bound open-vortex families (VD and 2V), Fig. 
<ref> includes the bifurcation of bound counter-rotating 
vortex rings in a dipole configuration (VRD), where for each element of 
vortex line in one of the rings there is another element in the parallel ring 
which is its mirror image.


Representative examples of bound solitary waves are shown in Figs. 
<ref> and  <ref>, 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 
two-straight-vortex states of Fig. <ref>, and in the 
two-vortex-rings shown in Fig. <ref>,
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. <ref> for the families of co-rotating open vortices (2V) 
and co-rotating 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. 
<ref>) 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. <ref>) 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. 
<ref>) shwows 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-straight vortices (VD). 
Stability is expected just after the bifurcation, inheriting 
this property from the parent bound-soliton 
states <cit.>. 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. <ref> is a 
stable example. In general, the 
typical decay of bound states show a 
complex vortex dynamics of rebounds or reconnections
<cit.>. We show an example of bound 
vortex decay in Fig. <ref>, 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. 
<ref>, 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. <ref>). As a result, the particle density distribution 
becomes strongly distorted and the two solitonic vortices move apart.




§ 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 two-counter-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 support dynamically 
stable states, thus feasible to experimental realization. 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 vortex-chain 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.



M. G. and R. M. 
acknowledge financial support from Ministerio de Economía y 
Competitividad (Spain), Agencia Estatal de Investigación (AEI) and Fondo Europeo de Desarrollo Regional (FEDER, EU) 
under Grants No. FIS2014-52285-C2-1-P and FIS2017-87801-P, and from Generalitat de Catalunya Grant 
No.  2017SGR533. 





