Stability of quantum degenerate Fermi gases of tilted polar molecules

TL;DR

This paper develops a mean-field variational approach to analyze the stability and Fermi surface deformation of quantum degenerate dipolar Fermi gases of polar molecules, revealing universal stability conditions and geometry-dependent effects.

Contribution

It introduces a new theoretical framework based on the Wigner function to study ground-state properties of dipolar Fermi gases, highlighting the role of trap geometry and dipole orientation.

Findings

Stability depends only on trap aspect ratios and dipole orientation.

The Fermi surface deforms significantly in molecular systems compared to atomic species.

The stability diagram is species-independent and universal.

Abstract

A recent experimental realization of quantum degenerate gas of $^{40}$K$^{87}$Rb molecules opens up prospects of exploring strong dipolar Fermi gases and many-body phenomena arising in that regime. Here we derive a mean-field variational approach based on the Wigner function for the description of ground-state properties of such systems. We show that the stability of dipolar fermions in a general harmonic trap is universal as it only depends on the trap aspect ratios and the dipoles' orientation. We calculate the species-independent stability diagram and the deformation of the Fermi surface (FS) for polarized molecules, whose electric dipoles are oriented along a preferential direction. Compared to atomic magnetic species, the stability of a molecular electric system turns out to strongly depend on its geometry and the FS deformation significantly increases.

Equations

$$E_{\mathrm{tot}}=\frac{N}{8}\left(\sum _i\frac{{\hslash }^2{K}_i^2}{2M}+\sum _{i,j}\frac{M{\omega }_i^2{\mathbb{R}}_{ij}^{\prime 2}{R}_j^2}{2}\right)-\frac{6N^2{c}_0}{R_x{R}_y{R}_z}$$

$$\times \left[F_A\left(\frac{R_x}{R_z},\frac{R_y}{R_z},\theta ,\phi ,{\theta }^{\prime },{\phi }^{\prime }\right)-F_A\left(\frac{K_z}{K_x},\frac{K_z}{K_y},\theta ,\phi ,{\theta }^{\prime \prime },{\phi }^{\prime \prime }\right)\right].$$

$${\epsilon }_{\mathrm{dd}}=\frac{d^2}{4\pi {\epsilon }_0}\sqrt{\frac{M^3}{{\hslash }^5}}({\omega }_x{\omega }_y{\omega }_{z}N{)}^{1/6},$$

$$A_R(t)=\sqrt{\frac{⟨r_{\mathrm{v}}^2⟩}{⟨r_{\mathrm{h}}^2⟩}}=\frac{R_z{b}_z(t)}{\sqrt{R_x^2{b}_x^2(t){\sin }^2\gamma +R_y^2{b}_y^2(t){\cos }^2\gamma }},$$

$$\nu (\mathbf{r},\mathbf{k})=\Theta \left(1-\sum _{i,j}{r}_i{\mathbb{A}}_{ij}{r}_j-\sum _{i,j}{k}_i{\mathbb{B}}_{ij}{k}_j\right),$$

$$E_{\mathrm{tot}}=\frac{N}{8}\left(\sum _i\frac{{\hslash }^2{K}_i^2}{2M}+\sum _{i,j}\frac{M{\omega }_i^2{\mathbb{R}}_{ij}^{\prime 2}{R}_j^2}{2}\right)-\frac{6N^2{c}_0}{R_x{R}_y{R}_z}\left[F_A\left(\frac{R_x}{R_z},\frac{R_y}{R_z},\theta ,\phi ,{\theta }^{\prime },{\phi }^{\prime }\right)-F_A\left(\frac{K_z}{K_x},\frac{K_z}{K_y},\theta ,\phi ,{\theta }^{\prime \prime },{\phi }^{\prime \prime }\right)\right],$$

$$\mathbb{R}(\alpha,\beta)=\left(\begin{array}[]{ccc}\cos\alpha\cos\beta&-\sin\beta&\sin\alpha\cos\beta\\ \cos\alpha\sin\varphi&\cos\beta&\sin\alpha\sin\beta\\ -\sin\alpha&0&\cos\alpha\end{array}\right)\,.$$

$$F_A(x,y,\theta ,\phi ,\tilde{\theta },\tilde{\varphi })={\left(\sum _i{\mathbb{R}}_{iz}{\tilde{\mathbb{R}}}_{ix}\right)}^{2}f\left(\frac{y}{x},\frac{1}{x}\right)+{\left(\sum _i{\mathbb{R}}_{iz}{\tilde{\mathbb{R}}}_{iy}\right)}^{2}f\left(\frac{x}{y},\frac{1}{y}\right)+{\left(\sum _i{\mathbb{R}}_{iz}{\tilde{\mathbb{R}}}_{iz}\right)}^{2}f(x,y),$$

$$f\left(x,y\right)=1-3xy\frac{E(\varphi ,\kappa )-F(\varphi ,\kappa )}{(1-y^2)\sqrt{1-x^2}},$$

$$F_A(x,y,0,0,0,0)=f(x,y),F_A(x,y,\alpha ,\beta ,\alpha ,\beta )=f(x,y).$$

$$N-\frac{1}{48}{R}_x{R}_y{R}_z{K}_x{K}_y{K}_z$$

$$K_x-K_y$$

$$\frac{2{\hslash }^2{K}_x^2}{3}-\frac{{\hslash }^2{K}_y^2}{3}-\frac{{\hslash }^2{K}_z^2}{3}+\frac{48MN{c}_0}{R_x{R}_y{R}_z}\frac{K_z}{K_x}{\partial }_{K_x}{F}_A\left(\frac{K_z}{K_x},\frac{K_z}{K_y},\theta ,\phi ,{\theta }^{\prime \prime },{\phi }^{\prime \prime }\right)$$

