Engineering Bright Matter-Wave Solitons of Dipolar Condensates

We present a comprehensive analysis of the form and interaction of dipolar bright solitons across the full parameter space afforded by dipolar Bose-Einstein condensates, revealing the rich behaviour introduced by the non-local nonlinearity. Working within an effective one-dimensional description, we map out the existence of the soliton solutions and show three collisional regimes: free collisions, bound state formation and soliton fusion. Finally, we examine the solitons in their full three-dimensional form through a variational approach; along with regimes of instability to collapse and runaway expansion, we identify regimes of stability which are accessible to current experiments.


Introduction
Solitons have been observed in a profusion of physical systems, and fundamentally exhibit a dual nature possessing both wave and particle like qualities. All solitons owe their existence to the balance between their kinetic and interaction energies, which gives them their defining characteristic of maintaining their form over great distances. Solitonʼs dual wave-particle nature has led to them being extolled as information carriers in optical systems [1], as well as being of fundamental interest as solutions to nonlinear systems [2].
The realization of Bose-Einstein condensates (BECs) formed from ultracold atomic gases has, over the last two decades, opened a new chapter in nonlinear physics, allowing unprecedented experimental and theoretical insight into quantum mechanical effects of these many-particle systems [3]. In the weakly-interacting, lowtemperature limit the intrinsic nonlinear dynamics of the condensate are well-described by a mean-field wavefunction that encompasses the behavior of the coherent matter-wave. At the mean-field level, the van der Waals interactions between particles gives rise to a local nonlinearity that can support soliton solutions in quasi-one-dimensional waveguides. Under repulsive van der Waals interactions, dark solitons (non-dispersive waves of depleted density) are supported in the condensate [4]. Meanwhile, under attractive van der Waals interactions, bright solitons (droplets which are self-bound along the axis of the waveguide) have been observed [5][6][7][8][9][10]. Subsequent experiments have explored the effect of collisions with barriers [11,12], and also the role of the relative phase for two trapped solitons [13]. Matter-wave interferometry using bright solitons has also been the focus of recent experimental [9] and theoretical [14] interest.
Powerful analytical tools such as the inverse scattering transform have allowed the investigation of higherorder solitons [15] and the derivation of a particle model for the soliton dynamics and interactions, based on the knowledge of the scattering phase shifts [16]. This, in turn, has lead to the identification of regimes of chaotic dynamics for three trapped bright solitons [17] and the observed frequency shifts of trapped bright solitons in a recent experiment [18]. Examining collisions between bright solitons introduces an extra parameter, namely the relative phase. While these solitons collide freely in one dimension, numerical simulations [19][20][21] have demonstrated that in-phase collisions promote the collapse of the condensate in three dimensions, while a relative phase of π suppresses the collapse. Indeed, such π phase differences are believed to have existed between experimental bright solitons [6,7,22], critical to the observed stability of these states. A strategy to control the relative phase of the solitons has been proposed [23]. Attractive nonlinearities can also facilitate molecule-like bound states of two bright solitons [15,18,24], with these states being sensitive to both the relative phase and velocity of the solitons [24].
The creation of condensates with atoms possessing significant magnetic dipole moments- 52 Cr [25,26], 164 Dy [27] as well as 160 Dy, 162 [28] and 168 Er [29]-afford a new opportunity to explore the interplay of magnetic effects with the coherent nature of the condensate. The dipole-dipole (DD) interaction is anisotropic and longranged, contributing a nonlocal nonlinearity to the mean-field equation of motion for the condensate [30]. Importantly, the relative strength of the local to nonlocal interactions can be precisely tuned using Feshbach resonances, allowing the creation of condensates possessing a dominantly dipolar character [31]. The recent observation of the Rosensweig instability in a condensate of 164 Dy atoms [32] has focused interest towards the manifestation of self-bound droplet phases [33][34][35], where conventionally one would expect the condensate to undergo collapse. The unexpected stability of these states has been attributed to many-body effects [36][37][38][39]. There has also recently been the creation of a spin-orbit coupled dipolar degenerate Fermi gas, opening a new route to the study of complex phases in analogy with condensed matter physics [40].
Solitons whose existence depends on a non-local rather that a purely local nonlinearity represent a burgeoning area of interest in many disciplines of physics [41][42][43], and dipolar condensates provide a highlytunable platform to study such solitons. The effect of varying the relative strength of the local and nonlocal interactions has revealed novel bright [44][45][46] and dark dipolar matter wave solitons [47][48][49][50] in one dimension. Two-dimensional bright solitons are also predicted to be supported by the anisotropy of the DD interaction [51][52][53][54]. These works highlight important prevailing physical characteristics, including the existence of bound states of multiple dipolar solitons due to the intrinsic long-ranged nature of the dipolar interaction.
We begin by deriving the one-dimensional mean-field equation of motion for the dipolar condensate confined to a quasi-one-dimensional waveguide (section 2). In section 3, we procure the dipolar bright soliton solutions as a function of the DD strength and polarization angle, for attractive and repulsive van der Waals interactions. Subsequently in section 4 we focus on the interplay of the relative phase with the DD interaction strength on soliton-soliton collisions, highlighting regimes of free collisions, bound state formation and soliton fusion, as well as the effect of noise on such collisions. In section 5 we characterize the stability of the bright solitons in three dimensions, revealing the parameter regimes where the solitons are stable to collapse and expansion. Our findings are then summarized in the conclusion (section 6).

