Engineering dark solitary waves in ring-trap Bose–Einstein condensates

The emerging field of atomtronics [1, 2] is associated with the creation of atomic circuit architectures based on ultracold atoms. A promising candidate for a closed prototype atomtronic circuit is based on laser-beam manipulation of ultracold atoms confined in toroidal geometries [3], a situation readily available in numerous laboratories [3–21]. Harnessing such circuits for technological applications (e.g. rotation sensors) requires a detailed understanding of the dynamics induced in such geometries through controlled perturbations, which has recently become very timely. Parallel to this, nonlinear excitations in the form of solitons could be useful for potential applications, e.g. due to their repetitive motion in a closed circuit and robustness against collisions.

The aim of this work is to demonstrate that although there are no known stable azimuthal solitonic solutions in toroidal geometries (somewhat related radial excitations in the form of 'ring dark solitons' have been discussed in [22–27]), soliton-like structures propagating at a fraction of the speed of sound and largely maintaining their shape after numerous collisions can nonetheless be generated through density engineering, even in the presence of thermal fluctuations.

In the context of ultracold atoms, solitonic nonlinear excitations arise spontaneously at the phase transition, as a consequence of quenching the system from the thermal to the condensed regimes through the Kibble–Zurek mechanism [28–32], or can be engineered by means of well known techniques such as phase imprinting [33–35], density engineering [33, 36–38], or a combination of both [33, 39]. The creation of solitons by engineering the density of the gas is typically performed by using a blue-detuned laser beam focussed in a narrow region of the system, on the scale of the healing length. The density distribution of the gas adapts to the presence of this perturbation, and the atoms are repelled from the region where the laser field is applied. Imposing a sharp density feature in a Bose–Einstein condensate determines a localised dip in the distribution which should then lead to the generation of solitons upon removal of the laser field. Solitons are one-dimensional (1D) objects originating from a unique balance between the kinetic energy, associated with spatial variations of the order parameter, and the atom–atom interaction energy; these waves largely preserve their shape after colliding with each other (undergoing only a phase shift). Although the initial engineered dark soliton experiments in harmonic traps led to both dynamical [40–47] and thermal [43, 48–50] instabilities, both long-lived [51] and stable [52, 53] solitons can now be routinely engineered in the lab.

In this work we demonstrate that long-lived structures resembling 1D dark solitons can also be engineered as counter-propagating pairs in ring-shaped traps within appropriate parameter regimes and excitation schemes, and thus study their stability, dynamics and interactions. More specifically, such structures are generated here numerically via the density engineering scheme, based on the (gradual) addition and removal of a Gaussian perturbation on a trapped Bose–Einstein condensate. To achieve optimal dynamical stability of such structures, one would need to restrict investigations to a very idealised regime of tight confinement in both transverse and radial directions, which effectively reduces the system to the 1D regime. Given the current significant experimental challenges in reaching this idealised 1D ring-trap regime, a pertinent question relates to how far from this regime one can deviate before dynamical instabilities dominate, severely limiting or even prohibiting solitonic behaviour, with a related question arising on the destabilising role of thermal fluctuations. Here, we demonstrate the existence of a broad experimentally relevant regime, where such engineered structures remain relatively robust both against dynamical and thermal instabilities, also surviving through multiple collisions. Our analysis highlights both the role of geometry and temperature in the evolution of such emerging solitary waves.

After discussing the system geometry and identifying relevant 'control parameters' (section 2), we focus on the question of optimisation of the generation of such solitary waves, by means of the Gross–Pitaevskii equation (section 3). Having identified optimum generation schemes, we then investigate the extent to which such structures could be obtained under realistic experimental conditions, in the presence of thermal fluctuations included here through numerical simulations of the stochastic Gross–Pitaevskii equation (SGPE) [54, 55] (section 4). The latter approach has already been demonstrated as an excellent model for ab initio equilibrium predictions of six independent quasi-2D [56] and quasi-1D [57, 58] Bose gas experiments, and has also been used to investigate condensate growth [54, 59] and dark soliton dynamics in quasi-1D geometries [60, 61], with the closely related stochastic projected Gross–Pitaevskii equation (SPGPE) [55, 62–64] used to study spontaneous defect formation following a quench [32, 65, 66], vortex dynamics [67], and decay of persistent currents [68]. Further details of dynamics following a nonoptimal choice of density engineering parameters and a comparison to the idealised 1D regime are discussed in two appendices.

We consider a trapped ultracold atomic gas (²³Na atoms, scattering length ${a}_{{\rm{s}}}=2.75\;{\rm{nm}}$) confined in a ring-shaped trap of the form (see figure 1(a) for a visual representation):

$$\begin{eqnarray}&&V({\bf{r}})={V}_{G}(1-{e}^{-2{({\rm{r}}-{{\rm{r}}}_{0})}^{2}/{w}^{2}})\end{eqnarray} \tag{ 1 }$$

where V_G, r₀ and w are respectively the depth, radius and $1/{e}^{2}$ half-width of the ring-gaussian potential. The radius, r₀, is the distance from the centre, where the potential reaches its minimum, whereas a length of $2w$ specifies the effective size of the ring channel wh ere the gas is confined.

Zoom In Zoom Out Reset image size

We assume tight transverse harmonic confinement with frequency ${\omega }_{\perp }$ in the direction perpendicular to the (x, y) plane, such that, for sufficiently small but realistic atom numbers, the gas can be brought into the transverse ground state, thus reducing all system dynamics to effectively two-dimensional.

We attempt to generate solitons via the density engineering technique, as this appears to be most relevant to recent experimental efforts [69]. Specifically we envisage perturbing an initial equilibrium density for a (quasi-2D) ring filled with a BEC by gradually ramping on the intensity of a blue-detuned laser sheet focussed in a localised region of the gas, and subsequently removing the laser. This method is modelled by adding (to the ring-trap of equation (1)) a narrow Gaussian potential of the form:

$$\begin{eqnarray}&&{V}_{{\rm{pert}}}({\bf{r}},t)={V}_{{\rm{L}}}(t){e}^{-{y}^{2}/2{\sigma }^{2}}(1-{\rm{\Theta }}(x)),\end{eqnarray} \tag{ 2 }$$

