Ground state of an ultracold Fermi gas of tilted dipoles in elongated traps

The physical properties of the above described system can be captured by means of the semiclassical Wigner function [50]. Indeed, the quantum-mechanical expectation values of the system observables, which are required for the calculation of the properties of nonrelativistic quantum systems based on exact diagonalization, can be obtained as their phase-space averages, weighted by the Wigner function. In the case where the dipolar orientation axis lies along one of the trapping axes, an accurate ansatz for the Wigner function takes the simple form [38–47]

$$\begin{eqnarray}&&{\rm{\Theta }}\left(1-\displaystyle \sum _{i}\displaystyle \frac{{r}_{i}^{2}}{{R}_{i}^{2}}-\displaystyle \sum _{i}\displaystyle \frac{{k}_{i}^{2}}{{K}_{i}^{2}}\right),\end{eqnarray} \tag{ 4 }$$

where Θ represents the Heaviside step function, while the variational parameters R_i and K_i stand for the Thomas–Fermi radius and the Fermi momentum in the direction $i=x,y,z$, respectively. This ansatz is motivated by the Fermi–Dirac distribution of a noninteracting Fermi gas at zero temperature, whose Wigner function has the above form. A theory based on the above ansatz [38, 47] was successfully used to determine the deformation of the FS, while its extension [46] enabled a detailed analysis of the ground state and the time-of-flight (TOF) expansion dynamics of the system for different collisional regimes. Furthermore, numerical comparisons [59, 60] have confirmed that, even in the case of polar molecules with masses of the order of 100 atomic units and an electric dipole moment as large as 1D, the above variational ansatz yields highly accurate results, within a fraction of per mille. This indicates that the general ansatz (4), first introduced in a slightly different manner in [38], is indeed very well suited to describe dipolar Fermi gases.

However, the experiment of [37] was performed for an arbitrary angle θ, and therefore the comparison of the theory [46, 47] was only possible for the special case of dipoles oriented along the z axis, i.e., for θ = 0°. Therefore, in order to model the global equilibrium distribution of the dipolar Fermi gas for arbitrarily oriented dipoles and to provide an accurate description of the experiment, it is necessary to generalise the ansatz (4). This is done in the present paper, where we apply an analogous reasoning and introduce the following ansatz for the Wigner distribution

$$\begin{eqnarray}&&{\nu }^{}({\bf{r}},{\bf{k}})={\rm{\Theta }}\left(1-\displaystyle \sum _{i,j}{r}_{i}{{\mathbb{A}}}_{{ij}}{r}_{j}-\displaystyle \sum _{i,j}{k}_{i}{{\mathbb{B}}}_{{ij}}{k}_{j}\right),\end{eqnarray} \tag{ 5 }$$

where ${{\mathbb{A}}}_{{ij}}$ and ${{\mathbb{B}}}_{{ij}}$ are matrix elements that account for the generalised geometry of the system and determine the shape of the cloud in real space and of the FS in momentum space, respectively.

The particle density distribution in real space is determined by the trapping potential and the Hartree direct energy. Therefore, one expects that the matrix ${\mathbb{A}}$ can be well approximated by a diagonal matrix in the coordinate system S, which is defined by the harmonic trap axes

$$\begin{eqnarray}{\mathbb{A}}=\left(\begin{array}{ccc}1/{R}_{x}^{2} & 0 & 0\\ 0 & 1/{R}_{y}^{2} & 0\\ 0 & 0 & 1/{R}_{z}^{2}\end{array}\right).\end{eqnarray} \tag{ 6 }$$

In this way, we neglect off-diagonal elements, which may arise due to the dipoles' arbitrary orientation. However, this is certainly justified for elongated traps, for which the cloud shape in the ground state is well determined by the trap. Therefore, we will consider here trap configurations that satisfy this condition.

On the other side, the momentum distribution is dominated by the interplay between the Pauli pressure, which is isotropic, and the Fock exchange energy, which is responsible for the deformation of the FS [38]. The experiment of [37] suggested that the FS follows the rotation of the external field, keeping the major axis of the FS always parallel to the direction of the maximum attraction of the DDI. Motivated by this, we will consider several possible scenarios for a detailed theoretical description, in order to verify the above hypothesis.

For the sake of completeness, we start with a simple spherical scenario, in which the FS remains a sphere, as displayed in figure 2(a). In that case all Fermi momenta are equal (K_i = K_F) and the matrix ${\mathbb{B}}$ is given by ${{\mathbb{B}}}_{1}={\mathbb{I}}/{K}_{{\rm{F}}}^{2}$, where ${\mathbb{I}}$ is the 3 × 3 identity matrix. We also consider a second, on-axis scenario, which includes the FS deformation such that it is an ellipsoid with fixed major axes coinciding with the trap axes, as shown in figure 2(b). Here the matrix ${\mathbb{B}}$ has a diagonal form

$$\begin{eqnarray}{{\mathbb{B}}}_{2}=\left(\begin{array}{ccc}1/{K}_{x}^{2} & 0 & 0\\ 0 & 1/{K}_{y}^{2} & 0\\ 0 & 0 & 1/{K}_{z}^{2}\end{array}\right),\end{eqnarray} \tag{ 7 }$$

in a similar way as the matrix ${\mathbb{A}}$, which recovers the old ansatz given by (4). We note that the first, spherically-symmetric scenario is a special case of the second ansatz, obtained by restricting the Fermi momenta to be equal. Finally, as a third and more general possibility, we consider the off-axis hypothesis of [37] and assume that the matrix ${\mathbb{B}}$ has a diagonal form ${{\mathbb{B}}}_{3}^{{\prime} }$ in a rotated coordinate system $S^{\prime}$, which is defined by the axes q_x, q_y and q_z, where the q_z axis remains parallel to the dipole moments, as depicted in figure 2(c)

$$\begin{eqnarray}{{\mathbb{B}}}_{3}^{{\prime} }=\left(\begin{array}{ccc}1/{K}_{x}^{{\prime} 2} & 0 & 0\\ 0 & 1/{K}_{y}^{{\prime} 2} & 0\\ 0 & 0 & 1/{K}_{z}^{{\prime} 2}\end{array}\right).\end{eqnarray} \tag{ 8 }$$

Here the parameters ${K}_{i}^{{\prime} }$ represent the Fermi momenta in the rotated coordinate system $S^{\prime}$. In order to describe the rotation from S to $S^{\prime}$, we introduce a rotation matrix ${\mathbb{R}}$,

$$\begin{eqnarray}{\mathbb{R}}=\left(\begin{array}{ccc}\cos \theta \cos \varphi & -\sin \varphi & \sin \theta \cos \varphi \\ \cos \theta \sin \varphi & \cos \varphi & \sin \theta \sin \varphi \\ -\sin \theta & 0 & \cos \theta \end{array}\right),\end{eqnarray} \tag{ 9 }$$

such that ${{\mathbb{B}}}_{3}^{{\prime} }={{\mathbb{R}}}^{T}{{\mathbb{B}}}_{3}{\mathbb{R}}$ and ${\bf{q}}={{\mathbb{R}}}^{T}{\bf{k}}$, where the angles θ and φ are defined in figure 1. We again note that the on-axis scenario is a special case of the off-axis one when the dipoles are parallel to one of the trap axes. For example, for θ = φ = 0° the matrix ${{\mathbb{B}}}_{3}$ is already diagonal, i.e., ${K}_{i}^{{\prime} }={K}_{i}$. We also note that, for all considered ansätze, the normalisation of the Wigner distribution $\nu ({\bf{r}},{\bf{k}})$ is given by

