Quantum droplets of quasi-one-dimensional dipolar Bose-Einstein condensates
QQuantum droplets of quasi-one-dimensional dipolar Bose-Einstein condensates
M. J. Edmonds, T. Bland, and N. G. Parker Department of Physics & Research and Education Center for Natural Sciences,Keio University, Hiyoshi 4-1-1, Yokohama, Kanagawa 223-8521, Japan Joint Quantum Centre (JQC) Durham–Newcastle,School of Mathematics, Statistics and Physics, Newcastle University,Newcastle upon Tyne, NE1 7RU, United Kingdom (Dated: November 20, 2020)Ultracold dipolar droplets have been realized in a series of ground-breaking experiments, where thestability of the droplet state is attributed to beyond-mean-field effects in the form of the celebratedLee-Huang-Yang (LHY) correction. We scrutinize the dipolar droplet states in a one-dimensionalcontext using a combination of analytical and numerical approaches, and identify experimentally vi-able parameters for accessing our findings for future experiments. In particular we identify regimesof stability in the restricted geometry, finding multiple roton instabilities as well as regions sup-porting quasi-one-dimensional droplet states. By applying an interaction quench to the droplet, amodulational instability is induced and multiple droplets are produced, along with bright solitonsand atomic radiation. We also assess the droplets robustness to collisions, revealing populationtransfer and droplet fission.
I. INTRODUCTION
Bose-Einstein condensates possessing long-ranged,anisotropic dipole-dipole interactions have been realizedin a series of ground-breaking experiments with highlymagnetic atoms. The first generation of experimentsachieved condensation of Cr [1, 2],
Dy [3, 4] and
Er [5]. Recently, a second series of experiments hasachieved a quantum analogue of the classical Rosensweiginstability [6], as well as the realization of droplet states[7, 8] – where the gas enters a high density phase whosestability has been attributed to the influence of quan-tum fluctuations [9–13]. Dipolar condensates constituteweakly correlated systems, and can exhibit propertiesand behaviour similar to that of a liquid in the beyond-mean-field limit [18–22]. Droplet states have also beenrealized in condensate mixtures [23, 24] – supported by abalance between attractive s -wave interactions betweenthe atoms and quantum fluctuations – and photonic sys-tems [25] – here the nonlinear medium provides a repul-sive d -wave contribution that stabilizes the light beam.The quantum liquid Helium II is well-known to ex-hibit a roton minimum in its excitation spectrum; this issupported by the strong interatomic interactions and cor-relations [26], where roton excitations typically occur atwavelengths comparable to the average inter-particle sep-aration, indicating that the superfluid is close to forminga crystalline structure [27, 28]. Although dipolar con-densates are weakly correlated, the nonlocal characterof the dipolar interaction supports roton-like excitations[14, 15], which have now been experimentally realizedin a gas of Er [16, 17]. Here the roton lengthscaleis dictated by the geometry of the gas. A plethora oftheoretical investigations have focussed on detailing thecorrespondence between rotons in the Helium II phaseand weakly interacting dipolar gases. Early work exam-ined the possibility of roton excitations in a quasi-one-dimensional setting [29], as well as the manifestation of rotons in a rotating dipolar condensate [30] and the iden-tification of a roton mode in trapped dipolar condensates– leading to the identification of ‘roton fingers’ [31], adiscrete manifestation of the instability found in the ho-mogeneous system.The existence of the roton mode indicates the proxim-ity of the quantum fluid to long-range crystalline order.If the quantum fluid can simultaneously support a super-fluid state, the system is a candidate for a supersolid –a phase of matter where these two forms of order coexist[32–35] (see also the review of Leggett, Ref. [36]). Thepossible unambiguous detection of a supersolid state inHelium has attracted long debate. Very recently, state-of-the-art experiments with dipolar condensates have re-ported supersolid phases in these systems in a quasi-one-dimensional setting [37, 38]. Complementary theoreticalinvestigations in lower-dimensional dipolar gases have ex-plored nonlinear wave structures in the form of dark soli-tons [39–42] and bright solitons [43–47]. Related the-oretical proposals have also examined the possibility ofsupersolids in binary condensates [48], whose existence issupported by a stripe-like phase.Within the weakly-interacting regime, dilute atomiccondensates are well described by the celebrated Gross-Pitaevskii description [49, 50]. A natural extensionto this is the Bogoliubov-de Gennes formalism, whichconstitutes the linear response theory of the nonlinearSchr¨odinger equation, and gives insight into the behav-ior of the elementary excitations and collective modesof the condensate [53]. Early work demonstrated thatthis approach could also be applied to dipolar conden-sates [54], revealing the effect of the dipolar interactionon the excitations. Such an approach is justified when thequantization of excitations is not required. For systemsin the beyond mean-field regime it is an open questionas to whether such an approach is justified. Nonethelessthere are several works that use this approach to studythe collective behavour of quantum droplets [12, 55, 56]. a r X i v : . [ c ond - m a t . qu a n t - g a s ] N ov r θ αz x ( a ) ( b ) + FIG. 1. Schematic representation of the atomic interac-tions (a). Each dipole (blue and red arrow) is separated bya distance r , while θ defines the angle between the dipolesand the polarization direction, per Eq. (1). The angle α defines the polarization direction of the dipoles in the x - z plane. The isotropic contact interactions are represented bythe gray spheres. (b) shows the two one-loop contributions tothe ground state energy, Eq. (2). Knowledge of the excitation spectrum gives insight intomany important properties of these systems – in manyphysical effects the dimensionality of the system playsa key role in the dipolar condensates overall behavior[57]. Dipolar condensates with spin degrees of free-dom have also been analyzed within the Bogoliubov-deGennes framework, exploring the interplay of the excita-tions with the magnetic phases inherent to these systems[58]. The anisotropy of the dipolar interaction leads tonovel ground states, including a concave (red blood cell)shaped solution in flattened trapping potentials [59–61],whose structure can be attributed to the excitation ofa roton mode in this system [59, 62]. Further work con-trasted the dipolar Bogoliubov-de Gennes equations witha variational approach in the pancake geometry [63]. Thesolutions to the Bogoliubov-de Gennes equations can alsobe used to study the depletion of the condensate, andin particular how the roton mode affects this importantquantity [64].Parallel to the ongoing experiments with atomic dipo-lar atoms, the realization of cold molecular gases has alsorevealed novel phenomenology. Different to their atomiccounterparts, atomic molecules possess additional vibra-tional and rotational degrees of freedom, which compli-cates their manipulation and cooling [65]. Many differ-ent molecules have now been cooled, including RbCs [66],fermionic NaK [67] and NaRb [68]. The rapid develop-ment of this field has seen the achievement of the molec-ular rovibrational ground states of KRb [69], controlledquantum chemical reactions [70], and molecular collisions[71]. The ability to control the dimensionality of ensem-bles of these molecules using optical lattices [72] has af-forded a new route towards two-dimensional systems [73],and quantum magnetism with molecular dipole-dipole in-teractions [74].Understanding the role of dimensionality is impor-tant for quantum fluids in general, since different trap-ping configurations can alter the stability and character of the nonlinear solutions to these systems. Moreover,since fluctuations and interactions are enhanced in lowerdimensional quantum systems, these regimes may pro-vide greater insight into quantum fluctuations in general.Several works have investigated droplets in a quasi-one-dimensional setting, including in binary mixtures [75, 76],and mixtures with coherent [77] and spin-orbit [78] cou-plings. The related crossover from a bright soliton to adroplet state was also investigated experimentally [79]. Itis the aim of this work to understand the regimes of sta-bility and the accompanying dynamics of dipolar dropletsin the quasi-one-dimensional setting.In this work we systematically investigate the solutionsto the quasi-one-dimensional dipolar Gross-Pitaevskiimodel in the beyond-mean-field regime. We begin in Sec-tion II by studying the extended dipolar Gross-Pitaevskiimodel, from which we derive the Bogoliubov-de Gennesequations including the beyond-mean-field contribution,and use this to identify regimes of roton stability inthe full parameter space of the model. Following thisin Section IV we solve numerically the extended Gross-Pitaevskii equation, examining the form of the solutionsin the beyond-mean-field limit. The nature of the modu-lational instability is then scrutinised, as well as the be-haviour of droplet collisions in this system. We concludewith a summary and outlook of our findings in SectionV.
II. THEORETICAL MODELA. Dipole-dipole Interactions
We consider a gas of dipolar bosons of mass m inter-acting through short-range s -wave and long-range dipole-dipole interactions. Then, the total atomic interactionpotential has the form U ( r ) = gδ ( r ) + U dd ( r ) with U dd ( r ) = C dd π − θ | r | (1)where g = 4 π (cid:126) a s /m and a s defines the s -wave scat-tering length, while C dd characterizes the strength of thedipole-dipole interaction, and θ defines the angle betweenthe vector r joining two dipoles and the polarization di-rection of the dipoles. The atomic interactions are il-lustrated in Fig. 1(a). The dipole polarization angle θ m = cos − (1 / √
3) defines the ‘magic’ angle at whichthe dipolar interaction vanishes. If C dd is positive and θ < θ m the dipoles are orientated in an attractive head-to-tail configuration, while for θ > θ m the dipoles lieside-by-side and are repulsive. The strength of the dipo-lar interaction is typically characterized in terms of thedimensionless parameter ε dd = C dd / g . The ‘anti-dipole’regime, where C dd is negative, has also been proposed inRef. [80] by performing a rapid rotation of the dipoles,such that the attractive and repulsive regimes are re-versed. Recent experimental work [81] indicated thatsuch a scenario can be achieved; however, the condensatelifetime is hampered by a dynamical instability [82, 83].Since we consider parameter regimes where the sign of ε dd can be either positive or negative, we consider stateswith the possibility of differing signs of the parameters C dd and a s , which can be accessed by rapid rotation ofthe dipoles and with the powerful tool of optical Feshbachresonances. B. Beyond-Mean-Field Dipolar Bogoliubov-deGennes Equations
The realization of stable droplet phases with highlymagnetic
Dy [6] has been attributed to quantum fluc-tuations . Theoretically, quantum fluctuations are formu-lated in terms of the Lee-Huang-Yang (LHY) correction[84], which within the local density approximation ap-pears as a term proportional to a non-integer power ofthe atomic density in the appropriate generalized Gross-Pitaevskii equation.To obtain the appropriate correction to the dipolarGross-Pitaevskii equation, one begins with the many-body Hamiltonian for a gas of homogeneous dipoles andfrom this computes the ground state energy of the sys-tem. Such a situation was originally studied by the au-thors of Refs. [85]. The technical details and analysisfor this are presented in Appendix A, the result of whichgives E QF V = 6415 g (cid:0) n (cid:1) (cid:114) n a s π (cid:18) ε (cid:19) . (2)Equation (2) is independent of the dipole polarizationangle α , a result which will be utilized in this work tounderstand the interplay of the polarization angle in thebeyond-mean-field regime. Here we adopt a quadratic ap-proximation for the beyond-mean-field contribution; thisis primarily motivated to allow us to perform a thoroughanalytical analysis of the dipolar system in the beyond-mean-field limit. Further details supporting our choiceof Eq. (2) are given in Appendix A.To obtain a finite ground state energy for the homo-geneous dipolar system (see Eq. (A5)), a renormaliza-tion procedure is required (see Appendix A for furtherdetails) in order to make the momentum integrals con-vergent, and hence physically meaningful. Then, the twoone-loop Feynman diagrams representing the renormal-ized ground state energy are shown in Fig. 1 (b). The firstdiagram (solid lines) shows the real part of the diagonalpropagator contribution from the fluctuations, while thesecond (dashed lines) shows the corresponding imaginarycontribution [88].The condensate of N atoms is parametrized bythe wave function Ψ( r , t ), normalized such that (cid:82) d r | Ψ( r , t ) | = N . Including the generalized dipo-lar LHY correction derived in Appendix A, the wavefunction is described by the generalized dipolar Gross- Pitaevskii equation, written as [13], i (cid:126) ∂ Ψ ∂t = (cid:20) ˆ p m + 12 mω ρ ρ + g | Ψ | + Φ dd [Ψ] + µ QF (cid:21) Ψ . (3)Here the mean field dipolar potential is defined asΦ dd [Ψ( r , t )] = (cid:82) d r (cid:48) U dd ( r − r (cid:48) ) | Ψ( r (cid:48) , t ) | , µ QF defines thedensity-dependent contribution from the quantum fluc-tuations treated in the local density approximation, and ω ρ defines the harmonic trapping frequency in the radial ρ = y + z direction. The beyond-mean-field chem-ical potential is calculated from Eq. (2) using µ QF = ∂E QF /∂N = γ QF n ( r , t ) / where the effective strength γ QF is defined γ QF = 32 g (cid:114) a s π (cid:18) ε (cid:19) . (4)In this work we consider a quasi-one-dimensional dipo-lar condensate such that the transversal dynamics ofthe atomic cloud are effectively frozen [89]. Then,a good ansatz for the wave function is Ψ( r , t ) =( a ρ √ π ) − exp( − ρ / a ρ ) ψ ( x, t ) where a ρ = (cid:112) (cid:126) /mω ρ de-fines the transverse harmonic length scale. The dimen-sional reduction is performed by inserting the ansatz forΨ( r , t ) into Eq. (3) and integrating over the transversearea of the atomic cloud. After dropping trivial energyoffsets, this yields the quasi-one-dimensional generalizedGross-Pitaevskii equation i (cid:126) ∂ψ∂t = (cid:20) ˆ p x m + g πa ρ | ψ | + Φ + 2 γ QF π / a ρ | ψ | (cid:21) ψ. (5)By writing Eq. (5) a number of approximations have beenused. For the gas to be in the quasi-one-dimensionallimit we require a ρ /ξ < ξ is the appro-priate healing length, see Sec. III for further details)while the beyond mean field treatment requires a s n (cid:38) . | ε dd | ∼
1. The parameters we use in our model aresimilar to those used in two recent dipolar experiments[37, 38, 91], in these works the trapping geometry gives aratio of length scales a z /a x = (cid:112) ω x /ω z (cid:39) .
47, whilethe strength of the dipolar interaction for the atomicspecies
Dy takes a typical value of ε dd ∼ .
2, sim-ilar to the values we consider through the course ofthis work. The dimensionally reduced dipolar inter-action appears as Φ [Ψ] = (cid:82) d x (cid:48) U ( x − x (cid:48) ) | ψ ( x (cid:48) ) | ,where the real-space form of U is given by U ( x ) = U U ( | x | /a ρ ) with U = C dd (1 + 3 cos 2 α ) / (32 πa ρ ) and U ( u ) = [2 u − √ π (1 + u ) e u / erfc( u/ √ / δ ( u )[92, 93]. Other works have even considered a generalthree-dimensional dipole polarization [94]. The form ofthe quantum fluctuation term appearing in Eq. (5) isconsistent with the derivation and analysis presentedin Ref. [20], who explored how the beyond mean-fieldterm changes as the density a s n is changed. An attrac-tive regime − n LHY was found at low densities, while for n a s (cid:38) . n / term is recovered. In ourwork, the typical value of n a s ∼
10 is associated withthe results presented throughout our work. To linearizeEq. (5) around the stationary state ψ ( x ) we introducethe ansatz ψ ( x, t ) = [ ψ ( x )+ δψ ( x, t )] exp( − iµt/ (cid:126) ) where δψ ( x, t ) = [ u ( x ) + v ( x )] e − iωt − [ u ( x ) − v ( x )] ∗ e iωt (6)describes the small-amplitude fluctuations of ψ ( x ) and u ( x ) and v ( x ) represent the mode functions, ω is theassociated excitation frequency and µ is the quasi-one-dimensional chemical potential. Then, by inserting theexpansion of ψ ( x, t ), including Eq. (6) into Eq. (5) andassuming the ground state ψ is real, we obtain the cou-pled equations (cid:20) H GP1D − µ + 2 MH GP1D − µ (cid:21) (cid:20) u ( x ) v ( x ) (cid:21) = (cid:126) ω (cid:20) u ( x ) v ( x ) (cid:21) (7)where H GP1D = ˆ p x m + g πa ρ | ψ | + Φ + 2 γ QF π / a ρ | ψ | (8)defines the quasi-one-dimensional Hamiltonian appearingin Eq. (5), while M = g πa ρ | ψ | + 3 γ QF π / a ρ | ψ | + χ (9)gives the exchange operator. The additional nonlocaloperator is defined as χ [ f ( x )] ≡ ψ ( x ) (cid:82) d x (cid:48) U ( x − x (cid:48) ) ψ ( x (cid:48) ) f ( x (cid:48) ). Then, the pair of Bogoliubov-de Gennesequations given by Eq. (7) can be straight-forwardly de-coupled. We focus on solutions for the u ( x ) mode func-tion, which obeys (cid:20) H GP1D − µ + 2 M (cid:21)(cid:20) H GP1D − µ (cid:21) u ( x ) = ( (cid:126) ω ) u ( x ) . (10)A similar expression for the v ( x ) mode function can beobtained except with the bracketed terms switched inEq. (10). Further details concerning the Bogoliubov-deGennes equations are given in Appendix B. III. HOMOGENEOUS ANALYSISA. Roton Analysis
In general the eigenvalues of the Bogoluibov de-Gennesequation Eq. (10) must be obtained numerically. How-ever in the homogeneous (infinite) limit one can obtainthe spectrum analytically from Eq. (10), since the non-local operator χ reduces to the one-dimensional Fouriertransform of the dipolar interaction. In this limit we ob-tain (cid:15) = (cid:15) k +2 n (cid:15) k (cid:20) a ρ U V (cid:18) k x a ρ (cid:19) + g πa ρ + 3 γ QF √ n π / a ρ (cid:21) (11) where the excitation energy is (cid:15) = (cid:126) ω and the single-particle energy appearing in Eq. (11) is (cid:15) k = (cid:126) k x / m .The term describing the Fourier transform of the dipolarinteraction is defined as V ( u ) = ue u E ( u ) − / E ( u ) = (cid:82) ∞ u d t t − e − t defines the exponential integral.Accompanying the Bogoluibov de-Gennes energy is thechemical potential of the homogeneous system, given by µ [ n ] = n g πa ρ + Φ + 25 π / a ρ γ QF n / , (12)where the homogeneous dipolar potential is defined asΦ = − ε dd gn [1 + 3 cos 2 α ] / πa ρ . In what follows weexpress dimensions in terms of the natural units of thehomogeneous system: the unit of energy is the chemicalpotential µ , the unit of length is the healing length ξ = (cid:126) / (cid:112) m | µ | , and the unit of time is (cid:126) / | µ | . Note that as µ → ξ → ∞ , an inherent pathologyof this choice of units. Obviously we cannot (and shouldnot) simulate this point of the parameter space, but aswe shall see we can instead simulate small | µ | close tothe transition to the droplet state, obtaining physicallysensible results.We use the excitation spectrum defined by Eq.(11)along with the definition of the chemical potential,Eq.(12) to gain an understanding of the properties ofthe dipolar condensate in the beyond-mean-field regime.In the absence of the LHY correction, the excitation en-ergies defined by Eq. (11) exhibit different regimes ofphysical behaviour depending on the choice of parame-ters. In the limit of small k x , the dispersion ω k = k x c s isphonon-like, depending linearly on momentum such that ω k (cid:39) k x (cid:26) − a ρ n U n g πa ρ + 35 π / a ρ γ QF n / (cid:27) . (13)Here the LHY term contributes an additional density de-pendence to the speed of sound c s in Eq. (13), giving anincreased value of c s in the high density phase. The pointat which the roton minima touches the k x = 0 axis can becalculated in the following manner. First, the (squared)dispersion relation Eq. (11) is differentiated with respectto k x and set equal to zero such that ∂(cid:15) /∂k x = 0 toobtain the two family of extrema from the dispersion,the maxon and the roton. Since we are interested in theroton, we can remove the maxon by combining this ex-pression with the value of the dispersion set equal to zero,yielding a quartic equation in the square of the momenta k x . This can be solved analytically to obtain k c − B ξ (cid:26)(cid:20)
1+ 2 A B (cid:21) − (cid:115)(cid:20)
1+ 2 A B (cid:21) − B σ (cid:20) − A B (cid:21)(cid:27) (14)which is solved simultaneously with k c n ξ (cid:20) AV (cid:18) k c a ρ (cid:19) + B (cid:21) = 0 . (15) ξn ξn ξn γ QF=0 . . . . . − − α/π . . − . . . . ξk x ε dd = 0 . ε RIdd ε dd = ε RIdd ε dd = 1 . ε RIdd − − ξn ‘ = 2 × − ‘ = 1 × − ‘ = 5 × − ( a ) ( b ) ( c ) α =0 α = π ε dd ( (cid:15) / µ ) ε dd FIG. 2. Beyond-mean-field roton analysis. Panel (a) shows the roton unstable regimes in the ( ε dd , α ) parameter space, with (cid:96) = 5 × − . Colored regions indicate roton unstable regimes for fixed ξn = { } (blue, green, red respectively)obtained from Eq. (14) and (15). Panel (b) shows example dispersions plotted from Eq.(11) corresponding to the pointsindicated by the cross, circle and triangle on panel (a). Panel (c) shows how the critical roton ε dd changes with density ξn for (cid:96) = { , , } × − . The two functions A and B that carry the dependenceof the physical parameters in the problem appearing inEqs. (14) and (15) are defined as A ≡ mξ (cid:126) a ρ U , (16a) B ≡ mξ (cid:126) (cid:18) g πa ρ + 35 π / a ρ γ QF √ n (cid:19) . (16b)Then for a given set of physical parameters we can com-pute the solutions to Eqs. (14) and (15) numericallyusing an iterative procedure. Additionally, two other pa-rameters emerge from the dimensionless analysis, the ra-tio of the transverse harmonic length a ρ and the heal-ing length ξ , defined as σ = a ρ /ξ . The second is theratio of the scattering length a s and harmonic lengths a ρ which arrises from the LHY term, and is defined as (cid:96) = a s /a ρ . This second dimensionless parameter has atypical value of (cid:96) (cid:39) − for dipolar gases, where thescattering length a s (cid:39) a and a typical radial har-monic length is a ρ (cid:39) µm . In a previous work [41]it was shown that for the mean-field case ( γ QF = 0)the roton instability will only appear for a s > σ (cid:38) .
8, hence in what follows we assume the arbitraryvalue σ = 1. Alternatively to our approach one can alsoassess the quasi-one-dimensional character of the conden-sate by individually considering the healing lengths asso-ciated with the contact and dipolar interactions, ξ CI and ξ D respectively. If one ignores the dependence of n on ξ , and assumes an arbitrary droplet size (e.g. ∼ µm ),then there will be an inconsistency between the size ofthe droplet and the resulting length scale ξ . As we willsee in Sec. IV, the typical size of the droplet L d ( ξ ) = λξ where λ ∼
50 in the droplet phase. Then combining thedefinition of the dimensionless density n = ˜ n /L d ( ξ )with that of the healing lengths ξ j = (cid:126) / (cid:112) m | µ j | , one obtains the definitions ξ CI = λa ρ a s ˜ n , (17a) ξ D = 2 λa ρ | ε dd | a s ˜ n [1 + 3 cos 2 α ] , (17b) ξ = λa ρ a s ˜ n [1 − ε dd (1 + 3 cos 2 α )] . (17c)With a flexible choice of physical parameters afforded bythe cold atom toolbox, one can show that the ratio oflength scales are a ρ /ξ CI ∼ . a ρ /ξ D ∼ . a ρ /ξ ∼ . Dy has been used. One can also considerthe the ratio of energies g n j / (cid:126) ω ρ (cid:28)
1, in whichcase one finds that g n / (cid:126) ω ρ ∼ . × − for a ρ = 10 µm, a s = 10 a and a ρ = 100 µm, a s = 10 a s re-spectively. These considerations support our claim thatthis model and the results obtained from it are capableof capturing the quasi-one-dimensional limit.Figure 2 explores the behaviour of the roton instabilityin the beyond-mean-field regime. Panel (a) shows the nu-merical solutions obtained from Eqs. (14) and (15) in the( ε dd , α ) parameter space, and we take (cid:96) = 5 × − . For γ QF (cid:54) = 0, two roton instabilities are observed in this sys-tem, due to the underlying quadratic dependence on ε dd of the dispersion relation, Eq. (11). Each shaded areacorresponds to a different atomic density, with increas-ing density causing the unstable region to shrink in area.In this way the LHY term acts to stabilize the param-eter space as the density is increased. The gray shadedarea corresponds to the (single-valued) mean-field resultin the limit γ QF = 0. The roton unstable regions alsoappear for α = π/
2, although these are smaller in areathan their α = 0 counterparts, which is attributed tothe dipoles being in a side-by-side repulsive alignment, α = α m .
25 0 . − α/π ε dd α = α m .
25 0 . α/π Bright solitonsHomogeneous Droplets(BMF) Homogeneous(MF) Homogeneous( a ) ( b ) FIG. 3. Droplet phase diagram. Panel (a) shows the stableregions in the ( ε dd , α ) parameter space in the limit γ QF = 0.Panel (b) shows the regions where droplets are expected for ξn = 10 , (cid:96) = 10 − and σ = 0 .
2. The dashed lines in bothpanels indicate the position of the magic angle α = α m . which additionally reduces the roton unstable region ofthe parameter space. It is an open question as to whethersuch a ‘reentrant’ roton behaviour could be experimen-tally observed. Example dispersion relations plotted us-ing Eq. (11) are depicted in Fig. 2 (b). Individual curvescorrespond to the cross, circle and triangle markers show-ing the dispersion slightly below (cross), above (triangle)and at (circle) the roton instability respectively. In theseexamples the density ξn = 2 × , while the dipolepolarization α = 0. The last panel (c) of Fig. 2 plotssemi-logarithmically the value of the dipolar interactionstrength ε dd against the density ξn at which the rotonmanifests for several examples of fixed (cid:96) . Each curve ter-minates at a maximum value of n , due to the additionalrepulsive LHY term overwhelming the attractive part ofthe dipolar interaction. Solutions for both α = 0 and π/ α = π/ n ; attributed again to the repulsive nature ofthis side-by-side polarization. The presence of multipleroton instabilities enriches the prospects for supersolid-ity within this system, which is a necessary condition for such a state to occur. B. Quantum Depletion
Even at zero-temperature, it is well established thatsome of the atoms will not be in the zero-momentumcondensate mode due to interactions. This is caused bythe interactions effectively mixing atoms in different fi-nite momentum states, which is known as quantum de-pletion. This effect is of particular relevance for quan-tum dipolar droplets; since it is known to qualitativelyscale with the effective diluteness parameter na s . Wecan compute expressions for this quantity by consider-ing the number operator associated with the many-bodyHamiltonian Eq. (A1)ˆ N = N + (cid:88) k (cid:54) =0 v p + (cid:88) k (cid:54) =0 ( u k + v k )ˆ a † k ˆ a k , (18)here the u k and v k represent the momentum-space formof the three-dimensional amplitude functions. Note thata fully-consistent calculation of the quantum depletion inthe quasi-one-dimensional limit would require diagonal-ization of the three-dimensional Bogoliubov-de Gennesequations to obtain the mode functions v k , however wecan still gain some intuitive insight by considering thefree-space form of the depletion. Then the quantum de-pletion can be calculated in three dimensions from thesecond term of Eq. (18) by taking the continuum limitsuch that n ex = (cid:90) d p (2 π (cid:126) ) (cid:15) p + n U k Σ (cid:113) (cid:15) p + 2 (cid:15) p n U k Σ − , (19)here n is the homogeneous density in three-dimensions,while U k Σ = g + U dd ( k ) is the Fourier transform of thetotal real-space interaction pseudo-potential, and U dd ( k )is the Fourier transform of the dipole-dipole interaction,Eq. (A2). This quantity can be evaluated exactly for thecase of α = 0 by switching to dimensionless coordinatesand working in spherical polar coordinates. This yieldsthe expression n ex n = (cid:114) n a s π (5 ± | ε dd | ) (cid:112) ± | ε dd | + (cid:115) | ε dd | (1 ∓ | ε dd | ) × sinh − (cid:18)(cid:113) | ε dd | −| ε dd | (cid:19) sgn( ε dd ) > , sin − (cid:18)(cid:113) | ε dd | | ε dd | (cid:19) sgn( ε dd ) < . (20)in evaluating Eq. (19) one has to consider the sign of ε dd as this leads to different expressions for the dipolar de-pletion n ex . The expression given by Eq. (20) scales withthe diluteness parameter, which for the parameter regimestudied in this work is of the order n a s ∼ × − , and hence will not cause a significant atom loss. It isalso worth commenting at this point as to the meaningof ψ in this context with respect to the experiments onHelium droplets. In the context of dilute atomic gases,the quantum depletion is a small effect and hence the N [ × ] − −
20 0 20 4000.511.5 x/ξ N = 0 . N = 0 . N = 2 . N = 5 − − − − − ε dd N [ × ] N [ × ] E / | µ | α =0 ε dd =2 α = π/ ε dd = − x / ξ x / ξ | ψ ( x ) | [ × ] [ × ]( a ) ( b ) ( c )( d ) ( e ) ( f )max( | ψ | ) FIG. 4. Dipolar droplet ground state computation. Panels (a) and (b) show droplet ground states for fixed atom number N = 2 × (a) and fixed dipolar strength ε dd = 2 (b), with α = 0 in both cases. Likewise panels (d) and (e) correspond to N = 2 . × and ε dd = − α = π/ assumption that the system obeys a single mode equa-tion is well supported. Experiments with Helium dropletson the other hand represent strongly correlated systems– which do not necessarily justify such an assumption.We note that the authors of [85] have also computed thedipolar depletion. C. Droplet Phases
The existence of the droplet phases in the ( ε dd , α ) pa-rameter space depends on the balance between attractiveand repulsive forces in the system. We can understandthis by examining Eq. (12), the homogeneous chemicalpotential. In the limit n a s (cid:28) µ = 0,which gives ε dd = 4 / [1 + 3 cos(2 α )]. For n a s ∼ µ [ n ] = 0 as ε ± dd = C ± (cid:112) C − C ( C − C − , (21)where we have defined C = (128 / π ) (cid:96) / √ σ √ ξn + 1and C = [1+3 cos(2 α )] /
4. Figure 3 explores the dropletsexistence in the beyond-mean-field regime. Panel (a)shows the repulsive ( µ >
0) and attractive ( µ < n a s ∼ × − , close to values obtained from theexperiment of Ref. [7] which gave an approximate bound-ary for dipolar droplet of 10 − (cid:46) n a s (cid:46) − ) which for repulsive interactions gives a homogeneous groundstate. However, when the mean-field phonon instabil-ity is crossed, the system remains homogeneous (blue re-gions). Droplets are found beyond a critical value of ε dd obtained from Eq. (21). For α < α m the droplet regionis larger than for α > α m , attributed to the head-to-tail(attractive) polarization. Increasing the atomic density n has the effect of enlarging (shrinking) these regionsfor α < α m ( α > α m ). There is no roton instability inthis analysis due to the choice of σ = 0 .
2. In realizingthe droplet states the sign of C dd and a s can’t be chosenarbitrarily, since the overall interactions should be attrac-tive. Considering Fig. 3(a-b), the top left soliton/droplet‘pockets’ of both panels (a) and (b) should have C dd > a s >
0; since ε dd >
1. Likewise for the equivalentbottom right ‘pockets’ we must instead have C dd < a s > IV. NUMERICAL SIMULATIONSA. Single Droplets
To calculate the solutions of the generalized dipolarGross-Pitaevskii equation Eq. (5) we employ a psuedo-spectral approach, the Fourier split-step method, for theresults presented in this section. Due to the size of theparameter space of the extended dipolar model, whichincludes the strength of the dipolar interactions, the po-larization angle of the dipoles, as well as the numberof atoms, it is instructive to consider the ground stateof Eq. (5) by fixing two of these physical parameterswhilst varying the other. Further, in using the healingunits (see Sec. IIIA ) the numerical value of the den-sity n appearing in the healing units is chosen such that n = max( | ψ ( x, t ) | ).Figure 4 presents examples of the droplets spatialdensity | ψ ( x ) | as a function of the dipolar interactionstrength, ε dd in Fig. 4(a) and (d) for the choices of polar-ization angle α = 0 and π/ n a s (cid:28) ε ± dd , Eq. (21). This leadsto a pair of solutions for a given set of parameters, whichcorrespond to α < α m and α > α m respectively. Thewidth of the computed ground state is observed to besensitive to the dipolar strength, ε dd . To understand theeffect of changing the number of atoms in the droplet, wesolve Eq. (5) at fixed dipolar interaction strength. Then,the dotted lines in Fig. 4(a) and (d) correspond to thevalues ε dd = 2 and ε dd = − w drop increases linearly withatom number such that w drop ∝ N . The width of thedroplet is largest for the α = π/ N . For low atom num-bers (black data) the solutions resemble a bright soliton(sech-like profile) while for increasing atom number theprofiles widen, developing the characteristic flat top asso-ciated with the droplet state (red and blue data). Panel(f) meanwhile shows the energy calculated from the def-inition E = (cid:90) d x (cid:20) ˆ p x m + g πa ρ | ψ | + Φ( x )2 | ψ | + 4 γ QF | ψ | π / a ρ (cid:21) (22)as a function of N for the data presented in panel (b)and (e). The ground state energy decreases monotoni-cally with N and dE/dN < B. Modulation Instability
The strength of the different nonlinearities directly de-termines the nature of the state that the dipolar conden-sate is in. By changing the a s scattering length from aninitial value of a is to a final value a fs such that a fs < a is , amodulational instability can be induced. The instabilityoriginates from long-wavelength perturbations that causethe break-up of a waveform into pulses; coming from thegrowth of nonlinear excitations in the system. This effecthas been used previously to investigate experimentallythe stability of individual bright solitons [96, 97] as well .
05 1 . .
15 1 . x / ξ
00 0.2511 . . a s ( t ) / a f s -100 0 100012 · -120-60060120 x / ξ t/t f x / ξ t/t f a is /a fs t/t f x/ξ | ψ ( x , t ) | m a x ( | ψ | ) ( a )( b ) ( c )( d ) , a is /a fs =1 .
11 ( e ) , a is /a fs =1 . f ) , a is /a fs =1 .
16 ( g ) , a is /a fs =1 . a is =1 . a fs t = { } tf , . , . FIG. 5. Modulational instability. Panel (a) shows a heatmapof the post-quench density | ψ ( x, t f ) | , as a function of thequench strength a is /a fs . Panel (b) shows the time-dependentscattering length a s ( t ) (see Eq. (23)), while (c) shows examplequench density data for a is /a fs = 1 .
12. The three lower panels(d-g) show example dynamics for different quenches, while thedashed lines indicate the time at which the quench occurs, t Q = 0 . t f . as dipolar droplets in a trapped three-dimensional con-text [98]. Figure 5 investigates the effect of performingan interaction quench on an initial dipolar droplet with ε dd = 2, α = 0 and N = 10 . The scattering length takesthe time-dependent form a s ( t ) = a fs + ( a is − a fs ) H ( t − t Q ) (23)where t Q defines the time at which the quench is applied,while H ( t ) is the Heaviside function. In our simulationswe take t Q = 0 . t f . Figure 5(a) shows the atomic den-sity | ψ ( x, t f ) | of the dipolar gas as a function of thequench strength a is /a fs taken at the final point of the nu-merical integration t = t f . For modest quench strengths( a is /a fs (cid:46) .
1) the droplets shape remains intact, withthe exception of the excitation of surface waves. Above acritical quench strength the droplet undergoes the modu-lation instability, and in general breaks apart into smallerdroplets and dipolar bright solitons, as well as the emis-sion of low density radiation. Panel (b) and (c) show re-spectively the quench protocol of Eq. (23) and exampledynamics for a is /a fs = 1 .
12. The lower row of panels (d-g)of Fig. 5 show individual examples of the quench dynam-ics for increasing a is /a fs . Here we observe both even andodd numbers of droplets, while panel (f) shows a situa-tion where two droplets and two bright solitons are cre-ated after the quench. For larger quench strengths, suchas that presented in panel (g), a central droplet is pro-duced along with increasing amounts of radiation, in thisexample short-lived bright soliton bound states are alsoobserved. In general we find that the onset of the modu-lational instability (and the final state after the quench)depend strongly on the atom number. If for examplethe initial number of atoms in the droplet is reduced,then the modulational instability occurs at larger valuesof a is /a fs compared to a larger initial atom number. Thegeneration of multiple droplets as presented here usingthe modulational instability relies on being able to tunethe scattering length a s of the condensate which can inturn lead to significant atomic losses. To address this,there are proposals to produce condensates with attrac-tive interactions with reduced noise and greater controlover the final experimental state of the system [99, 100]. C. Collisional Population Transfer and DropletFission
The coherent nature of the superfluid state provides aconvenient tool to explore quantum mechanical phenom-ena at macroscopic length scales. One striking manifesta-tion of matter-waves coherence is the so-called Josephsoneffect. Here, two superconductors or superfluids whichare separated by an insulting barrier can experience acurrent, originating from atomic tunneling between thetwo superconductors/superfluids. The dipolar dropletsrepresent an interesting addition to the superfluid fam-ily, since they are effectively an isolated (finite) regionof homogeneous fluid, so understanding their binary dy-namics is expected to yield novel phenomena. To inves-tigate the basic physics of binary droplet dynamics, weperform simulations with an initial state of the form ψ ( x ) = (cid:88) n = ± ψ ( x − x n ) e imv n x/ (cid:126) + iδ n . (24)This constitutes a symmetric state comprising twodroplets whose centres are initially separated by a dis-tance x + − x − = 60 ξ , which are traveling towards eachother at constant velocity v ± = ± v . The initial phasedifference δ = δ + − δ − between the two droplets is δ ∈ [0 , π ]. In Figure 6 we present results of dropletcollisions using the initial state defined by Eq. (24). Wetake for the physical parameters ε dd = 2, N = 3 × , (cid:96) = 10 − , σ = 0 .
2, and the total length of each numer-ical integration is t f = 140 (cid:126) / | µ | . Figures 6 (a-f) show -60-3003060 x / ξ x / ξ t/t f x / ξ t/t f δ ( t ) /π ∆ N ( t f ) / N GPE[ v = 0 . h/mξ ] GPE[ v = 0 . h/mξ ]GPE[ v = 0 . h/mξ ] ∆ N ( t = t ) /N ∆ N ( t = t ) /N P j N j ( a ) , δ =0 ( b ) , δ = π/ c ) , δ = π/ d ) , δ =3 π/ e ) , δ = π ( f ) , δ =3 π/ g ) max( | ψ | ) FIG. 6. Population transfer. Panels (a-f) show space-timedynamics for different initial phase differences, δ . The totalintegration time is t f = 140 (cid:126) / | µ | . Panel (g) shows a com-parison of the final droplet populations at t = t f calculatedusing Eq. (25). space-time plots for different initial phase differences, δ and initial velocity mξv / (cid:126) = 0 .
25. We observe thatpost collision two droplets emerge - with different pop-ulations (and sizes) that depend on the choice of initialphase difference. For example, choosing δ = π (panel6 (e)) produces two droplets with an equal number ofatoms present in each droplet post collision. The pres-ence of excitations in the form of sound waves can beseen here post collision, reflecting back and forth insideeach droplet (viz. Eq. (13)), which is indicative of non-integrable dynamics. To quantify this change in popula-tion, we can calculate the population difference betweenthe droplets at t = t f as a function of the initial phasedifference δ ( t ). The population difference is defined as0∆ N ( t f ) = N ( t f ) − N ( t f ) where N j ( t ) = (cid:90) x j + x j − d x | ψ j ( x, t ) | . (25)Here each integral computes the number of atoms inthe individual droplets between the edges of the dropletgiven by x = x j ± . Panel 6(g) shows some example re-sults with different choices of initial velocity, mξv / (cid:126) =0 . , . , .
35 (triangle, circle and square markers).Only droplet collisions for mξv / (cid:126) = 0 .
25 represents asituation where two droplets emerge post collision for thefull range of phase differences between δ = 0 and δ = π ,due to the existence of bound states (droplet molecules)and droplet fission at smaller and larger initial velocitiesrespectively, breaking droplet number conservation. Assuch, a parameter window exists where one can comparethe collisional transfer to the Josephson effect. Then,the semi-classical Josephson equations for a superfluidare written as [102] ddt ∆ N ( t ) = J (cid:126) (cid:112) N ( t ) N ( t ) sin( δ ( t )) , (26a) ddt δ ( t ) = J (cid:126) (cid:20) N ( t ) N ( t ) − N ( t ) N ( t ) (cid:21) cos( δ ( t )) . (26b)Here, the parameter J describes the strength of the cou-pling between the two superfluids. We will use the fi-nal integration time t = t f to fit our model to the nu-merical simulations of the extended GPE (green circles),since the total integration time determines the nature ofthe observed Josephson oscillations. Panel 6 (g) showsa comparison between the populations of the dropletsat t = t f , calculated from the numerical solutions tothe Josephson equations. The blue and black dashedlines are obtained from the pair of Josephson relationsEq. (26), for t f = 0 . (cid:126) /J and t f = 1 . (cid:126) /J respectively.For small t f , (blue dashed) the oscillation is linear, anddoes not follow the extended GPE data (green circles).For larger t f , (black dashed) the oscillation exhibits astronger nonlinear character, and follows the extendedGPE data quite well. It would be interesting to studythis effect in more detail, to understand if this analogywith the Josephson effect can be extended, for exam-ple into the highly nonlinear regime. The related ex-periment of Ref. [51] studied the collisions of (trapped)matter-wave solitons in an attractive condensate of Liatoms. They investigated how the relative phase ∆ φ be-tween the solitons affected the collision dynamics of thesolitons, observing constructive and destructive matter-wave interference at the point of collision for in-phase(∆ φ = 0) and out-of-phase (∆ φ = π ) respectively. Theyalso observed the in-phase collapse of the wave functionabove a critical atom number. Our work reveals similarphase-sensitive dynamics, as well as population transferbetween the droplets post collision. Unlike the solitons[52] the droplets cannot undergo collapse, which is sup-pressed by the repulsive LHY contribution. Hence we areable to probe the dynamics for all relative phases. ( a ) | µ | t/ ¯ h x / ξ
00 35 70 105 140-60-3003060 x / ξ | µ | t/ ¯ h x / ξ | µ | t/ ¯ h ( f )0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.912345 mξv / ¯ h N o . D r o p l e t s δ = πδ = 0 δ = 0 δ = π ( b ) , mξv / ¯ h =0 . c ) , mξv / ¯ h =0 . d ) , mξv / ¯ h =0 . e ) , mξv / ¯ h =1 m a x ( | ψ | ) FIG. 7. Droplet collisions. A long-lived droplet dimer isshown in panel (a), while panels (b-e) show in and out-of-phase dynamics for different initial velocities. (f) computesthe number of droplets as a function of the initial velocity, v for in-phase δ = 0 collisions. To understand the effect of the droplets initial velocityon the dynamics, Figure 7 shows simulations of dropletcollisions with fixed initial phase difference, for δ = 0 , π .Panel (a) shows a long-lived droplet dimer formed fromtwo initially out-of-phase stationary droplets. Then, ex-ample dynamics are shown in panels (b-e) for differentfinite initial velocities (see individual captions). The in-phase collisions demonstrate that the droplets undergofission [103–105] with multiple droplet states producedpost-collision, shown in panels (b) and (d). For the out-of-phase collisions ( δ = π ) two out-going droplets areobserved, for low incoming velocity (panel (c)) wheresmall amounts of sound are produced post collision. Fora greater initial velocity (panel (e)), larger amounts ofsound and radiation are produced, still with two out-going droplets. Finally panel (f) computes the numberof droplets post collision for in and out-of-phase colli-sions, as a function of the initial velocity. For even larger1velocities ( mξv / (cid:126) ∼
1) the delicate balance of kineticand potential energy that maintains the droplet is vio-lated, and the droplet breaks apart into atomic radiation.The splitting of a single droplet into multiple smallerdroplets is also noted to occur with classical fluids. Ex-perimental investigations have also explored the collisionsof droplets in a binary system [101]. It is interesting tocompare the results of this experiment with our findings.Ferioli et al. found that the outcome of the droplets col-lisions depends critically on the relative velocity of thedroplet pair. For smaller incident velocities the dropletscould form a bound state, while at larger velocities thedroplets remain separated post collision. Our results con-cerning dipolar droplets reveal that the droplets can in-stead split into increasing numbers of smaller droplets asthe speed of collision increases. It should be noted thatas well as the beyond-mean-field LHY term the authorsof Ref. [101] also model three body effects, which we havenot considered in our work, and would be a natural andinteresting extension for a future study.
V. CONCLUSIONS AND OUTLOOK
In this work we have investigated the properties of aquasi-one-dimensional dipolar Bose-Einstein condensatein the presence of quantum fluctuations. By calculatingthe excitations of this system within the framework of theBogoliubov-de Gennes formalism, we identified regimesunstable to the roton instability. Interestingly, the rotonunstable regions are found in general to appear in pairsfor a given dipole polarization angle; this is due to theunderlying quadratic dependence on the dipolar strengthfrom the LHY contribution to the excitation energy.We examined the nature of the ground states of thissystem, observing the appearance of droplet phases in re-gions of the parameter space where the total interactionsare net attractive. By applying an interaction quenchto a single large droplet, the nature of the modulationinstability was investigated. It was found that for mod-erate quenches, both even and odd numbers of dropletscan be generated. For larger quenches, bright solitonsand increasing amounts of atomic 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 thedroplets, atomic population transfer was observed be-tween the droplets. This was interpreted and comparedwith the superfluid Josephson effect, finding good agree-ment in a window of the parameter space. Droplet fissionwas also observed as the initial velocity of the droplet wasincreased.It would be interesting in the future to understandhow a harmonic trap changes the physics of the one-dimensional droplet, and in particular how the interplayof quenching both the interactions and trapping strengthchanges the number of droplets produced. One couldalso use this model to understand supersolid phases in the one-dimensional context, as well as studying dropletsand their dynamics with models that incorporate higher-dimensional effects [106, 107]. Finally, it would alsobe beneficial to understand the lifetime of the one-dimensional droplet, which can be computed from three-body atomic recombinations [108].
VI. ACKNOWLEDGEMENTS
We thank Jing Li for discussions and Thomas Flynn forcomments on the manuscript. This work was supportedby the Ministry of Education, Culture, Sports, Science(MEXT)-Supported Program for the Strategic ResearchFoundation at Private Universities “Topological Science(Grant No. S1511006). NGP acknowledges support fromthe Engineering and Physical Sciences Research Council(Grant No. EP/M005127/1) for support. TB acknowl-edges support from the Engineering and Physical Sci-ences Research Council (Grant No. EP/R51309X/1) forsupport.
Appendix A: Quantum Fluctuations
In this appendix we compute the form of the dipolarLHY term [85, 86] including the effect of the polarizationangle of the dipoles. The LHY correction arises from cor-recting the otherwise divergent ground state energy of ahomogeneous gas of bosons. The many-body Hamilto-nian ˆ H for a gas of N interacting dipoles of volume V can be written as [109]ˆ H = n N U k Σ + (cid:88) k (cid:54) = (cid:113) ε ( ε + 2 n U k Σ )ˆ a † k ˆ a k − (cid:88) k (cid:54) = (cid:20) n U k Σ + ε − (cid:113) ε ( ε +2 n U k Σ ) (cid:21) , (A1)here ˆ a k (ˆ a † k ) defines the annihilation (creation) operatorfor a quasi-particle with momentum k and energy ε = (cid:126) k / m , and U k Σ = g + U dd ( k ) is used as short-handfor the total Fourier transform of the interaction pseudo-potential, which is defined as [87] U dd ( k ) = C dd (cid:20) k x sin α + k z cos α ) k x + k y + k z − (cid:21) (A2)where k i defines the cartesian components of the mo-mentum, and α is the dipole polarization angle in the x - z plane (see Fig. 1 (a)). The ground state energy is de-fined in the usual quantum mechanical way as E QF = (cid:104) Ψ k (cid:54) =0 | ˆ H | Ψ k (cid:54) =0 (cid:105) , and since ˆ a † k (cid:54) = ˆ a k (cid:54) = | Ψ k (cid:54) =0 (cid:105) = 0, thesummation term on the first line of Eq. (A1) does notcontribute to E QF . Converting the remaining summa-tion on the second line of Eq. (A1) to an integral using2 ( a )02004006008001 ,
000 Re( I / ( ε dd ,α =0))1 + ε ε + ε ( b ) − −
10 0 10 2002004006008001 , ε ddRe( I / ( ε dd ,α =0))Im( I / ( ε dd ,α =0)) FIG. 8. Comparison of approximations. In (a) the blue circlesshow the real part of the exact (numerical) form of Eq. (A7)for n = 5 /
2, while the green triangles and red crosses show thequadratic and cubic approximations respectively of Eq. (A10).(b) shows the real and imaginary parts of Eq. (A7). (cid:88) k (cid:54) =0 → V (cid:90) d k (2 π ) (A3)we find that this term diverges as k → ∞ . Since theground state energy E QF should be a real finite quan-tity, we can renormalize this term to give a finite result.Then, expanding this term in inverse power of the kineticenergy, one can show that n U kΣ + ε − (cid:113) ε ( ε +2 n U kΣ ) (cid:39) n U kΣ2 ε k . (A4)Then, a finite value for the ground state energy can beobtained by subtracting this final term in Eq. (A4) fromthe summation on the second line of Eq. (A1). Hence wewrite E QF V = − (cid:90) d k (2 π ) (cid:20) n U k Σ − (cid:113) ε ( ε +2 n U kΣ ) − n U kΣ2 ε k (cid:21) , (A5) Although the integral defined by Eq. (A5) is convergent,it does not in general exist in closed form. After inte-grating out the radial momentum from Eq. (A5), onecan show that this expression can be re-written as E QF V = 8 √ n g (2 π ) (cid:18) mn g (cid:126) (cid:19) / I / ( ε dd , α ) (A6)where I n ( ε dd , α ) = (cid:90) dΩ4 π (cid:26) ε dd [3 L ( θ, φ, α ) − (cid:27) n . (A7)Here the short-hand L ( θ, φ, α ) = sin θ cos φ sin α +cos θ cos α has been used. Although the integral definedby Eq. (A7) does not exist in closed form, on physicalgrounds we are interested in situations where the param-eter ε dd (cid:38) ε dd , which will allow us to identify anappropriate analytical approximation to Eq. (A6). Thentaking n = 5 /
2, Eq. (A7) becomes I / ( ε dd , α ) = (cid:90) dΩ4 π (cid:26)
1+ 5 ε dd L ( θ, φ, α ) − ε L ( θ, φ, α ) − + 15 ε
72 [3 L ( θ, φ, α ) − + . . . (cid:27) . (A8)Evaluation of the integrals appearing in Eq.(A8) can beaccomplished using the result (cid:90) dΩ4 π L ( θ, φ, α ) n = 12 n + 1 , n ∈ Z ≥ , (A9)then after collating the results of Eqs. (A6), (A8) and(A9) we arrive at E QF V = 6415 g ( n ) (cid:114) n a s π (cid:18) ε + 211 ε (cid:19) . (A10)In Fig. 8(a) we compare the quadratic and cubic formsof Eq. (A10) with the exact numerically calculated valueof Eq. (A6). Figure 8(a) shows three sets of data. Theblue circles show the exact value of Eq. (A7) computednumerically, while terms up-to quadratic (green trian-gles) and cubic (red crosses) orders are plotted usingEq. (A10). For ε dd >
0, it is clear that the cubic form ofEq. (A10) gives closer agreement with the exact form of I / ( ε dd , α ). However for ε dd < ∼ ε dd , but gives closer agreement for negative ε dd , compared to the cubic form. As such, we adopta quadratic approximation, which overall allows us tostudy the role of the quantum fluctuations from a numeri-cal and analytical viewpoint, hence the quadratic natureof Eq. (2) of the main text. Panel (b) shows a com-parison of the real (black-dashed) and imaginary (blue-dashed) parts of I / . Note for − . ≤ ε dd ≤
1, one has3Im( I / ) = 0. For most of our results we only considervalues of | ε dd | (cid:46) | ε dd |∼
20, which accordingto Fig. 8(b) would be problematic since Im( I / ) (cid:29) σ (cid:46) a s ) <
0. Inthese regions the rotons are predicted closer to ε dd = 0,exactly where the full calculation of I / is purely real.This offers a possible opportunity for future studies ofthis interesting effect. Appendix B: Bogoliubov-de Gennes Equations
In this appendix we give further details of the deriva-tion of the beyond-mean-field Bogoliubov-de Gennesequations of the text, Eqs. (7). The derivation with themethodology we employed to solve Eq. (10) is qualita-tively the same for either mode function, here we focuson the u ( x ) mode function without loss of generality.To proceed, we note that the time-independent beyond-mean-field dipolar Gross-Pitaevskii equation [ H GP1D − µ ] φ j = (cid:15) GP j φ j (obtained from Eq. (5) by the substitution ψ ( x, t ) = φ j ( x ) exp( − i [ (cid:15) GP j + µ ] t/ (cid:126) )) possesses a spectralbasis with orthonormal modes φ j ( x ) and correspondingenergies (cid:15) GP j . Since the condensate mode has been re- moved, the resulting quasiparticle modes are automati-cally orthogonal to the condensate. Then by making theexpansion u ( x ) = (cid:88) λ c λ φ λ ( x ) (B1)and inserting Eq. (B1) into (10) whilst premultiplying by φ ∗ γ ( x ) and using the orthonormal property of the spectralbasis states (cid:90) d xφ ∗ γ ( x ) φ λ ( x ) = δ γλ , (B2)we obtain the matrix-valued equation for the ω eigenval-ues (cid:88) λ (cid:20) δ γλ (cid:15) GP λ + 2 M γλ (cid:21) (cid:15) GP λ c λ = ( (cid:126) ω ) c γ , (B3)with the exchange matrix elements M γλ = (cid:104) φ γ |M| φ λ (cid:105) defined as M γλ = (cid:90) d xφ ∗ γ ( x ) (cid:20) g πa ρ | ψ | + 3 γ QF π / a ρ | ψ | (cid:21) φ λ ( x )+ 12 π (cid:90) d kϕ ∗ γ ( − k ) U ( k ) ϕ λ ( k ) , (B4)with ϕ ν ( k ) = (cid:82) d xφ ν ( x ) ψ ( x ) exp( ikx ) [110]. [1] A. Griesmaier, J. Werner, S. Hensler, J. Stuhler, and T.Pfau, Phys. Rev. Lett. 94, 160401 (2005).[2] Q. Beaufils, R. Chicireanu, T. Zanon, B. Laburthe-Tolra, E. Mar´echal, L. Vernac, J. C. Keller, and O.Gorceix, Phys. Rev. A 77, 061601(R) (2008).[3] M. Lu, N. Q. Burdick, S. H. Youn, and B. L. Lev, Phys.Rev. Lett. 107, 190401 (2011).[4] Y. Tang, N. Q. Burdick, K. Baumann, and B. L. Lev,New J. Phys. 17, 045006 (2015).[5] K. Aikawa, A. Frisch, M. Mark, S. Baier, A. Rietzler, R.Grimm, and F. Ferlaino, Phys. Rev. Lett. 108, 210401(2012).[6] H. Kadau,M. Schmitt, M.Wenzel, C.Wink, T. Maier, I.Ferrier-Barbut, and T. Pfau, Nature 530, 194 (2016).[7] I. Ferrier-Barbut, H. Kadau, M. Schmitt, M.Wenzel,and T. Pfau, Phys. Rev. Lett. 116, 215301 (2016).[8] L. Chomaz, S. Baier, D. Petter, M. J. Mark, F.W¨achtler, L. Santos, and F. Ferlaino, Phys. Rev. X 6,041039 (2016).[9] D. Baillie, R. M. Wilson, R. N. Bisset, P. B. Blakie,Phys. Rev. A 94, 021602(R) (2016).[10] R. N. Bisset, R. M. Wilson, D. Baillie, P. B. Blakie,Phys. Rev. A 94, 033619 (2016).[11] D. Baillie and P. B. Blakie, Phys. Rev. Lett. 121, 195301(2018).[12] D. Baillie, R. M. Wilson, and P. B. Blakie, Phys. Rev.Lett. 119, 255302 (2017). [13] F. W¨achtler and L. Santos, Phys. Rev. A 93, 061603(R)(2016).[14] D. H. J. O’Dell, S. Giovanazzi, and G. Kurizki, Phys.Rev. Lett. 90, 110402 (2003).[15] L. Santos, G. V. Shlyapnikov, and M. Lewenstein, Phys.Rev. Lett. 90, 250403 (2003).[16] L. Chomaz, R. M. W. van Bijnen, D. Petter, G. Faraoni,S. Baier, J. H. Becher, M. J. Mark, F. W¨acher, L. San-tos, and F. Ferlaino, Nat. Phys. 14, 442 (2018).[17] D. Petter, G. Natale, R. M. W. van Bijnen, A. Patschei-der, M. J. Mark, L. Chomaz and F. Ferlaino, Phys. Rev.Lett. 122, 183401 (2019).[18] A. Bulgac, Phys. Rev. Lett. 89, 050402 (2002).[19] D. S. Petrov and G. E. Astrakharchik, Phys. Rev. Lett.117, 100401 (2016).[20] D. Edler, C. Mishra, F. W¨achtler, R. Nath, S. Sinha,and L. Santos, Phys. Rev. Lett. 119, 050403 (2017).[21] P. Zin, M. Pylak, T. Wasak, M. Gajda, and Z. Idziaszek,Phys. Rev. A 98, 051603(R) 2018.[22] R. O(cid:32)ldziejewski, W. G´orecki, K. Paw(cid:32)lowski, and K.Rz¸a˙zewski, Phys. Rev. Lett. 124, 090401 (2020).[23] G. Semeghini, G. Ferioli, L. Masi, C. Mazzinghi, L. Wol-swijk, F. Minardi, M. Modugno, G. Modugno, M. In-guscio, and M. Fattori, Phys. Rev. Lett. 120, 235301(2018).[24] C. R. Cabrera, L. Tanzi, J. Sanz, B. Naylor, P. Thomas,P. Cheiney, and L. Tarruell, Science 359, 301 (2018). [25] K. E. Wilson, N. Westerberg, M. Valiente, C. W. Dun-can, E. M. Wright, P. ¨Ohberg, and D. Faccio, Phys.Rev. Lett. 121, 133903 (2018).[26] A. J. Leggett. Quantum Liquids.
Oxford UniversityPress, 2006.[27] T. Schneider, C.P. Enz, Phys. Rev. Lett. 27, 1186(1971).[28] L.P. Pitaevskii, Zh. Eksp. Teor. Fiz. 39, 423 (1984).[29] S. Giovanazzi and D. H. J. O’Dell, Eur. Phys. J. D 31,439 (2004).[30] M. Jona-Lasinio, K. (cid:32)Lakomy, and L. Santos, Phys. Rev.A 88, 013619 (2013).[31] R. N. Bisset, D. Baillie, and P. B. Blakie, Phys. Rev. A88, 043606 (2013).[32] E. P. Gross, Ann. Phys. 4, 57 (1958).[33] A. F. Andreev and I. M. Lifshitz, Sov. Phys. JETP 29,1107 (1969)[34] G. V. Chester, Phys. Rev. A 2, 256 (1970).[35] A. J. Leggett, Phys. Rev. Lett. 25, 1543 (1970).[36] A. J. Leggett, J. Stat. Phys. 93, 927 (1998).[37] F. B¨ottcher, J.-N. Schmidt, M. Wenzel, J. Hertkorn, M.Guo, T. Langen, and T. Pfau, Phys. Rev. X 9, 011051(2019).[38] L. Tanzi, E. Lucioni, F. Fam´a, J. Catani, A. Fioretti, C.Gabbanini, R. N. Bisset, L. Santos, and G. Modugno,Phys. Rev. Lett. 122, 130405 (2019).[39] T. Bland, M. J. Edmonds, N. P. Proukakis, A. M. Mar-tin, D. H. J. O’Dell, and N. G. Parker, Phys. Rev. A 92,063601 (2015).[40] K. Paw(cid:32)lowski and K. Rz¸a˙zewski, New J. Phys. 17,105006 (2015).[41] M. J. Edmonds, T. Bland, D. H. J. O’Dell, and N. G.Parker, Phys. Rev. A 93, 063617 (2016).[42] T. Bland, K. Paw(cid:32)lowski, M. J. Edmonds, K. Rz¸a˙zewski,and N. G. Parker, Phys. Rev. A 95, 063622 (2017).[43] B. B. Baizakov, S. M. Al-Marzoug, and H. Bahlouli,Phys. Rev. A 92, 033605 (2015).[44] M. J. Edmonds, T. Bland, R. Doran, and N. G. Parker,New J. Phys. 19, 023019 (2017).[45] P. Pedri and L. Santos, Phys. Rev. Lett. 95, 200404(2005).[46] I. Tikhonenkov, B. A. Malomed, and A. Vardi, Phys.Rev. Lett. 100, 090406 (2008).[47] M. Raghunandan, C. Mishra, K. (cid:32)Lakomy, P. Pedri, L.Santos, and R. Nath, Phys. Rev. A 92, 013637 (2015).[48] R. Sachdeva, M. N. Tengstrand, S. M. Reimann, Phys.Rev. A 102, 043304 (2020).[49] C. F. Barenghi and N. G. Parker,
A Primer on QuantumFluids (Springer, Berlin, 2016).[50] C. J. Pethick and H. Smith,
Bose-Einstein Conden-sation in Dilute Gases (Cambridge University Press,Cambridge, 2002)[51] J. H. V. Nguyen, P. Dyke, D. Luo, B. A. Malomed, andRandall G. Hulet, Nat. Phys. 10, 918 (2014).[52] N. G. Parker, A. M. Martin, S. L. Cornish, and C. S.Adams, J. Phys. B: At. Mol. Opt. Phys. 41, 045303(2008).[53] F. Dalfovo, S. Giorgini, L. P. Pitaevskii, and S.Stringari, Rev. Mod. Phys. 71, 463 (1999).[54] S. Ronen, D. C. E. Bortolotti, and J. L. Bohn, Phys.Rev. A 74, 013623 (2006).[55] M. Tylutki, G. E. Astrakharchik, B. A. Malomed, andD. S. Petrov, Phys. Rev. A 101, 051601(R) (2020). [56] H. Hu and X.-Ji Liu, Phys. Rev. A 102, 053303 (2020).[57] D. Baillie and P. B. Blakie, New J. Phys. 17 033028(2015).[58] J. A. M. Huhtam¨aki and P. Kuopanportti, Phys. Rev.A 84, 043638 (2011).[59] S. Ronen, D. C. E. Bortolotti and J. L. Bohn, Phys.Rev. Lett. 98, 030406 (2007).[60] R. M. Wilson and J. L. Bohn, Phys. Rev. A 83, 023623(2011).[61] R. N. Bisset, D. Baillie and P. B. Blakie, Phys. Rev. A86, 033609 (2012).[62] A. D. Martin and P. B. Blakie, Phys. Rev. A 86, 053623(2012).[63] M. Kreibich, J. Main and G. Wunner, J. Phys. B: At.Mol. Opt. Phys. 46 045302 (2013).[64] P. B. Blakie, D. Baillie, and R. N. Bisset, Phys. Rev. A88, 013638 (2013).[65] J. L. Bohn, A. M. Rey, and J. Ye, Science 357, 1002(2017).[66] P. K. Molony, P. D. Gregory, Z. Ji, B. Lu, M. P.K¨oppinger, C. Ruth Le Sueur, C. L. Blackley, J. M.Hutson, and S. L. Cornish, Phys. Rev. Lett. 113, 255301(2014).[67] J. Woo Park, S. A. Will, and M. W. Zwierlein, Phys.Rev. Lett. 114, 205302 (2015).[68] M. Guo, B. Zhu, B. Lu, X. Ye, F. Wang, R. Vexiau,N. B.-Maafa, G. Qu´em´ener, O. Dulieu, and D. Wang,Phys. Rev. Lett. 116, 205303 (2016).[69] K.-K. Ni, S. Ospelkaus, M. H. G. de Miranda, A. Pe‘er,B. Neyenhuis, J. J. Zirbel, S. Kotochigova, P. S. Juli-enne, D. S. Jin, and J. Ye, Science 322, 231 (2008).[70] S. Ospelkaus, K.-K. Ni, D. Wang, M. H. G. de Miranda,B. Neyenhuis,1 G. Qu´em´ener, P. S. Julienne, J. L. Bohn,D. S. Jin, and J. Ye, Science 527, 853 (2010).[71] K.-K. Ni, S. Ospelkaus, D. Wang, G. Qu´em´ener, B.Neyenhuis, M. H. G. de Miranda, J. L. Bohn, J. Ye,and D. S. Jin, Nature 464, 1324 (2010).[72] J. G. Danzl, M. J. Mark, E. Haller, M. Gustavsson, R.Hart, J. Aldegunde, J. M. Hutson, and H.-C. N¨agerl,Nat. Phys. 6, 265 (2010).[73] M. H. G. de Miranda, A. Chotia, B. Neyenhuis, D.Wang, G. Qu´em´ener, S. Ospelkaus, J. L. Bohn, J. Ye,and D. S. Jin, Nat. Phys. 7 502, (2011).[74] B. Yan, S. A. Moses, B. Gadway, J. P. Covey, K. R. A.Hazzard, A. M. Rey, D. S. Jin, and J. Ye, Nature 501,521 (2013).[75] G. E. Astrakharchik and B. A. Malomed, Phys. Rev. A98, 013631 (2018).[76] T. Mithun, A. Maluckov, K. Kasamatsu, B. A. Mal-omed, and A. Khare, Symmetry 12, 174 (2020).[77] E. Chiquillo, Phys. Rev. A 99, 051601(R) (2019).[78] A. Tononi, Y. Wang, and L. Salasnich, Phys. Rev. A 99,063618 (2019).[79] P. Cheiney, C. R. Cabrera, J. Sanz, B. Naylor, L. Tanzi,and L. Tarruell, Phys. Rev. Lett. 120, 135301 (2018).[80] S. Giovanazzi, A. Gorlitz, and T. Pfau, Phys. Rev. Lett.89, 130401 (2002).[81] Y. Tang, W. Kao, K.-Y. Li, and B. L. Lev, Phys. Rev.Lett. 120, 230401 (2018).[82] D. Baillie and P. B. Blakie, Phys. Rev. A 101, 043606(2020).[83] S. B. Prasad, T. Bland, B. C. Mulkerin, N. G. Parker,and A. M. Martin, Phys. Rev. Lett. 122, 050401 (2019).5