$$\sum _i{\omega }_i^2{\mathbb{R}}_{ix}^{\prime 2}{R}_x^2-\sum _i\frac{{\hslash }^2{K_i}^2}{3M^2}-\frac{48N{c}_0}{M{R}_x{R}_y{R}_z}{F}_A\left(\frac{K_z}{K_x},\frac{K_z}{K_y},\theta ,\phi ,{\theta }^{\prime \prime },{\phi }^{\prime \prime }\right)$$

$$+\frac{48N{c}_0}{M{R}_x{R}_y{R}_z}\left[F_A\left(\frac{R_x}{R_z},\frac{R_y}{R_z},\theta ,\phi ,{\theta }^{\prime },{\phi }^{\prime }\right)-R_x{\partial }_{R_x}{F}_A\left(\frac{R_x}{R_z},\frac{R_y}{R_z},\theta ,\phi ,{\theta }^{\prime },{\phi }^{\prime }\right)\right]$$

$$\sum _i{\omega }_i^2{\mathbb{R}}_{iy}^{\prime 2}{R}_y^2-\sum _i\frac{{\hslash }^2{K_i}^2}{3M^2}-\frac{48N{c}_0}{M{R}_x{R}_y{R}_z}{F}_A\left(\frac{K_z}{K_x},\frac{K_z}{K_y},\theta ,\phi ,{\theta }^{\prime \prime },{\phi }^{\prime \prime }\right)$$

$$+\frac{48N{c}_0}{M{R}_x{R}_y{R}_z}\left[F_A\left(\frac{R_x}{R_z},\frac{R_y}{R_z},\theta ,\phi ,{\theta }^{\prime },{\phi }^{\prime }\right)-R_y{\partial }_{R_y}{F}_A\left(\frac{R_x}{R_z},\frac{R_y}{R_z},\theta ,\phi ,{\theta }^{\prime },{\phi }^{\prime }\right)\right]$$

$$\sum _i{\omega }_i^2{\mathbb{R}}_{iz}^{\prime 2}{R}_z^2-\sum _i\frac{{\hslash }^2{K_i}^2}{3M^2}-\frac{48N{c}_0}{M{R}_x{R}_y{R}_z}{F}_A\left(\frac{K_z}{K_x},\frac{K_z}{K_y},\theta ,\phi ,{\theta }^{\prime \prime },{\phi }^{\prime \prime }\right)$$

$$+\frac{48N{c}_0}{M{R}_x{R}_y{R}_z}\left[F_A\left(\frac{R_x}{R_z},\frac{R_y}{R_z},\theta ,\phi ,{\theta }^{\prime },{\phi }^{\prime }\right)-R_z{\partial }_{R_z}{F}_A\left(\frac{R_x}{R_z},\frac{R_y}{R_z},\theta ,\phi ,{\theta }^{\prime },{\phi }^{\prime }\right)\right]$$

$$\sum _{i,j}{\omega }_i^2{\mathbb{R}}_{ij}^{\prime }{\partial }_{{\theta }^{\prime }}{\mathbb{R}}_{ij}^{\prime }{R}_j^2-\frac{48N{c}_0}{M{R}_x{R}_y{R}_z}{\partial }_{{\theta }^{\prime }}{F}_A\left(\frac{R_x}{R_z},\frac{R_y}{R_z},\theta ,\phi ,{\theta }^{\prime },{\phi }^{\prime }\right)$$

$$\sum _{i,j}{\omega }_i^2{\mathbb{R}}_{ij}^{\prime }{\partial }_{{\phi }^{\prime }}{\mathbb{R}}_{ij}^{\prime }{R}_j^2-\frac{48N{c}_0}{M{R}_x{R}_y{R}_z}{\partial }_{{\phi }^{\prime }}{F}_A\left(\frac{R_x}{R_z},\frac{R_y}{R_z},\theta ,\phi ,{\theta }^{\prime },{\phi }^{\prime }\right)$$

$${\partial }_{{\theta }^{\prime \prime }}{F}_A\left(\frac{K_z}{K_x},\frac{K_z}{K_y},\theta ,\phi ,{\theta }^{\prime \prime },{\phi }^{\prime \prime }\right)$$

$${\partial }_{{\phi }^{\prime \prime }}{F}_A\left(\frac{K_z}{K_x},\frac{K_z}{K_y},\theta ,\phi ,{\theta }^{\prime \prime },{\phi }^{\prime \prime }\right)$$

$$1-R_x{R}_y{R}_z{K}_x{K}_y{K}_z$$

$$K_x-K_y$$

$$2K_x^2-K_y^2-K_z^2+\frac{3{\epsilon }_{\mathrm{dd}}{c}_{\mathrm{d}}}{R_x{R}_y{R}_z}\frac{K_z}{K_z}{\partial }_{K_x}f\left(\frac{K_z}{K_x},\frac{K_z}{K_y}\right)$$

$$\sum _i\frac{{\omega }_i^2}{{\omega }_x^2}{\mathbb{R}}_{ix}^{\prime 2}{R}_x^2-\frac{1}{3}\sum _i{K}_i^2-\frac{{\epsilon }_{\mathrm{dd}}{c}_{\mathrm{d}}}{R_x{R}_y{R}_z}f\left(\frac{K_z}{K_x},\frac{K_z}{K_y}\right)$$

$$+\frac{{\epsilon }_{\mathrm{dd}}{c}_{\mathrm{d}}}{R_x{R}_y{R}_z}\left[F_A\left(\frac{R_x{\omega }_z}{R_z{\omega }_x},\frac{R_y{\omega }_z}{R_z{\omega }_y},\theta ,\phi ,{\theta }^{\prime },{\phi }^{\prime }\right)-R_x{\partial }_{R_x}{F}_A\left(\frac{R_x{\omega }_z}{R_z{\omega }_x},\frac{R_y{\omega }_z}{R_z{\omega }_y},\theta ,\phi ,{\theta }^{\prime },{\phi }^{\prime }\right)\right]$$