$$\begin{eqnarray}&&N=\int \int \displaystyle \frac{{{\rm{d}}}^{3}{{r}{\rm{d}}}^{3}k}{{(2\pi )}^{3}}\nu ({\bf{r}},{\bf{k}})=\displaystyle \frac{{\bar{R}}^{3}\bar{K}{{\prime} }^{3}}{48},\end{eqnarray} \tag{ 10 }$$

where the bar denotes the geometric averaging: $\bar{R}={({R}_{x}{R}_{y}{R}_{z})}^{1/3}$ and $\bar{K}^{\prime} ={({K}_{x}^{{\prime} }{K}_{y}^{{\prime} }{K}_{z}^{{\prime} })}^{1/3}$.

Zoom In Zoom Out Reset image size

To determine the values of the variational parameters for each scenario, as usual, we require that they minimise the total Hartree–Fock energy of the system. This leads, together with the particle number conservation (10), to algebraically self-consistent equations that determine the Thomas–Fermi radii and momenta. In section 2.2 we calculate the total energy of the system for each of the outlined scenarios and then in section 3.1 we proceed to determine the one that yields a minimal energy and that will be used in the rest of the paper. We note that one can certainly consider even more general ansätze, however, as we will see from the comparison with the experimental data, the proposed approach is fully suitable for describing our system not only qualitatively, but also quantitatively.

Now that we have identified several relevant ansätze for describing the Wigner function of a dipolar Fermi gas of tilted dipoles, we proceed to determine the optimal values of the variational parameters. In order to do so, we have to minimise the total energy of the many-body Fermi system, which is in the Hartree–Fock mean-field theory given by the sum of the kinetic energy E_kin, the trapping energy E_trap, the Hartree direct energy ${E}_{\mathrm{dd}}^{{\rm{D}}}$, and the Fock exchange energy ${E}_{\mathrm{dd}}^{{\rm{E}}}$. Within a semiclassical theory, they can be written in terms of the Wigner function according to [50]

$$\begin{eqnarray}&&{E}_{\mathrm{kin}}=\int \int \displaystyle \frac{{{\rm{d}}}^{3}{{r}{\rm{d}}}^{3}k}{{(2\pi )}^{3}}\displaystyle \frac{{{\hslash }}^{2}{{\bf{k}}}^{2}}{2M}\nu ({\bf{r}},{\bf{k}}),\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&{E}_{\mathrm{trap}}=\int \int \displaystyle \frac{{{\rm{d}}}^{3}{{r}{\rm{d}}}^{3}k}{{(2\pi )}^{3}}{V}_{\mathrm{trap}}({\bf{r}})\nu ({\bf{r}},{\bf{k}}),\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&{E}_{\mathrm{dd}}^{{\rm{D}}}=\displaystyle \frac{1}{2}\int \int \int \int \displaystyle \frac{{{\rm{d}}}^{3}{{r}{\rm{d}}}^{3}r^{\prime} {{\rm{d}}}^{3}{{k}{\rm{d}}}^{3}k^{\prime} }{{(2\pi )}^{6}}{V}_{\mathrm{dd}}({\bf{r}}-{\bf{r}}^{\prime} )\nu ({\bf{r}},{\bf{k}})\nu ({\bf{r}}^{\prime} ,{\bf{k}}^{\prime} ),\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&{E}_{\mathrm{dd}}^{{\rm{E}}}=-\displaystyle \frac{1}{2}\int \int \int \int \displaystyle \frac{{{\rm{d}}}^{3}{{r}{\rm{d}}}^{3}r^{\prime} {{\rm{d}}}^{3}{{k}{\rm{d}}}^{3}k^{\prime} }{{(2\pi )}^{6}}{V}_{\mathrm{dd}}({\bf{r}}^{\prime} ){{\rm{e}}}^{{\rm{i}}({\bf{k}}-{\bf{k}}^{\prime} )\cdot {\bf{r}}^{\prime} }\nu ({\bf{r}},{\bf{k}})\nu ({\bf{r}},{\bf{k}}^{\prime} ),\end{eqnarray} \tag{ 14 }$$

and have already been calculated before with an ansatz (4) for the case when the dipoles are parallel to one of the trap axes [38, 44–47]. Whereas both the kinetic energy (11) and the trapping energy (12) yield simple integrals, the computation of the integrals in the Hartree term (13) and the Fock term (14) is nontrivial and is therefore presented for the most general case in appendices A and B, respectively.

In the spherical scenario, depicted in figure 2(a), the total energy of the system can be calculated using K_i = K_F in ansatz (5), where the Fock exchange term turns out to give no contribution, yielding

$$\begin{eqnarray}&&{E}_{\mathrm{tot}}^{(1)}=\displaystyle \frac{N}{8}\left(\displaystyle \frac{3{{\hslash }}^{2}{K}_{{\rm{F}}}^{2}}{2M}+\displaystyle \sum _{j}\displaystyle \frac{M{\omega }_{j}^{2}{R}_{j}^{2}}{2}\right)-\displaystyle \frac{6{N}^{2}{c}_{0}}{{\bar{R}}^{3}}{f}_{{\rm{A}}}\left(\displaystyle \frac{{R}_{x}}{{R}_{z}},\displaystyle \frac{{R}_{y}}{{R}_{z}},\theta ,\varphi \right).\end{eqnarray} \tag{ 15 }$$

Here ${c}_{0}={2}^{10}{C}_{\mathrm{dd}}/({3}^{4}\cdot 5\cdot 7\cdot {\pi }^{3})$ is a constant related to the dipolar strength, while the features of the DDI are embodied into the generalised anisotropy function ${f}_{{\rm{A}}}(x,y,\theta ,\phi )$, which includes explicitly the angular dependence of the DDI. It is defined as

$$\begin{eqnarray}&&{f}_{{\rm{A}}}(x,y,\theta ,\varphi )={\sin }^{2}\theta {\cos }^{2}\varphi \,f\left(\displaystyle \frac{y}{x},\displaystyle \frac{1}{x}\right)+{\sin }^{2}\theta {\sin }^{2}\varphi \,f\left(\displaystyle \frac{x}{y},\displaystyle \frac{1}{y}\right)+{\cos }^{2}\theta \,f(x,y),\end{eqnarray} \tag{ 16 }$$

where f(x, y) stands for the well-known anisotropy function derived, at first, for dipolar BECs [61]. Note that $f(x,y)={f}_{A}(x,y,0,0)$. This function has been encountered also in previous studies of fermionic dipolar systems [45] in the hydrodynamic collisional regime, as well as in the transition from the collisionless to the hydrodynamic regime in both the TOF expansion dynamics [46] and collective excitations [47]. More details on the anisotropy and the generalised anisotropy function are given in appendix C.

In the on-axis scenario, the FS is deformed to an ellipsoid whose axes are taken to be parallel to the trap axes, as shown in figure 2(b). This ansatz leads to the total energy of the system given by

