Bound solitonic states in trapped multidimensional Bose-Einstein condensates

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

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 Burger et al. (1999), by engineering the hyperfine states of the same atomic species Becker et al. (2008), by producing the interference of two condensates Weller et al. (2008), or by the Kibble-Zurek mechanism after crossing the phase transition to Bose-Einstein condensation Lamporesi et al. (2013). 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 Tsuzuki (1971); Zakharov and Shabat (1972), and, on the other hand, to go beyond in order to explore multidimensional solitary waves Anderson et al. (2001); Dutton et al. (2001); Ginsberg et al. (2005); Becker et al. (2013); Donadello et al. (2014); Ku et al. (2016); Aycock et al. (2017).

Dark solitons are unstable structures in systems with dimensions higher than one Derrick (1964); Kuznetsov and Turitsyn (1988). 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 Muryshev et al. (1999), and therefore where they can be experimentally realized Weller et al. (2008). Inside the trap, the spatial extension of a dark soliton 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 Ginsberg et al. (2005); Becker et al. (2013); Donadello et al. (2014); Ku et al. (2014). 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 Komineas and Papanicolaou (2003a); Muñoz Mateo and Brand (2014). In this way, vortex rings have been observed after the decay of dark solitons in three-dimensional systems with isotropic trapping Anderson et al. (2001). In elongated systems instead, the unstable planar solitons lead eventually to a single vortex line, a solitonic vortex Brand and Reinhardt (2001).

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 Alfimov and Zezyulin (2007); Zezyulin et al. (2008) and quasi-one-dimensional Theocharis et al. (2010) settings. In these works it is assumed that the soliton decay has been suppressed by selecting a tight enough transverse trapping Kuznetsov and Turitsyn (1988); Muryshev et al. (1999); Muñoz Mateo and Brand (2014). Under these conditions, it has been experimentally demonstrated Weller et al. (2008) 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 Zhao and Bourkoff (1989); Pawlowski and Rzazewski (2015); Bland et al. (2015). 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 Akhmediev and Ankiewicz (1993); Gagnon (1993). 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 Theocharis et al. (2010). 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 Morera et al. (2018).

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. 1Theocharis et al. (2010). 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 Fetter (1965); Guilleumas and Graham (2001); Jezek et al. (2008); Middelkamp et al. (2011); Navarro et al. (2013); Calderaro et al. (2017)). 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. Calderaro et al. (2017)). 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 Neely et al. (2010); Freilich et al. (2010). 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) does. A more complex dynamics is expected in three-dimensional (3D) systems, where vortex lines can bend in response to external forces Aftalion and Riviere (2001); García-Ripoll and Pérez-García (2001); Crasovan et al. (2003); Aftalion and Danaila (2003); Serafini et al. (2015, 2017) .

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 Ginsberg et al. (2005); Becker et al. (2013); Serafini et al. (2015); Ku et al. (2016); Serafini et al. (2017). 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. Contrary to the case of solitary waves emerging from single axial excitations Muñoz Mateo and Brand (2014), we have not found dynamically stable states made out of two static vortices along the axial trap.

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.1 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 :


where is the interparticle interaction strength characterized by the scattering length and the particle mass . We consider a BEC confined by a cylindrically symmetric harmonic trap; in 3D, and in 2D, with aspect ratio . The stationary states of Eq. (1) are , where is the chemical potential. The number of particles is fixed by normalization .

The dynamical stability of the stationary states is checked by introducing linear modes with energy around the equilibrium state, i.e. . 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 , 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 Rb atoms in the hyperfine state . The corresponding scattering length is m. The numerical results are obtained for a transverse harmonic trap of frequency  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 . The nonlinear continuation (for ) 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 . 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 .

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 as Zhao and Bourkoff (1989); Theocharis et al. (2010), where 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 from the trap center find an equilibrium configuration. From a particle-like approach for the soliton dynamics Scott et al. (2011); Muñoz Mateo et al. (2015), the soliton separation is determined at equilibrium by the balance between the inter-soliton force and the buoyancy-like force due to the inhomogeneous density background induced by the trap, where is evaluated at maximum density, and (0) is the number of particles depleted by the soliton. The equilibrium leads to a transcendental equation for


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, ) units. In healing length units, however, the inter-soliton distance increases with the chemical potential, which reflects a lower underlying force of soliton-soliton interaction.

