Quantum droplets of quasi-one-dimensional dipolar Bose-Einstein condensates

Ultracold dipolar droplets have been realized in a series of ground-breaking experiments, where the stability of the droplet state is attributed to beyond-mean-field effects in the form of the celebrated Lee-Huang-Yang (LHY) correction. We scrutinize the dipolar droplet states in a one-dimensional context using a combination of analytical and numerical approaches. In particular we identify regimes of stability in the restricted geometry, finding multiple roton instabilities as well as regions supporting one-dimensional droplet states. By applying an interaction quench to the droplet, a modulational instability is induced and multiple droplets are produced, along with bright solitons and atomic radiation. We also assess the droplets robustness to collisions, revealing population transfer and droplet fission.


I. INTRODUCTION
Bose-Einstein condensates possessing long-ranged, anisotropic dipole-dipole interactions have been realized in a series of ground-breaking experiments with highly magnetic atoms. The first generation of experiments achieved condensation of 52 Cr [1,2], 164 Dy [3,4] and 168 Er [5]. Recently, a second series of experiments has achieved a quantum analogue of the classical Rosensweig instability [6], as well as the realization of droplet states [7,8] -where the gas enters a high density phase whose stability has been attributed to the influence of quantum fluctuations [9][10][11][12][13]. Dipolar condensates constitute weakly correlated systems, and can exhibit properties and behaviour similar to that of a liquid in the beyondmean-field limit [18][19][20][21][22]. Droplet states have also been realized in condensate mixtures [23] -supported by a balance between attractive s-wave interactions between the atoms and quantum fluctuations -and photonic systems [24] -here the nonlinear medium provides a repulsive d-wave contribution that stabilizes the light beam.
The quantum liquid Helium II is well-known to exhibit a roton minimum in its excitation spectrum; this is supported by the strong interatomic interactions and correlations [25], where roton excitations typically occur at wavelengths comparable to the average inter-particle separation, indicating that the superfluid is close to forming a crystalline structure [26,27]. Although dipolar condensates are weakly correlated, the nonlocal character of the dipolar interaction supports roton-like excitations [14,15], which have now been experimentally realized in a gas of 166 Er [16,17]. Here the roton lengthscale is dictated by the geometry of the gas. A plethora of theoretical investigations have focussed on detailing the correspondence between rotons in the Helium II phase and weakly interacting dipolar gases. Early work examined the possibility of roton excitations in a quasione-dimensional setting [28], as well as the manifestation of rotons in a rotating dipolar condensates [29] and the identification of a roton mode in trapped dipolar condensates -leading to the identification of 'roton fingers' [30], a discrete manifestation of the instability found in the homogeneous system.
The existence of the roton mode indicates the proximity of the quantum fluid to long-range crystalline order. If the quantum fluid can simultaneouly support a superfluid state, the system is a candidate for a supersolida phase of matter where these two forms of order coexist [31]. The possible unambiguous detection of a supersolid state in Helium has attracted long debate. Very recently, state-of-the-art experiments with dipolar condensates have reported supersolid phases in these systems in a quasi-one-dimensional setting [32,33]. Complementary theoretical investigations in lower-dimensional dipolar gases have explored nonlinear wave structures in the form of dark solitons [34][35][36][37] and bright solitons [38][39][40][41][42]. Related theoretical proposals have also examined the possibility of supersolids in binary condensates [43], whose existence is supported by a stripe-like phase.
Within the weakly-interacting regime, dilute atomic condensates are well described by the celebrated Gross-Pitaevskii description [44,45]. A natural extension to this is the Bogoliubov-de Gennes formalism, which constitutes the linear response theory of the nonlinear Schrödinger equation, and gives insight into the behavior of the elementary excitations and collective modes of the condensate [46]. Early work demonstrated that this approach could also be applied to dipolar condensates [47], revealing the effect of the dipolar interaction on the excitations. Knowledge of the excitation spectrum gives insight into many important properties of these systemsin many physical effects the dimensionality of the system plays a key role in the dipolar condensates overall behavior [48]. Dipolar condensates with spin degrees of freedom have also been analyzed within the Bogoliubov-de Gennes framework, exploring the interplay of the excitations with the magnetic phases inherent to these systems [49]. The anisotropy of the dipolar interaction leads to novel ground states, including a concave (red blood cell) shaped solution in flattened trapping potentials [50][51][52], whose structure can be attributed to the excitation of a roton mode in this system [50,53]. Further work contrasted the dipolar Bogoliubov-de Gennes equations with a variational approach in the pancake geometry [54]. The solutions to the Bogoliubov-de Gennes equations can also be used to study the depletion of the condensate, and in particular how the roton mode affects this important quantity [55]. Understanding the role of dimensionality is important for quantum fluids in general, since different trapping configurations can alter the stability and character of the nonlinear solutions to these systems. Moreover, since fluctuations and interactions are enhanced in lower dimensional quantum systems, these regimes may provide greater insight into quantum fluctuations in general. Several works have investigated droplets in a quasi-onedimensional setting, including in binary mixtures [56,57], and mixtures with coherent [59] and spin-orbit [60] couplings. The related crossover from a bright soliton to a droplet state was also investigated experimentally [61]. It is the aim of this work to understand the regimes of stability and the accompanying dynamics of dipolar droplets in the quasi-one-dimensional setting.
In this work we systematically investigate the solutions to the quasi-one-dimensional dipolar Gross-Pitaevskii model in the beyond-mean-field regime. We begin in Section II by deriving the beyond-mean-field correction to the dipolar Gross-Pitaevskii equation, generalizing previous works to include the effect of the dipole polarization angle. Following this, we study the Bogoliubov-de Gennes equations including the beyond-mean-field contribution, and use this to identify regimes of roton stability in the full parameter space of the model. Following this in Section IV we solve numerically the extended Gross-Pitaevskii equation, examining the form of the solutions in the beyond-mean-field limit. The nature of the modulational instability is then scrutinised, as well as the behaviour of droplet collisions in this system. We conclude with a summary and outlook of our findings in Section V.