Theoretical model
We consider an ensemble of weakly-interacting atoms forming a BEC at zero temperature. Each atom has a mass m and permanent magnetic dipole moment d, polarized in a common direction by an external magnetic field. In the dilute limit described by the Gross-Pitaevskii theory, interactions are described by the two-body pseudopotential The first term in equation (1) describes the conventional short-range van der Waals interactions, with where a s defines the s-wave scattering length. The second term describes the long-ranged anisotropic DD interaction, given by where p = C d 4 dd 2 characterizes the DD interaction strength. Hereê j defines the unit vector in the direction of the coordinater j . Equivalently one can also show that equation (2) can be written in the form where θ defines the angle between the vector joining two dipoles and the direction of polarization. At the magic angle  q  54 m the DD interaction vanishes. When q q < m the DD interaction is attractive with dipoles lying head-to-tail. When q q > m the interaction is instead repulsive and the dipoles are orientated in a side-by-side configuration. Rather than defining individually the dipolar and van der Waals interaction strengths, it is convenient to work in terms the ratio e = C g 3 dd dd . It is further possible to consider the parameter region for < C 0 dd [55], where under rapid rotation the dipoles can be thought of as anti-dipoles, effectively reversing the attractive and repulsive regimes of interaction. The form of the DD interaction given by equation (2) will be used to perform the dimensional reduction in momentum space in what follows. The dynamics of the quantum gas are encompassed by the mean-field wave function ( ) Y t r, , which defines the atomic density distribution as The trapping potential appearing in equation (4) has a transverse trapping frequency w^, where = + r x y 2 2 2 defines the coordinate in the radial direction. As such the potential effectively forms a waveguide, which is uniform along the axial direction, z. Meanwhile, the nonlocal mean-field potential generated by the dipoles appearing in equation (4) takes the form We consider the zero temperature limit of the quasi-one-dimensional dipolar condensate, under which the three-dimensional GPE can be reduced to an equation of motion for the axial degree of freedom. Rigorous analysis of this dimensional reduction has been previously reported [57,58]. The transverse degrees of freedom are assumed to be in their harmonic ground state, which is encapsulated by the condition  w m , where μ is the 3D chemical potential (i.e. that associated with equation (4) above). The ansatz for the real-space wave function is assumed to be ( ) , where k i is the momentum associated with coordinate i and ( ) n k z z is the momentum space density in the axial direction. To perform the dimensional reduction, we first make use of the mathematical identity [59] (ˆˆ) Equation (6) introduces the quantity ( ) d^r jk which gives the transverse part of the δ-function, defined as for unit vector in momentum space,k j . Then, along with the Fourier representation of the δ-function, one can write down the Fourier transform of equation (2) as The one-dimensional analog of equation (5) can then be found using the definition given by equation (8)  . This gives where the parameter encapsulates the strength of the dimensionallyreduced DD interaction [60]. The momentum associated with the z-axis is defined as k z , while defines the exponential integral. Then, the one-dimensional mean-field model is encompassed by the pseudo-potential˜( ) , which includes the van der Waals as well as dipolar contribution. Equation (9) will be utilized to find the family of bright soliton solutions to the quasione-dimensional dipolar Gross-Pitaevskii equation. The effective 1D dipolar GPE is then written as where the one-dimensional dipolar potential is obtained via the convolution theorem as denotes the inverse Fourier transform. To obtain the solutions to equation (10) we work in momentum space using a split operator method. In what follows we adopt the so-called 'soliton' units [17], where length, time and energy are measured in terms of 2 defines the units of velocity and onedimensional van der Waals interaction strength respectively. We quantify how one-dimensional the system is through the ratio s =l l z . To be in the true one-dimensional limit one must have  s 1.

