Comptes Rendus
Physique

Wilhelm Zwerger

Reflections on dipolar quantum fluids

Volume 25 (2024), p. 389-413

Online since: 6 December 2024

https://doi.org/10.5802/crphys.222

<figure><img src="image_1.png" /></figure>

This article is licensed under the

CREATIVE COMMONS ATTRIBUTION 4.0 INTERNATIONAL LICENSE.

http://creativecommons.org/licenses/by/4.0/
Intervention in a conference / *Intervention en colloque*

# Reflections on dipolar quantum fluids

## Réflexions sur les fluides quantiques dipolaires

Wilhelm Zwerger<sup>●,a</sup>

<sup>a</sup> Physik Department TU München, James Franck Strasse, 85748 Garching, Germany

E-mail: zwerger@tum.de

**Abstract.** We present a thermodynamic description of ultracold gases with dipolar interactions which properly accounts for the long-range nature and broken rotation invariance of the interactions. It involves an additional thermodynamic field conjugate to the linear extension of the gas along the direction of the dipoles. The associated uniaxial pressure shows up as a deviation from the Gibbs-Duhem relation in the density profile of a trapped gas. It has to vanish in self-bound droplets, a condition which determines the observed dependence of the aspect ratio on particle number. A tensorial generalization of the virial theorem and a number of further exact thermodynamic relations are derived. Finally, extending a model due to Nozières, a simple criterion for the freezing transition to a superfluid mass density wave is given.

**Résumé.** Nous donnons une description thermodynamique des gaz froids en interaction dipolaire qui tient compte de la longue portée de ces interactions et de leur absence d'invariance par rotation. Elle fait intervenir un champ thermodynamique supplémentaire conjugué à l'extension linéaire du gaz selon la direction des dipôles. La pression uniaxiale associée se manifeste par un écart à la relation de Gibbs-Duhem dans le profil de densité du gaz piégé. Elle doit s'annuler dans les gouttelettes (des états liés du système), une condition qui détermine la dépendance observée du rapport d'aspect avec le nombre de particules. Une généralisation tensorielle du théorème du viriel et un certain nombre d'autres relations thermodynamiques exactes sont obtenues. Enfin, en étendant un modèle dû à Nozières, on aboutit à un critère simple pour la transition de solidification vers une onde de densité de masse superfluide.

**Keywords.** Ultracold atoms, Dipolar fluids, Supersolids, Thermodynamic limit, Long-range interactions.

**Mots-clés.** Atomes froids, Fluides dipolaires, Supersolides, Limite thermodynamique, Interactions à longue portée.

Manuscript received 15 November 2024, accepted 18 November 2024.

## 1. Introduction

The realization of a Bose Einstein condensate with Chromium [1] has opened a new field of research in ultracold gases. It allows to explore a wide range of phenomena which are uniquely tied to the long-range and partially attractive nature of dipolar interactions. Following the extension to condensates with Erbium or Dysprosium where the strength of these interactions may exceed the short-range repulsion, the field has grown immensely in recent years. This is based, in particular, on the discovery of self-bound droplets [2] and of supersolid phases [3–5], where superfluidity coexists with broken translation invariance: see [6] for a recent review of magnetic dipolar gases and the Lectures [7, 8] for an introduction to the underlying concepts. On the theory side, a successful qualitative description of the observations is provided by an extended Gross–Pitaevskii equation, where a non-analytic contribution proportional to |ψ|⁵ is added to the energy functional. As will be discussed below, a proper microscopic derivation of
this procedure in the relevant regime of strong dipolar interactions, where the chemical potential may assume negative values, is lacking. It is therefore of interest to develop a description which only relies on general thermodynamic relations, properly accounting for the long-range and anisotropic nature of the interactions. A step in this direction will be taken in the present contribution for quantum fluids in thermal equilibrium, extending concepts from classical polar fluids [9] and liquid crystals [10]. In particular, it is shown that the anisotropy and long-range nature of the interactions require the introduction of an additional contribution $h_1 dL_z$ in the differential of the free energy, with $L_z$ the system length along the direction of the dipoles. It gives rise to the shape-dependence of thermodynamic properties and an anisotropy in the momentum current tensor $\Pi(x)$ and thus effectively to pressure. The appearance of an additional extensive variable $L_z$ in the free energy $F(T, V, L_z, N)$ beyond particle number and volume leads to a violation of the Gibbs-Duhem relation which is reflected in the density profile of harmonically trapped gases. Moreover, for the case of self-bound droplets where the internal forces have to balance locally, the vanishing of the uniaxial contribution implies a scale-dependent aspect ratio $\kappa(N) \sim N^{-1/4}$, approaching zero in the thermodynamic limit where the droplets evolve into a needle-like shape. Finally, we address the nature of the transition from a homogeneous superfluid to a supersolid and discuss a simple model due to Nozières [11] which provides a criterion for the associated critical value of the roton gap beyond mean-field.

In general, the aim of the present contribution is to clarify a number of conceptual points in the theory of dipolar gases, thus providing a better understanding of some fundamental issues and open problems.

## 2. Thermodynamics of uniaxial quantum fluids

In close analogy to the standard description of classical polar fluids [9], the two-body interaction

$$V(x_{12}) = V_{sr}(r_{12}) - \frac{\mu_m^2}{4\pi} \frac{2P_2(\hat{z} \cdot \hat{x}_{12})}{r_{12}^3} \chi(r_{12}) \quad (1)$$

in Bose quantum fluids whose permanent dipoles $\mu_m$ are all oriented along the z-direction may be separated into a rotation invariant short-range plus the long-range, anisotropic dipolar contribution (we use magnetic dipoles in the following and units where $\mu_0 \equiv 1$ but our results also hold for electric ones with minor changes). Here, $P_2(x) = (3x^2 - 1)/2$ is the standard Legendre polynomial and $\hat{z}, \hat{x}_{12}$ are unit vectors along, respectively, the dipole orientation or the separation vector $x_{12} = x_1 - x_2$. In the ultracold limit, the two-body scattering amplitude associated with $V_{sr}(r_{12})$ only involves s-wave scattering and thus is fully characterized by the scattering length $a$, which is assumed to be positive. The short-range interaction can therefore be replaced by a zero-range pseudopotential with strength parameter $g = 4\pi\hbar^2 a/m$. A similar simplification is possible for the dipolar interaction, whose strength defines an additional effective length $a_{dd} > 0$ by $g_{dd} = \mu_m^2/3 = 4\pi\hbar^2 a_{dd}/m$. Indeed, dipolar scattering at low energies arises from large separations, where $V_{dd}(x_{12}) \sim \mu_m^2/r_{12}^3$ is weak. Asymptotically, the two-body scattering amplitude

$$f_{dd}(\hat{k}_i \rightarrow \hat{k}_f, |k|) \underset{|k|a_{dd} \ll 1}{\longrightarrow} f_{dd}^{(B)}(\hat{k}_i \rightarrow \hat{k}_f) = -a_{dd} \cdot 2P_2(\hat{z} \cdot \hat{q}) \quad \text{for } |q| = 2|k| \sin(\theta/2) \neq 0 \quad (2)$$

is therefore given by the Born approximation, proportional to the Fourier transform $V_{dd}(q) = g_{dd}2P_2(\hat{z} \cdot \hat{q})$ of the bare dipolar interaction at the momentum transfer wave vector $q = k_f - k_i$. Note that $V_{dd}(q)$ is independent of $|q|$ only for $|q| r_0 \ll 1$, where the short distance cutoff function $\chi(r)$ in Equation (1) may be replaced by unity and also that the forward scattering limit $\theta = 0$ is ill-defined. Moreover, since $V_{dd}(x_{12})$ is not rotationally symmetric, the scattering amplitude cannot be decomposed into separate angular momenta and also depends on the initial direction $\hat{k}_i$<fnref n="1" />.

<footnote n="1">¹A discussion of how this affects rethermalization in dipolar gases has been given by Bohn and Jin [12].</footnote>
The validity of the result (2) has been tested by Bohn et al. [13] based on a numerical solution of the full dipolar two-body scattering problem. Using a pure power law dependence of $V_{dd}(x_{12})$ with a sharp cutoff at a distance $r_0 \ll a_{dd}$, the cross section $\bar{\sigma}(|k|)$ averaged over all incident directions $\hat{k}_i$ is found to essentially coincide with the Born approximation value obtained for $|k|a_{dd} \to 0$ up to energies of order $\hbar^2/(ma_{dd}^2)$. For typical values $a_{dd} \approx a \approx 100 a_0$, dipolar scattering in ultracold gases is therefore fully described in terms of a single parameter $a_{dd}$ and deviations of the actual dipolar interaction from a pure power law at short distances $r_{12} < r_0$ are not relevant. This standard assumption, where point-like dipoles are combined with a zero-range pseudopotential has to be used with caution, however. Indeed, the cutoff scale $r_0$ explicitly enters the corrections beyond the Born approximation. Specifically, while both the real and the imaginary part $\text{Im} f_{dd}(0) = |k|\sigma_{dd}(\hat{k}_i)/(4\pi) \sim |k|a_{dd}^2$ of the forward scattering amplitude vanish at low energies and are independent of $r_0$, the real part $\text{Re} f_{dd}^{(2B)}(0) \sim a_{dd}^2/r_0$ is finite in second order Born approximation and linearly sensitive to the short distance cutoff. A more serious problem arises at the many-body level where the expectation value of the dipolar interaction

$$ \langle \hat{H}_{dd} \rangle = \int_{x_1} \int_{x_2} \left[ -\frac{\mu_m^2}{8\pi} \frac{2P_2(\hat{z} \cdot \hat{x}_{12})}{r_{12}^3} \right] \rho^{(2)}(x_1, x_2) \quad \text{with} \quad \rho^{(2)}(x_1, x_2) = \rho^{(1)}(x_1)\rho^{(1)}(x_2)g^{(2)}(x_1, x_2) \quad (3) $$

depends on the exact two-particle density $\rho^{(2)}(x_1, x_2)$. It thus requires knowledge not only of the one-particle density $\rho^{(1)}(x)$ of the atomic cloud but also of the associated two-particle distribution function $g^{(2)}(x_1, x_2)$ (we consider a general inhomogeneous situation, following the notation for classical fluids in [9]). Expressed in terms of center-of-mass and relative coordinates $X$ and $x_{12}$, the integral in (3) is well defined despite the short-distance singularity $1/r_{12}^3$ provided $\rho^{(2)}(X, x_{12})$ is regular in the limit $x_{12} \to 0$ (for a mathematically precise statement see e.g. [14]). The disappearance of the singularity due to the vanishing angular average of $P_2(\hat{z} \cdot \hat{x}_{12})$ is spoiled, however, if the effective range of $V_{sr}(r_{12})$ is set to zero. Indeed, provided the actual interaction becomes isotropic at short distances, the two-particle density is asymptotically determined by the one obtained from a pseudopotential approach where

$$ \rho^{(2)}(X, x_{12}) \xrightarrow{x_{12} \to 0} \frac{\mathcal{C}_2(X)}{(4\pi)^2} \left( \frac{1}{r_{12}^2} - \frac{2}{ar_{12}} + \dots \right) \quad (4) $$

is singular at short distances [15, 16]. Here, $\mathcal{C}_2(X) \to [4\pi\rho^{(1)}(X)a]^2$ is the Tan contact density in the weak interaction limit, which replaces the ill-defined local pair distribution function $g^{(2)}(0)$ for zero-range interactions (see [7] for an introduction to the Tan relations and the associated singular behavior of short distance correlations for Bosons). In practice, this problem—which does not show up at the mean-field level where $g^{(2)}(x_1, x_2) \to 1$ is replaced by one—is avoided if the expectation value of the dipolar interaction energy is evaluated with a two-particle density $\rho^{(2)}(X, x_{12})$ which is averaged over separations large compared to the actual interaction range. In particular, for $r_{12} \gg a$, the most singular term in Equation (4) is negligible and the next-to-leading term only gives rise to a small correction $\delta\rho^{(2)}(X, x_{12}) = -[\rho^{(1)}(X)]^2 2a/r_{12}$ to the mean-field result. As will be shown below, a cutoff at short distances is also required to properly account for the angular dependence of the pair distribution function.

Beyond the issue of the singularity of the dipolar interaction at short distances for Bosons, which generically have a finite probability to be at the same point in space, a major problem that needs to be addressed in dipolar fluids concerns the long-range and anisotropic, partially attractive nature of the interaction. Indeed, the standard proofs by Ruelle [17] and Fisher [18] on the existence of a proper thermodynamic limit require interactions which decay faster than $1/r^3$ for large distances, leaving dipolar fluids as a marginal case. This problem has been discussed first by Griffiths [19] for dipoles on a lattice and then, more generally, also for continuum fluids
in [20, 21]. What has been shown by these authors is that—in the presence of a stabilizing short distance regularization e.g. via a lattice or a hard-sphere potential—the long-range and partially attractive dipolar interaction gives rise to an extensive free energy $F(T, V, N) = N \cdot \tilde{f}(T, v = V/N)$ with a free energy per particle $\tilde{f}$ which is finite and independent of the boundary conditions provided the expectation value $\mu_m \sum_j \langle \tilde{d}_j \rangle = 0$ of the total dipole moment vanishes. The proof thus covers standard polar fluids like water whose electric dipoles point in an arbitrary direction or magnetic systems in the absence of an ordering external field. The situation is fundamentally different, however, for dipolar gases where the dipole directions $\tilde{d}_j$ are all identical and their sum $\sum_j \tilde{d}_j = N\hat{z}$ is extensive. As a consequence, even if the issue of the existence of a well defined limit $N \to \infty$ is ignored for typical particle numbers $N \approx 10^4$ [6], standard results like the Gibbs-Duhem relation $\mu = \tilde{f} + pv$ are not expected to hold and, moreover, thermodynamic quantities will become shape-dependent. These points will be addressed in some detail in the following.

## 2.1. Magneto-chemical potential and demagnetization tensor

A well known feature in the thermodynamics of particles with charge $q$ is that the full chemical potential $\mu(x) + q\phi(x)$ whose gradient determines the current in an inhomogenous situation contains the interaction energy with the local electrostatic potential $\phi(x)$ in addition to the contribution $\mu(x)$ associated with other interactions. In a completely analogous manner, for particles with permanent magnetic dipole moments pointing along $\hat{z}$, the relevant magnetochemical potential contains a contribution $\mu_m B_z(x)$, where $B_z(x)$ is the exact magnetic field at the position of the particle. In a formal manner, this can be derived by noting that the local chemical potential

