Supersolidity in an elongated dipolar condensate

We present a theory for the emergence of a supersolid state in a cigar-shaped dipolar quantum Bose gas. Our approach is based on a reduced three-dimensional (3D) theory, where the condensate wavefunction is decomposed into an axial field and a transverse part described variationally. This provides an accurate fully 3D description that is specific to the regime of current experiments and efficient to compute. We apply this theory to understand the phase diagram for a gas in an infinite tube potential. We find that the supersolid transition has continuous and discontinuous regions as the averaged density varies. We develop two simplified analytic models to characterize the phase diagram and elucidate the roles of quantum droplets and of the roton excitation.

Calculations of the ground states of this system using the extended Gross-Pitaevskii equation (eGPE) have shown good quantitative agreement with the observations of the experiments. The eGPE theory differs from the usual Gross-Pitaevskii equation for dipolar condensates by including the leading order effect of quantum fluctuations [26][27][28][29]. Currently little is known about the phase diagram, such as the nature of the transition between phases and how this depends on the confining potential or atomic density. Previous general studies of supersolidity have predicted that the transition is continuous in the one-dimensional (1D) case, while higher dimensional cases are generally discontinuous [30] (also see [13]). However, while the supersolid realized in dipolar gas experiments exhibits crystalline order (density modulation) along one spatial dimension, the transverse degrees of freedom are not frozen out. This intrinsic 3D character is not captured in earlier theories of 1D supersolids, and developing an appropriate theory for this regime is the focus of this paper.
Here we develop formalism for a dipolar gas in an infinite tube, i.e. with transverse harmonic confinement but free in the z-direction. We develop a simplified 3D theory in which the transverse wavefunction is described by two variational parameters and the axial field is treated numerically (cf. the planar theory of Ref. [19]). Our main results are the phase diagrams in Fig. 1(b) and (c) revealing the contrast of the density modulations (i.e. crystalline order) and the persistence of superfluidity, respectively. These results show that the con- Contrast of density modulation and (c) superfluid fraction as a function of the average linear density and the s-wave scattering length. The s-wave scattering length for the transition to the modulated state (a * , white line) and for the roton instability of the uniform BEC (a * rot , red line) are shown. The dashed line (where a * = a * rot ) between the circle markers at the densities nlow = 1.3 × 10 3 /µm and nhigh = 4.6 × 10 3 /µm indicates a continuous transition. Results for 164 Dy using a dd = 130.8a0, with ωx,y = 2π × 150Hz.
densate to supersolid transition is continuous for a range of intermediate densities, but is otherwise discontinuous. Reduced 3D theory for the elongated system -We consider a zero temperature dipolar Bose gas described by the field Ψ(x). The transverse confinement is harmonic with angu-arXiv:2004.12577v1 [cond-mat.quant-gas] 27 Apr 2020 lar frequencies (ω x , ω y ) and the atomic magnetic dipoles are aligned along the y-direction by an external magnetic field [see Fig. 1(a)]. Following [31] we decompose the field as Ψ(x) = ψ(z)χ(x, y) where ψ(z) describes the axial field and χ is a variational treatment of the transverse directions. We take χ(x, y) = 1 √ πl e −(ηx 2 +y 2 /η)/2l 2 to be a Gaussian function with variational parameters {l, η} describing its mean width and anisotropy. Our interest is in stationary solutions of specified average linear density n along z.
The uniform ground state is the form Ψ BEC = √ nχ, which we refer to as the BEC state. Here the energy per particle, computed from the eGPE energy functional, is given by the nonlinear function where x η + ω 2 y η) is the singleparticle energy of the transverse degrees of freedom. Note that U and g QF also depend on l and η [26][27][28][29]32]. The two-body interaction for the axial field in k z -space is with Ei being the exponential integral, and Q ≡ 1 2 √ ηk 2 z l 2 [31]. In Eq. (2) a s is the s-wave scattering length and a dd ≡ mµ 0 µ 2 m /12π 2 is the dipole length. The quantum fluctuations are described in Eq. (1) by the higher order nonlinearity with coefficient g QF = 256 2 15ml 3 a s a 3 s (1 + 3 2 2 dd ), where dd = a dd /a s . In Ref. [31] the accuracy of this variational theory has been established with detailed comparisons to full numerical solutions of the 3D eGPE.
Of most interest here is when the ground state of the system spontaneously breaks the translational symmetry along z and develops crystalline order, i.e. |ψ| 2 is periodically modulated in space. In this case we can define ψ on a unit cell of length L, i.e. uc = {− 1 2 L ≤ z < 1 2 L}, subject to periodic boundary conditions, and the normalization constraint uc dz |ψ| 2 = nL. The energy per particle is given by where the two-body interactions are described by the effective potential To obtain stationary solutions we vary {l, η, L} and ψ(z) to find local minima of (3). Phase diagram -In Fig. 1(b) and (c) we present phase diagrams for the development of crystalline and superfluid order found by solving the reduced 3D eGPE for the ground state by minimising Eq. (3) as a function of n and a s . Our results are for Dy atoms in a radially isotropic tube potential with ω x,y = 2π×150Hz, similar to the transverse trapping used in experiments [20][21][22][23][24][25].
We characterize the crystalline order by the linear density contrast C = (|Ψ| 2 max − |Ψ| 2 min |)/(|Ψ| 2 max + |Ψ| 2 min ), where |Ψ| max and |Ψ| min are the maximum and minimum of |Ψ| on the z axis, respectively. When C = 0 the ground state is uniform and is identical to Ψ BEC . The results in Fig. 1(b) show that density-modulated ground states occur below a certain scattering length a * that depends on n. For an intermediate range of densities n ∈ [n low , n high ] the transition between uniform and modulated states is continuous, i.e. the contrast develops continuously and the ground state wavefunction evolves smoothly as a s changes [also see Fig. 2 Outside of this density range the transition is discontinuous, i.e. with a discontinuity in contrast and a sudden change in the ground state wavefunction (also see Figs. 3 and 4).
We also quantify the superfluid order in the system using the bound proposed by Leggett [4,30] For states that modulate along one spatial dimension this measure is equivalent to the superfluid fraction based on the nonclassical translation inertia [18,33]. The results in Fig. 1(c) show that the BEC has f s = 1. This reduces when the density modulates, but remains appreciable close to the transition line a * and for n 1 × 10 3 /µm. We regard the modulated ground state with an appreciable superfluid fraction as being in the supersolid phase. As the scattering length is reduced the superfluid fraction rapidly reduces and the system becomes an insulating droplet crystal without any appreciable superfluid transport. In this regime quantum and thermal fluctuations (not included in the current theory) will be important and may cause the crystal to insulate at higher a s values.
Continuous transition region -Weakly density-modulated states are well-described by the cosine-modulated (CM) (also see [13,19,30]). The CM ansatz replaces the axial wavefunction by the variational parameters θ, describing the amplitude of the density modulation, and L, specifying the wavelength of the modulation [see Fig. 2(a)]. Thus Ψ CM depends on the four variational parameters {θ, L, l, η}. Note that the factor of √ 2 ensures that this ansatz has an average density n per unit cell (independent of θ). We restrict θ to the range 1 We note that θ directly relates to the density contrast as which decreases with increasing θ, until f s (ϕ) = 0. Using this ansatz we can analytically evaluate the energy per particle as where Λ(θ) = 1 32 (90 cos θ − 55 cos 3θ − 3 cos 5θ). We determine solutions Ψ CM by numerically minimising E CM with respect to its variational parameters. We expect this ansatz to accurately describe the system's behavior in the vicinity of the continuous transition.
In Fig. 2 6)]. Using the reduced 3D theory we find that the system undergoes a continuous transition from the uniform BEC state into a modulated state at a * ≈ 91.6a 0 . Both theories are in excellent agreement close to the transition and show that for a s < a * the contrast initially develops as C ∼ √ a * − a s . For a s values further below a * , the theories begin to deviate as a strong density modulation develops.
Our results also show the importance of the full 3D description of the system [see Figs. 2(c) and (d)]. While the confinement is radially isotropic magnetostriction causes χ to highly elongate in the dipole direction with an anisotropy of η ∼ 4. Interactions also cause l to be significantly greater than the harmonic oscillator length l ho = 0.64µm.
Using the Ψ CM ansatz we can identify the conditions for the continuous transition by looking for stationary points of E CM for small θ. To leading order in θ we have ∂ECM ∂θ = θ∆( 2π L ), where ∆(k) = 0 (k) + 2nŨ (k) + 3g QF n 3/2 , with 0 ≡ 2 k 2 /2m. Thus the stationary points of E CM in this regime are either: (i) a uniform solution with θ = 0 (in which case L is irrelevant); (ii) the modulated solutions with θ = 0 and L determined by ∆( 2π L ) = 0. We note that the Bogoliubov spectrum for the excitations of the uniform condensate θ = 0 is given by (k) = 0 (k)∆(k) [31], thus the condition ∆( 2π L ) = 0 means that an excitation of wavelength L has zero energy, i.e. a roton-like excitation in the system goes soft [34,35]. In Fig. 2(e) we show the uniform BEC excitation spectrum for various a s values, observing the formation of a roton at a s ≈ 95a 0 that softens to zero energy at 91.6a 0 . This marks the dynamic instability of the uniform BEC state and defines the roton critical value a * rot . At this critical point the modulated state develops with a wavelength corresponding to the roton wavevector [ Fig. 2(d)]. In the regime where the transition is continuous a * coincides with a * rot . We see that this holds for the results in Figs. 1(b) and (c), with a * rot = a * for n ∈ [n low , n high ]. Outside of this density range (where the transition is discontinuous) we see that a * rot < a * . length scales as as varies. Solid (dotted) lines are used for each theory when its energy is lower (higher) than ΨBEC. In (b) the CM ansatz results are also shown (red line). In (d) the critical roton wavelength is indicated at a * rot . In (e) σz is the 1/e halfwidth of |Ψ| 2 along the z axis. Parameters as in Fig. 1 with n = 0.7 × 10 3 /µm.
Low-density discontinuous transition region -For densities n < n low the contrast has a finite jump as the transition to the modulated state is crossed. In Figs. 3(b)-(e) we focus in on this regime of the phase diagrams in Fig. 1. We consider the particular case of n = 0.7×10 3 /µm with a transition point at a * = 83.3a 0 , identified as where the energy per particle E of the modulated stationary solution Ψ crosses E BEC [ Fig. 3(b)]. This transition point is appreciably higher than a * rot = 81.2a 0 , where the BEC becomes dynamically unstable, meaning that hysteresis can occur in a s ramps across the transition with the uniform or modulated state persisting as a metastable state. Results of the reduced 3D theory shows that the transition occurs with a sudden jump in the transverse properties of the ground state [Figs. 3(c) and (e)], and that the unit cell size is larger than (and disconnected from) the roton wavelength [ Fig. 3(d)]. The contrast (not shown) remains close to unity up until the modulated solution branch terminates as a metastable state.
The Ψ CM ansatz fails to describe this regime [e.g. Fig. 3(b)] and incorrectly predicts the transition to occur continuously at the point of roton softening. Here the reduced 3D model indicates that the system prefers to organize into localized well separated droplets, exhibiting a high contrast modulation of the density. This motivates us to introduce the droplet array represents a Gaussian droplet of z-width l z and containing N D atoms. The droplets repeat every L [i.e. one per unit cell, see Fig. 3(a)], with the relation N D = nL ensuring that the average density is fixed to n. We require well separated droplets (L l z ) for ψ DA to avoid a discontinuity at the unit cell boundary. We can evaluate the energy per particle of this ansatz where κ = η 1/4 l/l z , ζ is the Riemann zeta function, and f = 1+2κ 2 −3κ 2 atanh (also see Ref. [36]). All terms in Eq. (9) are evaluated exactly except the last term describing the long-range interaction between droplets. That term is obtained by approximating each droplet as a point dipole of N D atoms, which is a good approximation for our regime with l z L. The Ψ DA solutions are determined by numerically minimising E DA with respect to {L, l z , l, η}.
The results in Figs. 3(b)-(e) show that Ψ DA is in good agreement with Ψ. In Fig. 3(e) we show the droplet width σ z (l z ) along z from Ψ (Ψ DA ), and observe that these remain much smaller than their spacing L [ Fig. 3(d)] up until their solution branch terminates. At the termination point the modulated solution is metastable (E BEC is lower), and occurs because the droplets unbind (l z , σ z → ∞). High-density discontinuous transition region -In Figs. 4(a)-(d) we focus in on the high-density (n > n high ) discontinuous transition of the phase diagrams in Fig. 1. We In (c) we also indicate the roton wavelength at a * rot for reference. (d) Energy per particle comparison of the theories. The solid lines are used for each theory when its energy is lower than ΨBEC. Parameters as in Fig. 1 with n = 6.25 × 10 3 /µm. consider the particular case of n = 6.25 × 10 3 /µm, where the transition occurs at a * = 89.6a 0 (cf. a * rot = 88.6a 0 ). As we cross the transition the ground state Ψ suddenly changes its transverse profile (η and l) and develops a strong density modulation. In this regime the discontinuous transition arises from an interplay between the transverse degrees of freedom and the modulation along z mediated by the dominant role of interactions at higher densities. The behavior of Ψ is well described by the Ψ CM ansatz, however the small-θ stationary point of Ψ CM that we analyzed earlier to understand the continuous transition is an unstable saddle point in this regime. Here the (meta)-stable CM solution found [see Figs. 4(a)-(d)] has a large θ value (i.e. high contrast) and a longer wavelength than the roton wavelength [ Fig. 4(c)]. This solution also predicts transverse properties (l, η) and a * > a * rot , similar to the reduced 3D theory. We also note that at this density the DA solution terminates at a s ∼ 80a 0 , well below the scattering length range considered in Fig. 4. Summary and outlook -We have developed theory for a tube confined dipolar quantum gas to understand its phase diagram. In Fig. 5(a) we summarize the phase diagram quantitatively presented in Fig. 1, indicating the three phases. At low densities there is a direct discontinuous transition between the insulating droplet 2 and BEC phases. At higher densities a  [18]. The inset shows the rescaled data. supersolid phase emerges, separating the BEC and insulating droplet phases, and the BEC-supersolid transition can be continuous within a certain density range.
In Fig. 5(b) we characterize how the phase diagram changes with system parameters. We show the continuous transition range [n low , n high ] as vertical lines obtained from phase dia-grams like Fig. 5(a) computed for the two relevant experimental atomic species (Er and Dy) and various (isotropic) transverse confinement strengths. We also show n max , defined as the density where a * rot (and a * ) is maximised. These results indicate that the characteristic densities {n low , n high , n max } decrease with increasing radial trapping. Also that a continuous transition in Er requires a higher density than Dy, which will shorten the system lifetime due to three-body loss. A recent 3D eGPE calculation by Roccuzzo et al. [18] was performed for the tube-dipolar system and found a small jump in f s in the BEC-supersolid transition [for parameter set marked "x" in Fig. 5(b)]. This may indicate the weak first order transition, and that the continuous region is narrower in the full theory.
By rescaling the eGPE using a dd and ω dd = /ma 2 dd as units of length and frequency, respectively, the resulting equation only depends on the scaled axial density na dd , trap frequencies ω x,y /ω dd , and dd . The inset to Fig. 5(b) shows the collapse of the rescaled Er and Dy results in these units. We also observe that the characteristic densities scale with the transverse confinement frequency as ω −1/2 x,y . The phase diagram and scaling we have predicted could be explored in future experiments, e.g. by changing the atom number to pass through the transition at different densities, or by changing the transverse confinement to shift the location of the continuous transition region. Additionally, we note that having tighter confinement along the dipole direction relative to other transverse direction (i.e. ω y > ω x ) can introduce additional discontinuous transitions into the phase diagram, similar to those observed in oblate pancake shaped traps [37,38]. The CM and DA analytic models we have presented provide simple tools to map out the phase diagrams and should aid in exploring and better understanding this fascinating system. Another important future step will be to include the axial trapping potential which will introduce finite size effects and a spatially dependent mean axial density.