A general theory of flattened dipolar condensates

We develop theory for a flattened dipolar Bose-Einstein condensate (BEC) produced by harmonic confinement along one direction. The role of both short-ranged contact interactions and long-ranged dipole-dipole interactions (DDIs) is considered, and the dipoles are allowed to be polarised along an arbitrary direction. We discuss the symmetry properties of the condensate and the part of the excitation spectrum determining stability, and introduce two effective interaction parameters that allow us to provide a general description of the condensate properties, rotons, and stability. We diagonalize the full theory to obtain benchmark results for the condensate and quasiparticle excitations, and characterize the exact mean field stability of the system. We provide a unified formulation for a number of approximate schemes to describe the condensate and quasiparticles, including the standard quasi-two-dimensional (quasi-2D) approximation, two kinds of variational ansatz, and a Thomas-Fermi (TF) approximation. Some of these schemes have been widely used in the literature despite not being substantiated against the exact theory. We provide this validation and establish the regimes where the various theories perform well.

Quantitative calculations for the mean field ground states and quasiparticle excitations of flattened dipolar BECs are difficult, because fundamentally the problem is still fully three-dimensional (3D) and requires a number of specialized numerical techniques to deal with the DDI [30]. As such the most comprehensive calculations are restricted to cylindrically symmetric regimes, which requires the dipoles to be polarized along the symmetry axis of the trap. One simplification, commonly referred to as the quasi-two-dimensional (quasi-2D) approximation [31,32], is based on the assumption that the system is in the harmonic oscillator ground state in the tightly confined direction, which can then be integrated out leaving an effective 2D theory. This approximation is only valid when interactions are weak, and many regimes of interest (e.g. roton excitations or the collapse instability) occur where this condition is violated [31] (also see [33]). For systems with 3D confinement, a comparison of the quasi-2D approximation to other approaches was presented by Wilson and Bohn [34]. That work suggests that the quasi-2D approximation furnishes a qualitatively useful description, even though it is quantitatively inaccurate, e.g. the DDI strength at which collapse occurs is overestimated by a factor of about three (see figure 3 of [34]). An increasing number of studies of flattened dipolar BECs have been based on the quasi-2D approximation (e.g. see [32,[35][36][37][38][39]). It is therefore essential to develop a better quantitative theory for the flattened system in the regime where the quasi-2D approximation is only qualitatively correct, particularly because of the interesting physics that is predicted to occur in this regime.
In this paper we present a general treatment of a flattened dipolar BEC. We focus on the situation where the confinement is harmonic along one direction, with free in-planar motion. We allow for the dipoles to be tilted at an arbitrary angle into the plane. We formulate mean field equations for the condensate and quasiparticle excitations. We observe that the condensate possesses a symmetry, whereby its properties are determined by a single effective interaction parameter accounting for the contact interaction, DDI, and the effect of tilt. Additionally, we show that the part of the excitation spectrum that determines the onset of rotons and collapse instability possesses another rigorous symmetry. This motivates the introduction of a second effective interaction parameter, and allows us to provide a universal description of the stability boundary. We numerically solve the exact mean field equations for condensate and quasiparticles. This allows us to obtain benchmark results for the universal properties of the condensate and the stability phase diagram, which are provided as data files with this paper 1 . We then turn to considering approximate descriptions of the system which involve a non-trivial reduction of the system to a 2D problem by assigning a specified profile for the condensate and quasiparticle profiles in the confined direction. In detail we consider the standard quasi-2D approximation, a Gaussian (i.e. single mode) and two-mode variational ansatz, and a Thomas-Fermi (TF) approximation. We unify the formulation of these approximate descriptions and consider how well they perform relative to our exact results at predicting the properties of the condensate, the excitation spectrum and stability.

Gross-Pitaevskii equation
We consider here a dipolar BEC that is harmonically confined along the z direction and unconfined in the radial plane. The condensate wave function ψ 0 satisfies the non-local Gross-Pitaevskii equation (GPE) where ω z is the trap frequency and μ is the chemical potential. We consider DDIs for the case of dipoles polarized by an external field along the direction α α = + e z xˆcosˆsin , lying in the xz-plane. We refer to α as the tilt angle of the dipoles, with the untilted case α = 0 corresponding to dipoles polarized along z. The associated interaction potential is is the DDI coupling constant and = | | r r r . The particles can also interact by a short ranged contact interaction with coupling constant , where a s is the scattering length, so that the full interaction we use is [40,41]. The condensate solution takes the form ψ χ = ν n z r ( ) ( ) , obtained by taking the expectation of  GP . Their values, obtained from our numerical calculations for χ ν , are plotted in figure 2(b) and are made available as raw data files (see footnote 1). Noting that for the non-interacting case ν = 0, we have μ = = Ē¯z 1 2 and γ π = 1 2 .

Bogoliubov excitations
The excitations of a condensate are Bogoliubov quasiparticles. For our geometry these can be obtained by linearizing the time-dependent GPE 3 about the condensate using an ansatz (e.g. see [30,45]) of the form In the above expansion we have introduced the planar coordinate vector ρ = x y ( , )and wave vector = ρ k k k ( , ) x y . Because the system is homogeneous in the ρ-plane, the quasiparticles are of well defined momenta  ρ k . The additional index j labels the z vibrational mode. The Bogoliubov-de Gennes (BdG) equations that must be solved to determine denotes the Fourier transform in z, −  1 denotes its inverse, and (14) is the Fourier transform in r of the interactionU r ( ). We numerically solve the BdG equations (11) using the same discretization we employed to solve for χ ν . The discretized equations are then diagonalized using standard linear algebra packages [46]. For accurate calculation of the excitations it is important to employ a cutoff interaction, see appendix C (also [30]). An example of the spectrum and ground band quasiparticle amplitudes is shown in figure 3. We note that the ground band spectrum in figure 3(a) shows some key characteristics of a planar dipolar condensate: phonon (linear) behavior at low momentum (up to ∼ k l 0.3 z ), a softened spectrum at ∼ k l 1 z arising from the attractiveness of the DDI, and free particle (quadratic) behavior at high k where the kinetic energy dominates. For stronger dipolar interactions the spectrum can develop a roton minimum at ∼ k l 1 z . The role of interactions in higher bands ( > j 0) is less important (see figure 3(a)). 3 Equation (1) The backaction of quasiparticles could be included in more advanced theories [44].

k y -spectrum symmetry
In general the quasiparticles must be solved numerically, however the spectrum along the k y axis (i.e. , possesses an interesting symmetry property 4 . To see this we note that for excitations along k y , which we denote by setting → ρ k k y in the subscripts, equation (13) simplifies to The overall action of  remains constant if ν and ν d remain constant. These two constraints define a curve in α g g { , , } s d parameter space where the k y spectrum is invariant (cf the GPE solution χ ν z ( )which is invariant on ν constant surfaces). E.g., the k y -spectra for the parameter sets This result is useful because for > g 0 d the stability of the system, and the development of roton excitations, is addressed by examining excitations on the k y -axis, and in the j = 0 band. We note that if ν ν < 3 d then the interaction along k y is always repulsive.
In figure 4 we show the condensate mode and excitation spectrum for a dipolar system with α = = g 0 s and ν = 2.7434 d (i.e. ν = 5.4868). The spectrum reveals that the system has a roton (local minima) at ≈ ρ k l 1 z , and with approximately zero energy. Increasing g d further will cause the roton quasiparticles to become dynamically unstable, leading to local collapse of the condensate (e.g. see [7,10,[47][48][49][50][51]). For this untilted case the spectrum is isotropic, however this spectrum also applies along k y for other interaction parameter sets with α ≠ 0, according to the mapping (17a), (17b). As an example we present a result for α = π 3 in figure 4(b). Along the k y axis the spectrum is the same as the untilted case, but along the k x axis the spectrum differs 5 .
As noted earlier, the k y -spectrum is only a function of ν and ν d . Thus, using these as coordinates we can map out the stability and roton phase diagram for the system, with results shown in figure 5(a). In the shaded region the excitation spectrum is stable, but has a roton at ≠ k 0 y . We identify the stability boundary as where a quasiparticle mode has zero energy. That is, to the left of the stability boundary the mode softens further and its eigenvalue E k y becomes imaginary, signalling that the condensate is dynamically unstable. We emphasise that this phase diagram is universal for any tilt angle and the raw data for the stability boundary is provided (see footnote 1). As an example, in figure 5(c) we show the phase diagram mapped to the untilted case as a function of g s and g d .

General approximate schemes
In what follows we investigate different approximate schemes for obtaining the condensate and excitations. These all have in common the assumption that the axial profile of the excitations is identical to the condensate, with results independent of direction in the ρ k plane. This solution also applies to α′ = π 3 using the mapping (17a), (17b) with spectrum along the k y axis (solid, same as α = 0) and along the k x axis (dashed, only shown for BdG and the same shape approximation). 5 We note that u z ( ) can be taken real (e.g. see figures 3(b), (c)). However, for α ≠ 0 and ≠ k 0 i.e.
where ρ u k and ρ v k are constants and we have dropped the superscript j, restricting ourselves to the ground mode in z. Our full numerical solutions for the k y -quasiparticles (see figures 3(b) and (c), particularly the insets) reveal that this is typically a good approximation for the u-modes, but is a poorer approximation for the v-modes 6 . This deformation of the v-modes tends to get worse for large DDIs, notably when a roton forms (hence near instability), but we note that where the v-modes differ significantly in shape they tend to be relatively small in absolute norm, and thus it is not clear how important neglecting this feature will be. For this reason we pay attention to how well the stability boundary is predicted by the various theories we derive based upon this assumption.
The various approximate schemes we consider differ in the choice of χ ν . Making assumption (18) allows us to integrate out the z dimension, and reduce the problem to an effective 2D uniform system. Key to this simplification is the introduction of the part of the effective 2D k-space interaction that does not contribute to the GPE (3), i.e. 7 2 andU k ( )was given in equation (14). We have also introduced (8) and, since the integrand is always positive, ν G q ( )is monotonically increasing in q. The spectrum is then obtained from the BdG equations (11) For k x = 0 or α = 0 the spectrum simplifies to

Same shape approximation
The most immediate application is to use χ ν obtained from the numerical solution of the GPE mode to describe the axial shape of the quasiparticles. We refer to this as the same shape approximation. This approach is hence fully numerical, but does not require numerical diagonalization, and provides a useful comparison to the other approximate schemes we now introduce. For numerical calculations (i.e. when not using the analytical results from the approximate schemes below), it is important to use ν G q ( )subject to a cutoff interaction which is shown in appendix C. We emphasise that because the same shape approximation involves the full GPE solution, all condensate quantities (e.g. μ γ Ē, ,¯z) are exact, and we use the term same shape to refer to the approximate description of the excitations.

Variational and quasi-2D approximations
An additional variational approximation is to take χ ν to be an ansatz depending on one or more variational parameters to be determined by minimizing the energy: The variational Gaussian approximation is of the form with length scale σl z , where σ is a variational parameter, giving ν σ σ = + Ē ( ) 1 (4 ) 4 z 2 2 and γ π σ = 1 ( 2 ). Using χ σ , the k-space interaction (20) can be evaluated with q 0 2 with erfc the complementary error function. In general the optimal value of σ is most easily determined numerically, however for small ν The quasi-2D approximation is when interactions are sufficiently weak (ν → 0) that χ can be approximated as the non-interacting harmonic oscillator ground state with σ = 1 (e.g. see [27]). A more complex but more accurate two-mode variational approximation is based on the ground and second excited harmonic oscillator states is the second Hermite polynomial and a 2 is an additional variational parameter. This variational ansatz was first introduced by Wilson and Bohn [34] for application to multilayered condensates (i.e. sets of planar systems as produced by a one-dimensional optical lattice). Expressions for γ, ) using the two-mode approximation are given in appendix D.

Axial TF approximation
In the opposite regime to the quasi-2D case, i.e. for ν ≫ 1, interactions dominate the confinement energy in the z direction. In this regime a TF approximation is applicable (see [16]) This ansatz gives

Comparison of approximate schemes
The approximate theories can be compared to each other by assessing the differences in their predictions for the functions μ ν ( )(see figure 6), γ ν ( )(see figure 7) and ν G (see figure 8), which are all presented over a broad range of ν values. These results make it clear that the quasi-2D approximation is only applicable for ν ≪ 1. The twomode variational results provides the best approximation for the range of interactions considered (e.g. see insets to figures 6 and 7), although the TF approximation eventually becomes a better approximation at ν ∼ 25. While the variational Gaussian is everywhere worse than the two-mode scheme, it is more accurate than the TF approximation until ν ∼ 10.
The stability boundary is determined by the excitation spectrum, and the predictions of the various theories are compared against the full numerical result in figure 5. Similar to our findings for μ and γ, we see that the variational Gaussian is better than the axial TF description for interactions up to about ν ∼ 10. However, the same shape approximation always furnishes a better description than either the axial TF or variational Gaussian at all values of ν considered. We also consider the differences between ν G (21), which determines the spectrum via (22), for the approximate theories in figure 8. Notably, figure 8(b) shows that the difference between the same shape and variational Gaussian result increases with increasing ν. The two-mode variational result (figure 8(c)) has similar behavior, but the size of the difference is significantly smaller. In figure 8(d) the TF approximation is considered, and we observe that the difference from the same shape results decreases with increasing ν.
In figure 4 we show the condensate solutions and excitation spectrum using the full and various approximate theories for an untilted purely dipolar case at instability (i.e. at ν = 2.743 d and ν ν = 2 d ). The spectrum shows that the full solution has a zero energy roton, whereas the various approximate theories have the roton minima at higher energies, except for the quasi-2D theory which does not have a roton for the interaction parameters considered. Increasing the dipole strength (i.e. moving along the dashed-dotted line in figures 5(a) or (c)) until each approximate theory predicts a zero energy roton we find their respective critical values ν d * . These are tabulated in table 1 for comparison. The quasi-2D approximation is in worst agreement with the full theory and overestimates the stability value by a factor of about 1.6.
To put the results of this paper into context it is useful to consider parameters that are achievable in current experiments. Consider a condensate of 164 Dy atoms [3], which have a magnetic moment of μ μ = 10 B (where μ B is the Bohr magneton). For the case of = × N 1 10 5 condensed atoms confined in a pancake harmonic trap of  , and for untilted dipoles, we have ν = 2.06 d and ν = 5.69 at trap center (ρ = 0). 8

Conclusions and outlook
In this paper we have studied the properties of a dipolar BEC flattened by confinement in the z direction. We have pointed out the rigorous symmetries of the condensate and its excitation spectrum, which allow us to reduce the condensate solution to its dependence on a single effective interaction parameter ν that includes all aspects of the interactions. Similarly we find that the k y -excitation spectrum, which is key to determining the emergence of rotons and stability, depends on two effective interaction parameters ν ν { , } d , allowing us (aided by exact numerical calculations) to provide a universal stability phase diagram in this parameter space.
A second major focus of this paper has been to introduce and unify a number of analytic approximations for the condensate and excitations. The basis of these approaches consists of approximating the condensate profile in the confined direction by a specified profile. We introduced four analytic forms for this profile, some of which for ν from 0 (blue) to 20 (green) in steps of 2. Inset:G q ( ) 0 .  8 To obtain this estimate we have taken the radial density profile to be described by a TF approximation (see [33]), and we have assumed a background s-wave scattering length of 100 a 0 (a 0 is the Bohr radius). 9 i.e. g s = g d sin 2 α. This includes the untilted purely dipolar gas.
have been extensively used in the literature, particularly as a basis to simplify simulating the Gross-Pitaevskii dynamics of the condensate. These approximations are validated against our full numerical calculations (without approximation), and the same shape approximation. We validate the accuracy of these approximations by considering both condensate properties and the quasiparticle excitation spectrum. The best approximation is shown to be dependent on the strength of the interactions, and importantly that the widely used quasi-2D approximation is inaccurate except for very weak interactions where dipolar effects are small. The detailed understanding of the uniform planar case provided in this paper will underpin future work on tractable calculations of fully 3D trapped dipolar condensates without requiring cylindrical symmetry. This theory could be used to model the 3D trapped system with tight axial confinement and weaker in-plane trapping, by taking a wave function ansatz of the form where ν depends on ρ via the areal density ρ n ( ) (see equation (6)). With in-plane trapping, the radial Fourier transform of the density has a finite spread in ρ k (cf the delta function in equation (A.1) for uniform in-plane density). The GPE then has two interaction terms, a contact term proportional to ρ ν ( ) locally (see equation (4)) and a long-range term depending on ρ n ( ) everywhere. This long-range term has been discussed in the context of the quasi-2D approximation in [32], and can be important in weakly anisotropic traps (e.g. ω ω ≲ 20 z x y , , with ω x and ω y the radial trap frequencies) where it gives rise to density oscillating ground states [9,10]. For higher aspect ratio traps the radial density varies more slowly, and the long-range term becomes less significant.
Stability of dipolar BECs has been of significant experimental interest [7], and while dipole tilt is an easily accessed control parameter in experiments, accurate theory has been limited to situations of untilted dipoles (e.g. see [9,10]) where the system has cylindrical symmetry. Our universal stability results fully quantify the role of tilt and present an opportunity for validation in current experiments. With radial trapping, the density changes with the interactions, notably the density will tend to decrease as the interactions increase because the condensate radially spreads out. The instability condition will occur when the maximum value of ν reaches the critical value predicted in the uniform case, although there will be corrections due to the effect of the radial trapping on the quasi-particle excitations.
Finally, we note that our work can be also extended down a number of other avenues, e.g., as the basis of descriptions for local collapse [52] and condensate dynamics (e.g. [53]), the calculation of static [54] and dynamic correlation functions, and finite temperature formalism for a trapped partially condensed dipolar condensate [55].  Figure D1. Variational parameters (a) σ for Gaussian (blue) with the large ν limit σ ν π ≈ (2 ) 1 3 1 6 (blue dashed) and two-mode (cyan) with the large ν limit σ ν ≈ 0.4916 1 3 (cyan dashed) and (b) a 2 and the large ν limit ≈ a 0.3186 2 (dashed).