Exploring the Stability and Dynamics of Dipolar Matter-Wave Dark Solitons

We study the stability, form and interaction of single and multiple dark solitons in quasi-one-dimensional dipolar Bose-Einstein condensates. The solitons are found numerically as stationary solutions in the moving frame of a non-local Gross Pitaevskii equation, and characterized as a function of the key experimental parameters, namely the ratio of the dipolar atomic interactions to the van der Waals interactions, the polarization angle and the condensate width. The solutions and their integrals of motion are strongly affected by the phonon and roton instabilities of the system. Dipolar matter-wave dark solitons propagate without dispersion, and collide elastically away from these instabilities, with the dipolar interactions contributing an additional repulsion or attraction to the soliton-soliton interaction. However, close to the instabilities, the collisions are weakly dissipative.


I. INTRODUCTION
Solitary waves, or solitons, are excitations of nonlinear systems that possess both wave-like and particle-like qualities. They obey wave equations and yet do not disperse, maintaining their shape and speed by balancing dispersion with nonlinearity. Solitons appear across a wide range of physical systems that include water, light, plasmas and liquid crystals [1], and have been touted as playing a fundamental role in the fabric of our universe [2]. A more recent addition to this list is the atomic Bose-Einstein condensate (BEC) formed in an ultracold gases which are described by a nonlinear Schrödinger equation known as the Gross-Pitaevskii equation (GPE). Experimental demonstrations of solitons in BECs include both the generation of dark solitons [3][4][5][6][7][8][9][10][11], and bright solitons [12][13][14][15][16][17]. The interaction between the atoms in these experiments was short range and isotropic (predominantly of the van der Waals type) giving a local cubic nonlinearity in the GPE, with dark and bright solitons supported for repulsive and attractive nonlinearity, respectively. In binary condensates, dark-bright soliton complexes have also been probed experimentally [8,18]. The experiments confirm that solitons in BECs and in classical systems such as water or light are for many purposes the same phenomenon and share the same three particle-like defining properties, namely: permanent form, localization within a region, and emergence from collisions with other solitons unaltered, except for a phase shift [19]. It is important to bear in mind, however, that in BECs the soliton relies on quantum mechanical coherence across the sample, and is at heart a probability wave [20].
Ultracold Bose gases provide an appealing system in which to explore soliton physics because of the almost complete absence of dissipation (due to the superfluid nature of the gas) and the high degree of experimental control that can be exerted over the atoms and their interactions using lasers as well as magnetic and electric fields. This, for example, has led to proposals to access exotic solitons such as in spin-orbit coupled condensates [21][22][23], and chiral solitons in 'interacting' gauge theories [24], and non-local solitons in dipolar condensates [25][26][27][28] which are the subject of this paper. Furthermore, matter-wave solitons have been proposed for applications in precision interferometry [12,17,29,30] and surface interrogation [31].
The experimental achievement of condensation of bosonic elements possessing large magnetic dipole moments -52 Cr [32,33], 164 Dy [34,35] and 168 Er [36] has opened yet another chapter in BEC physics. Dipoles introduce long-range anisotropic interactions falling off as 1/r 3 , in contrast to the usual short-range isotropic interactions, and hence give rise to an additional non-local nonlinearity [37]. This has striking consequences, as observed experimentally in the form of magnetostriction of the condensate [38] and shape-dependent stability [39], anisotropic collapse and explosion [40], and droplet formation analogous to the Rosensweig instability in classical ferrofluids [41,42]. This modulational instability is a direct result of the roton dip the dipolar interactions introduce into the excitation spectrum [43,44]. Another prominent physical system that supports dark solitons with nonlocal interactions are nonlinear optical We consider an elongated Bose-Einstein condensate of atoms (shown as a density isosurface) with dipole moments (shown as arrows) which are copolarized in a common direction. For suitable interaction parameters, a dark soliton can be supported, characterized by a 1D density notch and a non-trivial 1D phase slip. media, where the interaction of the electric field of light with the material gives rise to a defocussing local nonlinearity [45], and a nonlocal nonlinearity can arise due to thermal conduction. This is typically modeled via a response function which decays exponentially with distance [46,47]. A strong mathematical analogy can be drawn between these optical systems and atomic gas condensates, since both are studied with a similar underlying model, at least in the purely local case [48]. Studies of the optical systems with nonlocal nonlinearity have focused on their stability [49], and the arising interaction forces between dark solitons [50,51]. This has in turn lead to the observation [52] of both repulsion and attraction between dark solitons, with the latter supporting bound states in the optical context. The generation of dark solitons from shocks in these systems has also been experimentally studied [53].
The local cubic defocussing Schrödinger equation represents a solvable model within the framework of the inverse scattering method [54,55]. The inclusion of a cubic non-local potential to this model greatly complicates its analytical treatment within this method, there currently being no-known exact solutions. Nevertheless, approximate methods including series approximations for the non-local potential [56] and variational calculations [57] have been successfully employed within this context to elucidate the physics of these models, with the latter capturing the existence of dark soliton bound-states.
A series of theoretical investigations have indicated that dipolar interactions in BECs also considerably enrich the properties of solitons. In quasi-one-dimensional geometries, bright [25] and dark solitons [27,28] can be supported when the net (van der Waals and dipolar) interactions are attractive and repulsive, respectively. Dark-in-bright and bright-in-dark dipolar solitons have also been predicted [58]. Yet perhaps the most interesting facet of solitons in general are their interactions and collisions [59], and in dipolar condensates the solitons themselves, considered as individual particle-like entities, inherit non-local soliton-soliton interactions in addition to the usual short-range soliton-soliton interaction [25,28]. The play off between these two contributions can lead to the formation of unconventional bound states of bright [26] and dark solitons [27,28]. Dipole-dipole interactions are also predicted to support two-dimensional bright solitons in quasi-2D geometries [60][61][62], and suppress the well-known transverse "snaking" instability of dark solitons in three-dimensional geometries [63].
In a recent work [28], we predicted the existence of dark soliton solutions in a homogeneous quasi-one-dimensional dipolar condensate (shown schematically in Fig. 1), and studied the non-local interaction between two such solitons. Here we expand upon this topic by presenting a comprehensive analysis of the family of dark soliton solutions and their interactions, across the main system parameters, namely the angle of polarization, the relative strength of the dipolar interactions, the soliton speed and the width of the quasi-1D system. In order to un-derstand the regimes of soliton stability, we establish the stability properties of the soliton-free ground state of the system; this allows us to relate the soliton stability to the phonon and roton instabilities of the system. We examine the family of single soliton solutions, including their integrals of motion, and finally explore the soliton-soliton collisions through simulations of the dipolar GPE.
The main body of the paper is organized into four sections. In Sec. II we derive the mean-field equation of motion for the dipolar condensate. Following this in Sec. III we analyze the homogeneous system, obtaining analytical expressions for the position of the phonon and roton instabilities. Section IV describes single dipolar dark soliton properties and solutions across the full parameter space of the problem, while Sec. V explores their collision dynamics. Our findings are summarized in the conclusion, Sec. VI. The body of the paper is supported by a technical appendix explaining the numerical method used to obtain the dark soliton solutions.