i.e. applied in the left half-plane (region across y = 0, for $x\lt 0$, see figure 1(a) for reference), where ${\rm{\Theta }}(x)$ is the step function, σ is the half-width of the Gaussian barrier and ${V}_{{\rm{L}}}(t)$ the time-dependent laser amplitude. The ramp needs to be kept on for a time long enough (few ${\rm{ms}}$—tens of ${\rm{ms}}$) to create a sufficiently deep notch in the density distribution, which should then lead to the generation of a (single) pair of counter-propagating solitary waves.

In order to minimise the linear (sound wave) excitations emerging from sudden perturbations, we follow a scheme similar to that used in the experimental work of [47], such that the perturbing potential is linearly ramped up over a time ${\tau }_{{\rm{on}}}$ to its maximum value, V₀, and then ramped down to 0 according to (perturbation on for times $0\leqslant t\leqslant {\tau }_{{\rm{pert}}}$ ):

$$\begin{eqnarray}{V}_{{\rm{L}}}(t) & = & {V}_{0}\left(\displaystyle \frac{t}{{\tau }_{{\rm{on}}}}\right)[{\rm{\Theta }}(t)-{\rm{\Theta }}(t-{\tau }_{{\rm{on}}})]\\ & & +{V}_{0}\left(\displaystyle \frac{1}{{\tau }_{{\rm{pert}}}-{\tau }_{{\rm{on}}}}\right)({\tau }_{{\rm{pert}}}-t)[{\rm{\Theta }}(t-{\tau }_{{\rm{on}}})-{\rm{\Theta }}(t-{\tau }_{{\rm{pert}}})],\end{eqnarray} \tag{ 3 }$$

where ${\tau }_{{\rm{on}}}$ the duration of the ramping-on sequence and $({\tau }_{{\rm{pert}}}-{\tau }_{{\rm{on}}})$ the ramping off timescale, which we have chosen to be relatively short.

Consistent with [47] we find that adiabatically ramping up the perturbation limits the sound emitted, thus leading to 'clean' profiles, without compromising the depth (or, equivalently, speed) of the emerging macroscopic excitations, which depends on the maximum value of ${V}_{0}(t)$. Although one could have used a slightly smoother (e.g. parabolic) turning on/off ramp to give less perturbation, our linear scheme seems to offer a sufficient reduction in background density noise facilitating our subsequent analysis. Thus, throughout this work, we show results for ${\tau }_{{\rm{on}}}=35.5$ ms and a ramping down time $({\tau }_{{\rm{pert}}}-{\tau }_{{\rm{on}}})=0.5$ ms. The particular shape of the perturbation is shown for the chosen parameters in figure 1(b).

The underlying system geometry and induced perturbation lead to three physically distinct sets of controllable parameters, respectively characterising: (i) the unperturbed ring trap potential $V({\bf{r}})$ (depth V_G, location of minimum of confined potential r₀, width of ring w) with the transverse frequency ${\omega }_{\perp }$ entering our analysis implicitly, through its contribution to the effective two-dimensional interaction parameter, ${g}_{{}_{{\rm{2D}}}}\propto \sqrt{{\omega }_{\perp }};$ (ii) the density-engineering perturbation ${V}_{{\rm{pert}}}$ (amplitude V₀, width σ, duration and slope of its application); and (iii) the details of the confined gas (atomic species, atom number N, characterised through the two-dimensional chemical potential μ, s-wave scattering length a_s, and temperature T). Although this cumulatively leads to a very broad parameter diagram to be probed, given a particular physical configuration and excitation scheme, the main physics is actually set only by a few parameters (or rather, their ratios).

To demonstrate this we choose to work with the particular experimental ring trap geometry of [70], namely: ${\omega }_{\perp }=2\pi \times 600\;{\rm{Hz}}$, ${r}_{0}=18.5\;\mu {\rm{m}}$, $w=9.45\;\mu {\rm{m}}$ and ${V}_{G}/{k}_{{\rm{B}}}=31.5\;{\rm{nK}}$ (Notice that in section 4 we work in the temperature range $1\;{\rm{nK}}\lt T\lt 10\;{\rm{nK}}$, i.e. $0.04\;{\rm{nK}}\lt {k}_{{\rm{B}}}T/{\hslash }{\omega }_{\perp }\lt 0.35\;{\rm{nK}}$). This, in turn, fixes the radial harmonic oscillator length, ${l}_{{\rm{r}}}=\sqrt{{\hslash }/m{\omega }_{r}}=2\;\mu {\rm{m}}$ (where ${\omega }_{r}=\sqrt{4{V}_{G}/{{mw}}^{2}}$) and the transverse spatial extent, ${l}_{\perp }=\sqrt{{\hslash }/m{\omega }_{\perp }}=0.86\;\mu {\rm{m}}$. Following [70], we introduce here a reference 'length unit' ${l}_{0}=\sqrt{({\hslash }/m{\omega }_{0})}=10\;\mu {\rm{m}}$ (using ${\omega }_{0}/(2\pi )=4.4\;{\rm{Hz}}$), to which all our results are scaled.

To maintain a quasi-2D geometry, suppressing instabilities due to coupling to dynamics outside the (x, y) plane, we work here with an atom number $N\approx 15\;000$ (which leads to a corresponding peak density of $2500{l}_{0}^{-2}$ corresponding roughly to 25 atoms$/\mu {{\rm{m}}}^{2}$). This choice ensures that the 2D condition $\mu \lt {\hslash }{\omega }_{\perp }$ holds: specifically, we choose a chemical potential $\mu \sim {\hslash }{\omega }_{\perp }/3$. For such an atom number we are typically in the 2D regime defined by ${l}_{\perp }\lt \xi \lt {l}_{{\rm{r}}}$, although a further dimensional reduction to a 1D regime (${l}_{\perp }\lt {l}_{{\rm{r}}}\lt \xi$) is theoretically feasible through a slight reduction (by a factor of 2–3) of the scattering length by means of Feshbach resonances.