$$\begin{eqnarray}\begin{array}{rcl}{E}_{\mathrm{tot}}^{(2)} & = & \displaystyle \frac{N}{8}\displaystyle \sum _{j}\left(\displaystyle \frac{{{\hslash }}^{2}{K}_{j}^{2}}{2M}+\displaystyle \frac{M{\omega }_{j}^{2}{R}_{j}^{2}}{2}\right)-\displaystyle \frac{6{N}^{2}{c}_{0}}{{\bar{R}}^{3}}\left[{f}_{{\rm{A}}}\left(\displaystyle \frac{{R}_{x}}{{R}_{z}},\displaystyle \frac{{R}_{y}}{{R}_{z}},\theta ,\varphi \right)-{f}_{{\rm{A}}}\left(\displaystyle \frac{{K}_{z}}{{K}_{x}},\displaystyle \frac{{K}_{z}}{{K}_{y}},\theta ,\varphi \right)\right].\\ & & \end{array}\end{eqnarray} \tag{ 17 }$$

Note that (17) reduces, indeed, to (15) for the special case of ${K}_{x}={K}_{y}={K}_{z}$, since ${f}_{{\rm{A}}}(1,1,\theta ,\varphi )=0$.

Finally, in the most-general considered off-axis scenario displayed in figure 2(c), we allow for both the FS deformation and its rotation so that one of its axes is parallel to the external field. In this case, the total energy of the system reads as

$$\begin{eqnarray}\begin{array}{rcl}{E}_{\mathrm{tot}}^{(3)} & = & \displaystyle \frac{N}{8}\displaystyle \sum _{j}\left(\displaystyle \frac{{{\hslash }}^{2}{K}_{j}^{{\prime} 2}}{2M}+\displaystyle \frac{M{\omega }_{j}^{2}{R}_{j}^{2}}{2}\right)-\displaystyle \frac{6{N}^{2}{c}_{0}}{{\bar{R}}^{3}}\left[{f}_{{\rm{A}}}\left(\displaystyle \frac{{R}_{x}}{{R}_{z}},\displaystyle \frac{{R}_{y}}{{R}_{z}},\theta ,\varphi \right)-f\left(\displaystyle \frac{{K}_{z}^{{\prime} }}{{K}_{x}^{{\prime} }},\displaystyle \frac{{K}_{z}^{{\prime} }}{{K}_{y}^{{\prime} }}\right)\right].\\ & & \end{array}\end{eqnarray} \tag{ 18 }$$

Note that the above form of the Fock energy in the last term is not surprising if we bear in mind the form of equation (16). Namely, in the rotated coordinate system $S^{\prime}$, the axis q_z coincides with the direction of the external field, so the generalised anisotropy function f_A reduces to the standard anisotropy function f, just with the arguments ${K}_{i}^{{\prime} }$ from the rotated system.

In order to compare the three ansätze, we solve the corresponding sets of equations and calculate the total energy of the system in each case. As a model system, we consider the case of a dipolar Fermi gas of atomic ¹⁶⁷Er, using the typical values from the Innsbruck experiments (see below and also [37]), N = 6.6 × 10⁴, (ω_x, ω_y, ω_z) = (579, 91, 611) × 2π Hz, unless otherwise specified. The underlying geometry of the experimental setup is depicted in figure 1.

In figure 3 we compare the total energy of the system as a function of angles θ and φ for all three different scenarios. The comparison is done in terms of the relative total energy shift

$$\begin{eqnarray}&&\delta E=\displaystyle \frac{{E}_{\mathrm{tot}}}{{E}_{0}}-1,\end{eqnarray} \tag{ 19 }$$

where ${E}_{0}=\tfrac{3}{4}{{NE}}_{{\rm{F}}}$ stands for the total energy of the ideal Fermi gas confined into a harmonic trap (3), and ${E}_{{\rm{F}}}={\hslash }\bar{\omega }{(6N)}^{1/3}$ denotes its Fermi energy, where $\bar{\omega }={({\omega }_{x}{\omega }_{y}{\omega }_{z})}^{1/3}$. Figure 3(a) presents the relative total energy shifts as functions of the angle θ for a fixed value of the angle φ = 14°, corresponding to the typical experimental configuration (see below). The three curves, from top to bottom, correspond to $\delta {E}_{}^{(1)}$, $\delta {E}_{}^{(2)}$, and $\delta {E}_{}^{(3)}$, respectively. As a cross-check, we note that the total energies ${E}_{\mathrm{tot}}^{(2)}$ and ${E}_{\mathrm{tot}}^{(3)}$ coincide for θ = 0°. This is expected, since the on-axis scenario is a special case of the off-axis one for θ = 0°. From this figure we immediately see that there are no intersections between the curves, and that it always holds ${E}_{\mathrm{tot}}^{(1)}\geqslant {E}_{\mathrm{tot}}^{(2)}\geqslant {E}_{\mathrm{tot}}^{(3)}$. As a consequence, we conclude that the off-axis scenario, in which the FS is deformed into an ellipsoid that follows the orientation of the dipoles, is favoured among the considered cases as it has the minimal energy. The same conclusion is obtained if we consider the φ-dependence of the relative total energy shifts, depicted in figure 3(b) for a fixed value of the angle θ = 10°. More exhaustive numerical calculations show that this remains to be true even for arbitrary values of the angles θ and φ.

Zoom In Zoom Out Reset image size

Comparing figures 3(a) and (b) we see that the relative total energy shift always remains small, of the order of 0.4%–0.7%, due to a relatively weak DDI between the erbium atoms compared to the energy scale set by the Fermi energy. We also see that the θ-dependence of the total energy is much stronger than the corresponding φ-dependence. We note that the shift would certainly be more significant for atomic and molecular species with a stronger DDI.

The above conclusion is valid not only for the parameters used in figure 3, but, in fact, we have numerically verified that the off-axis scenario for the ansatz (5) for the Wigner function in global equilibrium always yields a minimal energy given by (18) and, thus, we will use it throughout the rest of the paper. The corresponding equations determining all seven variational parameters are given in appendix D as (D.13)–(D.19). A closer examination of those equations reveals that the FS is always a cylindrically symmetric ellipsoid, where ${K}_{x}^{{\prime} }={K}_{y}^{{\prime} }$ holds. This is expected, since the orientation of the dipoles in the rotated coordinate system coincides with the q_z axis and singles this particular direction out, leaving the perpendicular plane perfectly symmetric in momentum space. Therefore, as shown in appendix D, the equations for the seven variational parameters can be rewritten in a more convenient form as (D.16)–(D.22).

Now that we have shown that the FS is, indeed, deformed by the DDI into an ellipsoid, which follows the orientation of the dipoles, we study the angular dependence of this deformation in more detail. To that end, and taking into account the cylindrical symmetry of the FS, we define the FS deformation as the difference between the momentum-space aspect ratio for the dipolar and the noninteracting Fermi gas in the rotated system $S^{\prime}$ according to

$$\begin{eqnarray}&&{\rm{\Delta }}=\displaystyle \frac{{K}_{z}^{{\prime} }}{{K}_{x}^{{\prime} }}-1.\end{eqnarray} \tag{ 20 }$$

This quantity measures the degree of deformation, which emerges purely due to the DDI. In particular, we investigate how the deformation Δ depends on the trap geometry, the orientation of the dipoles, the DDI strength, and the number of particles. The DDI strength is usually expressed in terms of a dimensionless relative strength ε_dd, defined as

$$\begin{eqnarray}&&{\varepsilon }_{\mathrm{dd}}=\displaystyle \frac{{C}_{\mathrm{dd}}}{4\pi }\sqrt{\displaystyle \frac{{M}^{3}\bar{\omega }}{{{\hslash }}^{5}}}{N}^{1/6},\end{eqnarray} \tag{ 21 }$$

which gives a rough estimate of the ratio between the mean dipolar interaction energy and the Fermi energy. We will use it to characterise the strength of the DDI when comparing its effects for different species.