Figure 1: Features of stationary two-bound-soliton states in harmonically trapped 1D BECs. Top panel: Probability density (lines) and phase (symbols) for two values of the chemical potential . The inset shows the imaginary frequencies of unstable excitation modes in the two-bound-soliton family. Bottom panel: (half) soliton separation as a function of , measured relative to the energy of the second excited linear eigenstate . The analytical expression (3) (solid line) follows in very good approximation the numerical data from the solution of the GP equation (open symbols). The healing length (dashed line) is shown for comparison.

iii.1 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 Theocharis et al. (2010). 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 Muryshev et al. (1999). 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. 2 shows a robust 2D state with in a trap with aspect ratio . After seeding a perturbative small amount of white noise on the stationary state, the equilibrium configuration is preserved for a long time evolution.

Figure 2: Time evolution of bound dark solitons across the direction in 2D trapped systems with aspect ratio . Sections of the density profile at are represented versus time. Top panel: robust state against perturbations in a BEC with . Bottom panel: soliton-vortex decay in a BEC with . The three insets represent snapshots of the density profile on the - plane at three different stages of the evolution: (a) initial two-soliton state, (b) decay into two solitonic vortices and, (c) recurrence of dark solitons.

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 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 oscillate around their center of mass. During the motion along the weak trap -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 two stationary states is another instance of a nonlinear recurrence in the GP equation that reflects the dynamical instability of both stationary states Kuznetsov (2017).

Figure 3: Bound solitary waves in elongated 2D BECs. Left panels: half axial separation between solitary waves of different type (top) and between dark solitons in different dimensional settings (bottom) as a function of the chemical potential (measured relative to the energy of the second excited non-interacting eigenstate along the weak-trap -axis , where is the number of dimensions). Right panel: density in the - plane of the 2D bound solitary waves (see text) marked by open symbols in the top left panel.

In elongated condensates, it has been shown that static states made of solitonic vortex lines bifurcate from a single dark soliton excitation Muñoz Mateo and Brand (2014). 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 dimensionality, 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 Calderaro et al. (2017)). However, in elongated settings the solitonic character of the vortices makes their features localized Brand and Reinhardt (2001) and the interaction turns essentially local Serafini et al. (2015). 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 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. 3), denoting higher interaction forces.

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

Figure 4: Unstable frequencies of 3D two-bound-dark-soliton states classified by the azimuthal polarity (see text) of the corresponding excitation mode. Solid lines correspond to pure imaginary frequencies, whereas dashed lines correspond to complex frequencies. The instabilities appear by pairs for given according to the axial () parity of the modes. For larger aspect ratio the stability window increases. The labels at the top panel indicate the nonlinear bifurcations described in the text.

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 , with energy , where the axial part corresponds to the second Hermite function 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 , 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 () out-of-phase instability Theocharis et al. (2010), but it is strongly suppressed at high aspect ratios . 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 Kuznetsov and Turitsyn (1988); Muryshev et al. (1999). 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 and presents a single transverse nodal line at each soliton plane Muñoz Mateo and Brand (2014). It can be viewed as the linear superposition (with equal weight) of a transverse vortex and an antivortex with angular momentum and respectively. The subsequent growth of this unstable mode produces a solitonic vortex Brand and Reinhardt (2001) 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 Komineas and Papanicolaou (2003b); Muñoz Mateo and Brand (2014).

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.

Figure 5: Bifurcation of first 3D bound solitary waves in an elongated harmonic trap with aspect ratio . The chemical potential (in energy units of the transverse trap ) is represented against the interaction parameter . The bifurcation points (open symbols) coincide with the corresponding emergence of a linear excitation mode (see Fig. 4), with equal azimuthal polarity and pure imaginary frequency , on the soliton planes.
Figure 6: Density isocontours (at 5 of the maximum density and colored by phase) of bound solitary waves in a harmonic trap with aspect ration . From top to bottom: two dark solitons (DS), two co-rotating vortices (2V) and two counter-rotating vortices (or vortex dipole VD). The three states are stationary solutions of the full GP equation for fixed interaction parameter . The bottom panel depicts the non-dimensional axial density (see text) of these states.

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 Muñoz Mateo and Brand (2014). 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 Strogatz (2018). Some of the unstable modes with (dashed lines) depicted in Fig. 4 belong to this set. Second, for given 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 , 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.

Figure 7: Same as Fig. 6 for bound-vortex-ring states. The upper pictures (from left to right) show a semi-transparent perspective view, and the lower pictures (from top to bottom) show a side view of co-rotating vortex rings (2VR) and counter-rotating vortex rings (VRD), respectively, both at interaction . The bottom panel, showing the axial densities , includes the dark-soliton and co-rotating-vortex profiles (at equal interaction) for comparison.