$$\sum _i\frac{{\omega }_i^2}{{\omega }_y^2}{\mathbb{R}}_{iy}^{\prime 2}{R}_y^2-\frac{1}{3}\sum _i{K}_i^2-\frac{{\epsilon }_{\mathrm{dd}}{c}_{\mathrm{d}}}{R_x{R}_y{R}_z}f\left(\frac{K_z}{K_x},\frac{K_z}{K_y}\right)$$

$$+\frac{{\epsilon }_{\mathrm{dd}}{c}_{\mathrm{d}}}{R_x{R}_y{R}_z}\left[F_A\left(\frac{R_x{\omega }_z}{R_z{\omega }_x},\frac{R_y{\omega }_z}{R_z{\omega }_y},\theta ,\phi ,{\theta }^{\prime },{\phi }^{\prime }\right)-R_y{\partial }_{R_y}{F}_A\left(\frac{R_x{\omega }_z}{R_z{\omega }_x},\frac{R_y{\omega }_z}{R_z{\omega }_y},\theta ,\phi ,{\theta }^{\prime },{\phi }^{\prime }\right)\right]$$

$$\sum _i\frac{{\omega }_i^2}{{\omega }_z^2}{\mathbb{R}}_{iz}^{\prime 2}{R}_z^2-\frac{1}{3}\sum _i{K}_i^2-\frac{{\epsilon }_{\mathrm{dd}}{c}_{\mathrm{d}}}{R_x{R}_y{R}_z}f\left(\frac{K_z}{K_x},\frac{K_z}{K_y}\right)$$

$$+\frac{{\epsilon }_{\mathrm{dd}}{c}_{\mathrm{d}}}{R_x{R}_y{R}_z}\left[F_A\left(\frac{R_x{\omega }_z}{R_z{\omega }_x},\frac{R_y{\omega }_z}{R_z{\omega }_y},\theta ,\phi ,{\theta }^{\prime },{\phi }^{\prime }\right)-R_z{\partial }_{R_z}{F}_A\left(\frac{R_x{\omega }_z}{R_z{\omega }_x},\frac{R_y{\omega }_z}{R_z{\omega }_y},\theta ,\phi ,{\theta }^{\prime },{\phi }^{\prime }\right)\right]$$

$$\sum _{i,j}\frac{{\omega }_i^2}{{\omega }_j^2}{\mathbb{R}}_{ij}^{\prime }{\partial }_{{\theta }^{\prime }}{\mathbb{R}}_{ij}^{\prime }{R}_j^2-\frac{{\epsilon }_{\mathrm{dd}}{c}_{\mathrm{d}}}{R_x{R}_y{R}_z}{\partial }_{{\theta }^{\prime }}{F}_A\left(\frac{R_x{\omega }_z}{R_z{\omega }_x},\frac{R_y{\omega }_z}{R_z{\omega }_y},\theta ,\phi ,{\theta }^{\prime },{\phi }^{\prime }\right)$$

$$\sum _{i,j}\frac{{\omega }_i^2}{{\omega }_j^2}{\mathbb{R}}_{ij}^{\prime }{\partial }_{{\phi }^{\prime }}{\mathbb{R}}_{ij}^{\prime }{R}_j^2-\frac{{\epsilon }_{\mathrm{dd}}{c}_{\mathrm{d}}}{R_x{R}_y{R}_z}{\partial }_{{\phi }^{\prime }}{F}_A\left(\frac{R_x{\omega }_z}{R_z{\omega }_x},\frac{R_y{\omega }_z}{R_z{\omega }_y},\theta ,\phi ,{\theta }^{\prime },{\phi }^{\prime }\right)$$

Full text

The Fermi surface (FS) is one of the fundamental pillars of modern condensed matter physics Ashcroft and Mermin (1976). It represents the surface in reciprocal space, which separates occupied from unoccupied fermionic states at zero temperature, and is a direct consequence of the Pauli exclusion principle. For instance, interacting electrons in a normal metal can be described within the Landau Fermi-liquid theory Lifshitz and Pitaevskii (1980) as noninteracting fermionic quasi-particles with an effective mass, whose ground state forms such a FS. Due to the isotropy of the Coulomb repulsion between electrons in a uniform space, the FS turns out to be a sphere, whose radius is given by the Fermi momentum. The concept of the FS is crucial for understanding transport processes in metals Kontani (2014) and the Cooper pairing in superconductors Tinkham (1975); Bennemann and (Eds.). However, in case of complex interactions the FS can get modified. For example, in strongly-correlated electron systems the Fermi-liquid picture breaks down, giving rise to a spontaneous breaking of rotational invariance, which manifests itself in a deformation of the FS Fradkin et al. (2010).

Studying Fermi surfaces has now also become accessible within the realm of ultracold quantum gases Bloch et al. (2008); Giorgini et al. (2008); Lahaye et al. (2009); Blume (2012); Massignan et al. (2014) due to their high degree of tunability. In Fermi gases consisting of atoms or molecules with a permanent or induced magnetic or electric dipole moment the anisotropic and long-range dipole-dipole interaction (DDI) competes with the large kinetic energy close to the FS Baranov (2008). As a consequence, many theoretical papers predicted an anisotropic version of the Landau Fermi-liquid theory Chan et al. (2010); Sieberer and Baranov (2011); Baranov et al. (2012), which involves a deformation of the Fermi sphere Góral et al. (2001); Miyakawa et al. (2008); Sogo et al. (2009); Lima and Pelster (2010a, b); Baillie et al. (2015). A recent experiment Aikawa et al. (2014); Veljić et al. (2018) measured that for a fermionic gas of magnetic dipolar erbium atoms an ellipsoidal deformation of the Fermi sphere occurs, which is of the order of 2%. This is expected to lead to novel many-body phenomena, in particular in connection with fermionic superfluidity Baranov et al. (2004a, b); Zhao et al. (2010); Shi et al. (2010); Trautmann et al. (2018); Mulkerin et al. (2019). In a polarized one-component Fermi gas an intriguing interplay between an anisotropic order parameter with odd partial waves and the FS deformation enhances superfluid pairing via modifying the density of states Zhao et al. (2010). In contrast to that the more conventional type of Cooper pairing is predicted in a two-component dipolar Fermi gas, where the usual BCS theory together with the deformed FS leads to both spin-singlet even partial wave or spin-triplet odd partial wave Cooper pairs Shi et al. (2010). And it is suggested to obtain and observe a topological $p$-wave superfluid of microwave-dressed polar fermionic molecules in 2D lattices at temperatures of the order of tens of nanokelvins Fedorov et al. (2016).