II. MEAN-FIELD MODEL OF THE DIPOLAR CONDENSATE
We consider a gaseous BEC of ultracold weaklyinteracting atoms with mass m and permanent magnetic dipole moment d. Within the mean-field Gross-Pitaevskii theory, the atom-atom interaction in the low energy limit is described by the pseudo-potential The first term arises from the short range isotropic van der Waals-type interactions, where g = 4πh 2 a s /m and a s defines the s-wave scattering length. The long-range and anisotropic interaction appears as the bare dipole-dipole interaction between point dipoles [64,65] where C dd = 4πd 2 characterizes the strength of the dipole-dipole interaction andê j is the unit vector in the coordinate directionr j . Equation (2) can also be written as where the angle θ is defined between the vector joining the dipoles and the polarization direction (see Fig. 1). At the so called magic angle θ m 54 • , the dipole-dipole interaction reduces to zero. Assuming C dd > 0, then for θ < θ m (dipoles orientated predominantly head-to-tail) the dipolar interaction is attractive. For θ > θ m (dipoles dominantly side-by-side) the interaction is repulsive. It is also possible to consider a regime of "anti-dipoles", C dd < 0, as proposed in Ref. [66] by rapidly rotating the dipoles, for which the attractive and repulsive regimes are reversed. It is convenient to specify the relative strength of the dipole-dipole interaction to the van der Waals interaction via the parameter ε dd = C dd /3g. By means of Feshbach resonances to tune g [74], as well as the above mentioned method of generating C dd < 0, it is experimentally possible to access systems over the full range −∞ < ε dd < ∞, with negative or positive g or C dd . The quantum state of the dipolar condensate is described by its mean-field wavefunction Ψ(r, t); the condensate density distribution follows as n(r) = |Ψ(r, t)| 2 . In the limit of zero temperature, the wavefunction obeys the dipolar GPE [37] ih The external potential V (r) which confines the cloud can in general take many forms, but we assume it to be a harmonic waveguide given by V (r) = mω 2 ⊥ r 2 ⊥ /2, where r ⊥ = x 2 + y 2 defines the radial coordinate, and the transverse trapping frequency is ω ⊥ . We neglect any axial confinement since the exact soliton solutions we seek exist only for axially uniform systems. It is worth noting that such axially-uniform waveguides are accessible experimentally [12]. Meanwhile, Φ dd (r, t) is the non-local mean-field potential generated by the dipoles Φ(r, t) = dr U dd (r − r )n(r , t).
In this work we consider a quasi-one-dimensional dipolar condensate [68], under which the 3D GPE can be reduced to an effective 1D equation. The transverse harmonic trapping is sufficiently tight (hω ⊥ µ, where µ is the chemical potential of the threedimensional system) that no transverse degrees of freedom are excited. The wavefunction is taken to follow the ansatz Ψ(r, t) = ψ ⊥ (r ⊥ )ψ(z, t), where ψ ⊥ (r ⊥ ) = (l ⊥ √ π) −1 exp(−r 2 ⊥ /2l 2 ⊥ ), and l ⊥ = h/mω ⊥ defines the transverse harmonic length scale. Such a state constitutes a single mode approximation (SMA) for the axial dynamics of the condensate. Inserting this ansatz into the 3D GPE (4) and integrating out the transverse (x-y) dimensions leads to the effective 1D dipolar GPE The effective 1D dipolar mean-field potential is where the effective dipolar pseudo-potential is U 1D (z) = U 0Ũ1D (|z|/l ⊥ ) with [69,70] The term in square brackets in Eq. (8) gives a non-local contribution to the mean-field interactions, while the second contact-like term describes a short-ranged local contribution to the mean-field interactions. The strength and orientation of the dipole-dipole interaction is captured by U 0 = C dd (1 + 3 cos 2θ)/32πl 3 ⊥ . It is also useful to know the energy of the condensate, which is given by the integral The three terms here represent kinetic energy, interaction energy due to vdWs interactions and interaction energy due to dipolar interactions. From now on our analysis of dipolar dark solitons will be based on the effective 1D dipolar GPE (6). Our results will be presented in terms of the natural quantities of the homogeneous (soliton-free) condensate. Taking n 0 as the uniform density, the chemical potential (the energy eigenvalue of Eq. (6)) follows as µ 0 = n 0 g/(2πl 2 ⊥ ) + Φ 0 . Here the first term represents the van der Waals contribution, and Φ 0 = −C dd n 0 [1 + 3 cos 2θ]/24πl 2 ⊥ is the dipolar contribution. The natural units of length and speed are the dipolar healing length, ξ =h/ √ mµ 0 , and the speed of sound, c = µ 0 /m. A time unit follows as τ = ξ/c. To parameterize the cross-over from three to one dimension we define σ = l ⊥ /ξ; the quasi-one-dimensional limit is valid for σ < ∼ 1 [76].