Dipolar bright solitons
The local cubic nonlinear Schrödinger equation (equation (10) in the limit e = 0 dd ) is known to possess bright soliton solutions for < a 0 s (corresponding to when the chemical potential m < 0). The single soliton solution, which for simplicity we take to be initially positioned at the origin is where u defines the velocity of the bright soliton, while f is the phase. This solution describes a sech-shaped profile which propagates at constant velocity. Here, we explore the family of dipolar bright solitons across the full parameter space afforded by the quasi-one-dimensional model. For a single bright soliton, we can independently vary four key parameters: e dd , θ, sgn(a s ) and also σ. Meanwhile the normalization is given by . Despite the additional dipolar term present in equation (10), we will see that the allowed dipolar bright soliton solutions are still sech-shaped. Throughout this work we take s = 0.2. Figure 1 shows the ground state densities obtained by solving equation (10) numerically in imaginary time. Each individual plot is divided into two regions: a flat homogeneous region, corresponding to regimes of net repulsive interactions, and a second region showing the inhomogeneous dipolar bright soliton solutions. The solid white line in each figure shows the crossover between these two parts, which corresponds to when the 1D chemical potential of the ground state solution crosses from positive (homogeneous state) to negative (bright soliton). The top (bottom) rows in figure 1 correspond to sgn(a s )<0 (>0). Fixing both the sign and value of the van der Waals interactions reveals that altering the polarization angle of the dipoles between q = 0 and q p = 2 has the effect of shifting the soliton solutions from , corresponding to anti-dipoles are found to support soliton solutions as the net interactions are attractive for dipoles polarized perpendicular to the z-axis. We note that in each of the cases presented in figure 1 the borderline between the homogeneous and soliton solutions does not coincide with e = 0 dd . This can be understood from the form of the dimensionally-reduced pseudo-potential, equation (9) which comprises both a nonlocal and local contribution, whose net effect is to shift the value of e dd at which the chemical potential μ changes sign.

Collisions
In the absence of dipolar interactions and within the one-dimensional nonlinear Schrödinger equation, bright solitons are known to collide elastically, emerging from the collision with their original speed and form. The net effect of soliton-soliton interaction is a position and phase shift in the outgoing solitons. In this section we study the collisions of two dipolar bright solitons, exploring the role the relative phase plays in collisions as a function of the DD interaction and the initial kinetic energy of the solitons. In what follows we simulate two-counter-