Since the first experimental realization of a quantum degenerate dipolar Fermi gas of 161Dy in 2012 Lu et al. (2012), several more fermionic species, such as 167Er Aikawa et al. (2014) and 53Cr Naylor et al. (2015), were successfully cooled down to quantum degeneracy, which enabled studies of the effects of weak to medium-range DDI strength. However, the study of the strongly dipolar regime is still in its infancy, and awaits experimental availability of ultracold heteronuclear polar molecules with large dipole moments. In the last decade, significant efforts to produce chemically stable cold polar molecules Carr and Ye (2009); Gadway and Yan (2016) were based on photoassociation or the stimulated Raman adiabatic passage (STIRAP) Vitanov et al. (2017). As a result, samples of fermionic 40K87Rb Ni et al. (2008), 23Na40K Park et al. (2015a, b); Gempel et al. (2016); Seeßelberg et al. (2018), 23Na6Li Rvachov et al. (2017) and bosonic 7Li133Cs Deiglmayr et al. (2008, 2010), 87Rb133Cs Takekoshi et al. (2014); Molony et al. (2014) and 23Na87Rb Guo et al. (2016) were obtained in deeply bound molecular states. However, the quantum degeneracy was still not reached. Only very recently a quantum degenerate dipolar Fermi gas of 40K87Rb has been realized at JILA De Marco et al. (2019). This experimental protocol enabled to produce tens of thousands of unpolarized molecules at a temperature as low as 50 nK, which are well described by the Fermi-Dirac distribution. However, the molecules’ dipoles can be straight-forwardly polarized in a preferential direction by an external electric field De Marco et al. (2019), such that the DDI dominates the behavior of the system. This would be a long-awaited significant step forward, which would open up the realm for experimentally investigating strong dipolar Fermi gases.

The stability of quantum degenerate dipolar Fermi gases against mechanical collapse is defined by the positivity of the compressibility due to the Pomeranchuk criterion, which is a special case of the well-known criterion of thermodynamic stability Pomeranchuk (1958); Lifshitz and Pitaevskii (1980). It was previously considered in cylindrically symmetric harmonic traps Góral et al. (2001); Miyakawa et al. (2008); Sogo et al. (2009); Lima and Pelster (2010a, b); Baillie et al. (2015), as well as in homogeneous systems Ronen and Bohn (2010). Here we study the ground state stability of ultracold Fermi gases with tilted dipoles in triaxial harmonic traps and reveal a universal behavior of the critical DDI strength. In particular, we investigate the stability of a polarized 40K87Rb gas in an experimentally realistic parameter regime and calculate critical values of the electric dipole moment and the corresponding FS deformation. Note that, in contrast to non-dipolar systems, the 3-body recombination of dipolar atoms presumably does not play an important role in determining the stability of the system in the parameter range of current experiments De Marco et al. (2019). Finally, we demonstrate effects of the strong DDI on the time-of-flight expansion and show that a nonballistic expansion theory is essential to understand the dynamical behavior of strongly dipolar Fermi gases.

To achieve this, we use a variational phase-space approach Wächtler et al. (2017); Veljić et al. (2017, 2018) for the Wigner function $\nu({\mathbf{r}},\,{\mathbf{k}})=\int d^{3}r^{\prime}\,e^{-i{\mathbf{k}}\cdot{\mathbf{r}}^{\prime}}\,\rho\left({\mathbf{r}}+\frac{1}{2}{\mathbf{r}}^{\prime},\,{\mathbf{r}}-\frac{1}{2}{\mathbf{r}}^{\prime}\right)$, which relies on the Hartree-Fock mean-field approximation. Here $\rho({\mathbf{r}},\,{\mathbf{r}}^{\prime})=\langle\hat{\Psi}({\mathbf{r}})\hat{\Psi}^{\dagger}({\mathbf{r}}^{\prime})\rangle$ represents the one-body density matrix. Note that the second-order terms in the DDI in the theory beyond Hartree-Fock Krieg et al. (2015); Lange et al. (2016) yield only a small correction even for polar molecules, although the geometry may have an impact (see Supplemental Material Sup for more details). Furthermore, this beyond-mean-field correction turns out to destabilize the system Krieg et al. (2015); Lange et al. (2016), so our results on the stability represent proper upper boundaries. This is in stark contrast to bosonic systems, where the quantum fluctuations have turned out to stabilize the system and lead, for instance, to the formation of quantum droplets Lima and Pelster (2011); Kadau et al. (2016); Ferrier-Barbut et al. (2016); Wächtler and Santos (2016); Chomaz et al. (2016); Schmitt et al. (2016) and supersolids Chomaz et al. (2019); Tanzi et al. (2019a); Böttcher et al. (2019); Tanzi et al. (2019b); Guo et al. (2019); Natale et al. (2019) in Bose-Einstein condensates of dysprosium and erbium in Stuttgart, Innsbruck, and Pisa experiments.

