Engineering Dark Solitary waves in Ring-Trap Bose-Einstein condensates

We demonstrate the feasibility of generation of quasi-stable counter-propagating solitonic structures in an atomic Bose-Einstein condensate confined in a realistic toroidal geometry, and identify optimal parameter regimes for their experimental observation. Using density engineering we numerically identify distinct regimes of motion of the emerging macroscopic excitations, including both solitonic motion along the azimuthal ring direction, such that structures remain visible after multiple collisions even in the presence of thermal fluctuations, and snaking instabilities leading to the decay of the excitations into vortical structures. Our analysis, which considers both mean field effects and fluctuations, is based on the ring trap geometry of Murray et al. Phys. Rev. A 88 053615 2013.


Introduction
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][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20][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][23][24][25][26][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][29][30][31][32], or can be engineered by means of well known techniques such as phase imprinting [33][34][35], density engineering [33,[36][37][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 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][41][42][43][44][45][46][47] and thermal [43,[48][49][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 one-dimensional 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 one-dimensional (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 destabilizing 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" (Sec. 2), we focus on the question of optimisation of the generation of such solitary waves, by means of the Gross-Pitaevskii equation (Sec. 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 [54,55] (Sec. 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 [55,[62][63][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 non-optimal choice of density engineering parameters and a comparison to the idealised 1D regime are discussed in two Appendices.

Physical Set-up and Parameter regime
We consider a trapped ultracold atomic gas ( 23 Na atoms, scattering length a s = 2.75nm) confined in a ring-shaped trap of the form (see Fig. 1a for a visual representation): where V G , r 0 and w are respectively the depth, radius and 1/e 2 half-width of the ringgaussian potential. The radius, r 0 , is the distance from the center, where the potential reaches its minimum, whereas a length of 2w specifies the effective size of the ring channel where the gas is confined. We assume tight transverse harmonic confinement with frequency ω ⊥ 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 Eq. (1)) a narrow Gaussian potential of the form: i.e. applied in the left half-plane (region across y = 0, for x < 0, see Fig. 1a for reference), where Θ(x) is the step function, σ is the half-width of the Gaussian barrier and V L (t) the time-dependent laser amplitude. The ramp needs to be kept on for a time long enough (few ms -tens of 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 minimize the linear (sound wave) excitations emerging from sudden perturbations, we follow a scheme similar to that used in the experimental work of Ref. [47], such that the perturbing potential is linearly ramped up over a time τ on to its maximum value, V 0 , and then ramped down to 0 according to (perturbation on for times 0 ≤ t ≤ τ pert ): where τ on the duration of the ramping-on sequence and (τ pert − τ on ) the ramping off timescale, which we have chosen to be relatively short. Consistent with Ref. [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 τ on = 35.5ms and a ramping down time (τ pert − τ on ) = 0.5ms. The particular shape of the perturbation is shown for the chosen parameters in Fig. 1b.
The underlying system geometry and induced perturbation lead to three physicallydistinct sets of controllable parameters, respectively characterising: (i) the unperturbed ring trap potential V (r) (depth V G , location of minimum of confined potential r 0 , width of ring w) with the transverse frequency ω ⊥ entering our analysis implicitly, through its contribution to the effective two-dimensional interaction parameter, g 2D ∝ √ ω ⊥ ; (ii) the density-engineering perturbation V pert (amplitude V 0 , 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 showing also where the blue-detuned laser is applied (see main text). (b) Evolution of the amplitude of the imposed perturbing potential in units of the chemical potential µ; the perturbing potential is linearly ramped up to V 0 = 2 µ (solid line) and down to V 0 = 0 (dashed line) over a time period of 35.5 ms and 0.5 ms respectively, as defined by Eq. (3); a vertical dotted line is shown here for reference, in order to distinguish our procedure from a sudden turning off of the barrier. (c) Density (top), renormalized 'carpet' (middle) and phase (bottom) plots at times t = 36, 43, 72, 119, 750 ms (left to right) after switching on the perturbation, with t = 36ms corresponding to the time when the perturbation has just been turned off. The renormalized '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. The emergence of 'solitonic' excitations is evident from the combined density and phase information with t = 119 and 750 ms respectively corresponding to the cases after one and thirteen collisions, thus demonstrating that the generated 'solitonic' structures remain largely unaffected by multiple collisions. To hide spurious features in the phase plots, a mask has been used where the density is lower than 10% of the peak density at equilibrium.
excitation scheme, the main physics is actually set only by a few parameters (or rather, their ratios).
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 ≈ 15000 (which leads to a corresponding peak density of 2500 l −2 0 corresponding roughly to 25 atoms/µm 2 ). This choice ensures that the 2D condition µ < ω ⊥ holds: specifically, we choose a chemical potential µ ∼ ω ⊥ /3. For such an atom number we are typically in the 2D regime defined by l ⊥ < ξ < l r , although a further dimensional reduction to a 1D regime (l ⊥ < l r < ξ) 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 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 σ ≈ ξ where ξ is the minimum value of the healing length, as calculated at the peak density [33].
Increasing the width of the perturbation σ to values σ/ξ > 1 (or significantly increasing V 0 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 parametrizing the spatial extent of the emerging macroscopic excitation (e.g. dark soliton or vortex): this is defined here as ξ = / (mg 2D n 0 ), where g 2D = √ 8π(a s /l ⊥ )( 2 /m) is the 2D interaction strength, a s is the scattering length and n 0 refers to the maximum density. For a fixed transverse confinement ω ⊥ 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 2D n 0 which controls the effective soliton width. For numerical convenience, when probing different parameter regimes, we choose to fix n 0 and vary g 2D by increasing or decreasing the value of a s by up to three times its background value (see subsequent Fig. 3).
• The maximum amplitude, V 0 , of the density perturbation, which parametrizes 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 0 ) is parametrised by σ.
Our subsequent generation and dynamical stability analysis is thus primarily based on the chosen control parameters (V 0 /µ) fixing the depth of the emerging solitary wave excitations, and (l r /ξ) setting the effective dimensionality of the system.

Dynamics at T = 0
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: 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 Fig. 1c. More specifically, this figure shows (in situ) condensate density (top), renormalized ('carpet') density (middle) and phase (bottom) of the system at times t = 36, 43, 72, 119, 750 ms (left to right) after initiating the (ramped) perturbation, where 36 ms corresponds to the time that the perturbation is switched off, and 119 ms and 750ms the times after one and thirteen collisions. The renormalized '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 −2 0 , where l 0 = 10 µm is our reference 'length unit'. Figure 1c reveals ‡ the emergence of counter-propagating sound waves moving rapidly away from the region of the density perturbation, followed by two slower counterpropagating 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 Fig. 1c. More specifically, for the case considered here (with V 0 = 2µ), the generated structures travel with a 'mean' velocity v ≈ 0.5c (based on the time taken to cross half the ring, i.e. t ≃ 66 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 Fig. 5d) 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 > 0 case (Fig. 7b).
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 characterizes the maximum depth that the emerging nonlinear ‡ See also movie (2d.solitonic.long.evolution) in supplementary material. excitations can acquire. We anticipate requiring an intensity V 0 µ, 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 Ref. [33] for elongated quasi-1D BECs). In our present work, we span intensity perturbations ranging from V 0 = 0.5 µ to V 0 = 10 µ 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 /µ lead to deeper (slower) solitary wave generation. To characterise this, Fig. 2a shows the dependence of the effective azimuthal propagation speed (scaled to the local speed of sound measured at r 0 ) 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 s /c ∼ (V 0 /µ) −α , with a numerically-extracted exponent α ≈ 0.18. Using the standard expression for pure one-dimensional solitons in homogeneous settings [22], we can re-write this formula in terms of the depth, n sol , of the soliton from the peak of the unperturbed density, as n sol /n 0 ∼ 1 − (V 0 /µ) −2α , whose dependence is shown in the inset to Fig. 2a.
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 2b 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-one-dimensional BECs [46,73]. However, as evident from Fig. 2b, 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 /µ) is shown Corresponding plot in terms of the soliton depth, n sol , scaled to the maximum unperturbed density n 0 , obtained using the standard homogeneous relation v s /c = 1 − n sol /n 0 valid for purely 1D solitons [74]. (b) Dependence of v s /c on time for different maximum amplitudes of the perturbing potentials V 0 /µ = 1 (black squares), 2 (blue diamonds), 5 (red circles). All velocity ratios given here are based on the ratio of the instantaneous value of the soliton depth at r 0 (the radial distance specifying the location of the trap minimum) to the (peak) unperturbed density at that point. The identifiable oscillations are likely due to the fact that the minimum of the solitary wave structure is not always located at r 0 [see also Fig. 5a]. We have also verified that the approximate determination of the soliton velocity based on measuring its motion around the ring yields similar results. by vertical dotted lines in Fig. 2b, with this figure spanning 5-7 revolutions. This suggests that the observed anti-damping may not be a direct consequence of the (headon) 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µ (intermediate blue diamonds in Fig. 2b), their velocity has increased by a mere ∼ 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 µ < ω ⊥ , or equivalently l ⊥ < ξ) 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 r to the healing length ξ, whose effect is discussed next (with the idealised limit of l ⊥ < l r < ξ corresponding to an effective 1D regime of practically stable solitonic propagation over experimental timescales).

Solitonic Behaviour and Dynamical Instabilities
Figure 1c 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 r /ξ), we can identify 3 reasonably distinct dynamical regimes over the broad range 0 < V 0 /µ < 10 of density perturbation amplitudes probed, with the corresponding 'phase diagram' shown in Fig. 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 Fig. 3 illustrate characteristic snapshots identifying the key features of each of those regimes.
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 Fig. 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 (Fig. 3), we have chosen to probe the distinct regimes by varying the value of (l r /ξ) for fixed l 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 −2 0 = 25 atoms/µ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 < l r /ξ < 2.3, the scattering length spans the range ≈ [0.4a s , 3a s ] (i.e. 1.1 nm < a s < 8.3 nm for 23 Na used here), while the chemical potential still satisfies the 2D criterion through the condition 0.25 < µ/ ω ⊥ < 1 (with the number of atoms lying in the range ≈ [12000, 28000] respectively). We have also verified that changing ω ⊥ (instead of a s ) in the range ≈ [0.2ω ⊥ , 9.2ω ⊥ ], while still keeping the peak density where the structures are highly excited and emit at least one pronounced density depression in their attempt to eventually maintain, after internal re-arrangement, some solitonic features. The characteristic behaviour defining each regime is displayed in the 2D carpet snapshots reported in each of those cases (for the clockwise propagating wave): the dashed semi circles in each of these images indicate where the density drops to 10% of the peak equilibrium value. Crossover to the 3D regime (vertical dashed line) occurs roughly at l r ≈ 2.3ξ (corresponding l ⊥ ≈ ξ). The solitonic regimes features an internal 'subdivision' around l r = ξ, with the 1D regime exhibiting perfectly symmetric solitons also in the radial direction [see subsequent Fig. 5a] and enhanced stability, facilitated by the suppression of radial excitations. The horizontal dotted line at V 0 = µ indicates the regime below which only rather shallow structures appear following the density perturbation, in the sense that the depth of the soliton (measured from the top) does not exceed 15% of the peak unperturbed density, hence the soliton may not be pronounced enough to observe experimentally. n 2D (r 0 ) ∝ (ξ 2 a s √ ω ⊥ ) −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 Fig. 4 showing detailed snapshots of the evolution of set times for all cases §.    Fig. 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 r /ξ.
One can also internally sub-divide this regime, based on the value of (l r /ξ), since for values l r < ξ one arrives at an effectively 1D geometry satisfying l ⊥ < l r < ξ.
A comparison of density profiles between the case with l r /ξ = 0.8 and our reference case of Fig. 1c (l r /ξ = 1.3) is shown in Fig. 5a. Figure 5a shows 2D carpet plots of the solitonic structures when they are at the top of the ring, also plotting corresponding one-dimensional density cuts in the radial (Fig. 5b) and azimuthal (Fig. 5c-d) directions. We find that, although for l r < ξ the solitonic structures are symmetric along the radial direction about the middle of the ring, increasing the ratio l r /ξ, makes the structures less symmetric around the minimum of the ring trap ( Fig. 5a-b), with the location of the density minimum shifted towards the outer edge of the trap; more specifically, by comparing the radial profiles (Fig. 5b), while for l r /ξ = 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 r /ξ = 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.
The 1D azimuthal profiles for each case are also shown in Fig. 5c-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 ( Fig. 4b): 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 r /ξ), 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 r /ξ) 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 /µ), with the crossover region satisfactorily accounting for such behaviour.
As the values of l r /ξ 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 ( Fig. 4c): 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.