In-phase collisions
We consider the collision dynamics of bright solitons with zero initial phase difference, f D = 0. As we shall see, the collisions of the dipolar solitons can be inelastic. In order to quantify the elasticity of the soliton dynamics, we compute the coefficient of restitution, defined as where v j is the velocity of soliton j, and t i and t f are the initial and final times, respectively. The coefficient of restitution is a measure of the amount of kinetic energy before and after a collision event. If h = 1 the incoming and outgoing speeds are identical and a collision is perfectly elastic. For h < 1 the outgoing speeds are less than the incoming speeds, and the collision is deemed inelastic. As we shall see it is also possible to realize h > 1, corresponding to the outgoing speeds exceeding the incoming speeds, and which can occur, for example, when interaction (van der Waals plus dipolar) energy is transformed into kinetic energy during the collision. The coefficient of restitution η is mapped out in the (e dd ,v i ) plane in figure 2 (d). Each pixel represents an individual simulated collision between two bright solitons, with the color representing the value of the inverse of the coefficient of restitution (we plot h 1 as this quantityʼs scale evolves at a slower rate over the parameter range considered). Three different regimes of dynamics can then be identified. For relatively weak DD interactions  e 2.5 dd , the collisions are almost perfectly elastic (h~1), independent of the initial velocity. Figure 2(a) shows a typical set of collision trajectories in this limit. Here the incoming solitons scatter elastically with each other and escape at longer times. In the intermediate regime, short-lived (meta-stable) bound states are found, whose dynamics are inelastic. Here, the balance of the initial kinetic energy of the solitons to the DD strength is favorable to the formation of a short-lived bound state; again a simulation typical of this situation is shown in figure 2(b). We note that similar dynamics were reported by [61] for kink solutions to the sine-Gordon equation, which also exhibited short-lived bound-states. In contrast to this work, we observe a single transition from bound state to free solitons, as the speed of the incoming soliton is increased. In the limit where the DD interactions are large, we instead observe soliton fusion with  h 1 1. A trajectory plot indicative of this limit is shown in figure 2(c), showing the collision and eventual merging of the two solitons. Here the individual solitons do not re-emerge at long times.
We can gain an understanding of the short-lived bound state by considering the effect of the in-phase collision on the solitonʼs kinetic and potential energies. After the bound state has initially formed, each successive collision event causes a redistribution of some interaction (van der Waals and dipolar) energy into kinetic energy, causing the effective oscillation frequency of the soliton molecule to increase. Eventually, this is enough to cause the bind to break, releasing the two solitons. Note that the total energy remains constant throughout the simulations.

Out-of-phase collisions
It is also pertinent to consider the equivalent dynamics for out-of-phase dipolar bright solitons, with figure 2(e) plotting h 1 in the (e dd ,v i ) plane. For relatively large incoming speed, the dynamics are almost elastic ( ) h~1 , while for increasing dipole strength or decreasing incoming velocity, the dynamics are instead found to be increasingly inelastic. This system has a regime of long-lived bound states occurring for low incoming speeds, with example collision trajectories shown in the panel figure 2(f).
The binding of two out-of-phase dipolar bright solitons, has been studied previously by [45,46]. Unlike their in-phase counterparts, the π phase difference preserves long-lived bound states. One cannot assign a value of η to the collisions in this regime.

Noise
In order to quantify the collisional sensitivity of the dipolar bright solitons, we calculate the coefficient of restitution in the presence of noise. The noise is implemented by introducing a random term to the phase with mean zero whose amplitude  noise is a given fraction of the initial peak density of the soliton density such that , and max ( ( ) n z ) is the peak soliton density. As an example we consider h 1 versus e dd , for a fixed incoming speed and for in-phase collisions. Figure 3 shows a comparison of h 1 with  = 0 0 (no noise) and  = -10 0 2 . In the presence of noise, each data point was obtained by averaging over 10 individual simulations. The error bars represent the standard deviation. For low values of e dd , the value of h 1 follows very closely the value obtained in the absence of noise. This is attributed to the elastic dynamics being insensitive to the phase noise. On the other hand, for stronger dipolar interactions, there is a marked deviation from the nonoise case, resulting in a larger value of h 1 . In this regime, the presence of noise introduces an apparent repulsion between the solitons, which is large enough to make the collisions shown in figure 3 elastic. Such an effect has also also been noted in non-dipolar bright soliton collisions [62]. We see qualitatively the same behaviour for different strengths of noise, with increasing amounts of noise causing the bound-states to break sooner, until no bound-states are formed.
A corollary of the phase noise is that the bound state dynamics in the presence of noise show fewer oscillations before escaping their bind. This effect is illustrated in the soliton trajectories with/without noise in insets figures 3(a) and (b). The presence of only small amounts of noise demonstrates how sensitive the binding dynamics are: figure 3(b) illustrates that only one oscillation can occur in this example before the solitons escape. For larger amounts of noise the bound states are no longer present, even at these larger values of e dd . On the other hand, the collsions for out-of-phase solitons, in contrast, are insensitve to noise. The effective repulsion in these collisions serves to stabilize the collisions against the noise.
Although the analysis presented here is rudimentary, it nonetheless allows one to comment on the effect of dissipative processes that are present in a real system, especially those at finite temperature. For example, for small changes in temperature, it is expected that the formation of bound states would be unfavorable. The dipolar interactions would play an increasingly diminished role in the system dynamics, since this term fundamentally depends on the condensate density, which is reduced due to the presence of the noncondensate [63].