We consider the dipolar Fermi system to be at zero temperature. This is justified as the temperature in the experiments De Marco et al. (2019) is about $T/T_{F}\approx 0.3$ and as thermal corrections to the total energy are proportional to $(T/T_{F})^{2}$ Howe et al. (2009). Thus, we assume a Heaviside-shaped Wigner function in the ground state and obtain the total energy of $N$ identical fermions of mass $M$ in terms of the Thomas-Fermi (TF) momenta $K_{i}$ and radii $R_{i}$,

[TABLE]

Here the angles ($\theta,\varphi$) determine the dipoles’ orientation, ($\theta^{\prime},\varphi^{\prime}$) describe the orientation of the molecular cloud, and ($\theta^{\prime\prime},\varphi^{\prime\prime}$) determine the FS orientation, as illustrated in Fig. 1. We stress that the molecular cloud orientation does not have to coincide with the trap orientation, due to the DDI effects, and $\mathbb{R}^{\prime}$ stands for the corresponding rotation matrix, while $\omega_{i}$ denote the trap frequencies and $F_{\rm A}$ is a generalized anisotropy function (see Supplemental Material Sup for further details). Note that with this our theory not only allows a quantitative analysis of current experimental data. In addition, it also provides a framework for disentangling reliably the rigid rotation of the Fermi ellipsoid, which occurs for the weak DDI in magnetic dipolar atoms in Innsbruck Aikawa et al. (2014); Veljić et al. (2018), from its deformation in case of the strong DDI in electric dipolar molecules investigated at JILA De Marco et al. (2019).

By extremizing the above energy with respect to the variational parameters ($R_{i},K_{i},\theta^{\prime},\varphi^{\prime},\theta^{\prime\prime},\varphi^{\prime\prime}$), we obtain the corresponding equations for the ground state, which can be rewritten in a dimensionless, species-independent form Sup such that, for a given orientation of the dipoles, they only depend on three parameters: the two trap aspect ratios $\omega_{z}/\omega_{x}$ and $\omega_{z}/\omega_{y}$, and the relative DDI strength

[TABLE]

where $d$ denotes the electric dipole moment. This remarkable result reveals a universality governing the ground-state properties of dipolar Fermi gases. Furthermore, it allows us to determine the stability diagram of the system, shown in Fig. 2(a) for the case $\theta=\varphi=0$, in terms of the maximal DDI strength $\varepsilon_{\mathrm{dd}}^{\mathrm{crit}}$ for which the ground state exists. We see that large aspect ratios significantly increase the critical DDI strength, thus stabilizing the system in a much broader parameter range. As an immediate consequence we read off from Fig. 2(a), for instance, that a dipolar Fermi gas can be stabilized against mechanical collapse, which arises in 3D for a sufficiently strong interaction, by confining the polar molecules to 2D, i.e., to a monolayer. Note that using the tilting angle of the dipole orientation relative to the monolayer and the DDI as control parameters, one can find apart from the normal Fermi liquid and the collapse also a superfluid phase and a density-wave phase Sieberer and Baranov (2011); Baranov et al. (2012).

We also note that $\varepsilon_{\mathrm{dd}}^{\mathrm{crit}}$ turns out to be a symmetric function of its arguments $\omega_{z}/\omega_{x}$ and $\omega_{z}/\omega_{y}$ Sup . If we consider the experimentally available species 40K87Rb, the stability diagram from Fig. 2(a) can be used to obtain a species-specific stability diagram for a particular value of one of the trap frequencies, as shown in Figs. 2(b) and 2(c). Here we see how the critical value of the dipole moment $d^{\mathrm{crit}}$ depends on $\omega_{x}$ and $\omega_{y}$ for a fixed value of $\omega_{z}$. If we take into account that the permanent electric dipole moment of 40K87Rb has the value $d=0.574$ D, denoted by black lines in Figs. 2(b) and 2(c), we read off that for $\omega_{z}=2\pi\times 50$ Hz the instability can kick in already for frequencies $\omega_{x}$, $\omega_{y}$ of that order or larger. In the experiment of Ref. De Marco et al. (2019) the frequencies used are $(\omega_{x},\,\omega_{y},\,\omega_{z})=2\pi\times(63,\,36,\,200)$ Hz, and Fig. 2(c) reveals that the system may easily become unstable for slightly larger frequencies if the dipoles would be polarized along $z$ axis.

The most striking effect that can be demonstrated in the strong DDI regime is the FS deformation $\Delta=K_{z}/K_{x}-1$, defined in terms of the TF momenta aspect ratio. It was experimentally observed for the first time for magnetic dipolar 167Er atoms Aikawa et al. (2014), where $\Delta$ is of the order of a few percent. This effect is much larger in gases of polar molecules, as can be seen for 40K87Rb in Fig. 3(a), obtained by solving the equations presented in the Supplemental Material Sup . Here it is assumed that the electric dipole moment is tuned down to $d=0.22$ D, such that it is below the minimal value of $d^{\mathrm{crit}}=0.24$ D obtained in Fig. 2(c). For realistic values of the trapping frequencies we obtain that $\Delta$ varies between 5% and 30%.

Furthermore, the theory presented here makes it possible to calculate the stability properties for experimentally relevant dipolar Fermi systems, where even relatively small changes in the dipolar moment strength can significantly affect the system’s stability. This is demonstrated in Fig. 3(b), where for a slightly larger value of $d=0.26$ D we read off that the FS deformation becomes significantly larger than in Fig. 3(a), namely up to 45%, and that an unstable region appears, which does not support a stable ground state of the system.