III. STABILITY OF THE HOMOGENEOUS SYSTEM
Dark solitons are excitations of a background condensate, and so the stability of these states is heavily influenced by that of the background condensate. Here we perform an analysis of the stability of the homogeneous quasi-1D dipolar condensate. This is a generalization of the results of Refs. [44,70] which only considered dipoles polarized along the long axis. Although strictly speaking the dark soliton state will possess a different excitation spectrum to that calculated in this section, we will see that this simple analysis agrees remarkably well with the position of the instabilities of the dark soliton solutions to Eq. (6).
The stability of the condensate can be deduced from the fate of small perturbations whose energies are given by the Bogoliubov excitation spectrum. For a homogeneous system, the Bogoliubov spectrum depends on the momentum (Fourier) space version of Eq. (8). However, rather than directly transforming Eq. (8), we can proceed more easily by returning to the original 3D dipolar interaction given in Eq. (2) and transform it into Fourier space by using the useful identity [83] 1 where the quantity δ ⊥ jk (r) is the transverse part of the delta function, defined as . These show, respectively, the conventional phonon/free-particle spectrum, the flattening off of the dispersion relation at intermediate k, and the emergence of a roton minimum. The condensate has arbitrary width σ = 1.
The quantityk j appearing in Eq. (11) is the jth component of the unit vectork in Fourier space. Using these results together with the Fourier representation of the delta function, the Fourier transform of Eq. (2) can then be directly found as The dimensional reduction to 1D can be performed directly in Fourier space in an analogous way to the realspace reduction, again by assuming a harmonic ground state in the transverse directions. This yields the momentum space equivalent of Eq. (8) [44] wherehk z is the momentum associated with the axial zdirection, while E 1 (x) = ∞ x dt t −1 e −t is the exponential integral. The total one-dimensional pseudo-potential, including van der Waals and dipolar contributions, is then (8) and (13) can be split into a non-local and local contribution; the former gives a contribution to the total contact interactions while the latter forms the important long-ranged part of the dipolar interaction. We note that in general the scattering length, and hence g, can be modified both by the dipolar interactions and confinement induced resonances (due to tight 1D trapping) and so may not take on their 3D dipole-free values [70]. Nevertheless, the general form of this pseudo-potential is expected to hold.
With the above results in hand, the Bogoliubov excitation spectrum for perturbations of momentumhk z from the background can be written as where k =h 2 k 2 z /2m defines the free particle energy and V 1D (q) = q exp(q)E 1 (q)−1/3. The excitations associated with Eq. (14) are running waves of the form that constitute small amplitude fluctuations about the stationary condensate. Figure 2 illustrates the dispersion relation corresponding to Eq. (14). For certain parameters, the dispersion relation has the same structure as for the non-dipolar case (dot-dashed blue line): for low k it is linear in k (the phonon regime), changing to a quadratic form at higher k (the free-particle regime). The system is prone to a phonon (long wavelength k z → 0) instability. In nondipolar homogeneous condensates, this arises when the mean-field van der Waals interactions become attractive, g < 0. Examining the small k z behavior of the dispersion E k =hω k relation given by Eq. (14) yields where c s is the speed of sound associated with this longwavelength regime. We identify the phonon instability as occurring when the bracketed term (which corresponds to the homogeneous chemical potential µ 0 ) is less than zero. This leads to imaginary frequencies, signifying the unstable exponential growth of long-wavelength perturbations.
In other parameter regimes, the dipoles can change the form of the dispersion relation at intermediate k. In certain instances, it causes a flattening off of the dispersion spectrum, as shown in Fig. 2 (dashed black line), while in more extreme cases a minima can form in the dispersion relation at finite momentum, i.e. the roton. An example is shown in Fig. 2 (solid red line), with a pronounced dip appearing at k z ∼ ξ −1 . The associated local maximum in the dispersion relation is termed the maxon. If the roton minimum touches the zero-energy axis then the condensate undergoes the roton instability. The roton is predominantly driven by transverse (off-axis) effects of the dipole-dipole interaction, and becomes less pronounced in the one-dimensional limit (σ → 0) [37].
We identify the roton instability as follows. The expression (14) for E 2 k (k) is differentiated with respect to k z and set equal to zero so as to identify the stationary points (which may correspond to the roton or the maxon). This is then combined with the dispersion relation (14) equated to zero (the maxon is automatically excluded from this result as it can never touch the zero energy axis). With some manipulation, the critical wavenumber at which the roton touches zero energy is found to be While the above expression provides k c for known system parameters, we wish to predict the onset of the instability as a function of the dipole-dipole interaction strength ε dd . We can eliminate k c by combining the above expression with a rearranged version of the dispersion relation (14) equated to zero: These two equations can be solved iteratively to predict the critical ε dd for the roton instability to occur as a function of θ and σ, as will be presented below.
In Figure 3 we map out stability diagrams in the (θ, ε dd ) plane, showing the regions of roton instability (shaded red) and phonon instability (shaded blue). To give insight into the role of the transverse condensate width, three values are considered: (a) σ = 0.1, (b) 0.5 and (c) 1. For each value, we distinguish between the g > 0 (left column) and g < 0 (right column) cases. The phonon instability is independent of σ throughout, and intuitively follows from the play-off between the van der Waals interactions and the dipolar contribution to the contact interactions. Consider, for example, the case of repulsive vdW interactions. For θ = 0 the dipoles lie in the attractive end-to-end configuration, and when the dipoles are sufficiently strong (ε dd > 1) they can overwhelm the repulsive vdW interactions, inducing phonon instability. Conversely, for θ = π/2 the dipoles are sideby-side; conventionally this is a repulsive configuration, but in the regime of anti-dipoles (C dd < 0) this configuration is attractive, and induces phonon instability when the anti-dipoles are sufficiently strong (ε dd < −2).
The regions of roton instability are sensitive to σ. For low σ (cases (a) and (b) in Figure 3), the roton instability arises only for attractive vdWs interactions (g < 0). Deep in the 1D regime (case (a)) the roton instability arises only in a narrow band in the (θ, ε dd ) plane; as σ is increased this band expands (case (b)). However once a critical value (σ crit ∼ 0.8) is exceeded, as per row (c), the roton instability shifts instead to appearing only for repulsive vdW interactions (g > 0). The value of σ crit does not depend on the angle θ.
To further explore the shift of the roton instability from g < 0 to g > 0 for increasing σ, Fig. 4 depicts the roton instability in the (σ, ε dd ) plane for the extreme polarization angles, (a) θ = 0 and (b) θ = π/2. In both cases The band (grey) of instability is bounded by the onset of the roton (dashed) and phonon (solid) instabilities. Note that either side of this unstable band, the system is stable for only g > 0 or g < 0, as indicated.
it is seen that in the quasi-1D limit (σ 1) the roton instability band narrows, approaching the onset of the phonon instability (horizontal dashed line). However as σ is gradually increased, the roton instability undergoes a change of sign at the critical value σ = σ crit . Figure 4 (c) shows the critical wavelength defined as λ c = 2π/k c plotted as a function of σ. It is seen that λ c is monotonically increasing, and is always greater than σ, indicating that our one-dimensional analysis is valid. The inset to Fig. 4 (c) shows a zoomed in portion of this graph centered around σ = σ c , and clearly shows the discontinuity in λ c , indicated by the dashed blue line.