3D stability
A repercussion of attractive interactions between particles is that, for sufficiently large number of particles and/ or interaction strength, the mean-field wave function of a 3D condensate is unstable to collapse. For systems possessing only short-ranged isotropic van der Waals type interactions, the critical point of collapse has been extensively studied in [64][65][66], and gives insight into the parameter regimes where one can expect stable soliton dynamics [67]. The presence of DD interactions are expected to modify the collapse point significantly [68], which has recently been explored for lower dimensional systems in [69][70][71] as well as for two-dimensional dipolar bright solitons using 3D simulations [72,73]. As well as this, the presence of the dipolar interaction can lead to the spectacular d-wave collapse of the condensate in three dimensions [74]. Understanding when the dipolar bright soliton is unstable to collapse in turn allows us to identify regimes of stability applicable to the quasi-one-dimensional dynamics described earlier in section 4.

Gaussian variational approach
We employ a variational approach that approximates the wave function of the dipolar soliton as a 3D Gaussian packet with variable width in each dimension [64]. This approach has provided important insight into the stability of non-dipolar bright solitons, predicting thresholds for instability which agree closely with experiments [66]. Under general conditions an appropriate variational ansatz is given by where the lengthscale¯ w = l m is based on the geometric mean of the transverse trapping frequencies ( ) w w w = x y 1 2 and s x y z , , are the dimensionless variational widths of the packet. Equation (14) is normalized to the total number of atoms, We seek to calculate the total energy of the packet in terms of the above parameters. We write the total energy as = constitutes the remaining energy contributions arising from kinetic energy, potential energy (from the transverse trapping ( ) 2 2 ) and van der Waals interaction energy. These contributions to the energy are handled in real space. Meanwhile, the dipolar contribution E dd is evaluated in momentum space, using the convolution theorem. We consider the case where the atoms forming the condensate are polarized by an external magnetic field such that their individual dipole moments form an angle α with the z axis, as shown schematically in figure 4. This configuration leads to a momentum space pseudo-potential  We perform the integrations appearing in equation (17) in spherical polar coordinates, and assume the dipoles are polarized parallel to the z axis of the condensate so that a = 0. Then, using equations (14)- (17) 1 18 x y x y