Previously it was always assumed Góral et al. (2001); Miyakawa et al. (2008); Sogo et al. (2009); Lima and Pelster (2010a, b); Baillie et al. (2015); Wächtler et al. (2017); Veljić et al. (2017, 2018) that the cloud shape in real space follows the trap orientation, which is a reasonable approximation for a weak DDI and an elongated trap. However, here we provide a theory capable to describe dipolar Fermi systems in a general trap geometry for all DDI strengths. Therefore, we leave the orientation angles ($\theta^{\prime},\varphi^{\prime}$) of the TF ellipsoid in real space as free variational parameters, together with the TF radii $R_{i}$, as illustrated in Fig. 1(b). It was previously experimentally verified and theoretically always assumed Aikawa et al. (2014); Veljić et al. (2018) that the FS stretches into an ellipsoid along the orientation of the dipoles, as can be expected due to symmetry reasons. Here we use a more general ansatz, where the FS orientation angles ($\theta^{\prime\prime},\varphi^{\prime\prime}$) are also taken as free variational parameters, together with the TF momenta $K_{i}$, Fig. 1(b). However, we show here Sup that the principle of minimizing the energy leads to the solution $\theta^{\prime\prime}=\theta$, $\varphi^{\prime\prime}=\varphi$ Sup , i.e., that our theory confirms the notion that the FS follows the dipoles’ orientation and properly captures the physical behavior of the system. This result also allows us to reduce the number of variational parameters to eight, ($R_{i}$, $K_{i}$, $\theta^{\prime}$, $\varphi^{\prime}$), as well as the number of equations Sup .

The cloud orientation obtained within our theory strongly depends on both the DDI strength and the elongation of the trap. In the special case of a spherical trap the cloud is elongated along the dipoles’ direction, as the FS, but in a general case the cloud orientation can only be determined numerically. Figure 4(a) shows the angular stability diagram for a 40K87Rb system in terms of the critical dipole moment $d^{\mathrm{crit}}$, where all variational parameters ($R_{i}$, $K_{i}$, $\theta^{\prime}$, $\varphi^{\prime}$) are numerically calculated for each configuration. If one would assume that the molecular cloud follows the trap shape, i.e., $\theta^{\prime}=\varphi^{\prime}=0$, the obtained values of $d^{\mathrm{crit}}$ would be significantly different from those calculated in Fig. 4(a) (for a comparison, see Supplemental Material Sup ). This demonstrates that the theory developed here is important for an accurate qualitative and quantitative description of dipolar Fermi systems with moderate to strong DDI. Figures 4(b) and 4(c) illustrate that the trap geometry also strongly affects the system’s behavior, and that the FS deformation and its angular distribution can be tuned by changing the trap frequencies. Not only the range of the FS deformation values can be increased or decreased this way, but also its minima and maxima can be freely modified. The observed strong angular dependence of the FS deformation has an important consequence, namely, that the FS does not only follow the dipoles’ orientation, but its shape gets modified as well. This is a qualitatively different behavior compared to atomic magnetic species, where the angular dependence of the FS deformation is quite weak Veljić et al. (2018), thus the FS just rigidly follows the dipoles’ orientation.

Effects of the DDI and their interplay with the geometry also quite strongly influence the dynamics of the system, which is of particular importance for interpreting experimental time-of-flight (TOF) imaging data Veljić et al. (2017, 2018). For polar molecules with a strong DDI the difference between the usually assumed ballistic and the actual nonballistic expansion can be huge, which we show here. The TOF expansion imaging is commonly used for experimental measurements of the properties of ultracold Fermi gases, and the deformation of the cloud shape is described in terms of the cloud aspect ratio $A_{R}(t)$, which is defined by the ratio of the average sizes of the cloud in vertical $\sqrt{\langle r^{2}_{\mathrm{v}}\rangle}$ and horizontal $\sqrt{\langle r^{2}_{\mathrm{h}}\rangle}$ direction in the imaging plane. Since the imaging axis in the experiment of Ref. De Marco et al. (2019) lies in the $xy$ plane and forms an angle $\gamma=22.5^{\circ}$ with respect to the $x$ axis, according to Ref. Veljić et al. (2017, 2018) the aspect ratio is given by

[TABLE]

where the scaling parameters $b_{i}(t)$ represent variations from the global equilibrium values of the TF radii and momenta Pedri et al. (2003). A detailed derivation of equations of motion for the scaling parameters based on the quantum Boltzmann equation within the self-consistent relaxation-time approach for $\theta=\varphi=0$ is given in Ref. Veljić et al. (2017). Here we numerically solve these equations for the general triaxial trap geometry and the parameters corresponding to the polar molecules of Ref. De Marco et al. (2019).

It was previously shown that, even for magnetic atomic species such as erbium, the DDI effects could be experimentally observed in the TOF dynamics, and that a nonballistic expansion has to be used in order to properly describe the system’s behavior Veljić et al. (2017, 2018). For polar molecules with a strong DDI we expect that nonballistic effects are more pronounced, as can be read off from Fig. 5. Even more significant are large variations of nonballistic effects, which can be as small as 8% or as large as 60% for quite similar configurations, as is illustrated for the two examples in Fig. 5(a). Although the trap geometry plays a role here, Fig. 5(b) reveals that the ballistic behavior is roughly the same, as is expected based on the system parameters, while the DDI strength gives a major contribution. Furthermore, the inset in Fig. 5(a) shows that even the qualitative behavior of the system can be incorrectly predicted (monotonous vs. nonmonotonous behavior) when nonballistic effects are neglected. This demonstrates that the DDI has to be taken into account even during the TOF expansion, and that the interpretation of experimental data are hugely affected by the model used. Therefore, the general theory presented here also enables an accurate modeling of the dynamics of strongly interacting dipolar Fermi systems.

In conclusion, we have presented a general mean-field theory for the ground state of polarized, harmonically trapped dipolar Fermi gases at zero-temperature, with an arbitrary orientation of the dipoles. We have derived a universal, species-independent set of equations for the ground state and investigated the stability of systems of polar molecules. We have shown that the molecular cloud shape and the FS deformation strongly depend on the dipoles’ orientation. Our results are important for the study of the interplay between the FS deformation and superfluid pairing Baranov et al. (2004a, b); Zhao et al. (2010); Shi et al. (2010), in particular to address the open question of how the anisotropic order parameter of the emergent superfluidity and its critical temperature are tunable by both the trap geometry and the dipoles’ orientation. The presented theory paves the way towards new methods for quantum engineering of properties of dipolar Fermi gases that depend on the FS shape, such as the emergence of superfluidity. In the outlook we mention that possible fermionic analogues of the already observed bosonic quantum droplets may exist. Although the leading beyond-mean-field correction seems to destabilize fermionic systems, further studies might reveal an alternative stabilization mechanism.