$$ \mu(x) = \frac{\delta F[\rho^{(1)}]}{\delta \rho^{(1)}(x)} = \mu_{sr}(x) - \mu_m B_z(x) \quad \text{with} \quad -\mu_m B_z(x) = \int_{x'} V_{dd}(x-x') \rho^{(1)}(x') g^{(2)}(x, x') \quad (5) $$

is obtained by a functional derivative of the free energy (or simply $\langle \hat{H} \rangle$ at $T = 0$) with respect to the one-particle density [9]. Here, the explicit expression for the exact local magnetic field due to all other dipoles follows from differentiation of the associated interaction energy (3) while the contribution $\mu_{sr}(x)$ arises from the short-range interactions, including the change in kinetic energy due to $\hat{H}_{dd}$. As a consequence of the long-range nature of the dipolar interaction, $\mu(x)$ is not—as usual—fully determined by temperature and the local one-particle density $\rho^{(1)}(x)$ but depends on the overall shape of the sample. This is evident already at the mean-field level $g_{mf}^{(2)} = 1$ where the calculation of $B_z(x)$ reduces to a problem in classical magnetostatics with given magnetization $M(x) = \mu_m \rho^{(1)}(x) \hat{z}$. The appearance of shape-dependent thermodynamics is well known in this context and it is a standard textbook problem to show e.g. that the internal magnetic field in a sphere with constant magnetization is $B_{in} = (2/3)M$ [22]. An exact solution is available also for more general geometries provided the magnetization current density $j_M(x) = \text{rot } M(x)$ may be approximated by a pure surface current. This applies naturally in solid ferromagnets with a sharp boundary and it carries over to the inner region of trapped or self-bound dipolar gases where the spatial dependence of $M(x)$ is negligible. The relevant demagnetization field $H = -N_M$ in $B = H + M$ is then determined by a shape-dependent demagnetization tensor $N$ which obeys $\text{Tr } N = 1$ [23]. Specifically, for a spheroid with aspect ratio $\kappa = R_x/R_z$ and $M$ pointing along the $z$-axis, the tensor $N$ is diagonal. Symmetry fixes the relevant eigenvalue $n^{(z)}$ to be equal to 1/3 in the case of a sphere while $n^{(z)}$ approaches zero for a long cylindrical spheroid with $\kappa \to 0$, where the demagnetization effect is negligible. Now, there is a subtle point which needs to be taken into account in an application of these results to the interaction induced field $B_z(x)$ that enters the local chemical potential (5). This has to do with the fact that the contribution $(2/3)M$ to the inner field of a spherical configuration comes entirely from the singular self-interaction term due to an effective delta function in the magnetic field of a point dipole [22]. In contrast to the corresponding
zero-range contribution to the interaction between an electron and a nuclear magnetic moment which is the origin of hyperfine shifts in s-states [24], it is physically reasonable to assume that there is no direct overlap between the partially filled electronic shells in rare-earth atoms. The effective delta function contribution must therefore be subtracted, which leads to $\mathbf{B}_{dd} = \mathbf{H} + \mathbf{M}/3$ for the inner field due to the dipolar interaction. Using $\mu_m^2 = 3g_{dd}$, the resulting contribution to the chemical potential at the mean-field level

$$ \mu_m B_z^{(\text{mf})}(x) = g_{dd}(1 - 3n^{(z)})\rho^{(1)}(x) = g_{dd}f(\kappa)\rho^{(1)}(x) \quad (6) $$

displays a simple shape-dependence via the demagnetization tensor element $n^{(z)}$. The explicit result for $n^{(z)}$ is derived e.g. in [23] in the context of the equivalent problem of depolarization fields of perfect conductors. For an oblate spheroid with eccentricity $e = \sqrt{\kappa^2 - 1}$ it is given by

$$ n^{(z)}(e) = (1 + e^2) \frac{e - \arctan(e)}{e^3} \rightarrow (1 - \tilde{e}^2) \frac{\arctanh(\tilde{e}) - \tilde{e}}{\tilde{e}^3} = n^{(z)}(\tilde{e}) \in (0, 1/3] \quad (7) $$

which also determines the demagnetization factor $n^{(z)}(\tilde{e})$ in a prolate situation with eccentricity $\tilde{e} = \sqrt{1 - \kappa^2}$ by straightforward analytic continuation $e \rightarrow i\tilde{e}$. The expression (6) reproduces a result that has been derived early on in the field, based on a Gaussian Ansatz for the density profile of a dipolar gas in a harmonic trap [25, 26]. It has been shown to hold also in the exact solution of the Gross-Pitaevskii equation in the Thomas-Fermi limit [27] but—somewhat surprisingly—the simple connection $f(\kappa) = 1 - 3n^{(z)}$ of the shape function $f(\kappa)$ with the standard demagnetization tensor of magnetostatics does not seem to have been realized in the literature so far. Note that due to $n^{(z)}(\tilde{e}) < 1/3$ in a prolate configuration $\kappa < 1$, the internal magnetic field $B_z(x)$ in the bulk of the atomic cloud is pointing along the positive z-direction. This leads to a lowering of the chemical potential with respect to the value $\mu_{sr}(x)$ determined by the short-range interactions, reflecting the dominance of attractively interacting dipoles in such a configuration.