A. Dipole-dipole Interactions
We consider a gas of dipolar bosons of mass m interacting through short-range s-wave and long-range dipoledipole interactions. Then, the total atomic interaction potential has the form U (r) = gδ(r) + U dd (r) with where g = 4π 2 a s /m and a s defines the s-wave scattering length, while C dd characterizes the strength of the dipole-dipole interaction, and θ defines the angle between the vector r joining two dipoles and the polarization direction of the dipoles. The atomic interactions are illustrated in Fig. 1(a). The dipole polarization angle θ m = cos −1 (1/ √ 3) defines the 'magic' angle at which the dipolar interaction vanishes. If C dd is positive and θ < θ m the dipoles are orientated in an attractive headto-tail configuration, while for θ > θ m the dipoles lie side-by-side and are repulsive. The strength of the dipolar interaction is typically characterized in terms of the dimensionless parameter ε dd = C dd /3g. The 'anti-dipole' regime, where C dd is negative, has also been proposed in Ref. [62] by performing a rapid rotation of the dipoles, such that the attractive and repulsive regimes are reversed. Recent experimental work [63] indicated that such a scenario can be achieved; however, the condensate lifetime is hampered by a dynamical instability [64,65].

B. Beyond-Mean-Field Effects
The realization of stable droplet phases with highly magnetic 164 Dy [6] has been attributed to quantum fluctuations. Theoretically, quantum fluctuations are formulated in terms of the Lee-Huang-Yang (LHY) correction [66], which within the local density approximation appears as a term proportional to a non-integer power of the atomic density in the appropriate generalized Gross-Pitaevskii equation. In this section we generalize the dipolar LHY term [67] to include the effect of the polarization angle of the dipoles (previous works assumed polarization along the z-axis of the condensate) and find counter-intuitively that this correction remains isotropic, independent of the polarization angle of the dipoles. The LHY correction arises from correcting the otherwise divergent ground state energy of a homogeneous gas of bosons, which for the dipolar Bose gas with density n 3D 0 = N/V is given by Here 0 k = 2 k 2 /(2m) defines the single-particle quasiparticle energy and U dd (k) accounts for the Fourier transform of the dipole-dipole interaction (Eq. (1)), which is defined as [68] where k i defines the cartesian components of the momentum, and α is the dipole polarization angle in the x-z plane (see Fig. 1 (a)). Although the integral defined by Eq. (2) is convergent, it does not in general exist in closed form. After integrating out the radial momentum from Eq. (2), one can show that this expression can be re-written as where Here the short-hand L(θ, φ, α) = sin θ cos φ sin α + cos θ cos α has been used. Although the integral defined by Eq. (5) does not exist in closed form, on physical grounds the parameter ε dd ∼ 1 typically, so we can expand the integrand to quadratic order, which should provide sufficient accuracy for typical experimental parameters. Then taking n = 5/2, Eq. (5) becomes Evaluation of the integrals appearing in Eq.(6) can be accomplished using the result Using Eqs. (4), (6) and (7) one finds that the correction to the ground state energy due to quantum fluctuations becomes Equation (8) is independent of the dipole polarization angle α, a result which will be utilized in this work to understand the interplay of the polarization angle in the beyond-mean-field regime. The two one-loop Feynman diagrams representing the renormalized ground state energy are also shown in Fig. 1 (b). The first diagram (solid lines) shows the real part of the diagonal propagator contribution from the fluctuations, while the second (dashed lines) shows the corresponding imaginary contribution [69].