Acknowledgements. We acknowledge inspiring discussions with J. P. Covey, A. R. P. Lima, L. De Marco, K. Matsuda, W. Tobias, G. Valtolina, I. Vasić, and J. Ye. This work was supported in part by the Ministry of Education, Science, and Technological Development of the Republic of Serbia under projects ON171017 and QDDB, by the German Academic and Exchange Service (DAAD) under project QDDB, and by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) via the Collaborative Research Centers SFB/TR185 (project No. 277625399) and SFB/TR49, as well as the Research Unit FOR 2247 (project No. PE 530/6-1). Numerical simulations were run on the PARADOX supercomputing facility at the Scientific Computing Laboratory of the Institute of Physics Belgrade.

Appendix A Global equilibrium

We consider an ultracold quantum degenerate dipolar Fermi gas at zero temperature in global equilibrium. The system consists of $N$ identical spin-polarized single-component fermions of mass $M$ with an electric dipole moment ${\bf d}$. The system is confined into a three-dimensional harmonic trap with the frequencies $\omega_{i}$, whose axes coincide with the axes of the laboratory coordinate system $S$, as it is depicted in Fig. 1. The dipole moments are aligned in the direction defined by spherical angles $(\theta,\varphi)$, as shown in Fig. 1(b). For the ideal Fermi gas, the molecular cloud shape is determined purely by the trap, while the FS is a sphere, Fig. 1(a). In contrast to this, when the DDI is present, the molecular cloud shape will depend on both the trap geometry and the orientation of the dipoles. Here we assume that it has an ellipsoidal shape, oriented along the direction defined by the angles $(\theta^{\prime},\varphi^{\prime})$, as depicted in Fig. 1(b). Similarly, the FS is stretched into an ellipsoid whose longitudinal axis is expected to coincide with the dipoles’ orientation, but we do not assume it from the beginning and will instead show that this can be obtained within the presented theory. Therefore, as depicted in Fig. 1(b), the FS ellipsoid is oriented in the direction defined by the angles $(\theta^{\prime\prime},\varphi^{\prime\prime})$, which are free parameters.

We use the following variational ansatz for the Wigner distribution function:

[TABLE]

where $\Theta$ represents the Heaviside step function, while ${\mathbb{A}}_{ij}$ and $\mathbb{B}_{ij}$ are matrix elements that account for the geometry of the system and determine the shape of the cloud in real space and of the FS in momentum space, as illustrated in Fig. 1(b). With this ansatz we obtain the total energy of the system in the Hartree-Fock approximation:

[TABLE]

where $R_{i}$ and $K_{i}$ are the TF radii and momenta, respectively, $c_{0}=2^{10}d^{2}/(\varepsilon_{0}\cdot 3^{4}\cdot 5\cdot 7\cdot\pi^{3})$ is a constant related to the DDI strength, and ${\mathbb{R}}^{\prime}_{ij}$ are matrix elements of the rotation matrix $\mathbb{R}^{\prime}=\mathbb{R}(\theta^{\prime},\varphi^{\prime})$, with

[TABLE]

The features of the DDI are embodied into the generalized anisotropy function $F_{A}(x,y,\theta,\varphi,\tilde{\theta},\tilde{\phi})$, which includes the dependence on the general orientation of the dipoles $\bf d$ and the corresponding TF ellipsoid:

[TABLE]

where $\mathbb{R}_{ij}$ and $\tilde{\mathbb{R}}_{ij}$ are matrix elements of the rotation matrix $\mathbb{R}=\mathbb{R}(\theta,\varphi)$ and $\tilde{\mathbb{R}}=\mathbb{R}(\tilde{\theta},\tilde{\phi})$, respectively. In the above definition, $f(x,y)$ stands for the well-known anisotropy function Giovanazzi et al. (2006)

[TABLE]

where $\phi=\arccos x$ and $\kappa^{2}=(1-y^{2})/(1-x^{2})$, while $F(\phi,k)$ and $E(\phi,k)$ are the elliptic integrals of the first and of the second kind, respectively. Note that in the two relevant limiting cases the function $F_{A}$ satisfies

[TABLE]

Due to these identities and the fact that $f(x,y)$ is a symmetric function, the obtained distributions of $\varepsilon_{\mathrm{dd}}^{\mathrm{crit}}$ and $d^{\mathrm{crit}}$ in Fig. 2, as well as the distributions of $\Delta$ in Fig. 3 are symmetric with respect to their arguments for $\theta=\varphi=0$. Note that the above definition of the generalized anisotropy function enables symmetric treatment of both the Hartree and the Fock term in the expression for the total energy (S2).

The dipolar Fermi system is determined by the 10 variational parameters ($R_{i},K_{i},\theta^{\prime},\varphi^{\prime},\theta^{\prime\prime},\varphi^{\prime\prime}$), which are obtained by minimizing the total energy (S2), under the constraint that the total number of particles is equal to $N$. This leads to the following set of algebraic equations:

[TABLE]

From Eq. (S8) we see that the FS remains cylindrically symmetric even in the case of a general orientation of the dipoles with respect to the trap.

Appendix B Orientation of the Fermi surface

Equations (S15) and (S16) can be solved analytically, independently of other equations, yielding the physically expected result $\theta^{\prime\prime}=\theta$, $\varphi^{\prime\prime}=\varphi$. This means that the FS stretches along the dipoles’ orientation, as it was verified both experimentally and theoretically for the atomic erbium gas Veljić et al. (2018). Here we obtain this result self-consistently within our approach, which demonstrates that ansatz (S1) properly captures the ground-state properties of dipolar Fermi gases.