Dynamics at T > 0
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][80][81][82]. In the simplest approach (Sec. 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 [Eq. (4)]. This approach is typically referred to as the 'classical field' method [81,[83][84][85][86] (being closely related to the finite temperature truncated Wigner [87][88][89][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 fluctuationdissipation 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 stochastic Gross-Pitaevskii equation (SGPE) [54,94,95], which in our current 2D setting takes the  ) for the two cases, with density cuts taken at y 0 = 1.85 l 0 and 2.05 l 0 , for l r /ξ = 0.8 and l r /ξ = 1.3 respectively. Analytical soliton solutions are constructed from their measured speed through the relation v s /c = (n min (0, y 0 )/n 0 (0, y 0 )), where n min is the soliton depth and n 0 the unperturbed equilibrium density, with the healing length calculated at the peak unperturbed density.
form [56]: where φ(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 Stochastic Projected Gross-Pitaevskii Equation [55,81]). This should be directly contrasted to the usual Gross-Pitaevskii equation [Eq. (4)], where ψ(r, t) denotes simply the condensate wavefunction. In Eq. (5), (thermal) fluctuations are mimicked by the presence of the noise term η(r, t) which has Gaussian correlations of the form η * (r, t)η(r ′ , , 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 Fig. 1c, 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.

'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 [Eq. (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 Fig. 6a for different temperatures of the unperturbed thermal state. As before, we show densities (top), renormalized densities or 'carpet' plots (middle) and phase (bottom) plots for T = 1, 9, 10 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.
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 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 See also movie in accompanying material [2d.solitonic.T9nK]. . Note that these plots show the entire classical field density, |φ| 2 , rather than the condensate density, |ψ| 2 , of the ordinary Gross-Pitaevskii equation shown until now. One could in principle perform further analysis to extract the corresponding density images for the (quasi-)condensate, which would look smoother; however, the location and nature of the solitons in the condensate would closely mimic the effects seen in the classical field plots, adding no further insight into the soliton stability. Moreover, the displayed noisy profiles are closer in nature to what would be observed in an experiment.
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 Fig. 6b. 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 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 ≈ 1.85l 0 , x ≈ −1.85l 0 and again at x ≈ 1.85l 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 Fig. 6b]. 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 one-dimensional geometries [60,97]. This appears to be in partial disagreement to the findings of Refs. [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 Note that the resolution of these images is sub-µm, implying that an experimental study would actually reveal smoother profiles, in agreement with our earlier experimental-theoretical comparisons [56,57].
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-densityengineering) subplots of Fig. 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 Fig. 7a 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 quantum and thermal fluctuations are retained in this approach, and the quasi-condensate is obtained by using an extension of the Bogoliubov theory to treat low-dimensional Bose gases [98]. In our model the equilibrium solution is determined self-consistently via the SGPE, and contains information about both density and phase fluctuations; the phase-coherent condensate, or suppressed density-fluctuations quasi-condensate could then be extracted a posteriori from the SGPE classical field, and its density is expected to qualitatively resemble the plotted classical field density, but with the fluctuations largely suppressed. do not appear to critically affect the underlying solitonic shape even after a few collisions [see rightmost image in Fig. 7a]. To verify the solitonic nature of the emerging structures prior to any collisions, Fig. 7b plots their azimuthal one-dimensional 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. 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 non-equilbrium soliton dynamics on top of a fluctuating initial state.

Full stochastic evolution
To further improve on our earlier T > 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 φ(r, t) exhibit full dynamical coupling to the high-lying modes of the system, i.e. maintaining here both dynamical noise η(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 Stochastic Projected Gross-Pitaevskii Equation (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 γ ∼ 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 Fig. 8 show the corresponding case with a much larger (heuristically chosen) value of γ ∼ 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 = 164ms, corresponding to two collisional events), as shown in the third set of plots in Fig. 8.
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].

Conclusions
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 ω ⊥ ) and small atom number (or equivalently chemical potential µ) such that the two-dimensional condition µ < ω ⊥ is satisfied, thus suppressing threedimensional 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 ⊥ < ξ < l r ), unless the effective radial ring length l r satisfied l r 1.5ξ, 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 r < ξ, would lead to an effectively one-dimensional geometry, significantly stabilizing 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 23 Na by a factor of 2.5 from its background value was sufficient to generate stable one-dimensional 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 B T ω ⊥ , we performed an analysis based on two complementary models commonly used for non-equilibrium soliton dynamics (classical field simulations and stochastic Gross-Pitaevskii equation). 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.

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

Appendix A: Role of instantaneous and broad density perturbations
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  Figure A1: Evolution of emerging solitonic structures for different density engineering protocols: (a) As in Fig. 1c, but with the perturbation added suddenly (i.e. over a time equal to our time unit ∼ 36 µs), shown here at the moment the perturbation is switched off (t = 36ms) and subsequent times t = 43, 72, 119 ms (left to right) when the system evolves freely. (b) As in Fig. 1c (so with a ramped perturbation), but with σ/ξ ≈ 1.5; these are shown here at slightly different times t = 36, 50, 65, 137 ms (left to right) in order to best reveal the two ensuing solitonic pair dynamics.
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 A1a shows the situation analogous to our T = 0 reference case, when the Gaussian perturbation is suddenly turned on (over a physical timescale ∼ 36µs corresponding to our time discretization unit), depicting again the ensuing density and phase dynamics at the same times as in Fig. 1c. More specifically, after being turned on, the perturbing potential is here kept at the constant maximum value of V 0 = 2µ for 36ms, before being again 'instantaneously' removed. A detailed comparison of Fig. A1a and Fig. 1c 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. Figure A1b shows the effect of increasing σ/ξ to the value 1.5, which is here shown (for V 0 = 2µ) 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 ∼ 2ξ), 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.

Appendix B: Solitonic Propagation in 1D versus 2D Regimes
For completeness, we present in Fig. B1 a comparison of typical snapshots of the longterm evolution in the solitonic regime between the 1D and 2D limits, as characterised by the parameter l r /ξ. 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 ⊥ < l r < ξ), for which all Figure B1: Typical carpet (top rows) and phase (bottom rows) snapshots depicting the initial generation (leftmost column) and subsequent propagation (after the number of indicated collisions) of the counter-propagating solitonic structures in the 1D (l ⊥ < l r < ξ) and 2D (l ⊥ < ξ < l r ) solitonic regimes at the indicated times. As the generated solitons have different speeds, and the snapshots have been chosen to depict times when the solitons have returned to their initial position after a certain number of collisions (number of revolutions is half the number of collisions), the actual times of those snapshots do not coincide in the 1D and 2D cases. Parameters as in Fig. 1c, except in the 1D regime where l r = 0.8ξ (facilitated through the use of the modified 0.4a s scattering length.) . 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 r slightly exceeds the healing length.