To ensure all atoms remain confined within the ring trap, we also restrict the system temperature, T, to values sufficiently below ${V}_{G}/{k}_{{\rm{B}}}$. In general, the density engineering method can lead to the generation of one (or more) solitons [33]. To simplify the dynamics and avoid the generation of multiple pairs of counter-propagating structures, we thus choose the width of the perturbation $\sigma \approx \xi$ where ξ is the minimum value of the healing length, as calculated at the peak density [33].

Increasing the width of the perturbation σ to values $\sigma /\xi \gt 1$ (or significantly increasing V₀ for a fixed value of σ) leads to the generation of more than one pair of counter-propagating solitary waves, as discussed in appendix A.

Based on the choices described above, this effectively leaves us with the more manageable task of only 3 control parameters affecting the generation and subsequent propagation of the nonlinear excitations:

• The healing length of the system, ξ, broadly parametrising the spatial extent of the emerging macroscopic excitation (e.g. dark soliton or vortex): this is defined here as $\xi ={\hslash }/\sqrt{({{mg}}_{{}_{{\rm{2D}}}}{n}_{0})}$, where ${g}_{{}_{{\rm{2D}}}}=\sqrt{8\pi }({a}_{{\rm{s}}}/{l}_{\perp })({{\hslash }}^{2}/m)$ is the 2D interaction strength, ${a}_{{\rm{s}}}$ is the scattering length and n₀ refers to the maximum density. For a fixed transverse confinement ${\omega }_{\perp }$ investigated here, the value of ξ can be controlled either by changing the number of atoms in a given geometry, which affects the system density, or by varying the s-wave scattering length, e.g. by means of Feshbach resonances [71]: ultimately it is the product ${g}_{{}_{{\rm{2D}}}}{n}_{0}$ which controls the effective soliton width. For numerical convenience, when probing different parameter regimes, we choose to fix n₀ and vary ${g}_{{}_{{\rm{2D}}}}$ by increasing or decreasing the value of a_s by up to three times its background value (see subsequent figure 3).
• The maximum amplitude, V₀, of the density perturbation, which parametrises the overall depth of the imprinted density notch, scaled to the gas chemical potential μ.
• The effective width of the density perturbation which (for a given V₀) is parametrised by σ.

Our subsequent generation and dynamical stability analysis is thus primarily based on the chosen control parameters $({V}_{0}/\mu )$ fixing the depth of the emerging solitary wave excitations, and $({l}_{{\rm{r}}}/\xi )$ setting the effective dimensionality of the system.

In order to characterise the role of the relevant control parameters, and thus identify optimum regimes for solitonic generation in the idealised mean-field regime, we restrict our initial analysis to the (two-dimensional) Gross–Pitaevskii equation:

$$\begin{eqnarray}&&{\rm{i}}{\hslash }\displaystyle \frac{\partial \psi ({\bf{r}},t)}{\partial t}=\left(-\displaystyle \frac{{{\hslash }}^{2}}{2m}{{\rm{\nabla }}}_{x,y}^{2}+V({\bf{r}})+{g}_{{}_{{\rm{2D}}}}| \psi {| }^{2}-\mu \right)\psi ({\bf{r}},t)\end{eqnarray} \tag{ 4 }$$

The idealised proof-of-principle generation of quasi-stable solitary waves is best demonstrated through a series of density and phase snapshots following the removal of the perturbing potential. Characteristic images are shown in figure 1(c). More specifically, this figure shows (in situ) condensate density (top), renormalised ('carpet') density (middle) and phase (bottom) of the system at times $t=36,43,72,119,750\;{\rm{ms}}$ (left to right) after initiating the (ramped) perturbation, where $36\;{\rm{ms}}$ corresponds to the time that the perturbation is switched off, and 119 and 750 ms the times after one and thirteen collisions. The renormalised 'carpet' plots are obtained in the usual way, by subtracting from the perturbed instantaneous density the static density profile prior to the addition of the perturbation. Throughout this work, density is given in units of ${l}_{0}^{-2}$, where ${l}_{0}=10\;\mu {\rm{m}}$ is our reference 'length unit'.

Figure 1(c) reveals² the emergence of counter-propagating sound waves moving rapidly away from the region of the density perturbation, followed by two slower counter-propagating structures of reduced density, which additionally feature a pronounced phase slip across the density minima (bottom images). Such generated structures propagate in opposite directions within the ring, collide with each other at the far end of the ring and emerge largely unaffected after the collisions, as shown in the two rightmost frames of figure 1(c). More specifically, for the case considered here (with ${V}_{0}=2\mu$), the generated structures travel with a 'mean' velocity $v\approx 0.5c$ (based on the time taken to cross half the ring, i.e. $t\simeq 66\;{\rm{ms}}$), where c is the sound velocity in the medium, calculated at the peak density (i.e. at $r={r}_{0}$). We will thus infer that such structures are solitary waves. Further evidence for this is provided by their azimuthal 1D density cuts (subsequent figure 5) revealing excellent agreement with the anticipated analytical 1D soliton solutions, for a soliton propagating at the same speed, with this observation broadly extendable also to the $T\gt 0$ case (figure 7).

The profiles shown here are relatively 'clean', due to the gradual ramping of the perturbing potential, in stark contrast to equal duration potentials which are abruptly switched on and off, an example of which is shown in appendix A.

The intensity of the laser is an important control parameter (for a given ramping on/off sequence), as it characterises the maximum depth that the emerging nonlinear excitations can acquire. We anticipate requiring an intensity ${V}_{0}\gtrsim \mu$, although we note that much higher values would imply that the two BECs become effectively disjoint in the region of perturbation (an effect identified in [33] for elongated quasi-1D BECs). In our present work, we span intensity perturbations ranging from ${V}_{0}=0.5\;\mu$ to ${V}_{0}=10\;\mu$ within our ring trap geometry. As the symmetry of the Gross–Pitaevskii equation implies that the two emerging structures are mirror images of each other, we focus here on one of the emerging structures, arbitrarily chosen here as the one propagating clockwise.