We now calculate the FS deformation of ¹⁶⁷Er for the parameters of [37], yielding the relative interaction strength ε_dd = 0.15. In figure 4(a) we present the angular dependence of Δ on θ and φ, whose values turn out to be around 2.6%, consistent with earlier experimental results [37]. We observe that there is a maximum deformation of the FS at θ = φ = 90°, which corresponds to the direction of the smallest trapping frequency ω_y (y axis). This can be understood heuristically, if one recalls that the DDI is attractive for dipoles oriented head-to-tail. Thus, a weaker trapping frequency favours the stretching of the gas in that direction so that, in turn, this cigar-shaped configuration enhances the relative contribution of the DDI to the total energy.

Zoom In Zoom Out Reset image size

Another aspect relevant for experiments is the influence of the particle number N and the trap geometry on the deformation of the FS. Tuning these parameters and the direction of the dipoles might lead to an enhancement of the DDI effects, and therefore to a stronger deformation of the FS. This is investigated in figures 4(b)–(d), where the FS deformation is given as a function of parameters N, θ and the trap anisotropy $\lambda =\sqrt{{\omega }_{x}{\omega }_{z}}/{\omega }_{y}$ for a fixed value of the angle φ = 90°. Figures 4(c) and (d) explore the FS deformation as a function of the trap anisotropy λ, which was varied by changing the frequencies ω_x = ω_z and ω_y, while keeping the mean frequency $\bar{\omega }=300\times 2\pi$ Hz constant. From all these figures we conclude that the increase in the particle number yields a dominant increase in Δ compared to all other parameters. We note that, in fact, Δ also depends on $\bar{\omega }$, which we do not show here, since it can be directly connected to the particle number dependence. Indeed, the FS deformation depends on ε_dd [37], yielding a dependence of Δ on ${N}^{1/6}{\bar{\omega }}^{1/2}$. As the trap frequencies can be more easily tuned than the particle number, $\bar{\omega }$ can be considered as a predominant control knob in the experiment. However, a precise control of the angles and the anisotropy, which is experimentally easy to realise, may help to achieve an even larger increase in the deformation of the FS. We also note that the λ-dependence is the weakest one, and therefore the formalism for calculating the angular dependence presented here is significant for a systematic study of the influence of the relevant parameters.

Finally, we study the role of the DDI strength and explore whether qualitative changes of the system's behaviour emerge by increasing the value of the dipole moment. To this aim, we compare the erbium case with a molecular Fermi gas of ⁴⁰K⁸⁷Rb, assuming that the same gas characteristics can be achieved in the same trap with this species [54]. The latter possesses an electric dipole moment of strength m = 0.56 D, yielding a much larger relative interaction strength ε_dd = 7.76 for the same parameters. Since the critical value of ε_dd, for which the system is stable, amounts to ${\varepsilon }_{\mathrm{dd}}^{\mathrm{crit}}=2.5$ [46], the molecular ⁴⁰K⁸⁷Rb gas in such a geometry and with the maximal strength of the DDI would in fact not be stable and would collapse under the attractive action of the DDI. For the sake of simplicity and comparison between the systems, we consider a molecular sample of similar geometry and atom number but in which the dipole moment has been tuned to m = 0.25 D by means of an external field [34]. This leads to the relative DDI strength ε_dd = 1.51, which we use in the following.

As we see, the FS deformation Δ has a much stronger angular dependence in figure 5 than in figure 4(a) for the erbium case. Furthermore, in figure 5 we see that Δ has a local minimum for φ = 0° around θ = 40°, while no such minimum exists in figure 4(a). Only a detailed numerical study based on the formalism developed here can provide a precise landscape of the FS deformation behaviour for a concrete experimental setup.

Zoom In Zoom Out Reset image size

Although the shapes of both angular dependencies in figures 4(a) and 5 are quite similar, the main difference is that the deformation of the FS for polar molecules is an order of magnitude larger than for erbium and has a value of around 30%. However, we also observe that the variation in the values of Δ for different angles θ and φ is around 0.03% in the case of an atomic erbium gas, while for the molecules it amounts to around 5%, i.e., the variations of Δ are two orders of magnitude larger for the molecular case. The reason for this increase in both the maximal FS deformation and its angular variation is the same, namely the increase in the relative DDI strength ε_dd, which is one order of magnitude larger for our molecules compared to ¹⁶⁷Er. While the FS deformation is proportional to ε_dd, as expected [46] and as evidenced by our results above, our findings suggest that its maximal angular variation is proportional to ε_dd².

The calculated angular dependence of the FS deformation on the DDI strength has the following important physical consequence. For erbium atoms, where ${\varepsilon }_{\mathrm{dd}}$ is small, the angular variation of the FS deformation is even smaller, since it is proportional to ${\varepsilon }_{\mathrm{dd}}^{2}$, and it would be difficult to observe in experiments. Therefore, one could say that the FS behaves as a rigid ellipsoid, which just rotates following the orientation of the dipoles, without changing its shape [37], as illustrated in figure 6(a). This also implies that the atomic cloud shape in real space is practically disentangled from the FS, and is mainly determined by the trap shape. On the other hand, when ε_dd is large enough, as in the case of ⁴⁰K⁸⁷Rb, the FS not only rotates, but also significantly changes its shape, since the angular variation can be as high as 5%, which is experimentally observable. This is schematically shown in figure 6(b), where the FS behaves as a soft ellipsoid, whose axes are stretched as it rotates. Although we know that the phase-space volume is preserved, according to the particle number conservation (10), figure 6(b) illustrates that the FS, i.e., the momentum-space volume increases (${K}_{i}^{{\prime} }$ increase), while in real space the volume of the cloud shape decreases (R_i decrease). From this we see that the real-space atomic cloud shape is indeed coupled to the FS, and this effect can become measurable in future dipolar fermion experiments, with sufficiently large values of ${\varepsilon }_{\mathrm{dd}}$.

Zoom In Zoom Out Reset image size

The presence of the DDI in both bosonic [62] and fermionic [50] quantum gases has been predicted and evidenced in experiment by detailed TOF expansion measurements [63] to induce magnetostriction in real space, i.e., a stretching of the gas cloud along the direction of the dipoles. In this section we investigate the dependence of this effect on the orientation of the dipoles for the fermionic case. To this end, we first define real-space aspect ratios ${A}_{{ij}}={R}_{i}/{R}_{j}$ of the corresponding Thomas–Fermi radii, as well as their noninteracting counterparts ${A}_{{ij}}^{0}={R}_{i}^{0}/{R}_{j}^{0}={\omega }_{j}/{\omega }_{i}$. The magnetostriction can now be studied in terms of the relative cloud deformations:

$$\begin{eqnarray}&&{\delta }_{{xz}}={A}_{{xz}}/{A}_{{xz}}^{0}-1,\quad {\delta }_{{yz}}={A}_{{yz}}/{A}_{{yz}}^{0}-1.\end{eqnarray} \tag{ 22 }$$

Here the anisotropies due to the harmonic trap are already taken into account and eliminated, so that only effects of the DDI contribute to the nontrivial value of δ_xz and δ_yz. This is in close analogy to the definition of the relative total energy shift of the system in (19), or the FS deformation in (20).