Figure 5 shows the bifurcations of the first bound solitary waves from bound two solitons in a system with . The chemical potential of the solitonic states is shown versus interaction, which is parametrized by . The bifurcations coincide with the emergence of unstable modes of same azimuthal polarity that possess pure imaginary frequencies (as shown in Fig. 4), whereas the complex frequencies lead to oscillatory dynamics. Apart from the bound open-vortex families (VD and 2V), Fig. 5 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. 6 and  7, obtained from the numerical solution of the 3D GP equation with interaction parameters and , respectively. The top panels represent the density isocontours (colored by phase) of the BEC at of its maximum density, and the bottom panels depict the (non-dimensional) axial density profiles , 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. 6, and in the two-vortex-rings shown in Fig. 7, where the inter-vortex separation is clearly larger for the vortex dipole configurations (VD and VRD).

Figure 8: Same as Fig. 7 for states of co-rotating vortex lines at higher interaction. The open-vortex configuration (2V) correspond to (or , same as states in Fig. 7) and the vortex rings (2VR) to (or ).
Figure 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 and ms.

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 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. 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) shwows a three vortex-antivortex-vortex sequence of rings.

A relevant difference with respect to the case of single solitary waves is the absence of dynamically stable states made of bound vortex lines. The vortex-vortex interaction produces the decay of these stationary states into moving vortex lines in the trap that can show a complex dynamics of rebounds or reconnections Becker et al. (2013); Serafini et al. (2017). 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 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 two-counter-rotating vortex lines having both a straight and a ring configuration. Among the families of solitonic states, only bound dark solitons 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 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.


  1. S. Burger, K. Bongs, S. Dettmer, W. Ertmer, K. Sengstock, A. Sanpera, G. V. Shlyapnikov,  and M. Lewenstein, Phys. Rev. Lett. 83, 5198 (1999).
  2. C. Becker, S. Stellmer, P. Soltan-Panahi, S. Dörscher, M. Baumert, E.-m. Richter, J. Kronjäger, K. Bongs,  and K. Sengstock, Nat. Phys. 4, 496 (2008).
  3. A. Weller, J. P. Ronzheimer, C. Gross, J. Esteve, M. K. Oberthaler, G. Theocharis,  and P. G. Kevrekidis, Phys. Rev. Lett. 101, 130401 (2008).
  4. G. Lamporesi, S. Donadello, S. Serafini, F. Dalfovo,  and G. Ferrari, Nat. Phys. 9, 656 (2013).
  5. T. Tsuzuki, J. Low Temp. Phys. 4, 441 (1971).
  6. V. E. Zakharov and A. B. Shabat, Soviet Physics JETP 34(1), 62 (1972).
  7. B. P. Anderson, P. C. Haljan, C. A. Regal, D. L. Feder, L. A. Collins, C. W. Clark,  and E. A. Cornell, Phys. Rev. Lett. 86, 2926 (2001).
  8. Z. Dutton, M. Budde, C. Slowe,  and L. V. Hau, Science (New York, N.Y.) 293, 663 (2001).
  9. N. Ginsberg, J. Brand,  and L. Hau, Phys. Rev. Lett. 94, 40403 (2005).
  10. C. Becker, K. Sengstock, P. Schmelcher, P. G. Kevrekidis,  and R. Carretero-González, New J. Phys. 15, 113028 (2013).
  11. S. Donadello, S. Serafini, M. Tylutki, L. P. Pitaevskii, F. Dalfovo, G. Lamporesi,  and G. Ferrari, Phys. Rev. Lett. 113, 065302 (2014).
  12. M. J. Ku, B. Mukherjee, T. Yefsah,  and M. W. Zwierlein, Phys. Rev. Let , 045304 (2016).
  13. L. M. Aycock, H. M. Hurst, D. Genkina, H.-I. Lu, V. Galitski,  and I. B. Spielman, PNAS  (February, 2017), 10.1073/pnas.1615004114.
  14. G. Derrick, Journal of Mathematical Physics 5, 1252 (1964).
  15. E. A. Kuznetsov and S. K. Turitsyn, Sov. Phys. JEPT 67, 1583 (1988).
  16. A. E. Muryshev, H. B. van Linden van den Heuvell,  and G. V. Shlyapnikov, Phys. Rev. A 60, R2665 (1999).
  17. M. J. Ku, W. Ji, B. Mukherjee, E. Guardado-Sanchez, L. W. Cheuk, T. Yefsah,  and M. W. Zwierlein, Phys. Rev. Lett. 113, 065301 (2014).
  18. S. Komineas and N. Papanicolaou, Phys. Rev. A 67, 023615 (2003a).
  19. A. Muñoz Mateo and J. Brand, Phys. Rev. Lett. 113, 255302 (2014).
  20. J. Brand and W. P. Reinhardt, J. Phys. B: At. Mol. Opt. Phys. 34, 4 (2001).
  21. G. L. Alfimov and D. a. Zezyulin, Nonlinearity 20, 2075 (2007).
  22. D. A. Zezyulin, G. L. Alfimov,  and V. V. Konotop, Phys. Rev. A , 013606 (2008).
  23. G. Theocharis, a. Weller, J. P. Ronzheimer, C. Gross, M. K. Oberthaler, P. G. Kevrekidis,  and D. J. Frantzeskakis, Phys. Rev. A 81, 063604 (2010).
  24. W. Zhao and E. Bourkoff, Optics letters 14, 1371 (1989).
  25. K. Pawlowski and K. Rzazewski, New Journal of Physics 17, 105006 (2015).
  26. T. Bland, M. Edmonds, N. Proukakis, A. Martin, D. O’Dell,  and N. Parker, Phys. Rev. A 92, 063601 (2015).
  27. N. Akhmediev and A. Ankiewicz, Phys. Rev. A 47, 3213 (1993).
  28. L. Gagnon, J. Opt. Soc. Am. B 10, 469 (1993).
  29. I. Morera, A. M. Mateo, A. Polls,  and B. Juliá-Díaz, Phys. Rev. A 97, 043621 (2018).
  30. A. L. Fetter, Phys. Rev. 138, A429 (1965).
  31. M. Guilleumas and R. Graham, Phys. Rev. A 64, 033607 (2001).
  32. D. M. Jezek, P. Capuzzi, M. Guilleumas,  and R. Mayol, Phys. Rev. A 78, 053616 (2008).
  33. S. Middelkamp, P. J. Torres, P. G. Kevrekidis, D. J. Frantzeskakis, R. Carretero-González, P. Schmelcher, D. V. Freilich,  and D. S. Hall, Phys. Rev. A 84, 011605(R) (2011).
  34. R. Navarro, R. Carretero-González, P. J. Torres, P. G. Kevrekidis, D. J. Frantzeskakis, M. W. Ray, E. Altuntas,  and D. S. Hall, Phys. Rev. Lett. 110, 225301 (2013).
  35. L. Calderaro, A. L. Fetter, P. Massignan,  and P. Wittek, Phys. Rev. A 95, 023605 (2017).
  36. T. W. Neely, E. C. Samson, A. S. Bradley, M. J. Davis,  and B. P. Anderson, Phys. Rev. Lett. 104, 160401 (2010).
  37. D. V. Freilich, D. M. Bianchi, A. M. Kaufman, T. K. Langin,  and D. S. Hall, Science 329, 1182 (2010).
  38. A. Aftalion and T. Riviere, Physical Review A 64, 043611 (2001).
  39. J. García-Ripoll and V. Pérez-García, Physical Review A 64, 053611 (2001).
  40. L.-C. Crasovan, V. Vekslerchik, V. M. Pérez-García, J. P. Torres, D. Mihalache,  and L. Torner, Phys. Rev. A 68, 063609 (2003).
  41. A. Aftalion and I. Danaila, Phys. Rev. A 68, 023603 (2003).
  42. S. Serafini, M. Barbiero, M. Debortoli, S. Donadello, F. Larcher, F. Dalfovo, G. Lamporesi,  and G. Ferrari, Phys. Rev. Lett. 115, 170402 (2015).
  43. S. Serafini, L. Galantucci, E. Iseni, T. Bienaimé, R. N. Bisset, C. F. Barenghi, F. Dalfovo, G. Lamporesi,  and G. Ferrari, Physical Review X 7, 021031 (2017).
  44. R. G. Scott, F. Dalfovo, L. P. Pitaevskii,  and S. Stringari, Phys. Rev. Lett. 106, 185301 (2011).
  45. A. Muñoz Mateo, A. Gallemí, M. Guilleumas,  and R. Mayol, Phys. Rev. A 91, 063625 (2015).
  46. E. A. Kuznetsov, JETP letters 105, 125 (2017).
  47. S. Komineas and N. Papanicolaou, Physical Review A 68, 043617 (2003b).
  48. S. H. Strogatz, Nonlinear dynamics and chaos: with applications to physics, biology, chemistry, and engineering (CRC Press, 2018).
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description