Appendix C Dimensionless form of equations for the ground state

If we eliminate the angles $\theta^{\prime\prime},\varphi^{\prime\prime}$ as outlined above and set $\theta^{\prime\prime}=\theta$, $\varphi^{\prime\prime}=\varphi$ in all equations, the system is now determined by the 8 variational parameters ($R_{i},K_{i},\theta^{\prime},\varphi^{\prime}$), which are obtained by solving the set of equations (S7)–(S14). They can be transformed into a dimensionless form by expressing the TF radii $R_{i}$ and momenta $K_{i}$ in units of $R_{i}^{0}=\sqrt{\frac{2E_{\rm F}}{M\omega_{i}^{2}}}$ and $K_{\rm F}=\sqrt{\frac{2ME_{\rm F}}{\hbar^{2}}}$, respectively, where $E_{\rm F}=\hbar(6N\omega_{x}\omega_{y}\omega_{z})^{1/3}$ stands for the Fermi energy. The quantities $R_{i}^{0}$ and $K_{\rm F}$ represent the TF radii and Fermi momentum of the ideal Fermi gas, illustrated in Fig. 1(a). The dimensionless radii and momenta are defined by $\tilde{R}_{i}=R_{i}/R_{i}^{0}$ and $\tilde{K}_{i}=K_{i}/K_{\rm F}$, and if we drop, for simplicity, the tilde signs, the set of equations (S7)-(S14) reduces to

[TABLE]

where $\varepsilon_{\mathrm{dd}}=\frac{d^{2}}{4\pi\varepsilon_{0}}\sqrt{\frac{M^{3}}{\hbar^{5}}}(\omega_{x}\omega_{y}\omega_{z}N)^{1/6}$ is the dimensionless relative DDI strength and $c_{\rm d}=\frac{2^{\frac{38}{3}}}{3^{\frac{23}{6}}\cdot 5\cdot 7\cdot\pi^{2}}$ is a number.

Appendix D Beyond-mean-field corrections

Here we estimate beyond-mean-field effects in the calculation of the Fermi surface shape and the stability of the system for strong dipolar interaction. We follow Ref. Krieg et al. (2015), which derives beyond-mean-field corrections to both the Fermi surface deformation and compressibility of the system. Note that this reference considers a homogeneous system, and that the estimates based on these results might not be fully applicable to a trapped system. However, the corresponding results for a trapped system are not available, and therefore we use Ref. Krieg et al. (2015) to estimate beyond-mean-field corrections in our case. In order to do so, we identify the homogeneous density of Ref. Krieg et al. (2015) with the average density of the trapped system, calculated as $N/V$, where $V=\frac{4\pi}{3}R_{x}R_{y}R_{z}$ is a volume of the TF ellipsoid in real space for specific parameters in the experimental setup.

Having this in mind, we first use Eq. (44) from Ref. Krieg et al. (2015) to estimate beyond-mean-field corrections to the Fermi surface deformation. This is illustrated in Fig. S1(a) for experimental system parameters De Marco et al. (2019) with $d=0.25$ D, which are used to obtain Fig. 4(b). It turns out that corrections are just a fraction of one percent. In Fig. S1(b) we see how the beyond-mean-field correction depends on the dipole moment. It amounts to a few percent even for the strongest values of $d$ that can be achieved in current experiments with 40K87Rb De Marco et al. (2019).

However, the situation is more complex when we consider the bulk modulus, i.e., the inverse compressibility $K$ of the system, which is used to estimate the stability border according to the Pomeranchuk criterion Pomeranchuk (1958); Lifshitz and Pitaevskii (1980), and whose beyond-mean-field correction is given by Eq. (46) in Ref. Krieg et al. (2015). For instance, for system parameters used to obtain Fig. 4(b), the corresponding second-order correction to the inverse compressibility $\delta K^{(2)}$ is of the order of one percent, as can be seen in Fig. S2(a), where we plot its angular dependence. These corrections are calculated for the dipole moment value $d=0.25$ D. In Fig. S2(b) we see, however, that the correction can be much higher for larger values of $d$, and that it strongly depends on the orientation of the dipoles. If we use a 10% threshold for the inverse compressibility correction, we see that $d$ can be as high as 0.35 D in the worst-case scenario, when the dipoles lie within the pancake plane, while for other values of the angles one can use even larger values of $d$. Taking into account that this coincides with the maximal achievable dipole moment in the current experiment with 40K87Rb De Marco et al. (2019), for their trap configuration our mean-field theory is applicable with reasonable accuracy, as shown in Fig. S2. However, for other trap configurations the mean-field theory could break down for smaller values of $d$, as the inverse compressibility $K$ of the system can have a strong angular dependency, depending on the underlying trap geometry. One can use a similar calculation as the one presented here to make an appropriate estimate for any given trap configuration.

Appendix E Comparison of stability diagrams obtained with self-consistent and fixed cloud orientation

The orientation of the dipoles with respect to the harmonic trap affects not only the shape of the molecular cloud in real space, but also its orientation. While previously it was always assumed that the axes of the molecular cloud coincide with the axes of the trap ($\theta^{\prime}=\varphi^{\prime}=0$), our theory takes into account the effects of the DDI and determines the angles $\theta^{\prime},\varphi^{\prime}$ in a self-consistent manner. In Fig. S3 we compare stability diagrams, expressed in terms of the critical electric dipole moment $d^{\mathrm{crit}}$, obtained in panel (a) by our theory and in panel (b) by assuming $\theta^{\prime}=\varphi^{\prime}=0$. The numerically calculated angular distributions are markedly different and the stability region is reduced when the full theory is applied. This makes the approach presented here important for the design of new experiments with polar molecules, in particular in the strong DDI regime.