C. Beyond-Mean-Field Dipolar Bogoliubov-de Gennes Equations
The condensate of N atoms is parametrized by the wave function Ψ(r, t), normalized such that d 3 r|Ψ(r, t)| 2 = N . Including the generalized dipolar LHY correction derived in Section II B, the wave function is described by the generalized dipolar Gross-Pitaevskii equation, written as, Here the mean field dipolar potential is defined as the contribution from the quantum fluctuations treated in the local density approximation, and ω ρ defines the harmonic trapping frequency in the radial ρ 2 = y 2 + z 2 direction. The beyond-mean-field chemical potential is calculated from Eq. (8) using In this work we consider a quasi-one-dimensional dipolar condensate such that the transversal dynamics of the atomic cloud are effectively frozen [70]. Then, a good ansatz for the wave function is Ψ( where a ρ = /mω ρ defines the transverse harmonic length scale. The dimensional reduction is performed by inserting the ansatz for Ψ(r, t) into Eq. (9) and integrating over the transverse area of the atomic cloud. After dropping trivial energy offsets, this yields the quasi-one-dimensional generalized Gross-Pitaevskii equation The dimensionally reduced dipolar interaction appears [71]. The form of the quantum fluctuation term appearing in Eq. (11) is consistent with the derivation and analysis presented in Ref. [20] under the condition for the 1D density (color online) Beyond-mean-field roton analysis. Panel (a) shows the roton unstable regimes in the (ε dd , α) parameter space, with = 5×10 −4 . Colored regions indicate roton unstable regimes for fixed ξn0={750,1000,1250} obtained from Eq. (20) and (21). Panel (b) shows example dispersions plotted from Eq. (17)  n 0 a s 0.6. The typical value of n 0 a s ∼ 10 3 is chosen for the results presented here. To linearize Eq. (11) around the stationary state ψ 0 (x) we introduce the ansatz ψ( (12) describes the small-amplitude fluctuations of ψ(x) and u(x) and v(x) represent the mode functions, ω is the associated excitation frequency and µ is the quasi-onedimensional chemical potential. Then, by inserting the expansion of ψ(x, t), including Eq. (12) into Eq. (11) and assuming the ground state ψ 0 is real, we obtain the coupled equations defines the quasi-one-dimensional Hamiltonian appearing in Eq. (11), while gives the exchange operator. The additional nonlocal operator is defined as . Then, the pair of Bogoliubov-de Gennes equations given by Eq. (13) can be straight-forwardly decoupled. We focus on solutions for the u(x) mode function, which obeys A similar expression for the v(x) mode function can be obtained except with the bracketed terms switched in Eq. (16). Further details concerning the Bogoliubov-de Gennes equations are given in Appendix A.

A. Roton Analysis
In general the eigenvalues of the Bogoluibov de-Gennes equation Eq. (16) must be obtained numerically. However in the homogeneous (infinite) limit one can obtain the spectrum analytically from Eq. (16), since the nonlocal operator χ reduces to the one-dimensional Fourier transform of the dipolar interaction. In this limit we obtain where the excitation energy is = ω and the singleparticle energy appearing in Eq. (17) Accompanying the Bogoluibov de-Gennes energy is the chemical potential of the homogeneous system, given by where the homogeneous dipolar potential is defined as In what follows we express dimensions in terms of the natural units of the homogeneous system: the unit of energy is the chemical potential µ 0 , the unit of length is the healing length ξ = / m|µ 0 |, and the unit of time is /|µ 0 |. We use the excitation spectrum defined by Eq. (17) along with the definition of the chemical potential, Eq. (18) to gain an understanding of the properties of the dipolar condensate in the beyond-mean-field regime. In the absence of the LHY correction, the excitation energies defined by Eq. (17) exhibit different regimes of physical behaviour depending on the choice of parameters. In the limit of small k x , the dispersion ω k = k x c s is phonon-like, depending linearly on momentum such that Here the LHY term contributes an additional density dependence to the speed of sound c s in Eq. (19), giving an increased value of c s in the high density phase. The point at which the roton minima touches the k x = 0 axis can be calculated in the following manner. First, the (squared) dispersion relation Eq. (17) is differentiated with respect to k x and set equal to zero such that ∂ 2 /∂k x = 0 to obtain the two family of extrema from the dispersion, the maxon and the roton. Since we are interested in the roton, we can remove the maxon by combining this expression with the value of the dispersion set equal to zero, yielding a quartic equation in the square of the momenta k 2 x . This can be solved analytically to obtain The two functions A and B that carry the dependence of the physical parameters in the problem appearing in Eqs. (20) and (21) are defined as Then for a given set of physical parameters we can compute the solutions to Eqs. (20) and (21) numerically using an iterative procedure. Additionally, two other parameters emerge from the dimensionless analysis, the ratio of the transverse harmonic length a ρ and the healing length ξ, defined as σ = a ρ /ξ. The second is the ratio of the van der Waals a s and harmonic lengths a ρ which arrises from the LHY term, and is defined as = a s /a ρ . This second dimensionless parameter has a typical value of 10 −3 for dipolar gases, where the van der-Waals length a s 100a 0 and a typical radial harmonic length is a ρ 1µm. In a previous work [36] it was shown that for the mean-field case (γ QF = 0) the roton instability will only appear for a s > 0 when σ 0.8, hence in what follows we assume the arbitrary value σ = 1. Figure 2 explores the behaviour of the roton instability in the beyond-mean-field regime. Panel (a) shows the numerical solutions obtained from Eqs. (20) and (21) in the (ε dd , α) parameter space, and we take = 5 × 10 −4 . For γ QF = 0, two roton instabilities are observed in this system, due to the underlying quadratic dependence on ε dd of the dispersion relation, Eq. (17). Including higher orders of the expansion in Eq. (5), like those presented in Ref. [10], will lead to the appearance of more roton instabilities, however the boundary and nature of the solution for the smallest |ε dd | ∼ 1 roton instability boundary-the solution closest to experimentally realizable parameters-will be almost unchanged. Each shaded area corresponds to a different atomic density, with increasing density causing the unstable region to shrink in area. In this way the LHY term acts to stabilize the parameter space as the density is increased. The gray shaded area corresponds to the (single-valued) mean-field result in the limit γ QF = 0. The roton unstable regions also appear for α = π/2, although these are smaller in area than their α = 0 counterparts, which is attributed to the dipoles being in a side-by-side repulsive alignment, which additionally reduces the roton unstable region of the parameter space. Example dispersion relations plotted using Eq. (17) are depicted in Fig. 2 (b). Individual curves correspond to the cross, circle and triangle markers showing the dispersion slightly below (cross), above (triangle) and at (circle) the roton instability respectively. In these examples the density ξn 0 = 1250, and the dipole polarization α = 0. The last panel (c) of Fig. 2 plots semi-logarithmically the value of the dipolar interaction strength ε dd against the density ξn 0 at which the roton manifests for several examples of fixed . Each curve terminates at a maximum value of n 0 , due to the additional repulsive LHY term overwhelming the attractive part of the dipolar interaction. Solutions for both α = 0 and π/2 are displayed, with the α = π/2 solutions shifted to lower n 0 ; attributed again to the repulsive nature of this side-by-side polarization. The presence of multiple roton instabilities enriches the prospects for supersolidity within this system, which are a necessary condition for such a state to occur.

B. Droplet Phases
The existence of the droplet phases in the (ε dd , α) parameter space depends on the balance between attractive and repulsive forces in the system. We can understand this by examining Eq. (18), the homogeneous chemical potential. In the limit n 3D 0 a 3 s 1 the boundary between the repulsive and attractive regions of the parameter space are independent of the density when µ 0 = 0, which gives ε dd = 4/[1 + 3 cos(2α)]. For n 3D 0 a 3 s ∼ 1 one can instead obtain two beyond-mean-field solutions to µ[n 0 ] = 0 as where we have defined C 0 = (128/15π) 3/2 √ σ √ ξn 0 + 1 and C 1 = [1+3 cos(2α)]/4. Figure 3 explores the droplets existence in the beyond-mean-field regime. Panel (a) shows the repulsive (µ 0 > 0) and attractive (µ 0 < 0) shaded regions, which lead to homogeneous and bright solitons respectively in the quasi one-dimensional setting. Panel (b) shows in the beyond-mean-field regime (n 3D 0 a 3 s ∼1) which for repulsive interactions gives a homogeneous ground state. However, when the mean-field phonon instability is crossed, the system remains homogeneous (blue regions). Droplets are found beyond a critical value of ε dd obtained from Eq. (23). For α < α m the droplet region is larger than for α > α m , attributed to the head-to-tail (attractive) polarization. Increasing the atomic density n 0 has the effect of enlarging (shrinking) these regions for α < α m (α > α m ). There is no roton instability in this analysis due to the choice of σ = 0.2.

A. Single Droplets
To calculate the solutions of the generalized dipolar Gross-Pitaevskii equation Eq. (11) we employ a psuedospectral approach, the Fourier split-step method, for the results presented in this section. Due to the size of the parameter space of the extended dipolar model, which includes the strength of the dipolar interactions, the polarization angle of the dipoles, as well as the number of atoms, it is instructive to consider the ground state of Eq. (11) by fixing two of these physical parameters whilst varying the other. Further, in using the healing units (see Sec. IIIA ) the density n 0 appearing in the healing units is chosen such that n 0 = max(|ψ(x, t)| 2 ). Figure 4 presents examples of the droplets spatial density |ψ(x)| 2 as a function of the dipolar interaction strength, ε dd in Fig. 4(a) and (d) for the choices of polarization angle α = 0 and π/2 respectively. In both cases the appearance of the droplet state is not achieved for arbitrary dipole strength, but only manifests after the beyond-mean-field phonon instability is crossed, rather than the usual mean-field phonon instability corresponding to the point where the interactions become attractive in the low density limit (n 3D 0 a 3 s 1). In both Fig. 4(a) and (d) this point is indicated by a dashed red line obtained from the the solutions ε ± dd , Eq. (23). This leads to a pair of solutions for a given set of parameters, which correspond to α < α m and α > α m respectively. The width of the computed ground state is observed to be sensitive to the dipolar strength, ε dd . To understand the effect of changing the number of atoms in the droplet, we solve Eq. (11) at fixed dipolar interaction strength. Then, the dotted lines in Fig. 4(a) and (d) correspond to the values ε dd = 2 and ε dd = −4 used in in panels (b) and (e) respectively. For both polarization angles, it is found that the width of the droplet w drop increases linearly with atom number such that w drop ∝ N . The width of the droplet is largest for the α = π/2 polarization angle due to the increased influence of repulsive interactions. Panel Fig. 4(c) shows example ground states taken from panel (b) for increasing atom number N . For low atom numbers (black data) the solutions resemble a bright soliton (sech-like profile) while for increasing atom number the profiles widen, developing the characteristic flat top associated with the droplet state (red and blue data). Panel (f) meanwhile shows the energy calculated from the definition as a function of N for the data presented in panel (b) and (e). The ground state energy decreases monotonically with N and dE/dN < 0 throughout, adhering to the Vakhitov-Kolokolov stability critereon for stationary solutions of self-attractive nonlinear waves [72].

B. Modulation Instability
The strength of the different nonlinearities directly determines the nature of the state that the dipolar condensate is in. By changing the a s scattering length from an initial value of a i s to a final value a f s such that a f s < a i s , a modulational instability can be induced. The instability originates from long-wavelength perturbations that cause the break-up of a waveform into pulses; coming from the growth of nonlinear excitations in the system. This effect has been used previously to investigate experimentally the stability of individual bright solitons [73,74] as well as dipolar droplets in a trapped three-dimensional context [75]. Figure 5 investigates the effect of performing an interaction quench on an initial dipolar droplet with ε dd = 2, α = 0 and N = 10 7 . The scattering length takes the time-dependent form  where t Q defines the time at which the quench is applied, while H(t) is the Heaviside function. In our simulations we take t Q = 0.15t f . Figure 5(a) shows the atomic density |ψ(x, t f )| 2 of the dipolar gas as a function of the quench strength a i s /a f s taken at the final point of the numerical integration t = t f . For modest quench strengths (a i s /a f s 1.1) the droplets shape remains intact, with the exception of the excitation of surface waves. Above a critical quench strength the droplet undergoes the modulation instability, and in general breaks apart into smaller droplets and dipolar bright solitons, as well as the emission of low density radiation. Panel (b) and (c) show respectively the quench protocol of Eq. (25) and example dynamics for a i s /a f s = 1.12. The lower row of panels (d-g) of Fig. 5 show individual examples of the quench dynamics for increasing a i s /a f s . Here we observe both even and odd numbers of droplets, while panel (f) shows a situation where two droplets and two bright solitons are created after the quench. For larger quench strengths, such as that presented in panel (g), a central droplet is produced along with increasing amounts of radiation, in this example short-lived bright soliton bound states are also observed. In general we find that the onset of the modulational instability (and the final state after the quench) depend strongly on the atom number. If for example the initial number of atoms in the droplet is reduced, then the modulational instability occurs at larger values of a i s /a f s compared to a larger initial atom number. The generation of multiple droplets as presented here using the modulational instability relies on being able to tune the scattering length a s of the condensate which can in turn lead to significant atomic losses. To address this, there are proposals to produce condensates with attractive interactions with reduced noise and greater control over the final experimental state of the system [76,77].

C. Collisional Population Transfer and Droplet Fission
The coherent nature of the superfluid state provides a convenient tool to explore quantum mechanical phenomena at macroscopic length scales. One striking manifestation of matter-waves coherence is the so-called Josephson effect. Here, two superconductors or superfluids which are separated by an insulting barrier can experience a current, originating from atomic tunneling between the two superconductors/superfluids. The dipolar droplets represent an interesting addition to the superfluid family, since they are effectively an isolated (finite) region of homogeneous fluid, so understanding their binary dynamics is expected to yield novel phenomena. To investigate the basic physics of binary droplet dynamics, we perform simulations with an initial state of the form This constitutes a symmetric state comprising two droplets whose centres are initially separated by a distance x + − x − = 60ξ, which are traveling towards each other at constant velocity v ± = ±v 0 . The initial phase difference δ = δ + − δ − between the two droplets is δ ∈ [0, 2π]. In Figure 6 we present results of droplet collisions using the initial state defined by Eq. (26). We take for the physical parameters ε dd = 2, N = 3 × 10 6 , = 10 −3 , σ = 0.2, and the total length of each numerical integration is t f = 140 /|µ 0 |. Figures 6 (a-f) show space-time plots for different initial phase differences, δ and initial velocity mξv 0 / = 0.25. We observe that post collision two droplets emerge -with different populations (and sizes) that depend on the choice of initial phase difference. For example, choosing δ = π (panel 6 (e)) produces two droplets with an equal number of atoms present in each droplet post collision. The presence of excitations in the form of sound waves can be seen here post collision, reflecting back and forth inside each droplet (viz. Eq. (19)), which is indicative of nonintegrable dynamics. To quantify this change in population, we can calculate the population difference between the droplets at t = t f as a function of the initial phase difference δ(t 0 ). The population difference is defined as Here each integral computes the number of atoms in the individual droplets between the edges of the droplet given by x = x j± . Panel 6 (g) shows some example results with different choices of initial velocity, mξv 0 / = 0.15, 0.25, 0.35 (triangle, circle and square markers). Only droplet collisions for mξv 0 / = 0.25 represents a situation where two droplets emerge post collision for the full range of phase differences between δ = 0 and δ = π, due to the existence of bound states (droplet molecules) and droplet fission at smaller and larger initial velocities respectively, breaking droplet number conservation. As such, a parameter window exists where one can compare the collisional transfer to the Josephson effect. Then, the semi-classical Josephson equations for a superfluid are written as [78] d dt Here, the parameter J describes the strength of the coupling between the two superfluids. We will use the final integration time t = t f to fit our model to the numerical simulations of the extended GPE (green circles), since the total integration time determines the nature of the observed Josephson oscillations. Panel 6 (g) shows a comparison between the populations of the droplets at t = t f , calculated from the numerical solutions to the Josephson equations. The blue and black dashed lines are obtained from the pair of Josephson relations Eq. (28), for t f = 0.5 /J and t f = 1.75 /J respectively. For small t f , (blue dashed) the oscillation is linear, and does not follow the extended GPE data (green circles). For larger t f , (black dashed) the oscillation exhibits a stronger nonlinear character, and follows the extended GPE data quite well. It would be interesting to study this effect in more detail, to understand if this analogy with the Josephson effect can be extended, for example into the highly nonlinear regime.
To understand the effect of the droplets initial velocity on the dynamics, Figure 7 shows simulations of droplet collisions with fixed initial phase difference, for δ = 0, π. Panel (a) shows a long-lived droplet dimer formed from two initially out-of-phase stationary droplets. Then, example dynamics are shown in panels (b-e) for different finite initial velocities (see individual captions). The in-phase collisions demonstrate that the droplets undergo fission [79][80][81] with multiple droplet states produced post-collision, shown in panels (b) and (d). For the out-of-phase collisions (δ = π) two out-going droplets are observed, for low incoming velocity (panel (c)) where there are small amounts of sound produced post collision. For a greater initial velocity (panel (e)), larger amounts of sound and radiation are produced, still with two outgoing droplets. Finally panel (f) computes the number of droplets post collision for in and out-of-phase collisions, as a function of the initial velocity. For even larger velocities (mξv 0 / ∼ 1) the delicate balance of kinetic and potential energy that maintains the droplet is violated, and the droplet breaks apart into atomic radiation.

V. CONCLUSIONS AND OUTLOOK
In this work we have investigated the properties of a quasi-one-dimensional dipolar Bose-Einstein condensate in the presence of quantum fluctuations. It was shown that in the beyond-mean-field regime, the polarization angle does not contribute an angular dependence to the LHY term which stabilizes the gas against collapse. By calculating the excitations of this system within the framework of the Bogoliubov-de Gennes formalism, we identified regimes unstable to the roton instability. Interestingly, the roton unstable regions are found in general to appear in pairs for a given dipole polarization angle; this is due to the underlying quadratic dependence on the dipolar strength from the LHY contribution to the excitation energy.
We examined the nature of the ground states of this system, observing the appearance of droplet phases in regions of the parameter space where the total interactions are net attractive. By applying an interaction quench to a single large droplet, the nature of the modulation instability was investigated. It was found that for moderate quenches, both even and odd numbers of droplets can be generated. For larger quenches, bright solitons and increasing amounts of radiation are produced, suggesting a window of parameters for useful quenches. The collisional properties of droplets were also explored. By modulating the initial phase difference between the droplets, atomic population transfer was observed between the droplets. This was interpreted and compared with the superfluid Josephson effect, finding good agreement in a window of the parameter space. Droplet fission was also observed as the initial velocity of the droplet was increased.
It would be interesting in the future to understand how a harmonic trap changes the physics of the onedimensional droplet, and in particular how the interplay of quenching both the interactions and trapping strength changes the number of droplets produced. One could also use this model to understand supersolid phases in the one-dimensional context, as well as studying droplets and their dynamics with models that incorportate higherdimensional effects [82]. Finally, it would also be beneficial to understand the lifetime of the one-dimensional droplet, which can be computed from three-body atomic recombinations [83].