A. Non-Dipolar Dark Soliton Solutions
In the absence of dipolar interactions (ε dd = 0) and for repulsive s-wave interactions (g > 0) it is well-known that the 1D GPE is integrable and supports dark soliton solutions of the form [45,77], Here β = 1 − v 2 /c 2 , where v is the velocity of the soliton, c is the sound speed in the condensate, and the dark soliton's centre of mass is initially placed at the origin. The family of solitons commonly feature a density depression and a phase slip, with the depression density, n d and total phase slip S related to the soliton speed, v, via v/c = 1 − (n d /n 0 ) = cos(S/2). The v = 0 "black" soliton (z = 0) has zero density at its centre and a π phase slip, while the v = c soliton is indistinguishable from the background density. Since dark solitons deplete the density profile, they are analogous to particles with negative mass [78].

B. Dipolar Dark Soliton Solutions
We now set about exploring the dipolar dark soliton solutions across the important parameters of soliton speed v, polarization angle θ, dipolar strength ε dd , and condensate width σ. We consider the dark solitons as stationary states in the moving frame, and obtain these solutions by numerically solving the 1D dipolar GPE in the moving frame. The details of this approach can be found in Appendix A. Figure 5(a,b) plots the spatial density profile n(z) of the black (v = 0) soliton solutions as a function of ε dd , for the extreme polarization angles of (a) θ = 0 and (b) θ = π/2. The former represents the typical behavior for all solutions in the range θ < θ m , and the latter shows the converse. Figure 5(c,d) shows the corresponding plot for a finite speed case, v = 0.5c. The grey bands represent the range of ε dd for which no stable condensate exists, consistent with the corresponding stability diagrams presented in row (a) of Fig. 3. Away from the unstable region, the density profiles resemble the tanh 2 -density of the conventional dark solitons, with a width of the order of the healing length ξ (which should be noted is itself a function of ε dd and θ). Near to the phonon instability (solid lines) the density profile diverges in width; this is due to a cancellation between the local interactions arising from the explicit van der Waals interactions and the implicit local contribution to the DDIs, with a similar effect seen for vortices in 2D dipolar condensates [71]. Meanwhile, as the roton instability (dashed lines) is ap- proached, density ripples form symmetrically around the soliton, decaying as they recede from the core. For the cases shown in Fig. 5, the ripples can rise to twice the background density with the most prominent parts being the two dominant lobes either side of the dark soliton (see also Fig. 6). They arise due to the soliton state mixing with the roton, an effect akin to that predicted for vortices in 2D [65,71,80]. The ripples can be understood from an energetic point of view by noticing that they occur when the dipolar interactions are repulsive along the axial direction meaning that the system can lower its energy by placing more dipoles near to the empty core. The repulsive dipolar interaction due to the lobes in turn causes a density reduction next to them, then another peak, and so on. The ripples are thus a direct effect of the long range nature of the dipolar interactions. Despite the density modulations, the soliton depth still follows the relation n d /n 0 = 1 − v 2 /c 2 , familiar from non-dipolar dark solitons. We also find that the density modulations are slightly enhanced for slower solitons. As previously discussed in Sec. III, the roton and its instability is sensitive to the condensate width, σ. To determine the effect of σ on dark solitons, Fig. 6 compares the black soliton solution close to the (a) roton instability and (b) phonon instability, for different values of σ. Note that we maintain σ < 1 throughout so to satisfy the governing criteria for a quasi-1D condensate (see Section II). As σ increases the system becomes less "1D" in nature and the effects of dipolar interactions become more pronounced. At the roton instability, the density ripples grow rapidly with σ, becoming as large as 15n 0 for σ = 0.3 (dotted black line). The length scale of the ripples also increases with σ; this is consistent with the earlier homogeneous analysis, which showed the roton wavelength to increase with σ ( Fig. 4 (c)). Meanwhile, at the phonon instability the soliton has a funnel shaped density profile, which widens with σ. Having explored the dependency of σ we will employ σ = 0.1 for the remainder of our work (the relevant homogeneous stability diagrams being shown in row (a) of Fig. 3).
We briefly comment on the manifestation of the phonon and roton instabilities on the soliton solutions. Imagine crossing the phonon instability threshold, from the stable side to the unstable side. The net interactions switch from repulsive to attractive, such that dark solitons are no longer stable. At the same time the background condensate undergoes a modulational instability, as per the non-dipolar attractive condensate [81], and fragments into bright soliton-like structures (the stable structures under net attractive interactions). Next consider crossing the roton instability, again from stable to unstable. The ripples surrounding the dark soliton grow and we find the BEC eventually collapses. However, in both cases, the sharp growth in density means that higher-order effects such as three-body losses [40] and quantum fluctuations [42] become important. Such physics is not contained within our dipolar GPE and is beyond the scope of this work.