An obvious question is to which extent corrections beyond mean-field which are contained in the deviation of the pair distribution function $g^{(2)}(x, x')$ in Equation (5) from the trivial limit one, do affect the result (6). It turns out that there is indeed an important qualitative change which arises from the anisotropy of the dipolar interaction. For a rotation invariant $V_{sr}(r_{12})$, the two-particle density $\rho^{(2)}(X, x_{12})$ only depends on the magnitude $r_{12} = |x_{12}|$ of the separation vector. This is no longer true for anisotropic interactions, where also the orientation of $x_{12}$ with respect to the dipole direction enters. The associated physics has been described in quantitative terms in a landmark paper by Wertheim on a classical fluid of hard spheres with point dipoles whose directions $\hat{d}_j$ fluctuate statistically [28]. Following his notation$^2$, the pair correlation function of a homogeneous dipolar fluid with rotation invariance around the z-axis can be decomposed in the form

$$ h^{(2)}(x_{12}) = g^{(2)}(x_{12}) - 1 = h_S(r_{12}) + h_D(r_{12}) \cdot 2P_2(\hat{z} \cdot \hat{x}_{12}) + \dots \quad (8) $$

Here, the dots indicate terms involving Legendre polynomials of higher order which are neglected since their contribution to $\mu_m B_x(x)$ contains angular integrals $\int_{\Omega} P_2 P_l$ which tend to cancel (only even $l$ appear as long as $z \rightarrow -z$ is a symmetry). The two scalar functions $h_S(r_{12})$ and $h_D(r_{12})$ describe quite different physics and turn out to have opposite sign. Specifically, $h_S(r_{12}) < 0$ accounts for an effective short-range repulsion, reducing the probability to find two particles at separation $r_{12}$ below that of an uncorrelated state. In turn, at fixed $r_{12}$, a positive $h_D(r_{12}) > 0$ favors attractive head-to-tail configurations with $(\hat{z} \cdot \hat{x}_{12})^2 > 1/3$. In explicit form, the different behavior of separation and angular orientation correlations with respect to the z-axis may be derived for a uniform Bose gas in the limit of weak dipolar interactions $\epsilon_{dd} = a_{dd}/a \ll 1$ by

<footnote>$^2$In [28], an additional contribution $h_{\Lambda}(r_{12})\hat{d}_1 \cdot \hat{d}_2$ appears which couples to the relative orientation of the dipoles. In the present case with $\hat{d}_1 \cdot \hat{d}_2 \equiv 1$, this contribution can be incorporated into the rotation invariant function $h_S(r_{12})$.</footnote>
expanding the static structure factor within the Bogoliubov approximation in the regime of small momenta

$$S(q) = 1 + nh^{(2)}(q) \underset{\text{Bog.}}{\overset{\text{Bog.}}{\longrightarrow}} [1 + 2nV(q)/\epsilon_q]^{-1/2} \underset{|q|\xi, \epsilon_{dd}\ll 1}{\overset{|q|\xi}{\longrightarrow}} \frac{|q|\xi}{\sqrt{2}} [1 - \epsilon_{dd}P_2(\hat{z} \cdot \hat{q}) + \cdots]. \quad (9)$$

Here, $V(q) = g + V_{dd}(q)$ is the Fourier transform of the interaction, $\xi = 1/\sqrt{8\pi na}$ the healing length associated with its short-range part and $\epsilon_q = \hbar^2 q^2/(2m)$ the free particle dispersion. The expansion (9) of the static structure factor at low momenta precisely matches the corresponding one in real space of Equation (8). Moreover, the non-analytic behavior in $q$ implies that the associated pair correlation functions $nh_S(r) \to -\xi/(\sqrt{2\pi^2 r^4})$ and $nh_D(r) \to +\epsilon_{dd} \cdot \sqrt{2\xi/(\pi^2 r^4)}$ follow a power law decay at $T = 0$ and have opposite sign. A further point which distinguishes the correlations of the magnitude of separation vectors from those which describe their orientation with respect to the dipole direction shows up at finite temperature: while $h_S(r)$ decays to zero exponentially beyond the thermal wavelength $\lambda_T$, the angle dependent contribution

$$h_D(r)|_{T \neq 0} \to h_\infty/r^3 \quad \text{with} \quad nh_\infty = \frac{3}{4\pi} \frac{(\epsilon(0)-1)^2}{\epsilon(0)} \frac{k_B T}{g_{dd} n} \quad (10)$$

retains a power law decay, following that of the dipolar interaction. The asymptotic behavior (10) is an exact result first derived by Wertheim [28] for classical fluids with randomly oriented electric dipoles and $\epsilon(0)$ the associated static dielectric constant. It is consistent with a correlation inequality at finite temperature which implies that correlations cannot decay faster than the interactions. Using the connection $S(q) \to k_B T/mc^2(\hat{q})$ between the static structure factor at low momenta and the angle dependent sound velocity $c(\hat{q})$ of dipolar gases [29], it is straightforward to show that—to linear order in $\epsilon(0) - 1 \to \epsilon_{dd} \ll 1$, where the homogeneous system with average density $n$ is stable—the result (10) continues to hold for quantum fluids of magnetic dipoles which are fully polarized.

An important consequence of the results above for the thermodynamics of dipolar gases shows up by considering nearly spherical configurations, where mean-field predicts that the dipolar interactions average to zero. Obviously, this result remains valid more generally for an arbitrary pair correlation function which is rotation invariant. Including the angular dependent contribution $h_D(r) \neq 0$, however, there is a non-vanishing magnetic field

$$\mu_m B_z(x)|_{K=1} = (12/5)g_{dd} \rho^{(1)}(x) \cdot K \quad \text{with} \quad K = \int_0^\infty dr \frac{h_D(r)}{r} > 0 \quad (11)$$

near the center which points along the positive $z$-direction. The dimensionless parameter $K$, which plays a central role in a microscopic theory of the dielectric constant in classical polar fluids$^3$, apparently requires a cutoff at short distances. As mentioned above, a description of dipolar gases beyond mean-field must therefore take into account the finite range of interactions which, in practice, is of the order of the short-range scattering length $a$. An experimental estimate of the parameter $K$ may be obtained by observing the finite critical strength $\epsilon_{dd}^{(c)}(K=1) = 5/(12K)$ of the dipolar interaction where even a spherical cloud becomes unstable because the chemical potential $\mu(x) \approx g\rho^{(1)}(x) - \mu_m B_z(x)$ turns negative at leading order.

## 2.2. Pressure and uniaxial tension from dilatations

In the discussion so far, the consequences of the anisotropy and long-range nature of the dipolar interactions have been addressed by assuming a given shape, characterized e.g. by the

<footnote>$^3$Combining point dipoles with hard spheres, the dielectric constant is a function only of the volume fraction $\eta = (\pi/6)n\sigma^3$ and the ratio $y = g_{dd}n/(3k_B T)$. Remarkably, this extremely oversimplified model reproduces quite well the observed value $\epsilon(0) = 78.4$ for water under normal conditions, where $\eta \approx 0.42$ and $y \approx 3.35$, see Figure 11.5 in [9].</footnote>
aspect ratio $\kappa$. In reality, the specific shape of the atomic cloud is not an independent variable but is determined by the interaction strength and the external confinement, if present. In thermal equilibrium, it follows from the balance of forces which—in quite general form—may be expressed in terms of the associated momentum current tensor $\Pi(x)$ by

$$ \nabla \Pi(x) \underset{\text{equ.}}{\longrightarrow} \begin{pmatrix} \partial_x \Pi^{xx} \\ \partial_y \Pi^{yy} \\ \partial_z \Pi^{zz} \end{pmatrix} = \begin{cases} 0 & \text{self-bound} \\ -\rho_{\phi}^{(1)}(x) \operatorname{grad} \phi(x) & \text{trapped.} \end{cases} \qquad (12) $$

Here, $\phi(x)$ is the external trap potential and $\rho_{\phi}^{(1)}(x)$ the resulting equilibrium value of the inhomogeneous one-particle density. Moreover, we have used the fact that in thermal equilibrium the tensor $\Pi(x)$ has no off-diagonal elements as would appear e.g. in the presence of a finite Hall viscosity. Now, for rotation invariant interactions, all diagonal elements of $\Pi(x)$ are equal and they define a local pressure $p(x)$. The left hand side of Equation (12) is then just the gradient of pressure in a standard hydrostatic equilibrium. For a harmonic confinement, where the external contribution to the total Hamiltonian

$$ \hat{H}_{\omega} = \hat{H} + \frac{1}{2} \operatorname{Tr}(\underline{\omega}^2 \hat{\theta}) \quad \text{with} \quad (\hat{\theta})^{ij} = \int_x x^i x^j \hat{\rho}(x) \qquad (13) $$

can be expressed in terms of a tensor $\underline{\omega}^2$ of the trap frequencies and the inertia tensor $\hat{\theta}$ defined via the mass density operator $\hat{\rho}(x)$, the resulting cloud shape follows from the exact relation

$$ \langle \hat{\theta} \rangle_{\text{eq}} = (\underline{\omega}^2)^{-1} \int_x p(x) = E_{\text{rel}}(\underline{\omega}^2)^{-1} \quad \text{for rotation invariant } V(x_{12}) = V(r_{12}) \qquad (14) $$

obtained by taking the scalar product of Equation (12) with $\mathbf{x}$ and integrating over all space. Note that the relation (14) provides only an implicit result for the cloud shape as characterized by its inertia tensor since the release energy $E_{\text{rel}} = \int_x p(x)$ is itself shape-dependent. An exact consequence of Equation (14), however, is that the aspect ratio which—irrespective of the detailed density profile—is defined by $\kappa = |\theta^{xx}/\theta^{zz}|^{1/2}$ coincides with the ratio $\lambda = \omega_z/\omega_x$ of the trap frequencies (we assume $\omega_x = \omega_y \neq \omega_z$ in general, i.e. there is rotation invariance around the $z$-axis).

The anisotropy of the dipolar interaction invalidates the simple relation (14) and leads to shapes with $\kappa < \lambda$ which are elongated along the $z$-direction since the energy is lowered if the dipoles tend to arrange in a head-to-tail configuration. Dipolar gases are thus uniaxial quantum fluids which—in quite general terms—may be defined by $\Pi^{xx} = \Pi^{yy} \neq \Pi^{zz}$ in equilibrium. To characterize their size and shape one needs two instead of just a single scalar function. Specifically, we will use $\operatorname{Tr} \Pi$ and the anisotropy $h = \Pi^{xx} - \Pi^{zz}$, which acts like a uniaxial tension tending to stretch the system along $z$. As will be discussed below, the system length $L_z$ will then appear as an additional thermodynamic variable which is conjugate to the field $h$. On a microscopic level, an explicit expression for both $\operatorname{Tr} \Pi(x)$ and $h(x)$ may be derived by extending the results of Martin and Schwinger [30] for the momentum current tensor in a non-relativistic quantum many-body problem with rotation invariant interactions. Quite generally, pressure— as the field conjugate to a local expansion—is defined by the change in energy density under a uniform isotropic dilatation $x \to bx$ with scale factor $b \to 1$. This only fixes $\operatorname{Tr} \Pi(x)$, however. In order to specify the two independent components $\Pi^{xx}$ and $\Pi^{zz}$ it is necessary to consider, in addition, the volume preserving dilatation $x \to \theta[x]$ where the components $z \to \theta z$ and $(x, y) \to (x, y)/\sqrt{\theta}$ are scaled differently. The full momentum current tensor of uniaxial quantum fluids is thus obtained from the microscopic Hamiltonian density $\hat{\mathcal{H}}(x)$ in $\hat{H} = \int_x \hat{\mathcal{H}}(x)$ by evaluating the response to the two independent dilatations (note the different signs)

$$ \operatorname{Tr} \Pi(x) = 2\Pi^{xx}(x) + \Pi^{zz}(x) = - \left. \frac{\partial \langle \hat{\mathcal{H}}(bx) \rangle_{\text{eq}}}{\partial b} \right|_{b=1} \quad \text{and} \quad h(x) = \left. \frac{\partial \langle \hat{\mathcal{H}}(\theta|x) \rangle_{\text{eq}}}{\partial \theta} \right|_{\theta=1}. \qquad (15) $$
In explicit form, the anomalous field $h(x)$ receives contributions only from the kinetic energy and the dipolar interaction. They can be expressed in terms of the field operators $\hat{\psi}$ as

$$h(x) = \frac{\hbar^2}{m} \left[ \frac{1}{2} \langle \nabla_{\perp} \hat{\psi}^{\dagger} \cdot \nabla_{\perp} \hat{\psi} \rangle_{\text{eq}}(x) - \langle \partial_z \hat{\psi}^{\dagger} \partial_z \hat{\psi} \rangle_{\text{eq}}(x) \right] + \partial_{\theta} \mathcal{D}(\theta[x])|_{\theta=1}. \quad (16)$$

Here, $\nabla_{\perp}$ denotes the gradient with respect to the radial coordinates $(x, y)$ and $\mathcal{D}(x)$ is a dipolar analog of the contact density $\mathcal{C}_2(x)$ for zero-range interactions. It is defined by the expectation value of the dipolar interaction energy density

$$\mathcal{D}(x) = \langle \mathcal{H}_{\text{dd}}(x) \rangle_{\text{eq}} = -\frac{1}{2} \langle \hat{M}_z(x) \hat{B}_z(x) \rangle_{\text{eq}} \quad (17)$$

which involves the product of the operators for the magnetization $\hat{M}_z(x)$ and the field $\hat{B}_z(x)$ from all other dipoles at position $x$. Due to the long-range nature of the interaction, $\mathcal{D}(x)$ is in fact non-local and shape-dependent, a point which will be discussed further below. This is very different from the energy density $\varepsilon_{\text{sr}}(x)$ and the associated contribution [7, 31]

$$p_{\text{sr}}(x) = \frac{2}{3} \varepsilon_{\text{sr}}(x) + \frac{\hbar^2}{24\pi m a} \mathcal{C}_2(x) \xrightarrow{\text{mf}} \frac{g}{2} [\rho^{(1)}(x)]^2 + \dots \quad (18)$$

to the pressure arising from short-range interactions which—at the mean-field level—is quadratic in the local one-particle density $\rho^{(1)}(x)$ and independent of the shape. Using the microscopic definition (15), the full momentum current tensor of dipolar gases in turn contains the shape-dependent contributions $\mathcal{D}(x)$ and $\partial_{\theta}\mathcal{D}(\theta[x])|_{\theta=1}$ according to

$$\Pi^{xx}(x) = p_{\text{sr}}(x) + \mathcal{D}(x) + \frac{1}{3} h(x). \quad (19)$$

The formal relations above turn out to lead to observable and intuitively plausible consequences for the density profile of harmonically trapped dipolar gases. In a situation with $\omega_x = \omega_y$, where one has rotation invariance around the $z$-axis, a radial integration of the $(x, y)$-components of the equilibrium condition (12) shows that the integrated column density

$$\tilde{n}_{\phi}(z) = \int dx dy \rho_{\phi}^{(1)}(x, y, z) = \frac{2\pi}{m\omega_x^2} \Pi^{xx}(0, 0, z) \quad (20)$$

at a given $z$ is a direct measure of the $xx$-component of the momentum current tensor along the center of the atomic cloud. This relation is equivalent to one derived by Ho and Zhou for rotation invariant interactions in the context of imbalanced Fermi gases [32]. A different result appears, however, for the $z$-derivative of the density profile which follows from the third component of the momentum balance (12). The fact that $\Pi^{xx}$ and $\Pi^{zz}$ no longer coincide in the uniaxial case leads to a violation of the Gibbs-Duhem relation which shows up as a non-vanishing contribution on the right-hand side of the equation

$$\frac{d\tilde{n}_{\phi}(z)}{dz} + 2\pi \frac{\omega_z^2}{\omega_x^2} z \rho_{\phi}^{(1)}(0, 0, z) = \frac{2\pi}{m\omega_x^2} \partial_z h(0, 0, z) > 0. \quad (21)$$

Since the anomalous field $h(x)$ is negative and decays along the $z$-direction away from the center, the profile $\tilde{n}_{\phi}(z)$ falls off more slowly than what is expected for rotation invariant interactions. The uniaxial nature of dipolar gases and the associated field $h(x)$ can thus be inferred directly from in-situ density profiles.

## 2.3. Thermodynamic relations and the tensor virial theorem

The appearance of a finite uniaxial tension in dipolar gases leads to a fundamental change in thermodynamic relations. This turns out to provide an understanding of the shape-dependence mentioned above and also the peculiar form of self-bound droplets that will be discussed in the following subsection. In order to derive the extension to uniaxial fluids of the standard relation $dF = -p dV$ for the change in free energy associated with a shape deformation at a fixed value of
temperature and particle number it is convenient to introduce a local deformation tensor $\underline{u}(x)$ as in solids. Quite generally, $\underline{u}(x)$ is the conjugate field to the momentum current (or stress) tensor $\underline{\Pi}(x)$ [10] which, by the definition of an uniaxial fluid in thermal equilibrium, is diagonal with just two independent components. As a consequence, the associated differential of the free energy

$$dF[\text{Tr} \underline{u}, \underline{u}_{zz}]|_{T,N} = - \int_x \text{Tr} [\underline{\Pi}(x) \underline{d}\underline{u}(x)] \xrightarrow{\text{fluid}} \int_x [-\underline{\Pi}^{xx}(x) \text{d}\text{Tr} \underline{u}(x) + h(x) \text{d}u_{zz}(x)] \quad (22)$$

depends on the two dimensionless variables $\text{Tr} \underline{u}$ and $\underline{u}_{zz}$. In physical terms, they describe a local change $\text{d}\text{Tr} \underline{u}(x) = \text{d}\underline{v}(x)/\underline{v}$ in volume per particle at a fixed extension along the $z$-direction or a change of length $\text{d}u_{zz} = \text{d}L_z/L_z$ at fixed volume$^4$. The relation (22) shows that for a uniaxial fluid $\underline{\Pi}^{xx} \to p$ can be identified with the pressure as the variable conjugate to volume at fixed length. In addition, the anomalous field $h(x)$ gives rise to a line tension contribution in the Gibbs relation

$$dF(T, V, L_z, N) = -S dT - p dV + h_1 dL_z + \mu dN \quad \text{with} \quad h_1 L_z = \int_x h(x), \quad (23)$$

favoring an increase of length $L_z$ at given volume to lower the free energy ($h_1 < 0$). As will be shown below, the fact that $L_z$ appears as an additional extensive thermodynamic variable leads to a violation of the Gibbs-Duhem relation. Typically, such violations become negligible in the thermodynamic limit, which holds e.g. for a finite surface tension $\sigma$ in a self-bound liquid state described by $dF = \sigma dA$ at fixed $T, N, V$. For short-range interactions, this contribution scales with the surface area $A \sim N^{2/3}$ and thus eventually becomes subdominant compared with the extensive terms. The situation is different here, however, due to the long-range nature of the dipolar interaction. Indeed, since the interaction contribution $\partial_{\theta} \mathcal{D}(\theta[x])|_{\theta=1}$ in Equation (16) to the underlying variable $h(x)$ is intensive, $dF = h_1 dL_z \sim N$ remains extensive for an arbitrary large system. In physical terms, this can be understood from the fact that the dipolar contact density $\mathcal{D}(x)$ in Equation (17) arises from the magnetic field which, as discussed in Section 2.1, is effectively generated by surface currents. The demagnetization field $\mathbf{H}(x) = -\text{grad} \phi_M(x)$ is thus non-local, arising from magnetic surface “charges” with density $\sigma_M = \mathbf{n} \cdot \mathbf{M}$ ($\mathbf{n}$ is the unit normal vector) [22]. As emphasized in [21], the $1/r^2$ decay of the field is precisely cancelled by the increase $\sim r^2$ of the surface area. Both the magnetic field and the associated energy density $\mathcal{D}(x) = \langle \mathcal{H} \mathcal{E}_{dd}(x) \rangle_{\text{eq}}$ are therefore independent of the system size.

The extended Gibbs relation (23) implies that—at fixed temperature $T$—the pressure and chemical potential involve the strain tensor element $u_{zz}$ as a further thermodynamic variable beyond the volume per particle $\underline{v}$. The associated dependence is determined by the Maxwell relations

$$\underline{v} \frac{\partial p}{\partial \underline{v}} = \frac{\partial \mu}{\partial \underline{v}} \quad \text{at fixed } u_{zz}, \quad \frac{\partial \mu}{\partial u_{zz}} = \frac{\partial h}{\partial (1/\underline{v})} \quad \text{and} \quad \frac{\partial p}{\partial u_{zz}} = \underline{v} \frac{\partial h}{\partial \underline{v}}. \quad (24)$$

In an inhomogeneous situation, where the density $\rho^{(1)}(x)$ and the field $h(x)$ are spatially varying, the standard Gibbs-Duhem relation $d\underline{\mu}|_T = \underline{v} dp|_T$ is therefore violated. This leads to an anomalous contribution to the axial derivative of the density profile of a harmonically trapped dipolar gas derived in Equation (21). A further consequence of the dependence of pressure $p(\underline{v}, u_{zz}, T)$ on the uniaxial strain $u_{zz}$ is the angular dependence of the sound velocity in a uniform dipolar gas. In quite general terms, the difference $\Delta f = f - f_{\text{eq}}$ in the free energy density with respect its value $f_{\text{eq}}$ in a homogeneous situation with vanishing deformations $(\text{Tr} \underline{u})_{\text{eq}} = (u_{zz})_{\text{eq}} \equiv 0$ can be written in the form of an elastic free energy of a uniaxial solid [10]

$$\Delta f[\text{Tr} \underline{u}, u_{zz}] = \frac{1}{2}(u_{xx} + u_{yy}, u_{zz}) \begin{pmatrix} K_{11} & K_{13} \\ K_{13} & K_{33} \end{pmatrix} \begin{pmatrix} u_{xx} + u_{yy} \\ u_{zz} \end{pmatrix} = \frac{K_{11}}{2} [\text{Tr} \underline{u} + \tilde{\gamma} u_{zz}]^2 + \frac{\tilde{B}}{2} u_{zz}^2 \quad (25)$$

<footnote>$^4$While a cloud of atoms has no sharp boundary beyond the Thomas-Fermi limit where the kinetic energy is negligible, the associated Thomas-Fermi radius $R_z$ still provides an appropriate measure for the length $L_z \to 2R_z$.</footnote>
which is rotation invariant in the x, y-plane. The condition that Equation (25) describes a fluid rather than a solid requires that the associated tensor of elastic constants has vanishing determinant $\tilde{B}K_{11} = K_{33}K_{11} - K_{13}^2 = 0$. The last term is then absent and there is only a single sound mode<fnref n="5" />. Its velocity depends on the direction with respect to the orientation of the dipoles according to $\tilde{\rho}c^2(\hat{q}) = K_{11} + (K_{33} - K_{11}) \cdot (\hat{z} \cdot \hat{q})^2$ [10, 33]. This angular dependence is a result of the coupling between density and strain described by the dimensionless parameter $\tilde{\gamma}$ in

$$h = \left. \frac{\partial f}{\partial u_{zz}} \right|_{\text{Tr}\underline{u}} \rightarrow \tilde{\gamma}K_{11}(\text{Tr}\underline{u} + \tilde{\gamma}u_{zz}) \quad \text{with } \tilde{\gamma} = \frac{K_{13}}{K_{11}} - 1 \underset{\text{mf}}{\longrightarrow} \left( \frac{1 + 2\epsilon_{\text{dd}}}{1 - \epsilon_{\text{dd}}} \right)^{1/2} - 1. \quad (26)$$

A physical understanding of this coupling is provided by the observation that—according to Equation (25)—the energy cost for a compression $\text{Tr}\underline{u} < 0$ of a uniform dipolar gas can be exactly compensated by an expansion along $z$ with $u_{zz} = -\text{Tr}\underline{u}/\tilde{\gamma} > 0$. The mean-field result for $\tilde{\gamma}$ in terms of the dimensionless strength $\epsilon_{dd}$ of the dipolar interaction follows from the angular dependence $V(q) = g[1 + \epsilon_{dd} \cdot 2P_2(\hat{z} \cdot \hat{q})]$ of the two-body interaction. It indicates that the uniform gas with a homogeneous mass density $\tilde{\rho}$ is stable only as long as $\epsilon_{dd} < 1$, a point which will be discussed in detail in the following subsection.

An exact relation that holds for many-body problems in both classical and quantum physics is the virial theorem. It is a scalar identity which follows by considering the change in energy under a uniform dilatation of the coordinates $(x_1...x_N) \rightarrow (bx_1...bx_N)$ [37]. Now, as discussed in Equation (15) above, uniform dilatations are not sufficient to characterize the thermodynamically relevant set of shape deformations in the case of interactions which are not invariant under rotation. To deal with this more general situation, it is necessary to extend the virial theorem to a tensorial form, as was first derived by Parker [38] for a classical system of point particles in an astrophysical context. The extension to the quantum many-body problem is straightforward in principle. Indeed, the virial theorem in its most general form simply states that the expectation value of the commutator

$$\langle i[\hat{H} + \hat{H}_{\text{ex}}, \hat{D}] \rangle = 0 \underset{\text{harm.}}{\longrightarrow} \langle i[\hat{H}, \hat{D}] \rangle = \langle \underline{\omega}^2 \underline{\theta} \rangle \quad \text{with} \quad (\hat{D})^{ij} = \int_x x^i \hat{g}^j(x) / \hbar \quad (27)$$

of the total Hamiltonian with the dilatation tensor operator $\hat{D}$ vanishes in eigenstates of $\hat{H} + \hat{H}_{\text{ex}}$. The dimensionless dilatation tensor involves the components $\hat{g}^j(x)$ of the momentum density operator and generates dilatations along the $j$-direction proportional to separate components $x^i$ of the coordinates [39]. In the particular case of a harmonic confinement, the contribution of $\hat{H}_{\text{ex}}$ is reduced to the expectation value of the tensor contraction $\underline{\omega}^2 \underline{\theta}$. Its trace yields twice the trap energy $E_{\text{tr}}$ according to Equation (13). For rotation invariant interactions, the equilibrium expectation value $\langle i[\hat{H}, \hat{D}] \rangle_{\text{eq}} = \frac{1}{\hbar} \int_x p(x)$ is proportional to the unit tensor. The tensorial form of the virial theorem then reduces to the standard one obtained from a uniform dilatation $x \rightarrow bx$. In particular, the moment of inertia tensor is given by Equation (14). This is violated in the uniaxial situation, where the aspect ratio $\kappa = [\theta^{xx}/\theta^{zz}]^{1/2}$ no longer coincides with the ratio $\lambda = \omega_z/\omega_x$ of external trap frequencies. The difference is determined by the exact relation

$$-\int_x h(x) = \omega_z^2 \theta^{zz} - \omega_x^2 \theta^{xx} = 2E_{\text{tr}} \cdot \frac{\lambda^2 - \kappa^2}{\lambda^2 + 2\kappa^2} > 0 \quad (28)$$

which follows from the tensor virial theorem (27) in an equilibrium state by considering the volume preserving dilatation $x \rightarrow \theta[x]$ introduced in Equation (15). A detailed discussion of

<footnote n="5">⁵For finite values of the parameter $\tilde{B}$, the free energy (25) describes a smectic-A phase with broken translation invariance along the $z$-direction [10]. It exhibits two independent sound modes, both in the case of a classical fluid [33] and also—in the quite different form of a wave-like propagation of defects in supersolids predicted by Andreev and Lifshitz [34]—in the superfluid [35]. These two modes have very recently been observed in a 2D smectic phase of a Bose Einstein condensate in the presence of a periodic in time modulation of the scattering length [36], allowing to determine both the superfluid fraction $f_s$ and the effective layer compression modulus $\tilde{B}$ in quantitative terms.</footnote>
the connection between $\kappa$ and $\lambda$ as a function of the dipolar interaction strength $\varepsilon_{dd}$ has been given by Eberlein et al. [27] within the Thomas-Fermi approximation, where the density profile $\rho_{\phi}^{(1)}(x)$ has a sharp boundary. In particular, they have shown that, in the presence of a strong external confinement $\lambda \gg 1$ along the z-direction, a dipolar gas remains stable even in the limit $\varepsilon_{dd} \gg 1$ of negligible short-range repulsion. This has been verified experimentally by Koch et al. [40] for particle numbers $N \approx 10^4$. Effects beyond mean-field, which are contained in Equation (28), therefore do not give rise to qualitative changes in a tightly confined situation. As will be shown below, the situation is quite different in self-bound droplets whose form is completely determined by such corrections. A second independent relation which follows from the tensorial virial theorem is connected with the uniform dilatation $x \to bx$ generated by the scalar operator $\hat{D} = \text{Tr} \hat{D}$. It can be written in terms of the local contact densities $\mathscr{C}_2(x)$ and $\mathscr{D}(x)$ associated with the short-range and dipolar interactions in the form

$$ \langle i[\hat{H}, \hat{D}] \rangle = \langle \hat{H} \rangle + E_{\text{tr}} + \frac{\hbar^2}{16\pi m a} \int_x \mathscr{C}_2(x) + \frac{1}{2} \int_x \mathscr{D}(x) = 2E_{\text{tr}}. \quad (29) $$

A similar relation was derived in [41] for dipolar gases in 2D, however there are two important differences: since the dipolar interaction in 2D also leads to s-wave scattering, the two-body contact density $\mathscr{C}_2(x)$ is defined by the adiabatic derivative of the energy with respect to the full scattering length of the combined interaction $V_{sr} + V_{dd}$. As a consequence, the exact expression for the dipolar contact density $\mathscr{D}(x)$ requires a subtraction in the pair distribution function at short distances associated with the contribution of $V_{sr}(r_{12})$ [41]. A second point is that the angular independent dipolar interactions $\sim \mu_m^2/r_{12}^3$ in 2D are effectively of a short-range nature. In contrast to the situation here, the associated contact density $\mathscr{D}(x)$ is then shape-independent.

## 2.4. Self-bound droplets and the LHY correction

One of the major new developments associated with ultracold dipolar gases was the surprising observation of self-bound droplets by Ferrier-Barbut et al. [2]. These droplets are very different from what is expected for droplets of a standard liquid, whose equilibrium configuration is spherical due to the presence of a finite surface tension $\sigma$. Indeed, as will be shown below, the combination of anisotropy and the dominantly attractive nature of dipolar interactions in a non-confined situation leads to droplets whose aspect ratio $\kappa(N) \sim N^{-1/4}$ approaches zero for large particle numbers$^6$.

Based on the quite general hydrostatic equilibrium condition (12), the existence of stable self-bound droplets requires that the divergence $\nabla \Pi(x) = 0$ of the momentum current tensor vanishes. In a uniaxial fluid with long-range and angle-dependent attractive interactions, this condition is obeyed only if both $\text{Tr} \Pi(x)$ and $h(x)$ are zero. In contrast to the situation in stars, where a single scalar equation for the balance of the radial component of the internal and gravitational pressure is sufficient to determine the equilibrium shape (see e.g. [43] for an instructive discussion), two separate conditions need to be fulfilled here. In order to formulate these in a physically intuitive manner, it is convenient to rewrite the virial theorem (29) in a form where the external trap is eliminated

$$ \langle i[\hat{H}, \hat{D}] \rangle = \int_x \left[ 2\varepsilon_{\text{sr}}(x) + \frac{\hbar^2}{8\pi m a} \mathscr{C}_2(x) + 3\mathscr{D}(x) \right] = 3 \int_x [p_{\text{sr}}(x) + \mathscr{D}(x)] = \int_x \text{Tr} \Pi(x). \quad (30) $$

<footnote>$^6$This scaling law has been seen originally in numerical simulations by Baillie et al. [42] and has recently been derived in analytical form by Dalibard [8], based on a different but physically equivalent line of arguments.</footnote>
Here, following the Tan pressure relation in Equation (18), we have introduced the local pressure $p_{sr}(x)$ which arises from the short-range interactions$^7$. Beyond the associated energy density $\epsilon_{sr}(x)$, it involves the two-body contact density $C_2(x)$ which is defined by the adiabatic derivative of $\langle \mathcal{H}(x) \rangle_{eq}$ with respect to the inverse scattering length $1/a$ [16,31]. Equation (30) shows that the condition of a vanishing $\mathrm{Tr} \, \Pi(x) \equiv 0$ is equivalent to a local pressure balance $p_{sr}(x)+\mathcal{D}(x) = 0$, with $\mathcal{D}(x)$ playing the role of the pressure associated with the dipolar interactions. Prior to a detailed discussion of this exact and quite intuitive relation a number of comments should be added:

* The fact that the pressure arising from dipolar interactions is equal to the associated energy density $\mathcal{D}(x)$ is a consequence of the assumption that the interaction $V_{dd}(x_{12}) \sim \mu_m^2/r_{12}^3$ follows an inverse cube power law at all distances.

* Despite the formal notation, $\mathcal{D}(x)$ is not a local variable but depends on the overall shape of the droplet, as becomes evident in the mean-field result (31) below.

* The separation of pressure into a short-range and a dipolar contribution follows from the additivity $V = V_{sr} + V_{dd}$ of the interactions in the underlying Hamiltonian. It does not imply that pressure is simply additive, however, which is true only in mean-field. In general, the contributions $p_{sr}(x)$ and $\mathcal{D}(x)$ will depend on both $g$ and $g_{dd}$.

* The identification of $\mathrm{Tr} \, \Pi(x)$ with $3[p_{sr}(x) + \mathcal{D}(x)]$ on a local level neglects a possible divergence which integrates to zero in the virial equation. Indeed, quite generally, the momentum current tensor is not unique and the equilibrium condition (12) remains unchanged under $(\Pi)^{ij} \rightarrow (\Pi)^{ij} + (\partial^i \partial^j - \delta^{ij} \nabla^2)\Phi$ with an arbitrary scalar function $\Phi(x)$ which vanishes at infinity [39]. This freedom is removed by the concrete choice (31) below for $\mathcal{D}(x)$ in the central region of the self-bound droplet.

Explicit results for both the dipolar contact density $\mathcal{D}(x)$ and the local anomalous field $h(x)$ are available at the mean-field level. Specifically, we use Equation (16) which connects the interaction contribution to $h(x)$ to the derivative of the dipolar contact density (17) with respect to the parameter $\theta$. Within mean-field, this yields

$$ \mathcal{D}^{(\mathrm{mf})}(x) = -\frac{1}{2} \langle \hat{M}_z(x) \rangle B_z^{(\mathrm{mf})}(x) = -\frac{1}{2} g_{dd} [\rho^{(1)}(x)]^2 f(\kappa) \rightarrow h^{(\mathrm{mf})}(x) = -\frac{3}{4} g_{dd} [\rho^{(1)}(x)]^2 b(\kappa) \quad (31) $$

where we have used the result (6) for the magnetic field and the fact that the aspect ratio changes according to $d\kappa = -(3/2)\kappa d\theta$ under $x \rightarrow \theta[x]$. The positive function

$$ b(\kappa) = -\kappa \frac{d f(\kappa)}{d\kappa} \rightarrow \begin{cases} 6\kappa^2 \ln(2/\kappa) & \text{for } \kappa \ll 1 \\ 3\pi/(2\kappa) & \text{for } \kappa \gg 1 \end{cases} \quad (32) $$

is determined by Equation (7) for the demagnetization tensor element $n^{(z)}$. It decays to zero for both strongly prolate or oblate configurations $\kappa \ll 1$ or $\kappa \gg 1$ and exhibits a broad maximum $b(\kappa = 1.563...)=0.852...$ slightly above the value $b(\kappa = 1) = 4/5$ in the spherical limit. Now, the combination of Equation (31) with the leading order result (18) for the pressure due to short-range interactions shows that the two separate equilibrium conditions $p_{sr}(x) + \mathcal{D}(x) = 0$ and $h(x) = 0$ for a self-bound droplet cannot be obeyed within mean-field where both $p_{sr}(x) + \mathcal{D}(x)$ and $h(x)$ are negative for $\epsilon_{dd} f(\kappa) > 1$ or arbitrary values of $\kappa$, respectively. A resolution of the first problem, already proposed in the original publication on self-bound droplets [2], is provided by adding a contribution to the short-range pressure which scales with a higher power in density than the quadratic behavior stated as the low-density limit of the exact expression in Equation (18). On a purely empirical level, a self-bound liquid can then be stabilized against

<footnote>$^7$In the application of the Tan relations to Bosons, we consistently neglect three-body correlations which are connected with the Efimov effect and the underlying dependence of $\langle \mathcal{H}(x) \rangle_{eq}$ on the short-range parameter $\kappa_*$ [7, 16, 31].</footnote>
collapse by extending the mean-field effective potential by a local term proportional to $|\psi(x)|^2$ with $k > 2$ in the form

$$V_{\text{eff}}[\psi] = -\mu|\psi(x)|^2 + \frac{g}{2}|\psi(x)|^4 + \frac{1}{2}|\psi(x)|^2 \int_{x'} V_{\text{dd}}(x-x')|\psi(x')|^2 + \frac{\lambda_k}{k}|\psi(x)|^2 k \quad (?) \quad (33)$$

where $\lambda_k > 0$ is assumed to be positive. From a microscopic point of view, the true effective potential must, of course, be derived from the underlying Hamiltonian. In fact, due to the anisotropy of the dipolar interaction, the simple form (33) where the stabilizing contribution is rotation invariant, is not expected to be valid. In principle, a microscopic derivation of $V_{\text{eff}}[\psi]$ is possible by expressing the exact partition function in terms of a coherent state functional integral [44]. The first three terms in Equation (33) then arise from minimizing the bare associated action. At this mean-field (or tree-) level, only two-body interactions enter. It provides a proper description of the zero density limit near the vacuum state $\psi(x) \equiv 0$, where two-body scattering amplitudes are sufficient to deal with the many-body problem at finite density. For self-bound droplets, which arise from the vacuum via a first-order transition at a critical value $\mu_c < 0$ of the chemical potential where the density jumps from zero to a finite value $\bar{n} \neq 0$, this is no longer valid. A formal procedure to extract the form of $V_{\text{eff}}[\psi]$ near the vacuum-to-liquid transition requires to determine the generating functional $\Gamma[\psi_c]$ for the vertex functions from the exact partition function $Z[J]$ in the presence of an external field $J$ [45, 46]. In practice, this has been achieved in [47] for rotation invariant short-range interactions which become attractive beyond a zero crossing of the scattering length, where $g < 0$. The exact effective potential obtained from $\Gamma[\psi_c]$ then indeed contains a cubic term $k=3$ as in (33) whose strength $\lambda_3 = \hbar^2 D/(2m)$ is proportional the three-body scattering hypervolume $D$ introduced by Tan [48]. It has dimension (length)$^4$ and may be determined for a given two-body interaction $V_{\text{sr}}(r_{12})$ from a solution of the three-body scattering problem. Specifically, it is defined by the asymptotic behavior

$$\psi_{E=0}(x_1, x_2, x_3)|_{a=0} = 1 - \frac{\sqrt{3} D}{2\pi^3(r_{12}^2 + r_{13}^2 + r_{23}^2)^2} + \dots \quad \text{if } V(x_{12}) = V(r_{12}) \quad (34)$$

of the three-body wave function at zero energy and vanishing scattering length. The parameter $\lambda_3 \sim D$ characterizes the strength of effective three-body interactions which arise beyond mean-field even for a microscopic Hamiltonian which only contains two-body interactions. Explicit results for $D$ are available for simple model interactions like hard spheres [48], an attractive square well [49] or a Lennard-Jones potential [50]. In the two latter cases, the existence of two-body bound states gives rise to a finite imaginary part of the hypervolume. It determines the three-body loss coefficient $L_3 = -(\hbar/m) \text{Im} D = (\hbar/m) L_{\text{rec}}^4$ or the equivalent recombination length $L_{\text{rec}}$ [51]. Now, for a possible stabilization of a liquid state due to three-body interactions in the regime where $g_{\text{eff}}(\kappa) = g - g_{\text{dd}}f(\kappa) < 0$ is negative, it is the real part of the hypervolume $D$ which is relevant. In the realistic situation of finite three-body losses, the parameter which determines the liquid density $\bar{n} = 3|g_{\text{eff}}|/(4\bar{\lambda}_3)$ from the condition $p_{\text{sr}}(x) + \mathcal{D}(x) = 0$ in the presence of the additional pressure $\Delta p_{3\text{-body}}(x) = (2/3)\bar{\lambda}_3 [\rho^{(1)}(x)]^3$ due to three-body interactions is therefore $\bar{\lambda}_3 = \hbar^2 \text{Re} D/(2m)$. Evidently, $\text{Re} D$ must be positive for a stable liquid, a condition which can be verified only by a concrete solution of the three-body problem. In particular, $\text{Re} D$ assumes both positive and negative values, depending on the position of three-body bound states which is determined by the poles of the hypervolume [49–51]. In [2] and the literature beyond, a stabilization of self-bound droplets by repulsive three-body interactions in the regime $g_{\text{eff}} < 0$ has been excluded by the argument that the experimental value of the loss coefficient $L_3$ gives rise to an estimate for the magnitude of the complex parameter $\lambda_3$ which is far too small to be consistent with the observed central droplet densities $\bar{n} \approx 10^{-14} \text{ cm}^{-3}$. Since the real part $\bar{\lambda}_3 \sim \text{Re} D$ may be much larger than the estimate of $\lambda_3$ based on $|\text{Im} D|$, this argument is not compelling, however. For the moment, we will postpone a further discussion of the origin of the
beyond mean-field contribution $\Delta p \sim [\rho^{(1)}(x)]^k$ to the short-range pressure in Equation (18) and proceed by assuming that such a term is present on purely empirical grounds.

The shape-dependence of the effective coupling constant $g_{\text{eff}}(\kappa)$ that determines the position of the mean-field instability gives rise to a corresponding shape-dependence of the liquid density, e.g. via $\tilde{n}(\kappa) = 3|g_{\text{eff}}(\kappa)|/(4\tilde{\lambda}_3)$ in the case of repulsive three-body interactions. For a homogeneous fluid, this is a rather strange conclusion and indeed it is immediately changed if one takes into account that equilibrium in a uniaxial dipolar fluid requires in addition to $p_{\text{sr}}(x) + \mathcal{D}(x) = 0$ also that the anomalous field $h(x) = 0$ must vanish. Using the mean-field result (31), this can be obeyed only if $b(\kappa) = 0$ which implies that $\kappa$ is fixed to be either zero or infinity. The second possibility of an extremely oblate configuration is excluded since the dipolar interactions are then repulsive and the zero pressure condition cannot be fulfilled. By contrast, the limit $\kappa \to 0$ of a strongly prolate droplet turns out to describe the actual physical situation, consistent with the tendency of oriented dipoles to arrange in a head-to-tail configuration. More precisely, the condition $h(x) = 0$ gives rise to a finite value for the aspect ratio of a self-bound droplet rather than the singular result $\kappa^{(\text{mf})} \equiv 0$ within mean-field if one includes the contribution to $h(x)$ in Equation (16) which arises from the kinetic energy $\hat{H}_0$. This is a local term which—in contrast to the ill-defined kinetic contribution to $\text{Tr}[\Pi(x)]$—is finite even if the range of $V_{\text{sr}}(r_{12})$ is taken to zero. Indeed, if the interactions become rotation invariant at very short distances, the contribution of large momenta in $h^{(0)} \sim \int_k(k_x^2 + k_y^2 - 2k_z^2) n(k)$ cancels since the asymptotic momentum distribution $n(k) \to \mathcal{C}_2/|k|^4$ is then independent of the direction $\hat{k}$. A concrete result for $h^{(0)}$ can be obtained by assuming an anisotropic Gaussian profile with central density $\rho^{(1)}(x=0) = \tilde{n}^8$. In the relevant limit $\kappa \ll 1$, only the radial contribution in Equation (16) remains and the condition $h(x) = 0$ in the center of the droplet is reduced to

$$ \frac{\hbar^2}{2m} \pi \left( \frac{\tilde{n}}{N\kappa} \right)^{2/3} = \frac{3}{4} g_{\text{dd}} \tilde{n}^2 \cdot 6\kappa^2 \ln(2/\kappa) \rightarrow \kappa^{8/3} \ln(2/\kappa) = [36 \tilde{n}^{1/3} a_{\text{dd}} \cdot N^{2/3}]^{-1} = \epsilon. \quad (35) $$

With typical particle numbers $N \approx 10^4$ and central densities $\tilde{n}^{1/3} a_{\text{dd}} = \mathcal{O}(1)$, the smallness of the parameter $\epsilon \approx 6 \cdot 10^{-5}$ allows an asymptotically exact solution of the transcendental equation for $\kappa(\epsilon)$ in the form

$$ \kappa(\epsilon) = \left[ \frac{\epsilon}{\ln(2/\epsilon^{3/8})} \right]^{3/8} \rightarrow \kappa(N) = \frac{0.26}{[\tilde{n}^{1/3} a_{\text{dd}} \cdot \ln[2/\epsilon^{3/8}(N)]]^{3/8}} \cdot N^{-1/4}. \quad (36) $$

Apart from the logarithmic correction in the denominator which changes the result only by a numerical factor of order one, the aspect ratio thus approaches zero according to $\kappa(N) \sim N^{-1/4}$ for large particle numbers. In practice, with $N < 10^5$, typical aspect ratios are in the range $1/\kappa \approx 5-30$ [42, 52]. An important point to note is that the result (36) holds irrespective of the specific mechanism which stabilizes the central density $\tilde{n}$ at a finite value provided only that $\tilde{n}^{1/3} a_{\text{dd}}$ does not itself depend on $N$. The dependence on the deviation $\epsilon_{\text{dd}} - 1$ from the critical strength $\epsilon_{\text{dd}}^c = 1$ of the dipolar interactions, which is now shape-independent due to $\kappa \approx 0$, is sensitive to that, however. In particular, for a beyond mean-field contribution $\sim |\psi(x)|^{2k}$, the liquid density approaches zero like $\tilde{n}(\epsilon_{\text{dd}}) \sim (\epsilon_{\text{dd}} - 1)^{1/(k-2)}$. On physical grounds, the fact that the equilibrium configuration of self-bound dipolar fluids corresponds to an increasingly prolate object for large particle numbers is a consequence of the anisotropy of the interactions. In contrast to the situation of attractive short-range interactions stabilized by a repulsive three-body force discussed in [47, 53], there is no homogeneous liquid state at all. This quite unusual behavior is elucidated by a few further comments:

<footnote n="8">⁸This assumption has been made in this context by Dalibard [8]. More precisely, an exact result due to Triay [14] shows that—within an extended Gross-Pitaevskii description and for a negative chemical potential—the density profile of a non-confined dipolar gas is smooth and falls off exponentially.</footnote>
* The existence of a finite aspect ratio rather than the mean-field value zero relies on the kinetic energy contribution to the anomalous field. This must be carefully distinguished from a standard quantum pressure term, which may stabilize only self-bound droplets whose size vanishes as $N \gg 1$. Examples are the bright soliton in the Lieb-Liniger model with positive two-body scattering length $a_1 > 0$ or the droplets formed by attractive three-body interactions studied by Sekino and Nishida [54]. In both cases, the droplet extension approaches zero for large $N$ according to $R_N \to a_1/N$ or $R_N \sim \exp[-(4/\sqrt{3\pi})N^2]$ instead of increasing like $L_z(N) \sim \sqrt{N}$ in the case of 3D dipolar gases.

* There is a bound $N_0$ on the particle number below which the droplets evaporate. This has been investigated by Baillie et al. [52], where $N_0(\epsilon_{dd})$ has been determined from an extended Gross-Pitaevskii description by the condition of a finite binding energy despite the competing quantum pressure term. The analogous problem for a Bose droplet with short-range interactions has been studied in [47] and—with concrete predictions for the lifetime of metastable droplets in a finite window $N_1 < N < N_0$—by Son et al. [55]. As discussed in [47], the critical number $N_0 \sim \sqrt{D/a^2}$ in this context determines the scattering lengths at which $N$-body bound states detach from the continuum.

* Non-spherical self-bound objects have been studied a long time ago for stars with strong, frozen-in magnetic fields B by Chandrasekhar and Fermi [56]. The balance between the decrease of the attractive gravitational energy with a finite eccentricity and the gain in magnetic energy due an expansion in the direction perpendicular to the field results in an oblate deformation with aspect ratio $\kappa - 1 \approx (B/B_*)^2$, where $B_*$ is the field beyond which the star is no longer bound since magnetic pressure overwhelms gravitation.

Following the suggestion in the original publication on self-bound droplets [2], the mechanism for stabilizing a non-confined dipolar gas beyond the critical interaction strength $\epsilon_{dd} = 1$ where the homogeneous fluid is unstable is commonly believed [6] to be a generalized form of the Lee, Huang and Yang (LHY) correction, originally calculated for a dilute, homogeneous hard-sphere Bose gas [57]. Including the effect of anisotropy in the presence of an additional dipolar interaction [58], it gives rise to a beyond mean-field contribution to the effective potential of the form introduced in Equation (33) with $k = 5/2$ and a positive (as long as $\epsilon_{dd} \le 1$) coefficient

$$ \lambda_{5/2}(g, \epsilon_{dd}) = \frac{4g}{3\pi^2} (4\pi a)^{3/2} \cdot Q_5(\epsilon_{dd}) \quad \text{where } Q_5(\epsilon_{dd}) = \int_0^1 dx [1 + \epsilon_{dd} \cdot 2P_2(x)]^{5/2}. \quad (37) $$

For the standard case of repulsive short-range interactions, the LHY term provides the leading correction beyond mean-field in the equation of state, giving rise to an additional contribution of order $gn^2\sqrt{na^3}$ to the pressure. It is important here to carefully separate the dependence on density from that on the chemical potential as the thermodynamically conjugate variable. In particular, for a homogeneous fluid, there is an exact relation $-V_{\text{eff}}[n(\mu)] = p(\mu)$ which allows to infer the effective potential evaluated at the equilibrium value $n(\mu) = \partial p(\mu)/\partial\mu$ of the density from the pressure at given $\mu$ [47]. The increase of pressure with density which results from the positive LHY contribution in $V_{\text{eff}}[\psi]$ thus arises from a corresponding negative correction to $p(\mu)$. In explicit form, this emerges within a field-theoretic formulation where the beyond mean-field contribution determined by LHY is just the one-loop correction [45, 46]

$$ p^{(1)}(\mu) = -\frac{1}{2} \text{Tr} \ln \det [\delta^{(2)}S] \underset{\text{hom.}}{\longrightarrow} -\frac{1}{2} \int_q \int_\omega \ln[E_q^2 + (\hbar\omega)^2] \underset{\text{reg.}}{\longrightarrow} -\frac{8}{15\pi^2} \mu \left(\frac{m\mu}{\hbar^2}\right)^{3/2} \quad (38) $$

to the tree-level result $p^{(0)}(\mu) = \mu^2/(2g)$ for the pressure of a dilute Bose gas<fnref n="9" />. In physical terms, it describes the contribution to pressure due to zero-point fluctuations of the Bogoliubov

<footnote n="9">⁹Note that $p^{(1)}(\mu)$ does not explicitly depend on the scattering length. The contact density $\epsilon_2(\mu) \to (m\mu/\hbar^2)^2$ is therefore still given by its mean-field value and the LHY correction of relative order $(na^3)^{1/2}$, which contributes to a corresponding one in the pressure (18), only shows up as a function of density.</footnote>
excitations. Similar to the analogous Casimir effect, it is negative and it comes entirely from the low-energy part of the spectrum [46, 59]. In particular, after regularization of the divergent integration over momenta $q$, the remaining expression involves an integral $\int_q E_q^2 \sim \int_0 q^2 dq c^2(\hat{q})q^2$. Here—in an extension to anisotropic interactions—we include a possible dependence of the sound velocity $c(\hat{q})$ on direction. Now, it is a crucial point that the correction (38) to the pressure is fully determined by the excitations near $q=0$. As a consequence, an anisotropic sound velocity $c^2(\hat{q}) = c^2[1 + \tilde{\epsilon}f(\hat{q})]$ with strength $\tilde{\epsilon}$ and an arbitrary function $f(\Omega)$ of direction changes the LHY contribution (38) just by a numerical factor $Q(\tilde{\epsilon}) = \int(d\Omega)/(4\pi)[1 + \tilde{\epsilon}f(\Omega)]^{5/2}$. For the special case of dipolar interactions, this immediately explains the origin of the function $Q(\tilde{\epsilon}) \rightarrow Q_5(\epsilon_{dd})$ derived by Lima and Pelster [58]. It is obvious, however, that the LHY correction is well defined only as long as $c^2(\hat{q}) > 0$ remains positive. This requirement is violated in dipolar gases if $\epsilon_{dd} > 1$, where the homogeneous fluid is unstable. Formally, this shows up as a finite imaginary part of the function $Q_5(\epsilon_{dd})$ which is, however, neglected in the extended Gross–Pitaevskii description of the inhomogeneous configurations that are considered in practice.

For a proper understanding of the mechanism which underlies the stabilization of self-bound droplets and which allows to fulfill the zero pressure condition $p_{sr}(x) + \mathcal{D}(x) = 0$, it is of course necessary to provide a physical argument that justifies the addition of a contribution $\Delta p_{LHY}(x) = (3/5)\lambda_{5/2}[\rho^{(1)}(x)]^{5/2}$ of the LHY-form to the short-range pressure in a regime where the chemical potential is negative and the microscopic derivation along the lines in Equation (38) evidently fails. On a rather qualitative level, the inclusion of an LHY term even in the regime $\epsilon_{dd} > 1$ may be justified by noting that the characteristic radial extension $R_\perp \sim N^{1/4}$ of self-bound droplets is only of the order of the healing length $\xi$ [8]. As a result, there is no sound propagating in the $(x, y)$-plane in such a configuration and the fact that $c^2(\hat{q}) < 0$ for the associated wave vectors is irrelevant. The argument is not conclusive, however, because it does not address the problem that an exponent $5/2$ in the relation between pressure and density requires a corresponding power law in $p(\mu)$, which is excluded for negative values of $\mu$. To obtain some further insight, it is useful to investigate the microscopic origin of the LHY correction from a different point of view. The associated characteristic exponent $k = 5/2$ indicates that it can be attributed neither to two-body nor to three-body interactions which—as discussed above—lead to $k=3$. Now, as noted already by Lee, Huang and Yang [57] and expanded in detail later by Lüscher [60] and by Tan [48], the LHY term can be understood to arise from a finite size correction in the two-body problem. More precisely, the repulsive contribution $\sim gn(na^3)^{1/2}$ in the interaction energy per particle beyond mean-field emerges from a correction of order $a/L$ beyond the leading term $E_{N=2}^{(0)}(L) = g/L^3$ in the two-particle problem in a box with periodic boundary conditions. At finite density, the relevant size $L \approx \xi$ for this correction is set by the healing length beyond which the pair distribution is no longer determined by two-body physics. In the presence of the long-range and anisotropic dipolar interactions, both $E_{N=2}^{(0)}(L)$ and its finite size correction will depend on the aspect ratio $\kappa$ and the dimensionless strength parameter $\epsilon_{dd}$. A microscopic derivation of the LHY correction for dipolar gases in a configuration which mimics the elongated droplets thus requires the solution of the two-body problem in an anisotropic box of size $L \times L \times L/\kappa$ with periodic boundary conditions. Specifically, the parameter $\lambda_{5/2}(\kappa, \epsilon_{dd})$ is determined by the coefficient of the $a/L$—contribution to the energy in a finite size expansion analogous to that derived by Lüscher for rotation invariant short-range interactions [60]. While straightforward in principle, proving that $\lambda_{5/2}$ is finite and positive along these lines is nontrivial. A fully microscopic derivation of either the LHY correction or a possible repulsive three-body force that might contribute to the stabilization of self-bound droplets thus remains an open problem. It should be emphasized, however, that this is essentially a conceptual challenge, not a practical one. Indeed, as demonstrated recently by Bombín et al. [61], numerical approaches to dipolar gases in the regime $\epsilon_{dd} > 1$ based on quantum Monte Carlo provide results which agree quite
well with experiment, e.g. for the critical number $N_0$ where droplets unbind or for the supersolid transition.

### 3. Effective theories for the transition to a supersolid

The prediction that dipolar gases in a pancake geometry with $\lambda \gg 1$ will exhibit a roton-maxon character in their excitation spectrum has been made even before the experimental realization of dipolar condensates [62, 63]. In contrast to the situation in $^4$He, where the roton gap $\Delta_r$ is essentially independent of pressure up to the superfluid-to-solid transition at $p_c \approx 25$ atm, the depth of the roton minimum is now tunable. In fact, within mean-field, the roton gap in dipolar gases is found to vanish beyond a critical value of the interaction strength where the spectrum of Bogoliubov excitations reaches zero [62, 63]. A simple model capturing the essential physics is based on assuming a uniform quasi-2D situation with a Gaussian density profile in the transverse direction with characteristic length $\ell_z$. The resulting effective dipole-dipole interaction [63, 64]

$$V_{dd}(q) = \frac{\hbar^2}{m} \tilde{g}_{dd} \left[ 2 - 3 \sqrt{\frac{\pi}{2}} (q\ell_z) \exp(q^2 \ell_z^2/2) \operatorname{erfc}(q\ell_z/\sqrt{2}) \right] \quad (39)$$

in momentum space $q = (q_x, q_y)$ approaches a positive constant $V_{dd}(q) \to 2\hbar^2 \tilde{g}_{dd}/m$ in the limit $|q|\ell_z \ll 1$, where $\tilde{g}_{dd} = \sqrt{8\pi a_{dd}/\ell_z}$ is a dimensionless coupling constant. The fact that $V_{dd}(q=0)$ is positive guarantees that—in contrast to the situation in 3D—the homogeneous fluid is now stable for an arbitrary strength of the dipolar interaction. This does not exclude an instability at finite wave vector, however. Indeed, the effective interaction $V_{dd}(q)$ turns negative beyond $|q|\ell_z$ of order one and asymptotically approaches the constant value $-\hbar^2 \tilde{g}_{dd}/m$. In physical terms, this describes attractive head-to-tail collisions between aligned dipoles at distances less than $\ell_z$ with an effective scattering length $-a_{dd}$. For strong dipolar interaction strengths, the negative contributions to the total interaction $V(q) = \hbar^2 \tilde{g}_{dd}/m + V_{dd}(q)$ give rise to a roton minimum in the spectrum $E_q^2 = 2nV(q)\epsilon_q + \epsilon_q^2$ of Bogoliubov excitations which eventually touches zero at a characteristic wave vector $q_0$. A crucial point in this context is that the scale for $q_0$ is set by the inverse of the confinement length $\ell_z$, which is unrelated to and much larger than the average interparticle spacing in the transverse direction.

The predicted appearance and subsequent softening of a roton has been observed in dipolar condensates of Dysprosium in a cigar-shaped trap by Petter et al. [65] using Bragg spectroscopy. By tuning $\epsilon_{dd} = a_{dd}/a$ via a Feshbach resonance that allows to change the short-range scattering length at fixed $a_{dd}$, the measured dynamic structure factor shows a sharp increase near a critical value $\epsilon_{dd}^c$ of order one, consistent with theoretical results based on a solution of the Bogoliubov equations [66]. This signals a macroscopic occupation of the roton mode associated with a spontaneous density modulation along the weakly confined direction in the trap, as seen directly from in-situ density profiles [67]. The formation of a static density wave at a coupling strength where the roton dip is still not very pronounced and the interpretation of this state in terms of a supersolid has been confirmed in a number of experiments [6] which will not be discussed here in detail. Instead, we will present an elementary approach which allows to understand the nature of the associated phase transition and the underlying physics independent of specific details. In particular, based on a simple model due to Nozières [11], we will show that the mean-field roton instability is preempted into a first-order transition by an amount which depends on the strength of the short-range repulsion. As a first step, it is necessary to properly define the notion of a supersolid which, in a rather broad sense, may be characterized by:

**In a supersolid, superfluidity (defined by a finite superfluid fraction) appears together with a non-vanishing modulation of the density due to spontaneously broken translation invariance.**
Now, according to this definition, any spatially modulated superfluid like the vortex lattice is also a supersolid. To exclude such well known cases, the notion of a genuine supersolid should therefore be restricted to phases where spontaneously broken translation invariance and superfluidity are present simultaneously as *two independent* order parameters. The modulation in density is then still present even after superfluidity is lost. An important result due to Leggett [68] states that in any superfluid with a non-uniform density, there is an upper bound on the superfluid fraction $f_s$ strictly smaller than one. Specifically, the Leggett bound is of the form

$$f_s(T=0) \leq \left[ \frac{\bar{n}}{d} \int_0^d \frac{dx}{n_1(x)} \right]^{-1}, \quad \text{in mean-field } f_s^{(\text{mf})}(T=0) = m/m_B \quad (40)$$

where the inverse of the density $n_1(x)$ is integrated over a unit cell of the lattice whose length is denoted by $d$. Here, $n_1(x)$ is the average of the microscopic density $\langle \hat{\rho}^{(1)}(x) \rangle$ over the transverse directions of the unit cell and $\bar{n}$ its average value (without loss of generality, the $x$-direction has been singled out). Apparently, the bound is always finite unless the density vanishes identically in some region. Within mean-field, the ground state wave function factorizes into a product of single-particle ones. The Leggett bound then turns into an equality and relates the superfluid fraction to the ratio $m/m_B$ of the bare and the band mass $m_B$ in the given periodic potential. For a conventional solid, where the particle density is concentrated near a discrete set $\{\mathbf{R}\}$ of lattice sites with an exponentially suppressed value at interstitial positions, the upper bound on $f_s$ is much smaller than one. In the limit of a fluid with uniform density, in turn, Equation (40) reduces to the trivial identity $f_s \leq 1$. Obviously, supersolids with an appreciable superfluid fraction can only be found in situations where the density modulation is weak. A few important points should be noted in this context: First, the bound (40) does not provide a sufficient criterion for superfluidity in a state with broken translation invariance: a finite value of the bound is still compatible with no superfluidity at all. A case in point is the solid phase of $^4$He, where the microscopic density profile is known from path integral Monte Carlo and the Leggett bound gives $f_s \leq 0.16$ [69] while the true value vanishes$^{10}$. A second point is that the bound (40) makes no assumption about the physical origin of the density modulation. It may arise from a spontaneous breaking of translation invariance but it also holds if the density modulation is externally imposed. In fact, the latter method has been used in a recent experimental test of the Leggett bound [72]. Finally, a quite subtle point is that the bound holds independent of whether the number of particles within a unit cell happens to be an integer or not. Now, as emphasized in the classic paper by Andreev and Lifshitz [34], the generic realization of supersolids requires an incommensurate situation with a finite concentration of defects in the ground state. A detailed argument which shows that delocalized vacancies or interstitial atoms are indeed a necessary condition for a supersolid has been given by Prokof'ev and Svistunov [73]. This suggests that supersolids are in general just superfluid mass-density waves, with the homogeneous part of the density playing the role of delocalized “defects”. The ground state of generic solids like $^4$He, in turn, has an integer number of particles per unit cell and vanishing defect concentration [70]. The finite energy necessary for the creation of either vacancies or interstitials identifies such a state as a Mott-insulator, which cannot have a finite superfluid fraction.

### 3.1. Mean-field theory of freezing in classical and quantum systems

It is a classic argument due to Landau [74] that the transition from a uniform fluid to a state with a periodic modulation of the density is of first order. The argument relies on an expansion of the

<footnote>$^{10}$For a discussion of why a supersolid phase of bulk $^4$He can be ruled out, see [69, 70] and the review [71].</footnote>
free energy in terms of the Fourier components $n_G \neq 0$ of the density

$$ \rho^{(1)}(x) = \langle \rho^{(1)}(x) \rangle + \sum_{G \neq 0} n_G e^{iG \cdot x} = n_{\text{sol}} + \sum_{G \neq 0} n_G e^{iG \cdot x} \quad (41) $$

in the symmetry broken phase with average density $n_{sol}$, where $\{G\}$ denotes the set of associated reciprocal lattice vectors [10]. In the weak crystallization limit, the density jump $\Delta n = n_{sol} - n_{liq}$ is small and the instability is dominated by a single wave vector $q_0$ where the static structure factor of the homogeneous fluid exhibits a pronounced maximum $S(q_0)$. The magnitude $|G| = q_0$ of the reciprocal lattice vectors is then fixed and the problem reduces to finding the lattice with the lowest free energy. Within Landau theory, this is of the generic form [10]

$$ \Delta f_L(T) = f_{\text{sol}} - f_{\text{liq}} = \frac{r(T)}{2} n_G^2 - w n_G^3 + u n_G^4 \quad \text{with } r(T) = T - T^* \text{ and } u > 0. \quad (42) $$

Here, $n_G^2 \sim \sum'_G |n_G|^2$ is a sum of the magnitudes of the non-vanishing Fourier components with fixed $|G| = q_0$, rescaled in such a way that makes it independent of their number in the specific lattice and dimensionless<fnref n="11" />. Now, the simplified form (42) hides the dependence on the specific lattice structure which is contained in the detailed values of $w$ and $u$. For a finite $w \neq 0$, reciprocal lattices in which triads of different $G$'s add up to zero are favored. At fixed $|G| = q_0$, therefore, the planar hexagonal, the fcc and icosahedral lattices are the only possibilities [10], giving rise to a triangular lattice in real space as the unique option in 2D. Quite generally, the transition temperature $T_c$ is determined by the condition $\Delta f_L(T_c) = 0$ at a finite value $n_G(T_c) = w/(2u) \neq 0$, where the fluid with $n_G \equiv 0$ is degenerate with a state where the density develops a non-vanishing spatial modulation. Due to the presence of the third-order term, the critical temperature $T_c = T^* + w^2/(2u)$ lies above the temperature $T^*$ introduced in the phenomenological Ansatz for $r(T)$. Physically, $T^*$ is the lowest temperature down to which the fluid exists as a metastable configuration. Despite its purely phenomenological nature, the Landau theory of crystallization captures many aspects of real first-order fluid-to-solid transitions, like the preference for (real space) bcc lattices near the melting line [75] or the fact that a fluid can be undercooled by a far larger amount than solids can be overheated. Indeed, the temperature below which the solid exists at least as a metastable configuration is given by $T_1 = T^* + 9w^2/(16u)$ and the Ansatz in Equation (42) thus leads to $T_c - T^* = 8(T_1 - T_c)$. A further consequence of Landau theory, not mentioned usually, is that the static structure factor in the uniform fluid right at the transition has a quasi-universal value $S(q_0)|_{T_c} = T_c/r(T_c) = T_c/(T_c - T^*)$ of order one. A similar result was in fact found to hold for strong crystallization transitions e.g. of a Lennard-Jones fluid by Hansen and Verlet [76]. It has led to the empirical Hansen-Verlet criterion $S_{HV}(q_0)|_{T_c} = 2.85$ as a system-independent estimate for the position of the fluid-to-solid transition in classical systems with short-range interactions from properties within the fluid phase itself.

The Landau theory for the fluid-to-solid transition in a classical system can be extended into a quantum theory which describes the appearance of a mass-density wave in a superfluid. Based on an approach indicated already in the ground-breaking work of Gross on inhomogeneous Bose fluids [77], this has been developed by Pomeau et al. [78, 79]. Focussing again on density fluctuations with wave vectors of a fixed length $|G| = q_0$, the leading order expansion $\sum_q |\delta n_q|^2/2\chi(q)$ in the change of energy through density fluctuations $\delta n_q \neq 0$ in any fluid gives rise to a contribution $\tilde{r} n_G^2/2$, where $\tilde{r} = 1/\chi(q_0)$ is the inverse of the static density response function. Now, within the Bogoliubov approximation, this response is determined by the ratio $\chi(q) = 2\epsilon_q/E_q^2$ of the free particle energy $\epsilon_q$ and the dispersion $E_q$ of the collective excitations of the superfluid [80]. In particular, for interactions that lead to a roton minimum with energy $\Delta_r = E_{|q|=q_0}$, the parameter $\tilde{r}$ near the corresponding wave vector is given by $\tilde{r} = \Delta_r^2/(2\epsilon_r)$ with $\epsilon_r = \epsilon_{|q|=q_0}$. It is important

<footnote n="11">¹¹Note that the sign of $w$, assumed to positive here, is irrelevant since it can be changed trivially by $n_G \to -n_G$.</footnote>
to note that—within the Bogoliubov approximation—a roton minimum requires two-body interactions $V(q)$ which are negative in a certain range of wave vectors but not necessarily an attractive interaction in real space as in dipolar gases. As an example, a purely repulsive box potential $V(r_{12})$ of strength $V_0$ and size $\sigma$ gives rise to a negative $V(q)$ in a range $4.5 < q\sigma < 7.7$. The energy of the resulting roton minimum and the associated parameter $\tilde{r}$ can then be tuned by changing the dimensionless coupling strength $n\sigma^3 \cdot V_0 m\sigma^2/\hbar^2$. This leads to a roton instability at the point where $\tilde{r}$ vanishes and the uniform superfluid becomes unstable. The nontrivial question now is whether this instability can be cured by a nonlinear contribution to the energy analogous to the one which is introduced by hand in Landau's Ansatz (42). As shown by Pomeau and Rica [78] this is indeed the case and the resulting energy functional

$$ \Delta f_{\text{GP}} = f[n(x)] - f[\tilde{n}] = \frac{\tilde{r}}{2} n_G^2 - \tilde{w} n_G^3 + \tilde{u} n_G^4 \quad \text{with} \quad \tilde{r} = \Delta_r^2 / (2\varepsilon_r) \quad (43) $$

is in fact of the same form. Moreover, in contrast to the case of classical fluid-to-solid transitions where the parameters $r(T)$, $w$, $u$ are introduced in a purely phenomenological manner, the coefficients in Equation (43) can now be derived from a microscopic energy functional. In particular, the nonlinear terms involving $\tilde{w}$ and $\tilde{u}$ follow from the expansion of the denominator in the quantum pressure contribution $\hbar^2$ [grad $n(x)$]$^2$/$[8m n(x)]$ around the fluid state with uniform density $\tilde{n}$. For a given characteristic wave vector $q_0$ of the emerging density wave, they are on the order of the associated recoil energy $\tilde{w}, \tilde{u} \approx \varepsilon_r$, with detailed values again depending on the specific lattice. As in the standard Landau theory, the presence of a finite third-order term $\tilde{w} \neq 0$ implies that the uniform superfluid freezes with a jump in density determined by $n_G|_c = \tilde{w}/(2\tilde{u})$ at a critical value $\tilde{r}_c = \tilde{w}^2/2\tilde{u}$. The mean-field roton instability at $\tilde{r} = 0$, which characterizes the interaction strength up to which the homogeneous superfluid exists at least as a metastable configuration, is thus preempted by first-order transition at a finite critical value $\Delta_r/\varepsilon_r|_c = \mathcal{O}(1)$ of the roton gap. For the specific case of a triangular lattice in 2D, this value turns out to be $\Delta_r/\varepsilon_r|_c = 0.23$ [78].

## 3.2. Effects beyond mean-field and a two-mode model

The extension above of Landau's classical theory to the freezing of a uniform superfluid into a state with broken translation invariance describes a mean-field supersolid in the sense that its condensate fraction stays at the non-interacting value $f_0^{(\text{mf})} = 1$. By contrast—consistent with the Leggett bound (40)—the superfluid fraction is reduced below one depending on the magnitude of the density modulation. The model thus provides a qualitative description of the supersolid transition in dipolar gases where the emerging density wave is a phase-coherent superfluid. This is rather different from the situation in $^4$He, where the solid phase beyond the critical pressure $p_c$ is a commensurate Mott-insulator [70]. The result that the roton gap has a finite value $\Delta_r^c$ at the supersolid transition implies a quantum generalization of the Hansen–Verlet criterion. Indeed, the expression $\chi(q) = 2\varepsilon_q/E_q^2$ for the density response also holds for finite temperatures below the superfluid transition as long as the condensate fraction remains close to one [80]. As a result, the static structure factor of the homogeneous superfluid at finite temperature is related to its ground state value $S(q) = \varepsilon_q/E_q$ by a thermal enhancement factor coth($E_q/2T$). Right at the supersolid transition, it has thus a quasi-universal value

$$ S_c(q_0) = \left( \frac{\varepsilon_r}{\Delta_r^c} \right) \coth \left( \frac{\Delta_r^c}{2T} \right) \xrightarrow{T \gg \Delta_r^c} \frac{T}{\tilde{r}_c} \gg 1 \quad (44) $$

determined only by the ratios between the roton gap $\Delta_r^c$ and the recoil energy or temperature. As indicated, the value is large compared to one in the experimentally relevant limit $T \gg \Delta_r^c$. In particular, this applies to the measurements of the static structure factor in [67], where a pronounced peak in $S(q_0)$ near the ordering wave vector is found at the supersolid transition.
There are two major shortcomings of the approach sketched above: First of all, the description fails in cases where the emerging lattice does not give rise to a third-order invariant in the free energy. This applies e.g. to the situation realized in many experiments where the density modulation appears along a single direction preferred by the geometry in cigar-shaped traps. The functional (43) then leads to a continuous transition at the point where the roton gap vanishes which is not consistent with observation<fnref n="12" />, see e.g. [67]. A second point is that, even at finite values of $\tilde{w}$, the shift of the roton instability away from its mean-field value zero to a finite $\tilde{r}_c = \tilde{w}^2/2\tilde{u}$ predicted by Equation (43) gives rise to a critical roton gap of order $\epsilon_r$, independent of the strength of interactions. In the following, we will outline an approach in which both of these problems are absent. The crucial point to recognize is that beyond mean-field there is a fundamental difference between density fluctuations and the quasi-particles of the superfluid. It is only in the limit of weak interactions where they coincide and one may infer the density response $\chi(q) \to 2\epsilon_q/E_q^2$ from the quasi-particle dispersion $E_q$. A simple model for dealing with the interplay of density modes and superfluid quasi-particles was suggested by Nozières [11], whose aim was to shed light on the complex question about the role of the roton minimum in superfluid ⁴He for the transition to a solid. Based on a diagrammatic argument, Nozières showed that beyond the standard leading order term $\Sigma_q(\omega) \to n_0 V(q)$ of the self-energy, there is a contribution proportional to $V^2(q) \cdot \chi(q, \omega)$ which involves the dynamic density response function $\chi(q, \omega)$. Within a simple Ansatz $\chi(q, \omega) \approx 2\epsilon_q/(\Omega_q^2 - \omega^2)$ for this response in the absence of the coupling to superfluid quasi-particles, a two-mode model emerges whose excitation energies $E_{\pm}^2(q)$ follow from a bi-quadratic equation [11]

$$E_{\pm}^2(q) = E_q^2 - f_0 \frac{\Lambda_q}{\Omega_q^2 - E_{\pm}^2(q)} \rightarrow \tilde{r}^{\text{eff}}(q) = \frac{E_{-}^2(q)}{2\epsilon_q} \approx \frac{E_q^2}{2\epsilon_q} - f_0 \frac{\Lambda_q}{2\epsilon_q \cdot E_q^2} = \tilde{r}^{\text{mf}}(q) - \frac{f_0 \Lambda_q / (2\epsilon_q)^2}{\tilde{r}^{\text{mf}}(q)}. \quad (45)$$

The effective coupling contains the condensate fraction $f_0$ and the parameter $\Lambda_q$ which is at least quadratic in the interactions. With increasing coupling, the two modes $E_{\pm}^2(q)$ shift in an opposite direction and it is the lower one $E_{-}^2(q)$ whose softening signals the onset of a density wave at given $q$. In particular, the associated effective stiffness $\tilde{r}^{\text{eff}}(q)$ is renormalized down compared with the mean-field value by an amount which scales inversely with $\tilde{r}^{\text{mf}}(q)$ itself<fnref n="13" />. The fact that a mean-field approximation overestimates the stiffness for density fluctuations is consistent with the exact inequality $m_p^2(q) \le m_{p+1}(q) m_{p-1}(q)$ for the pth moment of the dynamic structure factor [80]. In fact, in the special case $p=0$, this reduces to

$$S^2(q) \le \epsilon_q \cdot \chi(q)/2 \rightarrow 1/\chi(q) \le \frac{\epsilon_q}{2S^2(q)} \xrightarrow{\text{mf}} \frac{\Delta_r^2}{2\epsilon_r} \quad (46)$$

in general, respectively evaluated at the roton minimum within mean-field. It is only in the limit where the dynamic structure factor has a single sharp peak that Equation (46) becomes an equality. A quantitative result for the renormalization of $\tilde{r}^{\text{mf}}(q)$ requires to determine the parameter $\Lambda_q$. This may be estimated by assuming that the relevant wave vector is large enough that $\Lambda_q$ can be inferred from the exact short-distance expansion

$$\frac{1}{\chi(q)} \rightarrow \frac{\epsilon_q}{2} \left[ 1 - \frac{\pi c_2}{8nq} + \dots \right] \rightarrow \frac{\Lambda_q}{\epsilon_q^4} = \frac{\pi c_2}{8nq} \quad (47)$$

for Bose gases with repulsive short-range interactions derived in [83]. The requirement that this is consistent with Equation (45) then fixes $\Lambda_q \sim \epsilon_q^4 c_2 / q$ to be proportional to the two-body

<footnote n="12">¹²This problem may be eliminated by incorporating a further wave vector $q'_0$ for the density wave, which breaks the symmetry $n_G \rightarrow -n_G$. The option comes, however, at the expense of introducing at least two additional free parameters.</footnote>

<footnote n="13">¹³This is reminiscent of the Brazovskii equation in the beyond mean-field description of classical fluid-to-solid transitions [81, 82] but note that the sign here is opposite and the effect scales with the condensate fraction $f_0$.</footnote>
contact density $\mathcal{C}_2$. Now, the characteristic wave vectors of the actually observed supersolids are considerably smaller than those where the expansion for high momenta in (47) is expected to apply and—moreover—effects of the long-range dipolar interactions have been ignored. The result is thus only of a rather qualitative nature. Nevertheless, it indicates how the problems of the mean-field approximation mentioned above can be resolved. Specifically, from the fact that the inverse static density response $\tilde{r}^{\text{eff}}(q_0)$ of the homogeneous fluid at the freezing transition is expected to be finite but very small, the critical value

$$ \left( \frac{\Delta_r}{\varepsilon_r} \right)_c^4 \approx f_0 \frac{\pi \mathcal{C}_2}{8nq_0} \qquad (48) $$

for the ratio between the roton gap $\Delta_r$ and the recoil energy $\varepsilon_r$ at the transition into a supersolid state can be estimated. Independent of the existence of a third-order invariant, it is always finite. Moreover, it depends explicitely on the interaction strength via the two-body contact density $\mathcal{C}_2$. In particular, the simple mean-field roton instability with a vanishing $\Delta_r^{c1}(\text{mf}) \equiv 0$ is obtained in the limit of vanishing short-range repulsion, where $\mathcal{C}_2 \sim a^2 \to 0$. It is obvious that the rather crude arguments above need to be replaced by a proper microscopic theory for effects beyond mean-field in the transition to a supersolid. This is an open but clearly quite challenging problem.

## 4. Conclusion and open problems

The thermodynamic approach developed in this work provides a description of dipolar quantum fluids which fully accounts for the anisotropy and the long-range nature of the interactions. It allows to derive a number of exact results, e.g. for the effective magnetic field in the center of a spherical cloud (11) and the resulting instability beyond a critical value of the interaction strength, the violation of the Gibbs-Duhem relation in the density profile of trapped gases (21) or the number-dependence (36) of the aspect ratio in a self-bound droplet. An interesting perspective to obtain further exact results is provided by a possible extension into an effective field theory, similar to the one that has been developed by Son and Wingate for the unitary Fermi gas [84]. Concerning the transition into a supersolid phase, the discussion in this work has been of a rather qualitative nature and a better understanding of the effects beyond mean-field is an open problem.

It is obvious, that a number of basic problems in the theory of dipolar gases have only been raised but have not been answered. This is true, in particular, for a microscopic derivation of the exact effective potential which replaces the simple ad hoc Ansatz assumed in Equation (33). Specifically, both the detailed form of an LHY-like contribution $\sim |\psi|^5$ with a proper positive strength $\lambda_{5/2}(g, \varepsilon_{dd}, \kappa) > 0$ and a quantitative theory for a possible contribution of three-body forces are missing. In addition, the transition from the supersolid phase to a collection of separate, incoherent droplets [6] is still not well understood.

Finally, the approach to uniaxial quantum fluids presented here applies to a situation where the breaking of rotation invariance is imposed externally by the fixed orientation of the dipoles. It is a challenge to see whether it can be extended to cover a situation where the symmetry breaking is spontaneous, e.g. in a quantum version of a nematic liquid crystal<fnref n="14" />.

## Declaration of interests

The authors do not work for, advise, own shares in, or receive funds from any organization that could benefit from this article, and have declared no affiliations other than their research organizations.

<footnote n="14">¹⁴I am grateful to Alan Dorsey for pointing out such a possibility.</footnote>
Acknowledgements

The present contribution is based on a Lecture on dipolar gases and supersolids given at the Institut Henri Poincaré in July 2024. I am very grateful to Yvan Castin for the invitation to this Lecture and for extensive discussions. It is also a pleasure to thank Jean Dalibard for a number of comments and a copy of the Notes for his Lectures on the subject at the Collège de France [8]. Finally, I am indebted to Markus Holzmann for clarifying some basic points about the freezing transition of quantum fluids.

References

[1] A. Griesmaier, J. Werner, S. Hensler, J. Stuhler, T. Pfau, “Bose–Einstein condensation of chromium”, Phys. Rev. Lett. 94 (2005), article no. 160401.

[2] I. Ferrier-Barbut, H. Kadau, M. Schmitt, M. Wenzel, T. Pfau, “Observation of quantum droplets in a strongly dipolar Bose gas”, Phys. Rev. Lett. 116 (2016), article no. 215301.

[3] F. Böttcher, J.-N. Schmidt, M. Wenzel, M. Hertkorn, M. Guo, T. Langen, T. Pfau, “Transient supersolid properties in an array of dipolar quantum droplets”, Phys. Rev. X 9 (2019), article no. 011051.

[4] L. Tanzi, E. Lucioni, F. Famà et al., “Observation of a dipolar quantum gas with metastable supersolid properties”, Phys. Rev. Lett. 122 (2019), article no. 130405.

[5] L. Chomaz, D. Petter, P. Ilzhöfer et al., “Long-lived and transient supersolid behaviors in dipolar quantum gases”, Phys. Rev. X 9 (2019), article no. 021012.

[6] L. Chomaz, I. Ferrier-Barbut, F. Ferlaino, B. Laburthe-Tolra, B. L. Lev, T. Pfau, “Dipolar physics: a review of experiments with magnetic quantum gases”, Rep. Prog. Phys. 86 (2023), article no. 026401.

[7] W. Zwerger, “Basic concepts and some current directions in ultracold gases”, in Lectures at the Collège de France, 2021, https://pro.college-de-france.fr/jean.dalibard/CdF/2021.

[8] J. Dalibard, “Magnetic interactions between cold atoms: Quantum droplets and supersolid states”, in Lectures at the Collège de France, 2024, https://pro.college-de-france.fr/jean.dalibard/CdF/2024.

[9] J.-P. Hansen, I. McDonald, Theory of Simple Liquids, 3rd ed., Elsevier, Academic Press, Amsterdam, 2006.

[10] P. M. Chaikin, T. C. Lubensky, Principles of Condensed Matter Physics, Cambridge University Press, Cambridge, 1995.

[11] P. Nozières, “Is the roton in superfluid helium 4 the ghost of a Bragg spot?”, J. Low Temp. Phys. 137 (2004), p. 45-67.

[12] J. L. Bohn, D. S. Jin, “Differential scattering and rethermalization in ultracold dipolar gases”, Phys. Rev. A 89 (2014), article no. 022702.

[13] J. L. Bohn, M. Cavagnero, D. Ticknor, “Quasi-universal dipolar scattering in cold and ultracold gases”, New J. Phys. 11 (2009), article no. 055039.

[14] A. Triay, “Existence of minimizers in generalized Gross–Pitaevskii theory with the Lee–Huang–Yang correction”, preprint, 2019, https://arxiv.org/abs/1904.10672.

[15] S. Tan, “Energetics of a strongly correlated Fermi gas”, Ann. Phys. 323 (2008), p. 2952-2970.

[16] E. Braaten, D. Kang, L. Platter, “Universal relations for identical Bosons from three-body physics”, Phys. Rev. Lett. 106 (2011), article no. 153005.

[17] D. Ruelle, “Classical statistical mechanics of a system of particles”, Helvetica Phys. Acta 36 (1963), p. 183-197.

[18] M. E. Fisher, “The free energy of a macroscopic system”, Arch. Rat. Mech. Anal. 17 (1964), p. 377-410.

[19] R. B. Griffiths, “Free energy of interacting magnetic dipoles”, Phys. Rev. 176 (1968), p. 655-659.

[20] J. Fröhlich, Y.-M. Park, “Correlation inequalities and the thermodynamic limit for classical and quantum continuous systems”, Comm. Math. Phys. 59 (1978), p. 235-266.

[21] S. Banerjee, R. B. Griffiths, M. Widom, “Thermodynamic limit for dipolar media”, J. Stat. Phys. 93 (1998), p. 109-141.

[22] J. D. Jackson, Classical Electrodynamics, 3rd ed., John Wiley, New York, 1999.

[23] L. Landau, E. Lifshitz, Electrodynamics of Continuous Media, Pergamon Press, New York, 1984.

[24] D. J. Griffiths, “Hyperfine splitting in the ground state of hydrogen”, Am. J. Phys. 50 (1982), p. 698-703.

[25] S. Yi, L. You, “Trapped condensates of atoms with dipole interactions”, Phys. Rev. A 63 (2001), article no. 053607.

[26] S. Giovanazzi, A. Görlitz, T. Pfau, “Ballistic expansion of a dipolar condensate”, J. Opt. B: Quantum Semiclass. Opt. 5 (2003), p. S208-S211.

[27] C. Eberlein, S. Giovanazzi, D. H. J. O'Dell, “Exact solution of the Thomas–Fermi equation for a trapped Bose–Einstein condensate with dipole–dipole interactions”, Phys. Rev. A 71 (2005), article no. 033618.

[28] M. S. Wertheim, “Exact solution of the mean spherical model for fluids of hard spheres with permanent electric dipole moments”, J. Chem. Phys. 55 (1971), p. 4291-4298.

[29] G. Bismut, B. Laburthe-Tolra, E. Maréchal, P. Pedri, O. Gorceix, L. Vernac, “Anisotropic excitation spectrum of a dipolar quantum Bose gas”, Phys. Rev. Lett. 109 (2012), article no. 155302.
[30] P. C. Martin, J. Schwinger, "Theory of many-particle systems. I", Phys. Rev. 115 (1959), p. 1342-1373.

[31] F. Werner, Y. Castin, "General relations for quantum gases in two and three dimensions. II. Bosons and mixtures", Phys. Rev. A 86 (2012), article no. 053633.

[32] T.-L. Ho, Q. Zhou, "Obtaining the phase diagram and thermodynamic quantities of bulk systems from the densities of trapped gases", Nat. Phys. 6 (2010), p. 131-134.

[33] P. C. Martin, O. Parodi, P. S. Pershan, "Unified hydrodynamic theory for crystals, liquid crystals, and normal fluids", Phys. Rev. A 6 (1972), p. 2401-2420.

[34] A. Andreev, I. Lifshitz, "Quantum theory of defects in crystals", Sov. Phys. JETP 29 (1969), p. 1107-1113.

[35] J. Hofmann, W. Zwerger, "Hydrodynamics of a superfluid smectic", J. Statist. Mech.: Theory Exp. (2021), no. 3, article no. 033104.

[36] N. Liebster, M. Sparn, E. Kath, J. Duchene, H. Strobel, M. K. Oberthaler, "Observation of supersolid sound modes in a driven quantum gas", preprint, 2024, University of Heidelberg.

[37] F. Werner, "Virial theorems for trapped cold atoms", Phys. Rev. A 78 (2008), article no. 025601.

[38] E. N. Parker, "Tensor virial equations", Phys. Rev. 96 (1954), p. 1686-1689.

[39] Y. Nishida, D. T. Son, "Nonrelativistic conformal field theories", Phys. Rev. D 76 (2007), article no. 086004.

[40] T. Koch, T. Lahaye, J. Metz, B. Fröhlich, A. Griesmaier, T. Pfau, "Stabilization of a purely dipolar quantum gas against collapse", Nat. Phys. 4 (2008), p. 218-222.

[41] J. Hofmann, W. Zwerger, "Universal relations for dipolar quantum gases", Phys. Rev. Res. 3 (2021), article no. 013088.

[42] D. Baillie, R. M. Wilson, P. B. Blakie, "Collective excitations of self-bound droplets of a dipolar quantum fluid", Phys. Rev. Lett. 119 (2017), article no. 255302.

[43] R. Balian, J. Blaizot, "Stars and statistical physics: A teaching experience", Am. J. Phys. 67 (1999), p. 1189-1206.

[44] J. W. Negele, H. Orland, *Quantum Many-Particle Systems, Advanced Book Classics*, Westview Press, Boulder, CO, 1998.

[45] C. Itzykson, J.-M. Drouffe, *Statistical Field Theory*, Academic Press, Springer, New York, 1991.

[46] A. Zee, *Quantum Field Theory in a Nutshell*, 2nd ed., Princeton University Press, Princeton, NJ, 2010.

[47] W. Zwerger, "Quantum-unbinding near a zero temperature liquid-gas transition", J. Statist. Mech.: Theory Exp. (2019), no. 10, article no. 103104.

[48] S. Tan, "Three-boson problem at low energy and implications for dilute Bose-Einstein condensates", Phys. Rev. A 78 (2008), article no. 013636.

[49] P. M. A. Mestrom, V. E. Colussi, T. Secker, S. J. J. M. F. Kokkelmans, "Scattering hypervolume for ultracold bosons from weak to strong interactions", Phys. Rev. A 100 (2019), article no. 050702.

[50] P. M. A. Mestrom, V. E. Colussi, T. Secker, G. P. Groeneveld, S. J. J. M. F. Kokkelmans, "van der Waals Universality near a quantum tricritical point", Phys. Rev. Lett. 124 (2020), article no. 143401.

[51] S. Zhu, S. Tan, "Three-body scattering hypervolumes of particles with short-range interactions", preprint, 2017, https://arxiv.org/abs/1710.04147.

[52] D. Baillie, R. M. Wilson, R. N. Bisset, P. B. Blakie, "Self-bound dipolar droplet: A localized matter wave in free space", Phys. Rev. A 94 (2016), article no. 021602.

[53] D. T. Son, M. Stephanov, H.-U. Yee, "The phase diagram of ultra quantum liquids", J. Statist. Mech.: Theory Exp. (2021), no. 1, article no. 013105.

[54] Y. Sekino, Y. Nishida, "Quantum droplet of one-dimensional bosons with a three-body attraction", Phys. Rev. A 97 (2018), article no. 011602.

[55] D. T. Son, M. Stephanov, H.-U. Yee, "Fate of multiparticle resonances: From Q-balls to $^3$He droplets", Phys. Rev. A 106 (2022), article no. L050801.

[56] S. Chandrasekhar, E. Fermi, "Problems of gravitational stability in the presence of a magnetic field", Astrophys. J. 118 (1953), p. 116-141.

[57] T. D. Lee, K. Huang, C. N. Yang, "Eigenvalues and eigenfunctions of a Bose system of hard spheres and its low-temperature properties", Phys. Rev. 106 (1957), p. 1135-1145.

[58] A. R. P. Lima, A. Pelster, "Quantum fluctuations in dipolar Bose gases", Phys. Rev. A 84 (2011), article no. 041604.

[59] M. D. Schwartz, *Quantum Field Theory and the Standard Model*, Cambridge University Press, Cambridge, 2014.

[60] M. Lüscher, "Volume dependence of the energy spectrum in massive quantum field theories", Commun. Math. Phys. 105 (1986), p. 153-188.

[61] R. Bombín, V. Cikojević, F. Mazzanti, J. Boronat, "Quantum-Monte-Carlo-based functional for dysprosium dipolar systems", Phys. Rev. A 109 (2024), article no. 033312.

[62] L. Santos, G. V. Shlyapnikov, M. Lewenstein, "Roton-maxon spectrum and stability of trapped dipolar Bose-Einstein condensates", Phys. Rev. Lett. 90 (2003), article no. 250403.

[63] D. H. J. O'Dell, S. Giovanazzi, G. Kurizki, "Rotons in gaseous Bose-Einstein condensates irradiated by a laser", Phys. Rev. Lett. 90 (2003), article no. 110402.

[64] U. R. Fischer, "Stability of quasi-two-dimensional Bose-Einstein condensates with dominant dipole-dipole interactions", Phys. Rev. A 73 (2006), article no. 031602.
[65] D. Petter, G. Natale, R. M. W. van Bijnen, A. Patscheider, M. J. Mark, L. Chomaz, F. Ferlaino, "Probing the roton excitation spectrum of a stable dipolar Bose gas", Phys. Rev. Lett. 122 (2019), article no. 183401.

[66] P. B. Blakie, D. Baillie, R. N. Bisset, "Roton spectroscopy in a harmonically trapped dipolar Bose-Einstein condensate", Phys. Rev. A 86 (2012), article no. 021604.

[67] J. Hertkorn, J.-N. Schmidt, F. Böttcher et al., "Density fluctuations across the superfluid-supersolid phase transition in a dipolar quantum gas", Phys. Rev. X 11 (2021), article no. 011037.

[68] A. J. Leggett, "Can a Solid Be "Superfluid"?", Phys. Rev. Lett. 25 (1970), p. 1543-1546.

[69] D. M. Ceperley, B. Bernu, "Ring exchanges and the supersolid phase of $^4$He", Phys. Rev. Lett. 93 (2004), article no. 155303.

[70] M. Boninsegni, A. B. Kuklov, L. Pollet, N. V. Prokof'ev, B. V. Svistunov, M. Troyer, "Fate of vacancy-induced supersolidity in $^4$He", Phys. Rev. Lett. 97 (2006), article no. 080401.

[71] M. Boninsegni, N. V. Prokof'ev, "Colloquium: Supersolids: What and where are they?", Rev. Mod. Phys. 84 (2012), p. 759-776.

[72] G. Chauveau, C. Maury, F. Rabec et al., "Superfluid fraction in an interacting spatially modulated Bose-Einstein condensate", Phys. Rev. Lett. 130 (2023), article no. 226003.

[73] N. Prokof'ev, B. Svistunov, "Supersolid state of matter", Phys. Rev. Lett. 94 (2005), article no. 155302.

[74] L. Landau, "Zur Theorie der Phasenumwandlungen II", Phys. Z. Sowjetunion 11 (1937), p. 26-35.

[75] S. Alexander, J. McTague, "Should all crystals be bcc? Landau theory of solidification and crystal nucleation", Phys. Rev. Lett. 41 (1978), no. 10, p. 702-705.

[76] J.-P. Hansen, L. Verlet, "Phase transitions of the Lennard-Jones system", Phys. Rev. 184 (1969), p. 151-161.

[77] E. P. Gross, "Classical theory of boson wave fields", Ann. Phys. 4 (1958), p. 57-74.

[78] Y. Pomeau, S. Rica, "Dynamics of a model of a supersolid", Phys. Rev. Lett. 72 (1994), no. 15, p. 2426-2429.

[79] C. Josserand, Y. Pomeau, S. Rica, "Coexistence of ordinary elasticity and superfluidity in a model of a defect-free supersolid", Phys. Rev. Lett. 98 (2007), article no. 195301.

[80] L. Pitaevskii, S. Stringari, Bose-Einstein Condensation and Superfluidity, Oxford University Press, Oxford, 2016.

[81] S. Brazovskii, "Phase transition of an isotropic system to a nonuniform state", Sov. Phys. JETP 41 (1975), p. 85-89.

[82] S. Brazovskii, I. Dzyaloshinskii, A. Muratov, "Theory of weak crystallization", Sov. Phys. JETP 66 (1987), p. 625-633.

[83] J. Hofmann, W. Zwerger, "Deep inelastic scattering on ultracold gases", Phys. Rev. X 7 (2017), article no. 011022.

[84] D. T. Son, M. Wingate, "General coordinate invariance and conformal invariance in nonrelativistic physics: Unitary Fermi gas", Ann. Phys. (NY) 321 (2006), p. 197-224.