As anticipated, higher values of ${V}_{0}/\mu$ lead to deeper (slower) solitary wave generation. To characterise this, figure 2(a) shows the dependence of the effective azimuthal propagation speed (scaled to the local speed of sound measured at r₀) on the maximum amplitude of the perturbing potential. As different structures travel at different speeds, and in order to avoid dependence on any initial excitations or related transient features (e.g. sound waves), we have chosen to characterise the propagation speed at the point when the emerging solitary waves reach the top of the ring, i.e. around x = 0. We find ${v}_{{\rm{s}}}/c\sim {({V}_{0}/\mu )}^{-\alpha }$, with a numerically extracted exponent $\alpha \approx 0.18$. Using the standard expression for pure 1D solitons in homogeneous settings [22], we can re-write this formula in terms of the depth, ${n}_{{\rm{sol}}}$, of the soliton from the peak of the unperturbed density, as ${n}_{{\rm{sol}}}/{n}_{0}\sim 1-{({V}_{0}/\mu )}^{-2\alpha }$, whose dependence is shown in the inset to figure 2(a).

Zoom In Zoom Out Reset image size

Our numerical analysis indicates that our excitation scheme leads to the initial generation of highly excited nonlinear structures, which gradually evolve towards more robust structures which we shall henceforth refer to as 'solitary waves'; nonetheless, such quasi-stable structures still feature some intrinsic dynamics. Thus, our subsequent analysis is further complicated by the fact that the emerging structures only approximately maintain their shape in time, in contrast to the case of a typical purely 1D soliton. Over longer timescales following the initial generation, we find that the curvature and closed geometry of the ring trap, which imply that the solitary waves which feature their own internal dynamics are continuously accelerated in their circular motion towards/against each other, actually leads to their gradual decay. Such decay manifests itself in the usual form of 'anti-damping' [72], i.e. growth of oscillation amplitude due to energy loss. Although this decay rate is relatively slow, it does imply that the apparent depth (or equivalently speed) of the solitary waves decreases (increases) with time in a (semi)-monotonic way. Figure 2(b) reveals some oscillations which could be attributed to a combination of the previously mentioned internal dynamics of the solitary waves, and their interactions with the propagating sound; the latter is somewhat reminiscent of (regular) oscillations induced by soliton-sound interactions in harmonically trapped quasi-1D BECs [46, 73]. However, as evident from figure 2(b), the decay of such structures in time is slow enough to allow for multiple collisions between the counter-propagating nonlinear structures, whose shape and speed appear to be only mildly perturbed by the collisions. The time of revolution for each pair of counter-propagating solitary waves (i.e. each value of ${V}_{0}/\mu$) is shown by vertical dotted lines in figure 2(b), with this figure spanning 5–7 revolutions. This suggests that the observed anti-damping may not be a direct consequence of the (head-on) collisions, but rather a complicated effect due to a combination of the internally excited state of the solitary wave and the related azimuthal-radial mode coupling, in conjunction with the accelerated circular motion through the ring and the interaction with the propagating background sound excitations.

To characterise the extent of such anti-damping numerically, we note that, after the emerging nonlinear excitations have completed ∼7 revolutions for the case ${V}_{0}=2\mu$ (intermediate blue diamonds in figure 2), their velocity has increased by a mere $\sim 15\%$ compared to the initial value. This can be taken as concrete evidence supporting our interpretation of such structures as (slowly decaying) solitary waves.

Having identified the potential for quasi-stable solitary-wave-like propagation in the ring, and the slow geometry- and interaction-induced underlying dissipation, there are two further main goals that we address in this work, namely the emergence of a regime where the solitary waves are reasonably stable and can be classified as 'solitonic', and the role of thermal fluctuations.

Throughout this work, the quasi-2D nature of the system (fixed by $\mu \lt {\hslash }{\omega }_{\perp }$, or equivalently ${l}_{\perp }\lt \xi$) implies that transverse excitations outside the (x, y) plane associated with 3D dynamical instabilities are suppressed. However, dynamical instabilities can also emerge in this two-dimensional geometry, depending on the ratio of the effective ring width ${l}_{{\rm{r}}}$ to the healing length ξ, whose effect is discussed next (with the idealised limit of ${l}_{\perp }\lt {l}_{{\rm{r}}}\lt \xi$ corresponding to an effective 1D regime of practically stable solitonic propagation over experimental timescales).

3.1. Solitonic behaviour and dynamical instabilities

Figure 1(c) clearly demonstrates that quasi-stable solitary wave propagation is possible around the ring; however it is important to further characterise such ensuing dynamics and identify regimes of rapid dynamical instabilities even in the 2D regime. By studying the dependence of the motion of the emerging solitary wave pairs on the ratio $({l}_{{\rm{r}}}/\xi )$, we can identify three reasonably distinct dynamical regimes over the broad range $0\lt {V}_{0}/\mu \lt 10$ of density perturbation amplitudes probed, with the corresponding 'phase diagram' shown in figure 3. Specifically, there are two very distinct regimes, respectively associated with (quasi-stable) 'solitonic' propagation and dynamically unstable 'snaking' behaviour, gradually mediated by an intermediate regime that we have termed 'shedding', due to the pronounced density emission from the stretched solitary wave. The insets to figure 3 illustrate characteristic snapshots identifying the key features of each of those regimes.

Zoom In Zoom Out Reset image size

Note that propagating structures resembling dark solitons have also been observed (but not analysed in detail) in parallel recent work [75], based on a similar experimental setup and perturbation scheme. In particular, numerical simulations reported in figure 8(left) of that paper show solitonic structures which appear to split into two, suggesting this could be somewhat analogous to the behaviour observed by us in the 'shedding' regime below.

In preparing the phase diagram (figure 3), we have chosen to probe the distinct regimes by varying the value of (${l}_{{\rm{r}}}/\xi$) for fixed ${l}_{{\rm{r}}}$. We choose to control the size of the healing length, ξ, by changing the value of the scattering length from its background value, a_s. In so doing, we have decided to work with a constant peak density of $2500\;{l}_{0}^{-2}=25$ atoms$/\mu {{\rm{m}}}^{2}$, which facilitates an easier comparison of the different emerging density profiles, rather than fixing the total atom number; in turn, this implies also adjusting the value of the chemical potential μ. More specifically, for the probed regime $0.8\lt {l}_{{\rm{r}}}/\xi \lt 2.3$, the scattering length spans the range $\approx [0.4{a}_{{\rm{s}}},3{a}_{{\rm{s}}}]$ (i.e. $1.1\;{\rm{nm}}\lt {a}_{{\rm{s}}}\lt 8.3\;{\rm{nm}}$ for ²³Na used here), while the chemical potential still satisfies the 2D criterion through the condition $0.25\lt \mu /{\hslash }{\omega }_{\perp }\lt 1$ (with the number of atoms lying in the range $\approx [12\;000,28\;000]$ respectively). We have also verified that changing ${\omega }_{\perp }$ (instead of a_s) in the range $\approx [0.2{\omega }_{\perp },9.2{\omega }_{\perp }]$, while still keeping the peak density ${n}_{{}_{{\rm{2D}}}}{({r}_{0})\propto ({\xi }^{2}{a}_{{\rm{s}}}\sqrt{{\omega }_{\perp }})}^{-1}$ fixed to the same constant value, yields the same physics.