C. Integrals of Motion
The family of non-dipolar dark solitons [Eq. (19)] possess an infinite number of integrals of motion (viz. conserved quantities). The first three of these have a clear physical meaning: the soliton normalization, momentum, and energy [45,72]. In order to calculate finite values for each of these quantities, we compute their renormalized versions, that is, the difference between the quantity in the presence and absence of the soliton. In this subsection we generalize these quantities for dipolar dark solitons and explore their dependence on the dipolar parameters.
The renormalized norm and momentum of the dark soliton are defined as where n 0 is the homogeneous background density. The renormalized energy must explicitly include the contribution from the dipoles. To calculate this we follow a similar approach to Ref. [73]. According to Eq. (9) the energy of a homogeneous system of length L is E 0 = (gn 2 0 /4πl 2 ⊥ + Φ 0 n 0 /2)L, where Φ 0 is the homogeneous dipolar potential. In order to make a direct comparison between a homogeneous system and a system containing a soliton the quantity E −µN must be considered, so as to account for the different particle number N between the two systems. Thus the renormalized soliton energy can be written as In the absence of a soliton, for which Φ dd = Φ 0 and |ψ| 2 = n 0 , this expression correctly reduces to zero. In the absence of dipoles and using the soliton solution (19), the renormalized norm, momentum and energy follow as where, recall, β = 1 − v 2 /c 2 . Meanwhile the effective mass of the non-dipolar dark soliton is found from the relation m 0 = ∂P 0 sol /∂v, and is given by We evaluate the norm, momentum and energy of the dipolar solitons numerically according to Eqs. (20)(21)(22). In the left rows of Fig. 7 (a-c) these conserved quantities are plotted as a function of the soliton speed for three values of ε dd , with the polarization angle fixed to θ = 0. Note that these three values of ε dd correspond to 168 Er parameters (ε dd = 0.4) and near to the phonon (ε dd = 0.8) and roton instabilities (ε dd = 5). Throughout, the dipolar results show the same qualitative structure as the non-dipolar result (solid black line), but can be quite different quantitatively. The soliton norm N sol [ Fig. 7 (a), left] decreases monotonically with speed due to the reduction in the soliton depth for larger speeds, becoming zero for v = c when the soliton is indistinguishable from the background. The soliton momentum P sol [ Fig. 7 (b), left] also decreases monotonically with v, due to the decreasing norm and decreasing phase gradient (the phase gradient determines the local fluid velocity according to u(z) = (h/m)∂ z S(z), where S(z) is the phase profile) across the soliton. It has a universal value for v = 0 of P sol = πhn 0 due to the π-phase-step profile of the v = 0 solitons. The momentum becomes zero for v = c when the norm and phase gradient reach zero. The dipoles have the greatest effect on P sol for intermediate velocities. Finally, the soliton energy E sol [ Fig. 7 (c), left] also decreases with v, associated with the decreasing interaction energy as the soliton gets increasingly shallow and the decreasing kinetic energy due to the reduced density and phase gradients, and at v = c E sol = 0.
Close to the phonon instability (blue dashed lines) the integrals of motion tend to be larger than the non-dipolar case. This can be related to the wide funnel-shaped profile which develops close to this instability. This vastly increases the effective volume of the soliton core, i.e. the norm. This in turn raises the momentum and energy (in the latter case, due primarily to the increased interaction energy associated with the larger density depletion). The momentum and energy are also modified by density gradients. Meanwhile, close to the roton instability (purple dot-dashed line) the integrals of motion tend to be smaller. The density ripples which form here act to reduce the effective volume of the soliton core, which reduces the momentum and energy relative to the nondipolar case.
Finally, in the right rows of Fig. 7 (a-c) the conserved quantities are plotted as a function of the polarization angle for the three values of ε dd , while keeping the soliton speed fixed (v = 0.5c). The non-dipolar result is constant in each plot. What is particularly prominent is that the θ dependence is the same for all 3 integrals of motion, with only the scale changing. At the magic angle, θ m ≈ 0.3π, the integrals equal the non-dipolar result, due to the vanishing of the dipolar potential here. Note that the gap in the curve for ε dd = 5 is due to the absence of stable solutions here, consistent with the stability diagram in Fig. 3 (a).
The effective mass of the soliton, defined as m = ∂P sol /∂v, is shown in Fig. 7 (d). The effective mass is negative throughout, tending towards zero effective mass when v = c, as expected. For most cases the effective mass increases monotonically with v. However, close to phonon instability it has a unusual form, being approximately constant for v/c < ∼ 0.4 and decreasing to a local minimum at v/c ≈ 0.75 m /m.