In figure 7 we present the angular dependence of the relative cloud deformations for the Fermi gas of erbium with the same parameters as in figure 4(a). We see that both deformations for a fixed angle $\theta \gt 0$ turn out to possess a minimum for φ = 90°, while along the φ direction the FS deformation monotonously increases. If we compare this to the behaviour in momentum space, see figure 4(b), we see that along the φ direction we have qualitatively the same type of increasing dependence, while the behaviour in the θ direction is markedly different. Indeed, in momentum space it exhibits a maximum for φ = 90°, while in real space we observe a minimum value in that direction for φ = 90°. A related effect has been previously found, showing that the Bose gas momentum becomes distorted in the opposite sense to that of the Fermi gas [40]. There, the effect can be traced back to the differences in the statistical nature of bosons and fermions. Here, however, the different behaviour is due to the anisotropic nature of the DDI and its interplay with the direction of the dipoles and the trap geometry.

Zoom In Zoom Out Reset image size

The TOF images are taken in the plane perpendicular to the imaging axis and the deformation of the atomic cloud is described in terms of the time-dependent cloud aspect ratio ${A}_{{\rm{R}}}$, which is defined as a ratio of vertical and horizontal radii of the cloud in the imaging plane. As depicted in figure 1, the imaging axis $y^{\prime}$ is in the xy plane, rotated by an angle α to the y axis, and the aspect ratio A_R is given by [46]

$$\begin{eqnarray}&&{A}_{{\rm{R}}}(t)=\sqrt{\displaystyle \frac{\langle {r}_{z}^{2}(t)\rangle }{\langle {r}_{x}^{2}(t)\rangle {\cos }^{2}\alpha +\langle {r}_{y}^{2}(t)\rangle {\sin }^{2}\alpha }},\end{eqnarray} \tag{ 23 }$$

where $\langle {r}_{i}^{2}(t)\rangle$ is the average size of the cloud in the direction i = x, y, z after TOF. These quantities are directly measurable in the experiment, and we use them to extract the value of the deformation of the FS, which is connected to the aspect ratio in momentum space. It is defined similarly as A_R [46], according to

$$\begin{eqnarray}&&{A}_{{\rm{K}}}=\sqrt{\displaystyle \frac{\langle {k}_{z}^{2}\rangle }{\langle {k}_{x}^{2}\rangle {\cos }^{2}\alpha +\langle {k}_{y}^{2}\rangle {\sin }^{2}\alpha }},\end{eqnarray} \tag{ 24 }$$

where $\langle {k}_{i}^{2}\rangle$ is the average size of the FS in the direction i = x, y, z in global equilibrium, before the trap is released. Using the definition (24), a straightforward but lengthy calculation yields the following expression for the aspect ratio in momentum space in terms of the Fermi momenta ${K}_{i}^{{\prime} }$ in the rotated coordinate system:

$$\begin{eqnarray}&&{A}_{{\rm{K}}}=\sqrt{\displaystyle \frac{{K}_{x}^{{\prime} 2}{\sin }^{2}\theta +{K}_{z}^{{\prime} 2}{\cos }^{2}\theta }{{K}_{x}^{{\prime} 2}[1-{\sin }^{2}\theta ({\cos }^{2}\varphi {\cos }^{2}\alpha +{\sin }^{2}\varphi {\sin }^{2}\alpha )]+{K}_{z}^{{\prime} 2}{\sin }^{2}\theta ({\cos }^{2}\varphi {\cos }^{2}\alpha +{\sin }^{2}\varphi {\sin }^{2}\alpha )}}.\end{eqnarray} \tag{ 25 }$$

Please note that only for θ = 0°, when the dipoles are parallel to the z axis, the above momentum-space aspect ratio coincides with the ratio between the Fermi momenta, ${A}_{{\rm{K}}}={K}_{z}^{{\prime} }/{K}_{x}^{{\prime} }=1+{\rm{\Delta }}$, where Δ denotes the previously introduced deformation of the FS. In general, however, the relation between A_K and Δ is nonlinear and can be obtained from (25), as follows:

$$\begin{eqnarray}&&{\rm{\Delta }}=\sqrt{\displaystyle \frac{{A}_{{\rm{K}}}^{2}[1-{\sin }^{2}\theta ({\cos }^{2}\varphi {\cos }^{2}\alpha +{\sin }^{2}\varphi {\sin }^{2}\alpha )]-{\sin }^{2}\theta }{{\cos }^{2}\theta -{A}_{{\rm{K}}}^{2}{\sin }^{2}\theta ({\cos }^{2}\varphi {\cos }^{2}\alpha +{\sin }^{2}\varphi {\sin }^{2}\alpha )}}-1.\end{eqnarray} \tag{ 26 }$$

In order to extract the value of the deformation of the FS from the experimental data using the above equation, we still need to calculate the momentum-space aspect ratio A_K. This is done by using the fact that the long-time expansion is mainly dominated by the velocity distribution right after the release from the trap. Here we rely on the ballistic approximation, which assumes that the TOF images, that show the shape of the atomic cloud in real space, purely reflect the momentum distribution in the global equilibrium, i.e.

$$\begin{eqnarray}&&{A}_{{\rm{K}}}\approx \mathop{\mathrm{lim}}\limits_{t\to \infty }{A}_{{\rm{R}}}^{\mathrm{bal}}(t).\end{eqnarray} \tag{ 27 }$$

In principle, this is true just in the case of ballistic expansion, when the effects of the DDI can safely be neglected during the TOF. However, since the DDI is long-range, it should be taken into account, rendering the TOF results always non-ballistic. A general theory that would allow such a treatment is not yet available and is beyond the scope of the present ground-state study. Nevertheless, if the DDI is weak enough, as in the case of erbium atomic gases, the difference between the ballistic (free) and non-ballistic expansion is small, as already shown in [46]. Thus, (27) can approximately be used in our case and the value of A_K in global equilibrium can be extracted from the long-time limit of A_R, which is available from the experimental data. We highlight that in some limiting cases it is still possible to take into account a non-ballistic expansion by using the previously developed dynamical theory [46]. This is expected to yield a more precise value of the aspect ratio, as will be illustrated in the next section.

With those cautionary remarks in mind, we complete the algorithm for analysing our experimental data by calculating the FS deformation from the extracted aspect ratio using (26), which enables its comparison with our numerical results.

Here we consider three different datasets corresponding to the experimental parameters listed in table 1. Case 1 corresponds to the experimental results published in [37], while Case 2 and Case 3 are new results reported here. Table 1 also gives the mean frequency of the trap $\bar{\omega }$ and the trap anisotropy λ for each case. While Cases 1 and 2 represent cigar-shaped traps, Case 2 is selected so that it has the same value of ω_y as Case 1, but a smaller anisotropy λ. On the other hand, Case 3 is chosen so that its mean frequency $\bar{\omega }$ is approximately the same as for Case 1, but its anisotropy λ is much reduced. For each dataset, we probe the FS deformation for various angles θ and a fixed angle φ = 14°. The measurement for each experimental configuration is repeated a large number of times, typically twenty, so that the mean value can be reliably estimated and the statistical error is reduced below 0.2%.

Table 1.  Number of atoms N, trap frequencies ω_i, mean frequencies $\bar{\omega }$ and anisotropies λ for three sets of experimental parameters used throughout this paper. Case 1 corresponds to [37], while Case 2 and Case 3 are new results.

167Er	$N\ (\times {10}^{4})$	ωx (Hz)	ωy (Hz)	ωz (Hz)	$\bar{\omega }\ (\mathrm{Hz})$	λ
Case 1	6.6	579 × 2π	91 × 2π	611 × 2π	318 × 2π	6.54
Case 2	6.3	428 × 2π	91 × 2π	459 × 2π	261 × 2π	4.87
Case 3	6.1	408 × 2π	212 × 2π	349 × 2π	311 × 2π	1.78