The distinct regimes identified, separated by a crossover rather than a sharp boundary, are discussed further below, with reference to figure 4 showing detailed snapshots of the evolution of set times for all cases³ .

Zoom In Zoom Out Reset image size

'Solitonic' regime (figure 4(a)): this regime, an example of which was discussed in figure 1(c), is characterised by the propagation of quasi-stable solitary waves which extend practically over the whole width of the ring, and are thus classified here as 'solitonic'. The example demonstrated in figure 2 revealed a gradual decay of the emerging solitons, through a slow decrease in their depth, with their overall structure remaining largely unaffected. Such structures do however still exhibit some internal dynamics, mainly associated with coupling between radial and azimuthal degrees of freedom with increasing values of ${l}_{{\rm{r}}}/\xi$.

One can also internally sub-divide this regime, based on the value of $({l}_{{\rm{r}}}/\xi )$, since for values ${l}_{{\rm{r}}}\lt \xi$ one arrives at an effectively 1D geometry satisfying ${l}_{\perp }\lt {l}_{{\rm{r}}}\lt \xi$.

A comparison of density profiles between the case with ${l}_{{\rm{r}}}/\xi =0.8$ and our reference case of figure 1(c) (${l}_{{\rm{r}}}/\xi =1.3$) is shown in figure 5(a). Figure 5(a) shows 2D carpet plots of the solitonic structures when they are at the top of the ring, also plotting corresponding 1D density cuts in the radial (figure 5(b)) and azimuthal (figures 5(c) and (d)) directions. We find that, although for ${l}_{{\rm{r}}}\lt \xi$ the solitonic structures are symmetric along the radial direction about the middle of the ring, increasing the ratio ${l}_{{\rm{r}}}/\xi$, makes the structures less symmetric around the minimum of the ring trap (figures 5(a) and (b)), with the location of the density minimum shifted towards the outer edge of the trap; more specifically, by comparing the radial profiles (figure 5(b)), while for ${l}_{{\rm{r}}}/\xi =0.8$ the density minimum lies exactly at the point where the trap reaches its minimum $(x=0,y=1.85\;{l}_{0})$, in the case of ${l}_{{\rm{r}}}/\xi =1.3$ the minimum occurs at $(x=0,y=2.05\;{l}_{0})$ instead, which is presumably related to the excited solitonic dynamics seen in its subsequent evolution.

Zoom In Zoom Out Reset image size

The 1D azimuthal profiles for each case are also shown in figures 5(c) and (d). Comparing these to the anticipated analytical 1D soliton profile [72, 76, 77] for the same speed, we find excellent agreement, thus fully supporting our claim that such structures can be termed 'solitonic'. The restriction of radial excitations significantly enhances the solitonic stability, as further illustrated in appendix B.

'Shedding' regime (figure 4(b)): this is an intermediate regime in which the emerging nonlinear structures display pronounced internal dynamics at early times. A defining characteristic in this regime is that, following the removal of the perturbing potential, the initial azimuthal stretching of the propagating density depressions is balanced by a pronounced density re-arrangement, which results in the gradual separation of a significant density wave from the main depression, with the emitted density wave eventually dispersing: in some cases (smaller values of ${l}_{{\rm{r}}}/\xi$), the remaining structure partially 'recovers' towards a more shallow 'solitonic' profile spanning a significant fraction of the width of the entire ring radially, which is however less stable than those in the identified 'solitonic' regime; in other cases (higher values of ${l}_{{\rm{r}}}/\xi$) the width of the remaining solitary-wave excitation remains clearly less than the width of the ring. In both cases, such structures continue moving azimuthally (even if they only span a fraction of the radial ring width), and such 'solitary waves' appear to still survive multiple collisions. The emitted density waves can also be thought of as secondary shallower solitary waves, and their respective initial depths (and thus survival lifetimes) are increased with increasing $({V}_{0}/\mu )$, with the crossover region satisfactorily accounting for such behaviour. As the values of ${l}_{{\rm{r}}}/\xi$ increase beyond a certain thereshold region, the emerging structure stretches so much that it actually bends and breaks due to the background density gradient, in a manner reminiscent of the snaking instability observed e.g. in optical media [78] or elongated BECs [42].

'Snaking' regime (figure 4(c)): the nonlinear structures emerging after the removal of the perturbation deform so substantially along the ring, becoming dynamically unstable, as the azimuthal width of stretched density depressions greatly exceeds the healing length ξ, implying that solitary wave solutions can no longer be the lowest energy states of the system. Each of the two counter-propagating nonlinear structures then breaks into two 2D vortices (representing the planar mapping of 3D vortex rings), located near the inner and outer edges of the ring trap. Such dynamics is highly reminiscent of the observed 'snaking instability' in which 3D dark solitary waves decay into vortex rings [22, 44, 45, 47].

Having investigated in reasonable detail the role of the various 'geometrical' control parameters for the optimal generation of solitonic structures in ring-trap BECs, we now briefly address the important role of temperature and fluctuations on the form and lifetime of the emerging solitonic structures.

Temperature can be introduced into the Gross–Pitaevskii model in two closely related ways, by the controlled addition of fluctuations into the numerical simulations [79–82]. In the simplest approach (section 4.1), we start with an appropriately thermalised initial state, at some temperature T, described by a fluctuating classical field which is then propagated by the usual Gross–Pitaevskii equation (equation (4)). This approach is typically referred to as the 'classical field' method [81, 83–86] (being closely related to the finite temperature truncated Wigner [87–90]), and relies on the ergodicity of the Gross–Pitaevskii equation. Such a model has been used to study, among other phenomena, spontaneous soliton generation [91] and dark soliton stability [92, 93].