V. DYNAMICS OF DIPOLAR DARK SOLITONS
In this section we explore the rudimentary dynamics of the dipolar dark solitons. In particular we seek to estab- lish their soliton-like nature. We will approach this by reference to the general definition of a soliton given by Johnson and Drazin [19], which is of three key properties: (i) permanent form, (ii) localized within a region of space, and (iii) emergence from collisions unchanged, barring a phase shift. The results presented here are based on numerical propagation of the 1D (lab-frame) dipolar GPE using the Crank-Nicolson method, using a suitable initial condition featuring soliton solutions obtained from the BCGM method.
A. Propagation Figure 8 shows the evolution of a v = 0.5c dipolar dark soliton (with fixed θ = 0) for three values of ε dd : 0.4 (corresponding to 168 Er), 0.8 (close to the phonon instability) and 5 (close to the roton instability). Insets show the density and phase profiles across the soliton. For all three cases, the soliton maintains a permanent and localized form, with no radiative losses. It also undergoes centre-of-mass translation at the expected speed. As such, these states satisfy the soliton criteria (i) and (ii) above. It is also worth observing the phase profile across the soliton. For ε dd = 0.4 [ Fig. 8(a)], the phase profile is practically identical to that of the non-dipolar dark soliton, with a tanh-shaped step which relaxes to the asymptotic value over a short length scale of the order of the healing length. Close to the instabilities, the phase relaxes over a much larger length scale, of around ∼ 50ξ close to the phonon instability [ Fig. 8(b)] and around 400ξ close to the roton instability [ Fig. 8(c)]. In all cases the total asymptotic phase slip is the same as for the non-dipolar soliton (this is not directly evident from the inset in (c) due to the limited length range of this plot). Close to the instabilities the phase profile also features distinctive prominences. At the phonon instability these are broad, while at the roton instability they are of order of the roton length scale.