Figure 8 shows a direct comparison between our experimental and theoretical results without any free parameters. Experimentally we measure the mean value of the aspect ratio A_R in free expansion using the TOF t = 12 ms, which is taken to be sufficiently long so that the approximation (27) can be used, and yet not too long so that the cloud does not get too dilute and a reliable fit of the density distribution from the absorption images is possible. In figure 8(a) we show the θ dependence of the measured quantity ${A}_{{\rm{R}}}(12\,\mathrm{ms})$ for the parameters of Case 1 (red circles), Case 2 (blue squares) and Case 3 (black triangles), as well as the corresponding theoretical curves (solid red, dashed blue and dotted black line, respectively) for A_K at global equilibrium, calculated according to (25). We see that the agreement is generally very good, and that the experimental data are closely matched by the shape predicted by theory. At the same time, this figure also presents an a posteriori justification for using the ballistic approximation in those three cases.

Zoom In Zoom Out Reset image size

The discrepancies observed in the figure can be accredited to the effects of the DDI, which are neglected during the TOF by using the ballistic approximation. Even better agreement between the experiment and the theory can be expected if a non-ballistic expansion would be taken into account. Although a theory for this is not yet available for an arbitrary orientation of the dipoles, [46] allows us to perform a non-ballistic expansion calculation for the special case θ = 0° in the collisionless regime. The comparison of the results is given in table 2, where we see that accounting for the DDI during the TOF yields theoretical values of the TOF real-space aspect ratio equal to the experimental ones, within the error bars of the order of 0.1%. Table 2 also shows that non-ballistic effects amount to 0.7% for Case 1, which has the largest anisotropy, and becomes smaller as the trap is closer to a spherical shape, i.e., as the trap anisotropy approaches the value of 1. Therefore, we conclude that the agreement of experimental data and theoretical results in figure 8(a) can be further improved by developing a theory for a non-ballistic expansion for a general experiment geometry, which is out of the scope of the present study.

Table 2.  Comparison of theoretical values of aspect ratios in momentum space A_K in global equilibrium and TOF aspect ratios in real space: theoretical value of ${A}_{{\rm{R}}}^{\mathrm{nbal}}$ and experimental value of ${A}_{{\rm{R}}}^{\exp }$, with corresponding statistical errors ${\rm{\Delta }}{A}_{{\rm{R}}}^{\exp }$. Real-space aspect ratios correspond to TOF of t = 12 ms and θ = 0°. Last two columns give trap mean frequency $\bar{\omega }$ and anisotropy λ for each case.

167Er	AK	${A}_{{\rm{R}}}^{\mathrm{nbal}}$	${A}_{{\rm{R}}}^{\exp }$	${\rm{\Delta }}{A}_{{\rm{R}}}^{\exp }$	$\bar{\omega }\ (\mathrm{Hz})$	λ
Case 1	1.0258	1.0324	1.0321	0.0012	318 × 2π	6.54
Case 2	1.0232	1.0282	1.0292	0.0015	261 × 2π	4.87
Case 3	1.0253	1.0270	1.0258	0.0020	311 × 2π	1.78

Figure 8(b) shows a comparison of our theoretical and experimental results for the deformation Δ of the FS for the three considered cases, where the experimental values are calculated according to (26), assuming ballistic expansion (27) and using the real-space aspect ratios shown in figure 8(a). Although the statistical error bars ${\rm{\Delta }}{A}_{{\rm{R}}}^{\exp }$ for the experimentally measured values of the real-space aspect ratios are small and almost constant, the corresponding errors for the FS deformation, calculated as

$$\begin{eqnarray}&&{\rm{\Delta }}{A}_{{\rm{R}}}^{\exp }{\left|\displaystyle \frac{\partial {\rm{\Delta }}}{\partial {A}_{{\rm{K}}}}\right|}_{{A}_{{\rm{K}}}={A}_{{\rm{R}}}^{\exp }},\end{eqnarray} \tag{ 28 }$$

show a strong angular dependence, due to the presence of a pole in the function $\partial {\rm{\Delta }}/\partial {A}_{{\rm{K}}}$. For the parameters of figure 8, the pole emerges at around θ = 50°. Therefore, the error bars appear significantly larger in the neighbouring region, which justifies to drop the data points around θ = 50° (shaded area in the graph).

As can be seen in figure 8(b), for all three cases the deformation of the FS is almost constant for all angles θ. Therefore, from the experimental point of view, it would be enough just to measure the aspect ratio for one value of θ, e.g., θ = 0° in order to determine the deformation of the FS. However, this is only true for a weak enough DDI. Nevertheless, even if this is the case, the measurement of the angular dependence of A_R is an indispensable tool for a full verification of the developed theory, as demonstrated in figures 8(a) and (b).

As already observed in figure 8(a), the A_K curves for all three considered cases intersect at a special point (θ*, ${A}_{{\rm{K}}}^{* }=1$). Figure 9 reveals that this is not just a coincidence. It shows the θ-dependence of the momentum-space aspect ratio A_K for several trapping geometries for erbium atomic gases, ranging from a cigar-shaped trap, through a spherical, to a pancake-shaped trap. The azimuthal angle is kept constant at the value φ = 0°, as well as the trapping frequencies ${\omega }_{x}={\omega }_{z}=500\times 2\pi$ Hz, while the frequency ${\omega }_{y}=n\times 100\times 2\pi$ Hz is varied by changing the value $n\in \{1,2,3,5,7,9\}$, which corresponds to the trap anisotropy λ = 5/n. The number of particles was fixed at N = 7 × 10⁴. We observe again that all curves intersect for ${A}_{{\rm{K}}}^{* }=1$, which suggests that this is a general rule. Indeed, if we take into account that ${K}_{z}^{{\prime} }\geqslant {K}_{x}^{{\prime} }\gt 0$, for ${A}_{{\rm{K}}}^{* }=1$ we can show from (25) that the following relation holds, which connects the intersection angles θ* and φ*:

$$\begin{eqnarray}&&{\sin }^{2}{\theta }^{* }=\displaystyle \frac{1}{1+{\cos }^{2}{\varphi }^{* }{\cos }^{2}\alpha +{\sin }^{2}{\varphi }^{* }{\sin }^{2}\alpha }.\end{eqnarray} \tag{ 29 }$$

This result is universal, i.e., it is independent on other system parameters as the trap geometry, the number of particles, and the DDI strength. In other words, this intersection point is purely a consequence of the geometry, and for any orientation of the dipoles there exists an imaging angle such that the aspect ratio is given by A_K = 1, while the FS deformation Δ can be nontrivial and even can have a significant value. We note that for larger ε_dd values additional parameter-specific intersection points may appear for some geometries, but the intersection point for A_K = 1 is universal and always present.

Zoom In Zoom Out Reset image size

To further illustrate this, in figure 10 we plot a diagram in the $({\theta }^{* },{\varphi }^{* })$-plane for α = 28°, where the regions with A_K > 1 and A_K < 1 are delineated by a solid line defined by (29). The two black dots correspond to intersection points from figure 9 for φ = 0° and from figure 8(a) for φ = 14°, respectively.

Zoom In Zoom Out Reset image size

The five variational parameters (K_F, R_i, μ) are determined by minimising the grand-canonical potential Ω⁽¹⁾, which leads to the following set of algebraic equations:

$$\begin{eqnarray}&&\mu =\displaystyle \frac{{{\hslash }}^{2}{K}_{{\rm{F}}}^{2}}{8M},\end{eqnarray} \tag{ D.1 }$$

$$\begin{eqnarray}&&{\omega }_{x}^{2}{R}_{x}^{2}+\displaystyle \frac{48{{Nc}}_{0}}{M{\bar{R}}^{3}}\left[{f}_{{\rm{A}}}\left(\displaystyle \frac{{R}_{x}}{{R}_{z}},\displaystyle \frac{{R}_{y}}{{R}_{z}},\theta ,\varphi \right)-{R}_{x}{\partial }_{{R}_{x}}\,{f}_{{\rm{A}}}\left(\displaystyle \frac{{R}_{x}}{{R}_{z}},\displaystyle \frac{{R}_{y}}{{R}_{z}},\theta ,\varphi \right)\right]-\displaystyle \frac{8\mu }{M}=0,\end{eqnarray} \tag{ D.2 }$$

$$\begin{eqnarray}&&{\omega }_{y}^{2}{R}_{y}^{2}+\displaystyle \frac{48{{Nc}}_{0}}{M{\bar{R}}^{3}}\left[{f}_{{\rm{A}}}\left(\displaystyle \frac{{R}_{x}}{{R}_{z}},\displaystyle \frac{{R}_{y}}{{R}_{z}},\theta ,\varphi \right)-{R}_{y}{\partial }_{{R}_{y}}\,{f}_{{\rm{A}}}\left(\displaystyle \frac{{R}_{x}}{{R}_{z}},\displaystyle \frac{{R}_{y}}{{R}_{z}},\theta ,\varphi \right)\right]-\displaystyle \frac{8\mu }{M}=0,\end{eqnarray} \tag{ D.3 }$$

$$\begin{eqnarray}&&{\omega }_{z}^{2}{R}_{z}^{2}+\displaystyle \frac{48{{Nc}}_{0}}{M{\bar{R}}^{3}}\left[{f}_{{\rm{A}}}\left(\displaystyle \frac{{R}_{x}}{{R}_{z}},\displaystyle \frac{{R}_{y}}{{R}_{z}},\theta ,\varphi \right)-{R}_{z}{\partial }_{{R}_{z}}\,{f}_{{\rm{A}}}\left(\displaystyle \frac{{R}_{x}}{{R}_{z}},\displaystyle \frac{{R}_{y}}{{R}_{z}},\theta ,\varphi \right)\right]-\displaystyle \frac{8\mu }{M}=0,\end{eqnarray} \tag{ D.4 }$$

$$\begin{eqnarray}&&N=\displaystyle \frac{1}{48}{\bar{R}}^{3}{K}_{{\rm{F}}}^{3}.\end{eqnarray} \tag{ D.5 }$$

Note that (D.5) represents the particle-number conservation constraint, which is the special case of (10) for ${K}_{x}={K}_{y}={K}_{z}={K}_{{\rm{F}}}$.

The seven variational parameters (K_i, R_i, μ) are determined by minimising the grand-canonical potential Ω⁽²⁾, which leads to the following set of algebraic equations:

$$\begin{eqnarray}&&\displaystyle \frac{{{\hslash }}^{2}{K}_{x}^{2}}{2M}+\displaystyle \frac{24{{Nc}}_{0}}{{\bar{R}}^{3}}{K}_{x}{\partial }_{{K}_{x}}\,{f}_{{\rm{A}}}\left(\displaystyle \frac{{K}_{z}}{{K}_{x}},\displaystyle \frac{{K}_{z}}{{K}_{y}},\theta ,\varphi \right)-4\mu =0,\end{eqnarray} \tag{ D.6 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{{\hslash }}^{2}{K}_{y}^{2}}{2M}+\displaystyle \frac{24{{Nc}}_{0}}{{\bar{R}}^{3}}{K}_{y}{\partial }_{{K}_{y}}\,{f}_{{\rm{A}}}\left(\displaystyle \frac{{K}_{z}}{{K}_{x}},\displaystyle \frac{{K}_{z}}{{K}_{y}},\theta ,\varphi \right)-4\mu =0,\end{eqnarray} \tag{ D.7 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{{\hslash }}^{2}{K}_{z}^{2}}{2M}+\displaystyle \frac{24{{Nc}}_{0}}{{\bar{R}}^{3}}{K}_{z}{\partial }_{{K}_{z}}\,{f}_{{\rm{A}}}\left(\displaystyle \frac{{K}_{z}}{{K}_{x}},\displaystyle \frac{{K}_{z}}{{K}_{y}},\theta ,\varphi \right)-4\mu =0,\end{eqnarray} \tag{ D.8 }$$

$$\begin{eqnarray}&&{\omega }_{x}^{2}{R}_{x}^{2}+\displaystyle \frac{48{{Nc}}_{0}}{M{\bar{R}}^{3}}\left[{f}_{{\rm{A}}}\left(\displaystyle \frac{{R}_{x}}{{R}_{z}},\displaystyle \frac{{R}_{y}}{{R}_{z}},\theta ,\varphi \right)-{f}_{{\rm{A}}}\left(\displaystyle \frac{{K}_{z}}{{K}_{x}},\displaystyle \frac{{K}_{z}}{{K}_{y}},\theta ,\varphi \right)-{R}_{x}{\partial }_{{R}_{x}}\,{f}_{{\rm{A}}}\left(\displaystyle \frac{{R}_{x}}{{R}_{z}},\displaystyle \frac{{R}_{y}}{{R}_{z}},\theta ,\varphi \right)\right]-\displaystyle \frac{8\mu }{M}=0,\end{eqnarray} \tag{ D.9 }$$

$$\begin{eqnarray}&&{\omega }_{y}^{2}{R}_{y}^{2}+\displaystyle \frac{48{{Nc}}_{0}}{M{\bar{R}}^{3}}\left[{f}_{{\rm{A}}}\left(\displaystyle \frac{{R}_{x}}{{R}_{z}},\displaystyle \frac{{R}_{y}}{{R}_{z}},\theta ,\varphi \right)-{f}_{{\rm{A}}}\left(\displaystyle \frac{{K}_{z}}{{K}_{x}},\displaystyle \frac{{K}_{z}}{{K}_{y}},\theta ,\varphi \right)-{R}_{y}{\partial }_{{R}_{y}}\,{f}_{{\rm{A}}}\left(\displaystyle \frac{{R}_{x}}{{R}_{z}},\displaystyle \frac{{R}_{y}}{{R}_{z}},\theta ,\varphi \right)\right]-\displaystyle \frac{8\mu }{M}=0,\end{eqnarray} \tag{ D.10 }$$

$$\begin{eqnarray}&&{\omega }_{z}^{2}{R}_{z}^{2}+\displaystyle \frac{48{{Nc}}_{0}}{M{\bar{R}}^{3}}\left[{f}_{{\rm{A}}}\left(\displaystyle \frac{{R}_{x}}{{R}_{z}},\displaystyle \frac{{R}_{y}}{{R}_{z}},\theta ,\varphi \right)-{f}_{{\rm{A}}}\left(\displaystyle \frac{{K}_{z}}{{K}_{x}},\displaystyle \frac{{K}_{z}}{{K}_{y}},\theta ,\varphi \right)-{R}_{z}{\partial }_{{R}_{z}}\,{f}_{{\rm{A}}}\left(\displaystyle \frac{{R}_{x}}{{R}_{z}},\displaystyle \frac{{R}_{y}}{{R}_{z}},\theta ,\varphi \right)\right]-\displaystyle \frac{8\mu }{M}=0,\end{eqnarray} \tag{ D.11 }$$

