Stability of quantum degenerate Fermi gases of tilted polar molecules

A recent experimental realization of quantum degenerate gas of $^{40}$K$^{87}$Rb molecules opens up prospects of exploring strongly 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.

A recent experimental realization of quantum degenerate gas of 40 K 87 Rb molecules opens up prospects of exploring strongly 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.
Ultracold gases at quantum degeneracy offer a wide playground for studying many-body physics [1][2][3][4][5][6][7], with important applications for the quantum simulation of a diverse range of systems and models [8][9][10][11][12][13][14][15][16], as well as for quantum information and quantum computing [17][18][19]. While usually the s-wave scattering is the main type of interaction in such systems, the presence of the anisotropic and long-range dipole-dipole interaction (DDI) in gases consisting of atoms or molecules with a permanent or induced magnetic or electric dipole moment leads to an even richer phenomena landscape. In particular, this is enabled by competing interaction effects and a high degree of their tunability. Furthermore, this remains true even for dipolar Fermi gases where, although the s-wave scattering is absent due to the Pauli exclusion principle, the anisotropic DDI competes with the large kinetic energy close to the FS, yielding a complex-enough energy landscape. This leads to novel many-body phenomena, including the recently observed deformation of the FS [20,21] and a predicted fermionic pairing of one-and two-components systems [22][23][24][25][26][27].
Since the first experimental realization of a quantum degenerate dipolar Fermi gas of 161 Dy in 2012 [28], several more fermionic species, such as 167 Er [29] and chromium 53 Cr [30], 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 [31,32] were based on photoassociation of the stimulated Raman adiabatic passage (STIRAP) [33]. As a result, samples of fermionic 40 [40,41], 87 Rb 133 Cs [42,43] and 23 Na 87 Rb [44] 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 40 K 87 Rb has been realized at JILA [45]. The new 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 [45], 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 normal and quantum degenerate dipolar fermions was previously studied in cylindrically symmetric harmonic traps with the dipoles along one of the trap axes [46][47][48][49][50][51] and in homogeneous systems [52]. 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 using a variational phase-space approach based on the Wigner function [21,53,54]. In particular, we investigate the stability of a polarized 40 K 87 Rb gas in an experimentally realistic parameter regime and calculate critical values of the electric dipole moment and the corresponding deformation of the FS. Furthermore, we obtain stability diagrams for an arbitrary orientation of the dipoles and for an angular dependence of the FS deformation. Finally, we demonstrate effects of the strong DDI on the time-of-flight (TOF) expansion and show that a nonballistic expansion theory is essential to understand the behavior of such strong dipolar Fermi gases and to interpret the experimental results. We consider the system to be at zero temperature, which is justified if we consider that T /T F ≈ 0.3 in the experiments [45], and that the thermal corrections to the total energy are proportional to (T /T F ) 2 [55].
The variational approach [21,53,54] for the Wigner function ν(r, k) = d 3 r e −ik·r ρ r + 1 2 r , r − 1 2 r , where ρ(r, r ) = Ψ (r)Ψ † (r ) represents the one-body density matrix, relies on the Hartree-Fock mean-field approximation. Note that the second-order terms in the DDI in the theory beyond Hartree-Fock [56][57][58] yield a negligible correction even for polar molecules. Assuming a Heaviside-shaped Wigner function in the ground state, we 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 , where the angles (θ, ϕ) determine the dipoles' orientation, (θ , ϕ ) describe the orientation of the molecular cloud, and (θ , ϕ ) determine the FS orientation, as illustrated in Figs. 3(a)-(c), respectively. We stress that the molecular cloud orientation does not have to coincide with the trap orientation, due to the DDI effects, and R stands for the corresponding rotation matrix, while ω i denote the trap frequencies and F A is a generalized anisotropy function (see Supplemental Material [59] for further details). By extremizing the above energy with respect to the variational parameters (R i , K i , θ , ϕ , θ , ϕ ), we obtain the corresponding equations for the ground state, which can be rewritten in a dimensionless, species-independent form [59] such that, for a given orientation of the dipoles, they only depend on three parameters: the two trap aspect ratios ω z /ω x and ω z /ω y , and the relative DDI strength 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. 1(a) for the case θ = ϕ = 0, in terms of the maximal DDI strength ε crit dd 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. We also note that ε crit dd turns out to be a symmetric function of its arguments ω z /ω x and ω z /ω y [59]. If we consider the experimentally available species 40 K 87 Rb, the stability diagram from Fig. 1(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. 1(b) and 1(c). Here we see how the critical value of the dipole moment d crit depends on ω x and ω y for a fixed value of ω z . If we take into account that the permanent electric dipole moment of 40 K 87 Rb has the value d = 0.574 D, denoted by black lines in Figs. 1(b) and 1(c), we read off that for ω z = 2π × 50 Hz the instability can kick in already for frequencies ω x , ω y of that order or larger. In the experiment of Ref. [45] the frequencies used are (ω x , ω y , ω z ) = 2π × (63, 36, 200) Hz, and Fig. 1(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 ∆ = 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 167 Er atoms [20], where ∆ is of the order of a few percent. This effect is much larger in gases of polar molecules, as can be seen for 40 K 87 Rb in Fig. 2(a), obtained by solving the equations presented in the Supplemental Material [59]. 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 crit = 0.24 D obtained in Fig. 1(c). For realistic values of the trapping frequencies we obtain that ∆ 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 rela- tively small changes in the dipolar moment strength can significantly affect the system's stability. This is demonstrated in Fig. 2(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. 2(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 [21, 46-51, 53, 54] 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 (θ , ϕ ) of the TF ellipsoid in real space as free variational parameters, together with the TF radii R i , as illustrated in Fig. 3(b). It was previously experimentally verified and theoretically always assumed [20,21] 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 (θ , ϕ ) are also taken as free variational parameters, together with the TF momenta K i , Fig. 3(c). However, we show here [59] that the principle of minimizing the energy leads to the solution θ = θ, ϕ = ϕ [59], 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 , θ , ϕ ), as well as the number of equations [59].
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 be only determined numerically. Figure 3(d) shows the angular stability diagram for a 40 K 87 Rb system in terms of the critical dipole moment d crit , where all variational parameters (R i , K i , θ , ϕ ) are numerically calculated for each configuration. If one would assume that the molecular cloud follows the trap shape, i.e., θ = ϕ = 0, the obtained values of d crit would be significantly different from those calculated in Fig. 3(d)   In both plots solid lines correspond to a nonballistic expansion, where the DDI is taken into account, while the dashed lines represent calculated results for a free (ballistic) expansion, N = 3 · 10 4 , θ = ϕ = 0. AR is calculated using the imaging angle 22.5 • of Ref. [45], in the geometry of Ref. [21].
systems with moderate to strong DDI. Figures 3(e) and 3(f) 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. In contrast to atomic magnetic species, where the angular dependence of the FS deformation is of the order of few per mill [21] and, thus, quite weak, the strong DDI in the samples of polar molecules leads to a much stronger angular dependence, as shown in Figs. 3(e) and 3(f). Namely, for stronger DDI we expect not only an increased critical temperature of Cooper pairing, but also a higher degree of tunability as the deformation of the FS depends on the dipoles' orientation relative to the trap geometry. Effects of the DDI and their interplay with the geometry quite strongly influence also the dynamics of the system. This is of particular importance for interpreting TOF imaging data, which are commonly used for experimental measurements of the properties of ultracold Fermi gases. Even for magnetic atomic species such as erbium this could be experimentally observed and a nonballistic expansion has to be used in order to describe the system's dynamics [54,60]. For polar molecules with a strong DDI we expect that nonballistic effects are more pronounced, as can be read off from Fig. 4

(see Supplemental Material
[59] for further details). 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 of Fig. 4(a). Although the trap geometry plays a role here, Fig. 4(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. 4(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. The general theory presented here 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 [22][23][24][25], 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.

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 d. The dipole moments are aligned in the direction defined by spherical angles (θ, ϕ), as shown in Fig. 3(a). The system is confined into a three-dimensional harmonic trap with the frequencies ω i , whose axes coincide with the axes of the laboratory coordinate system S, as it is depicted in Fig. S1. For the ideal Fermi gas, the molecular cloud shape is determined purely by the trap, Fig. S1(a), while the FS is a sphere, Fig. S1(b). 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 as depicted in Fig. S1(c). 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. S1(d), the FS ellipsoid is oriented in the direction defined by the angles (θ , ϕ ), which are free parameters.
FIG. S1. Schematic illustration of the ansatz for the shape of the molecular cloud and the FS.
We use the following variational ansatz for the Wigner distribution function: where Θ represents the Heaviside step function, while A ij and 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. S1. With this ansatz we obtain the total energy of the system in the Hartree-Fock approximation: where R i and K i are the TF radii and momenta, respectively, c 0 = 2 10 d 2 /(ε 0 · 3 4 · 5 · 7 · π 3 ) is a constant related to the DDI strength, and R ij are matrix elements of the rotation matrix R = R(θ , ϕ ), with R(α, β) =   cos α cos β − sin β sin α cos β cos α sin ϕ cos β sin α sin β − sin α 0 cos α   .
The features of the DDI are embodied into the generalized anisotropy function F A (x, y, θ, ϕ,θ,φ), which includes dependence on the general orientation of the dipoles d and the corresponding TF ellipsoid: where R ij andR ij are matrix elements of the rotation matrix R = R(θ, ϕ) andR = R(θ,φ), respectively. In the above definition, f (x, y) stands for the well-known anisotropy function [61] f (x, y) = 1 − 3xy where φ = arccos x and κ 2 = (1 − y 2 )/(1 − x 2 ), while F (φ, k) and E(φ, 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 Due to these identities and the fact that f (x, y) is a symmetric function, the obtained distributions of ε crit dd and d crit in Fig. 1, as well as the distributions of ∆ in Fig. 2 are symmetric with respect to their arguments for θ = ϕ = 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 now determined by the 10 variational parameters (R i , K i , θ , ϕ , θ , ϕ ), 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: 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.

Orientation of the Fermi surface
Equations (S15) and (S16) can be solved analytically, independently of other equations, yielding the physically expected result θ = θ, ϕ = ϕ. This means that the FS stretches along the dipoles' orientation, as it was verified both experimentally and theoretically for the atomic erbium gas [21]. 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. 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.