B. Collisions
In non-dipolar systems the interaction between multiple dark solitons has been experimentally observed and theoretically studied [75,82]. In a symmetric collision for solitons satisfying 0 < v/c < 1 2 the s-wave interactions create a repulsive force causing the solitons to appear to reflect at short distance. For velocities satisfying 1 2 ≤ v/c < 1 the solitons appear to pass through each other. In both cases the outgoing solitons are unchanged from the incoming form, barring a phase shift. In the presence of dipolar interactions this behavior is modified [27,28]. In Ref. [28] the effect of polarization angle on the collisions was explored, revealing an additional non-local repulsion or attraction on the soliton-soliton interaction due to the dipolar interactions. The latter case, in combination with the conventional short-range repulsion, was shown to support bound states. Here we explore the soliton collisions further, exploring the role of the experimentally-tunable interaction parameter ε dd and the soliton speed. We investigate the effect of the dipolar contribution near to the instabilities, for a system comprised of 168 Er atoms [36], with polarization θ = 0.  and C dd > 0 for this case, the dipoles lie end-to-end and attract each other. However, the dark soliton is a region of depleted dipoles, and can be interpreted as a positive density of anti-dipoles [28,86] which repel each other. This repulsive contribution to the dark soliton interaction will arise whenever the dipoles are net repulsive, i.e. when C dd > 0 with θ < θ m or when C dd < 0 with θ > θ m . The repulsive nature of the collision becomes washed out at higher incoming speeds (row (c), right panel). While the soliton collision is stable at low speed (left panel), at high speed the collision is inelastic, with energy lost from the solitons via the emission of sound waves (visible as bands propagating away from the collision at the speed of sound).
The case of ε dd = 5 (row (d)) instead has C dd < 0, i.e. repulsive dipoles. This in turn leads to an attractive contribution to the soliton interaction, and this is clearly observed in the corresponding soliton collisions. Note the distinctive sharp "pinching" of the solitons during their collision at low speed. More generally, this attractive contribution arises whenever the dipolar interactions are net attractive, i.e. when C dd < 0 with θ < θ m or when C dd > 0 with θ > θ m . Here, for both low and higher incoming speeds, the collisions are inelastic through sound emission.
Note that when the dipoles are polarized at the magic angle θ m the dynamics are equivalent to the non-dipolar case [28].
Away from the phonon and roton instabilities, the solitons collide elastically, and emerge unscathed from the collision. This satisfies the third soliton criteria outlined above. However, close to the instabilities, the collisions become dissipative, with sound being radiated away. This is particularly prevalent for higher speed collisions. We note, however, that the energy dissipated into sound waves during a single collision is typically very small, for example, in the maximally-dissipative case presented in Fig. 9 ((d), right panel), the energy loss is ∼ 10 −3 %E sol .