A more complete treatment of fluctuations requires both time-dependent stochastic (noise) fields and a dissipation term (with the two related through a fluctuation-dissipation relation [94]). In this case, both fluctuations and dissipation arise from the coupling of the stochastic classical field, representing the low-lying, highly occupied 'classical' modes of the system up to a cutoff to higher-lying (thermal) modes. The addition of the dissipation implies that the system relaxes (with a rate dictated by γ) to the equilibrium set by the heat bath parameters (temperature T and chemical potential μ).

Both approaches can atually appear as different limits of the SGPE [54, 94, 95], which in our current 2D setting takes the form [56]:

$$\begin{eqnarray}&&{\rm{i}}{\hslash }\displaystyle \frac{\partial \phi ({\bf{r}},t)}{\partial t}=[1-{\rm{i}}\gamma ]\left(-\displaystyle \frac{{{\hslash }}^{2}}{2m}{{\rm{\nabla }}}_{x,y}^{2}+V({\bf{r}})+{g}_{{}_{{\rm{2D}}}}| \phi {| }^{2}-\mu \right)\phi ({\bf{r}},t)+\eta ({\bf{r}},t),\quad \quad \end{eqnarray} \tag{ 5 }$$

where $\phi ({\bf{r}},t)$ now represents the multi-mode stochastic 'classical' field cumulatively describing the low-lying modes of the Bose gas (see also the closely related SPGPE [55, 81]). This should be directly contrasted to the usual Gross–Pitaevskii equation (equation (4)), where $\psi ({\bf{r}},t)$ denotes simply the condensate wavefunction. In equation (5), (thermal) fluctuations are mimicked by the presence of the noise term $\eta ({\bf{r}},t)$ which has Gaussian correlations of the form $\langle {\eta }^{*}({\bf{r}},t)\eta ({{\bf{r}}}^{\prime },{t}^{\prime })\rangle =2{\hslash }\gamma {k}_{{\rm{B}}}T\delta ({\bf{r}}-{{\bf{r}}}^{\prime })\delta (t-{t}^{\prime })$, where γ parametrises the strength of the noise and damping.

Although such an equation should be solved numerous times with different stochastic fields, with the results appropriately averaged, one can actually attribute an indirect physical interpretation to each numerical trajectory, as representing a plausible experimental run. For a discussion of the usefulness of single stochastic trajectory analysis and how to extract meaningful averaged parameters from this, see e.g. our earlier work on stochastic dark soliton dynamics in harmonic traps [60].

In order to investigate the effect of temperature on the soliton dynamics, we focus on our (largely sound free, 2D) reference case of figure 1(c), for which relatively deep solitons were clearly visible in the T = 0 limit. We now use those two different limits of the SGPE to discuss the ensuing soliton motion at finite temperatures, through indicative single-trajectory results.

4.1. 'Classical field' method

In this section we generate the 'initial' state, i.e. state prior to adding the density perturbation, as an appropriate thermal noisy equilibrium state, via dynamical equilibration of the Stochastic Gross–Pitaevskii equation (equation (5)). After equilibration, we switch off both dynamical noise and dissipation, which amounts to propagating our noisy thermalised initial state via the ordinary Gross–Pitaevskii equation. Representative images of equilibrium density profiles in the presence of fluctuations are shown⁴ in figure 6(a) for different temperatures of the unperturbed thermal state. As before, we show densities (top), renormalised densities or 'carpet' plots (middle) and phase (bottom) plots for $T=1,9,10\;{\rm{nK}}$ (left to right). The carpet plots (middle) are generated by subtracting from the single stochastic run perturbed density at a given time the corresponding T = 0 (mean field) unperturbed equilibrium result. As expected, the fluctuations in the background density increase with increasing temperature.

Zoom In Zoom Out Reset image size

A few comments are in order here:

• (i)  
  The stochastic numerical evolution leads to the generation of a different random phase in each numerical simulation, such that the underlying phase differs from run to run and temperature to temperature. To facilitate a more direct comparison of the soliton dynamics between the different temperature cases, we therefore numerically eliminate the initial random phase difference (i.e. the phase difference of the equilibrium configuration prior to turning on the perturbation) among the cases T = 1 and $9\;{\rm{nK}}$ shown here.
• (ii)  
  The T = 10 nK case we have chosen to show here is slightly different, as it contains a persistent current (here with a winding number 1) at our t = 0 time labelled as 'equilibrium', which is simply a reflection that the system has not yet actually fully equilibrated. This persistent current has appeared here spontaneously during our equilibration process (and will eventually decay after a sufficiently prolonged evolution). The reason for this appearance can be traced back to our method of generating the initial state, which is actually based on dynamical equilibration following an instantaneous numerical quench, a process which is known to support such spontaneous defect formation, in accordance with the Kibble–Zurek mechanism [96]. We have also checked that over numerous simulations we get a distribution of both flow-free solutions and persistent currents with positive and negative winding numbers, including also higher winding numbers, in qualitative agreement with experiments [13]. Persistent currents can also be generated in lower temperature cases, so our choice of displaying the 10 nK case here with a persistent current is because it yields a clean persistent current over strong background density fluctuations, thus providing clear evidence of a strikingly different motion around the ring, combining both persistent currents and density fluctuations. Note that 10 nK is also the highest (optimal) temperature we can realistically probe in our setup, to avoid atoms populating transversally excited modes, which are not accounted for in our purely 2D scheme. While instructive to show how the presence of the persistent current affects the generation/propagation of the solitary waves, and although we could ensure that the perturbing potential is added after the persistent current decays [68], we have chosen not to investigate this case further, since experiments aiming to generate solitons would also choose initial conditions without an intrinsic flow pattern.
• (iii)  
  As our simulations have been done at constant chemical potential, the atom numbers increase slightly with increasing temperature (up to 20%), but we do not expect this to have a significant effect on our presented analysis (other than, e.g., in making the speed of sound slightly temperature-dependent.)