$$\begin{eqnarray}&&N=\displaystyle \frac{1}{48}{\bar{R}}^{3}{\bar{K}}^{3}.\end{eqnarray} \tag{ D.12 }$$

Similarly as in the spherical scenario, (D.12) coincides with the particle-number conservation equation (10).

The seven variational parameters (${K}_{i}^{{\prime} }$, R_i, μ) are determined by minimising the grand-canonical potential Ω⁽³⁾, which leads to the following set of algebraic equations:

$$\begin{eqnarray}&&\displaystyle \frac{{{\hslash }}^{2}{K}_{x}^{{\prime} 2}}{2M}+\displaystyle \frac{24{{Nc}}_{0}}{{\bar{R}}^{3}}{K}_{x}^{{\prime} }{\partial }_{{K}_{x}^{{\prime} }}f\left(\displaystyle \frac{{K}_{z}^{{\prime} }}{{K}_{x}^{{\prime} }},\displaystyle \frac{{K}_{z}^{{\prime} }}{{K}_{y}^{{\prime} }}\right)-4\mu =0,\end{eqnarray} \tag{ D.13 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{{\hslash }}^{2}{K}_{y}^{{\prime} 2}}{2M}+\displaystyle \frac{24{{Nc}}_{0}}{{\bar{R}}^{3}}{K}_{y}^{{\prime} }{\partial }_{{K}_{y}^{{\prime} }}f\left(\displaystyle \frac{{K}_{z}^{{\prime} }}{{K}_{x}^{{\prime} }},\displaystyle \frac{{K}_{z}^{{\prime} }}{{K}_{y}^{{\prime} }}\right)-4\mu =0,\end{eqnarray} \tag{ D.14 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{{\hslash }}^{2}{K}_{z}^{{\prime} 2}}{2M}+\displaystyle \frac{24{{Nc}}_{0}}{{\bar{R}}^{3}}{K}_{z}^{{\prime} }{\partial }_{{K}_{z}^{{\prime} }}f\left(\displaystyle \frac{{K}_{z}^{{\prime} }}{{K}_{x}^{{\prime} }},\displaystyle \frac{{K}_{z}^{{\prime} }}{{K}_{y}^{{\prime} }}\right)-4\mu =0,\end{eqnarray} \tag{ D.15 }$$

$$\begin{eqnarray}&&{\omega }_{x}^{2}{R}_{x}^{2}+\displaystyle \frac{48{{Nc}}_{0}}{M{\bar{R}}^{3}}\left[{f}_{{\rm{A}}}\left(\displaystyle \frac{{R}_{x}}{{R}_{z}},\displaystyle \frac{{R}_{y}}{{R}_{z}},\theta ,\varphi \right)-f\left(\displaystyle \frac{{K}_{z}^{{\prime} }}{{K}_{x}^{{\prime} }},\displaystyle \frac{{K}_{z}^{{\prime} }}{{K}_{y}^{{\prime} }}\right)-{R}_{x}{\partial }_{{R}_{x}}\,{f}_{{\rm{A}}}\left(\displaystyle \frac{{R}_{x}}{{R}_{z}},\displaystyle \frac{{R}_{y}}{{R}_{z}},\theta ,\varphi \right)\right]-\displaystyle \frac{8\mu }{M}=0,\end{eqnarray} \tag{ D.16 }$$

$$\begin{eqnarray}&&{\omega }_{y}^{2}{R}_{y}^{2}+\displaystyle \frac{48{{Nc}}_{0}}{M{\bar{R}}^{3}}\left[{f}_{{\rm{A}}}\left(\displaystyle \frac{{R}_{x}}{{R}_{z}},\displaystyle \frac{{R}_{y}}{{R}_{z}},\theta ,\varphi \right)-f\left(\displaystyle \frac{{K}_{z}^{{\prime} }}{{K}_{x}^{{\prime} }},\displaystyle \frac{{K}_{z}^{{\prime} }}{{K}_{y}^{{\prime} }}\right)-{R}_{y}{\partial }_{{R}_{y}}\,{f}_{{\rm{A}}}\left(\displaystyle \frac{{R}_{x}}{{R}_{z}},\displaystyle \frac{{R}_{y}}{{R}_{z}},\theta ,\varphi \right)\right]-\displaystyle \frac{8\mu }{M}=0,\end{eqnarray} \tag{ D.17 }$$

$$\begin{eqnarray}&&{\omega }_{z}^{2}{R}_{z}^{2}+\displaystyle \frac{48{{Nc}}_{0}}{M{\bar{R}}^{3}}\left[{f}_{{\rm{A}}}\left(\displaystyle \frac{{R}_{x}}{{R}_{z}},\displaystyle \frac{{R}_{y}}{{R}_{z}},\theta ,\varphi \right)-f\left(\displaystyle \frac{{K}_{z}^{{\prime} }}{{K}_{x}^{{\prime} }},\displaystyle \frac{{K}_{z}^{{\prime} }}{{K}_{y}^{{\prime} }}\right)-{R}_{z}{\partial }_{{R}_{z}}\,{f}_{{\rm{A}}}\left(\displaystyle \frac{{R}_{x}}{{R}_{z}},\displaystyle \frac{{R}_{y}}{{R}_{z}},\theta ,\varphi \right)\right]-\displaystyle \frac{8\mu }{M}=0,\end{eqnarray} \tag{ D.18 }$$

$$\begin{eqnarray}&&N=\displaystyle \frac{1}{48}{\bar{R}}^{3}\bar{K}{{\prime} }^{3}.\end{eqnarray} \tag{ D.19 }$$

As before, (D.19) coincides with the particle-number conservation equation (10).

Due to the symmetry of the anisotropy function f(x, y) = f(y, x), it follows from (D.13) and (D.14) that ${K}_{x}^{{\prime} }={K}_{y}^{{\prime} }$, i.e., that the FS is cylindrically symmetric with respect to the dipoles' orientation. Additionally, in close analogy with the special case when the dipoles are aligned with one of the trapping axes [44–47], the three equations (D.13)–(D.15) can be rewritten in the following form:

$$\begin{eqnarray}&&{K}_{x}^{{\prime} }={K}_{y}^{{\prime} },\end{eqnarray} \tag{ D.20 }$$

$$\begin{eqnarray}&&{K}_{z}^{{\prime} 2}-{K}_{x}^{{\prime} 2}=\displaystyle \frac{144{{MNc}}_{0}}{{{\hslash }}^{2}{\bar{R}}^{3}}\left[1+\displaystyle \frac{(2{K}_{x}^{{\prime} 2}+{K}_{z}^{{\prime} 2}){f}_{{\rm{s}}}\left(\tfrac{{K}_{z}^{{\prime} }}{{K}_{x}^{{\prime} }}\right)}{2({{K}_{z}^{{\prime} }}^{2}-{{K}_{x}^{{\prime} }}^{2})}\right],\end{eqnarray} \tag{ D.21 }$$

$$\begin{eqnarray}&&\mu =\displaystyle \frac{1}{12}\displaystyle \sum _{j}\displaystyle \frac{{{\hslash }}^{2}{K}_{j}^{{\prime} 2}}{2M}.\end{eqnarray} \tag{ D.22 }$$