VI. CONCLUSIONS
In this work the family of dark solitons supported in a quasi-one-dimensional dipolar Bose-Einstein condensate were studied. A bi-conjugate gradient method was implemented to numerically obtain these non-trivial solitons as stationary solutions in the moving frame, as a function of the dipole-dipole interaction strength (ε dd ), the polarization angle θ and the soliton speed. The phonon and roton instabilities of the system play a key role in modifying the density and phase profiles of the solitons, which can deviate significantly from the non-dipolar form in these regimes. The dipolar dark solitons were characterized in terms of their integrals of motion (norm, momentum and energy). Due to the modified profiles in the presence of dipolar interactions, these quantities differ from their non-dipolar form, particularly so close to the instabilities. The prominent role of the phonon and roton instabilities in the soliton solutions motivated a detailed and general analytical treatment of these effects in the quasi-1D dipolar condensate. This, in particular, revealed the sensitivity of the roton to the transverse condensate size, σ, and the increase in the roton length scale as the dimensionality crossover, σ ∼ 1, is approached.
In isolation, the solitons propagate with unchanged form throughout the parameter space. Away from the instabilities their collisions are elastic, but become dissipative via emission of sound waves close to the instabilities. Thus, close to the instabilities these structures deviate from solitons in a strict sense, although it should be noted that the energy dissipated in a single collision is very small.
Data supporting this publication is openly available In this appendix, we describe how the soliton solutions in the main body of the paper were obtained. We consider the dark solitons as stationary solutions of the GPE in the moving frame. We obtain them by numerically solving the moving-frame time-independent GPE, (Ĥ + vp z − µ)φ = 0, wherep z = −ih∂/∂z defines the momentum operator in the z-direction. This is performed using the bi-conjugate gradient method [84] (BCGM). This technique has been used to obtain moving-frame vortex solutions in the 2D/3D GPE [79,85]. It is worth noting that this approach provides the true dipolar dark soliton solutions for arbitrary speed. In contrast, the approach to finding soliton solutions based on imaginary time propagation of the GPE [27] requires a priori knowledge of the soliton phase, and so is only capable of obtaining black (v = 0) soliton solutions, for which the phase profile is known to be a step function of amplitude π.
The BCGM is an iterative method based on the Newton-Raphson method for finding roots of equations.
Numerical implementation was handled in MATLAB with the function bicgstab. In practice, it is convenient to concatenate the real and imaginary components of f into a single vector of length 2N , and similarly for φ. J is then of size 2N × 2N . Taking the non-dipolar soliton solution, as defined by Eq. (19), as the initial guess for φ (centered at the origin), we find that the BCGM robustly converges to the required dipolar soliton solution (also centered at the origin). Note that a different choice for initial guess may lead instead to the homogeneous ground state.
The absolute value of the phase is arbitrary, and to aid convergence we fix the value of the phase at one end of the grid. The grid spacing ∆z is typically 0.1ξ. Away from the phonon/roton instabilities, we employ a box of typical length 100ξ. However, close to these instabilities, box sizes of up to 1600ξ were required to ensure good approximation to the infinite limit (that is, for the meanfield dipolar potential to reach its homogeneous value at the boundaries). The solutions were deemed converged when the relative residual, calculated as ||Jδφ + f ||/||f ||, had fallen below an arbitrary tolerance of 10 −5 .
The above numerical method is akin to solving the equation Jδφ = −f for δφ, then setting φ (p+1) = φ (p) + δφ at each step p. The advantage of using the BCGM is that we only require knowledge of the transpose of the Jacobian, which is numerically faster to obtain than matrix inversion.