The evolution of density and phase on these noisy initial states after the addition and removal of the perturbing potential is shown in figure 6(b). In all cases, analogously to the T = 0 case, we can detect two emerging structures which tend to propagate in opposite directions ($T=1,9\;{\rm{nK}}$) and appear to remain largely unchanged through their collisions. The profiles shown here are taken after the solitons have already undergone one and a half revolution (time t = 238 ms), such that they have met each other and interacted three times (at $x\approx 1.85{l}_{0}$, $x\approx -1.85{l}_{0}$ and again at $x\approx 1.85{l}_{0}$).

We have performed a study based on numerous individual classical field simulations, based on completely random initial conditions (generated through SGPE equilibration), and observe a range of features commented upon below:

As the temperature (represented by thermally induced background density inhomogeneity) increases, the motion of the two counter-propagating solitons reveals small differences (although the mean propagation speed remains approximately constant). We can attribute this to a combination of two effects (largely guided here by our earlier work on dark solitons [60]): on the one hand, the presence of random fluctuations in the initial state, implies that the emerging dark solitons are not identical; moreover, even though the average noise amplitude at each temperature is fixed, at any time each soliton is nonetheless propagating through a different random noisy background configuration, which introduces small random 'kicks' to the soliton motion around the ring. As a result, the two solitons do not collide exactly at y = 0, and at any given time their respective positions are not exactly mirror-images of each other (see e.g. T = 9 nK case in figure 6(b)). The fact that the generated structures appear in the same location after three consecutive interactions (and having done more than one full revolution in the ring) strongly suggests that the solitonic nature of such structures persists even in the presence of initial fluctuations.

Interestingly, we do not find a systematic net effect of temperature, i.e. the average position of the solitonic waves after few revolutions (averaged over ∼10 stochastic runs), is only mildly perturbed from the corresponding T = 0 case without displaying a clear dependence on temperature. This is a feature that we have also observed in our previous work on dark soliton dynamics in purely 1D geometries [60, 97]. This appears to be in partial disagreement to the findings of [92, 93], where it has been argued that dark solitons propagating on an initially fluctuating background exhibit some decay⁵ , whereas our work provides evidence of temperature-dependent modulations, but no net decay. While our analysis does not give (or intend to give) a conclusive answer to this issue, it does suggest that if classical field simulations correctly predict the soliton dynamics, then multiple soliton collisions should routinely arise in carefully engineered experiments.

The presence of the persistent current in the third (t = 0) and sixth (post-density-engineering) subplots of figure 6 imparts an additional flow velocity to the two solitons, thereby significantly speeding the motion of the co-flowing soliton, while simultaneously decelerating the soliton travelling against the flow. As a result, in the presence of a persistent current, the two solitons exhibit a net relative speed between them, and the motion in this case deviates significantly from that when no persistent current is present, where the mean soliton's x-coordinates were found to be approximately equal.

The role of fluctuations on the soliton density is illustrated in figure 7(a) which compares its form in the absence (T = 0) or presence (T = 1, 9 nK) of background fluctuations in the initial state. Although the actual position of the soliton in an individual numerical run jitters about the mean equilibrium position, and its profile becomes less well defined due to the underlying fluctuations, the fluctuations themselves do not appear to critically affect the underlying solitonic shape even after a few collisions (see rightmost image in figure 7(a)). To verify the solitonic nature of the emerging structures prior to any collisions, figure 7(b) plots their azimuthal 1D density cuts slightly after their generation (when located at the top of the ring) for finite temperatures, contrasting them to the pure T = 0 case. This figure clearly shows that although the fluctuations noticeably modify the density profiles, the underlying solitonic nature reflected by the central width of the density depression set by the healing length ξ remains clearly visible.

Zoom In Zoom Out Reset image size

Based on all above findings, we would thus argue that the solitonic nature appears to persist both in the initial and dynamical regimes, when modelling the nonequilbrium soliton dynamics on top of a fluctuating initial state.

4.2. Full stochastic evolution

To further improve on our earlier $T\gt 0$ predictions, we now consider the dynamics resulting from density engineering in the context of the dynamical SGPE, in which the classical modes of the system described by $\phi ({\bf{r}},t)$ exhibit full dynamical coupling to the high-lying modes of the system, i.e. maintaining here both dynamical noise $\eta ({\bf{r}},t)$ and dissipationγ.

The presence of γ ensures the system eventually relaxes to an equilibrium profile dictated by the heat bath T and μ, thus clearly leading to the gradual decay of any generated excitations, at a rate directly dependent on γ. Although γ is often considered a 'phenomenological parameter', an analytical prediction for its value does exist [80, 81, 97, 99, 100]. Importantly, recent work with the closely related SPGPE demonstrated excellent agreement between theoretical predictions based on the theoretically predicted γ value and experimental findings, in the context of persistent current decay in a ring trap [68], thus suggesting that such simple analytical estimates yield reasonably realistic values.

Using the predicted analytical expression [68, 99], leads in our system to an estimated value of $\gamma \sim {10}^{-5}$, which we use here simply as a guide. Figure 8 (left two images) shows snapshots of the post-density-engineering evolution of the solitonic structures, revealing the persistence of clearly identifiable solitonic structures even after 2 full revolutions (or 4 mutual collisions). Given the somewhat crude estimated values for the decay parameter γ, the two rightmost plots of figure 8 show the corresponding case with a much larger (heuristically chosen) value of $\gamma \sim {10}^{-2}$. Even in this case, which features enhanced soliton decay, we still find evidence of the (attenuated) solitonic structures surviving after at least one full revolution around the ring (t = 164 ms, corresponding to two collisional events), as shown in the third set of plots in figure 8.

Zoom In Zoom Out Reset image size

We thus conclude that although thermal excitations can significantly perturb the shape and reduce the lifetime of the solitonic excitations, their presence and collisions could be observable under realistic experimental conditions, provided the temperature is not too high. This is in qualitative agreement with previous discussions of soliton stability in elongated 3D harmonically trapped BECs [50].

We have investigated the conditions under which the addition of a carefully engineered density perturbation to an atomic Bose–Einstein condensate contained within a ring trap can controllably generate pairs of counter-propagating solitonic excitations, demonstrating that such structures should in fact survive (multiple) collisions and revolutions around the ring, even at finite temperatures.