Stability analysis
For a given set of system parameters (b e , dd and λ), equation (18a) defines an energy landscape as a function of the dimensionless length scales s s , x y and s z . A stable variational solution is an energy minimum (local or global minimum) in this landscape. It is instructive to consider the non-dipolar case as a pedagogical example, which has been studied numerically and using a variational approach [65,66] to determine the parameter regimes where the bright soliton is stable to collapse. For moderate attractive van der Waals interactions ( | |  N a l 0.7 s ), the energy landscape has a global minimum at the origin (s s = = 0 z ), representing a collapsed state, and a local minimum at finite s z and s^, representing the 3D bright soliton solution. This local energy minimum is preserved by a delicate balance between the kinetic and interaction energy. A saddle point separates the local and global minimum. For | | » N a l 0.7 s the saddle point and local minimum merge; this marks the threshold at which the 3D bright solitons become unstable to a runaway collapse. Meanwhile, for repulsive van der Waals interactions, the energy landscape decreases monotonically towards s  ¥; z any wavepacket will disperse axially and no stable solutions exist (unless axial trapping is applied, which we do not considered here). It is the goal of this section to obtain and analyze the nature of the critical points of the full 3D dipolar bright soliton as a function of the interaction parameters β and e dd .
We seek the points in the energy landscape defined by equation (18a) at which an instability manifests. This analysis in principle requires us to solve the set of equations defined by simultaneously with the determinant of the Hessian matrix set equal to zero [67], which are in general a set of four coupled algebraic equations in four unknowns. In order to simplify the general equations (18a) and (18b) but still gain useful insight, we consider the axially-symmetric case for which w w = x y and s s x y , . The integral over the azimuthal angle f appearing in equation (18a) can then be evaluated as the f dependency of ( ) k f 2 is removed in this limit. Although we have assumed the bright soliton is well-described by a Gaussian wave function (equation (14)), it is worth contrasting the value at which the instability manifests found from the equivalent analysis for e = 0 dd assuming either a sech or Gaussian like variational wave function. If * b denotes the dimensionless critical collapse parameter, then one finds that * * b b p = 3 s g , which means that the difference between the two approaches is ∼2%, supporting our choice of a Gaussian variational ansatz. Proceeding, we wish to solve , and the matrix elements of the Jacobian J are defined as We employ an iterative procedure to obtain numerical solutions for the collapse point in the ( ) b e , dd parameter space, using the exact analytical result for e = 0 dd at which the bright soliton collapses as the initial point to start the numerical calculation from. , dd parameter space is shown in figure 5(b). This particular plot shows a stable minima indicative of regions where the bright soliton is stable. Figure 5 For < C 0 dd however, one finds instead that  k 1, indicative of a very elongated (cigar) like cloud. Here the solitons exist closer to the one-dimensional limit. Finally, we note that close to the unstable boundary (red lines) the bright soliton can pass through a region where k > 1. Such states would be challenging to observe, as the the system would preferentially wish to collapse due to the presence of thermal or quantum fluctuations.
Although our presented results only consider a particular choice of dipole polarization, one can still comment on the stability of the dipolar bright soliton when the dipoles are, say, polarized perpendicular to the axis of the waveguide. As was noted in section 3, altering the polarization of the dipoles from q = 0 to q p = 2 has the effect of swapping the regimes of e dd where one obtains bright soliton solutions. We can speculate that a similar effect would occur when examining the stability of the dipolar system, except here we would see the regions associated with collapse and runaway expansion switch. However, this case breaks cylindrical symmetry, requiring a fully anisotropic ansatz to capture the stability of the system. This greatly complicates the analysis, as one has three variational width parameters to consider.

Conclusions
We have analyzed the solutions, quasi-one-dimensional dynamics and full three-dimensional stability of dipolar bright solitons. The bright soliton solutions obtained from the dipolar Gross-Pitaevskii equation exhibit a number of novel features, including collisions which have regimes of elastic behavior, bound state formation and soliton fusion. These regimes where shown to depend sensitively on the dipolar interactions and the presence of noise, which modify the phase shifts of the solitons. We quantified the collisional behaviour in terms of the coefficient of restitution. Analysis of soliton dynamics in terms of the coefficient of restitution could then provide important insight for systems where the full scattering phase shifts may be difficult to obtain analytically. The stability of the full three-dimensional dipolar system was explored; in particular it emerged that the dipolar interactions can destabilize the bright soliton in two distinct ways, either through a traditional collapse, or instead through a runaway expansion along the axis. For axially-polarized dipoles the former occurs when the DD interaction is positive, while the later is associated with regimes of anti-dipoles for which the DD interaction is instead negative.
Our results provide a benchmark for future experimental studies of nonlocal soliton in dipolar condensates. In turn, this system offers unique opportunities to explore the fundamental properties of nonlocal solitons in general with the immense tunability of atomic physics.
Data supporting this publication is openly available under an Open Data Commons Open Database License [76].