Optimum experimentally relevant conditions for their observation include tight transverse confinement along the direction orthogonal to the plane of the ring (denoted by the frequency ${\omega }_{\perp }$) and small atom number (or equivalently chemical potential μ) such that the two-dimensional condition $\mu \lt {\hslash }{\omega }_{\perp }$ is satisfied, thus suppressing three-dimensional dynamical instabilities. Nonetheless, the azimuthal and radial degrees of freedom can still couple with each other, and a form of dynamical 'snaking' instability was found to persist even in two-dimensional geometries (${l}_{\perp }\lt \xi \lt {l}_{{\rm{r}}}$), unless the effective radial ring length ${l}_{{\rm{r}}}$ satisfied ${l}_{{\rm{r}}}\lesssim 1.5\xi$, where ξ denotes the healing length of the gas determining the soliton width. Although experimentally challenging, a further reduction in the radial width of the ring trap, such that ${l}_{{\rm{r}}}\lt \xi$, would lead to an effectively 1D geometry, significantly stabilising the soliton against dynamical decay. An alternative, perhaps more easily accessible way to achieve the same 1D dimensional reduction, could be based on reducing the scattering length by means of a Feshbach resonance. In the particular realistic geometry discussed throughout this paper, a reduction in the scattering length of ²³Na by a factor of 2.5 from its background value was sufficient to generate stable 1D solitonic structures over the probed regime of numerous collisions.

To better distinguish the solitonic nature of the excitations over other (linear/background) excitations, we found it advantageous to use a density engineering protocol in which the intensity of the perturbing laser beam is gradually turned on over a period of few tens of ms, reaching a maximum intensity of few times the chemical potential. To simpify the ensuing dynamics, and the observation of the propagating solitonic structures, it is advantageous to only generate a single counter-propagating soliton pair, which requires the waist of the laser beam to be narrow, broadly comparable to the healing length.

Looking at the role of thermal effects in the quasi-two-dimensional regime ${k}_{{\rm{B}}}T\lesssim {\hslash }{\omega }_{\perp }$, we performed an analysis based on two complementary models commonly used for nonequilibrium soliton dynamics (classical field simulations and SGPE). Despite their somewhat distinct predictions, both models consistently indicated a high likelihood of observing solitonic generation, azimuthal propagation, and occurence of (possibly a few) solitonic collisions under realistic experimental conditions and temperatures.

We thus hope that our study will assist experimentalists in engineering quasi-stable solitonic propagation in closed ring-trap circuits, and that such nonlinear excitations could in the future prove useful for atomtronic applications.

We are grateful to Mark Edwards for originally suggesting this project to us, and to Bill Phillips for pointing us towards density engineering in the early stages of the development of this project. Funding was provided by the UK EPSRC (Grants No. EP/I019413/1 and EP/K03250X/1), and we also acknowledge discussions with Andrew Baggaley, Carlo Barenghi, Simon Gardiner, Mariella Loffredo and Nick Parker. Data supporting this publication is openly available under an 'Open Data Commons Open Database License'. Additional metadata are available at: 10.17634/122626-1. Please contact Newcastle Research Data Service at rdm@ncl.ac.uk for access instructions.

In the main text we have argued that efficient sound-free generation of a single counter-propagating solitonic pair requires, in addition to the other carefully considered parameters, a gradual excitation scheme, and a narrow laser beam. To highlight the importance of those additional control parameters, here we give (for fixed other parameters) evidence of the post-perturbation dynamics when either of those criteria is not satisfied.

Figure A1 shows the situation analogous to our T = 0 reference case, when the Gaussian perturbation is suddenly turned on (over a physical timescale $\sim 36\;\mu {\rm{s}}$ corresponding to our time discretisation unit), depicting again the ensuing density and phase dynamics at the same times as in figure 1(c). More specifically, after being turned on, the perturbing potential is here kept at the constant maximum value of ${V}_{0}=2\mu$ for 36 ms, before being again 'instantaneously' removed. A detailed comparison of figures A1(a) and (c) reveals that the turning on/off sequence therefore does not appear to significantly modify the details (depth/speed) of the emerging solitonic structure, but rather it controls the amount of emitted sound during the generation process, which in turn indirectly affects the long-term soliton evolution due to soliton-sound interactions.

Zoom In Zoom Out Reset image size

Figure A1(b) shows the effect of increasing $\sigma /\xi$ to the value 1.5, which is here shown (for ${V}_{0}=2\mu$) to lead to the eventual dynamical generation of more than one pair of counterpropagating dark solitons (of different depths). We have checked that for values of σ up to the value of ξ (such that the Gaussian perturbation half-width at $1/{e}^{2}$ is $\sim 2\xi$), a single pair of counterpropagating solitons is generated, placing an effective limit on experimental perturbing potentials able to generate only single (as opposed to multiple) dark soliton pairs.

For completeness, we present in figure B1 a comparison of typical snapshots of the long-term evolution in the solitonic regime between the 1D and 2D limits, as characterised by the parameter ${l}_{{\rm{r}}}/\xi$. Shown here are images shortly after the solitonic structures are first engineered (left columns), and then when the solitons have returned to the same position after undergoing two, six, and fourteen collisions (corresponding to one, three and seven revolutions around the ring respectively). This clearly demonstrates the robustness of the solitonic structures against collisions, while also showing the much more confined nature of the excitations in the 1D regime (${l}_{\perp }\lt {l}_{{\rm{r}}}\lt \xi$), for which all radial excitations are completely suppressed. As a result, any solitonic excitations in 1D happen along the ring (structures occasionally appear more 'oval-shaped' than circular), with our numerics indicating no noticeable change in the soliton speed/depth over the probed timescales (other than a small oscillation in their respective values). This is in contrast to the 2D regime (top images), where the solitons, although still reasonably robust to collisions, do exhibit changes in their profiles in time (exhibiting a coupling between azimuthal and radial degrees of freedom), and also gradually decay (albeit at a rather slow rate). While the 1D regime evidently provides optimal conditions for observing such an effect, our simulations indicate that the main effects should still be largely visible even when ${l}_{{\rm{r}}}$ slightly exceeds the healing length.

Zoom In Zoom Out Reset image size