Fourth- and fifth-order virial expansion of harmonically trapped fermions at unitarity

By generalizing our automated algebra approach from homogeneous space to harmonically trapped systems, we have calculated the fourth- and fifth-order virial coefficients of universal spin-1/2 fermions in the unitary limit, confined in an isotropic harmonic potential. We present results for said coefficients as a function of trapping frequency (or, equivalently, temperature), which compare favorably with previous Monte Carlo calculations (available only at fourth order) as well as with our previous estimates in the untrapped limit (high temperature, low frequency). We use our estimates of the virial expansion, together with resummation techniques, to calculate the compressibility and spin susceptibility.

Introduction.-At low enough temperatures, or high enough densities, matter invariably displays its quantum mechanical nature, first and foremost by virtue of quantum statistics (i.e. particles are ultimately bosonic or fermionic, at least in three spatial dimensions) but also due to interaction effects that may alter the nature of the equilibrium state. As the temperature is raised, these systems eventually undergo a quantum-classical crossover (QCC) in which interactions still play a role, but where quantum mechanical effects are slowly washed out by temperature fluctuations. This regime is especially interesting for strongly coupled matter, in particular in cases where a superfluid phase is present, as the behavior above the superfluid critical temperature (i.e. in the un-ordered phase) is still significantly affected by the interactions (e.g. inducing pairing correlations) but there are no obvious effective theory descriptions [1][2][3][4][5][6][7][8][9][10].
The QCC is governed by the so-called virial expansion (VE) [11], which breaks the quantum many-body problem into n-particle subspaces, captured in the so-called virial coefficients (see Ref. [12] for a review). For bulk thermodynamic quantities the virial coefficients are denoted by b n , and their change due to interactions is ∆b n . The calculation of ∆b n has a sparse history that started with ∆b 2 in 1937, by Beth and Uhlenbeck [13] and remained largely quiet until the early 21st century. On the theory side, this quiet period can be attributed to the well-known fact that the quantum two-body problem is considerably easier to solve (and to relate to two-body scattering properties) than its three-and higher-body counterparts. On the experimental side, these quantum virial coefficients became increasingly relevant in the early 2000's with the rise of ultracold atom experiments around the world and their ever-increasing ability to create, manipulate, and measure atomic clouds [14].
One of the most famous systems studied with ultracold atoms is the so-called unitary limit of the spin-1/2 Fermi gas [15], which represents a universal regime relevant for atomic and nuclear physics [16][17][18][19][20][21]. In this work we investigate the QCC of this universal regime using the VE up to fifth order for a system confined by a harmonic oscillator (HO) potential. Previous numerical work calcu-lated ∆b 3 [22][23][24][25][26][27] and ∆b 4 [28][29][30][31] (see also Ref. [32][33][34]). More recent work [35] studied analytic expressions in the so-called semiclassical approximation (previously implemented in a wide variety of situations [36][37][38][39][40][41]), which uses a coarse discretization of imaginary time. On the experimental side, there have also been attempts to determine ∆b 4 at unitarity in the untrapped limit, using measurements of the equation of state [42,43], However, those analyses are numerically challenging because one must fit a fourth-order polynomial assuming higher-order contributions are small (which is not necessarily the case, as shown in Ref. [40]).
In this work we generalize the above calculations to include ∆b 5 and go far beyond the semiclassical approximation, extrapolating to the continuous imaginary time limit. While we restrict ourselves to the unitary limit for the most part, we provide approximate analytic formulas that apply to arbitrary interaction strengths, trap frequency, and spatial dimension.
Hamiltonian and virial expansion.-In this work we focus on a system of harmonically trapped spin-1/2 fermions interacting via a short-range interaction. Thus, the Hamiltonian isĤ =T +V ext +V , wherê is the kinetic energy operator, is the external potential energy operator, and is the interaction. Above, m is the mass of the particles, ω is the isotropic harmonic trapping frequency, g is the bare coupling,n s (r) =ψ † s (r)ψ s (r) is the particle density operator for spin-s particles, andψ † s (r) andψ s (r) are, respectively, the creation and annihilation operators for particles of spin s at position r. We use units such that = k B = m = 1 from this point on. Naturally, the noninteracting pieceT +V ext can be diagonalized exactly in the single-particle subspace of the Fock space, which leads to the HO basis we will refer to below. The contact interaction of Eq. (3) is singular in three spatial dimensions and must therefore be regularized and renormalized. To that end, we place the system on a spatial lattice of spacing and implicitly take the continuum limit by transforming spatial sums into integrals at the end. In the process, we renormalize by tuning the coupling so that the known two-body answer for the second-order virial coefficient is reproduced (see below).
The VE accesses thermodynamics by breaking down the calculation by particle number. Specifically, one expands the grand thermodynamic potential Ω in powers of the fugacity z = exp(βµ) as where β is the inverse temperature, Q 1 is the singleparticle partition function, and b n is the n-th order virial coefficient. The b n capture, in a nonperturbative fashion, the contribution of the n-body problem to the full Ω. Plugging in the definition of the grand-canonical partition function Z, namely into Eq. (4) and expanding ln Z in powers of z, the b n can be written in terms of the N -particle canonical partition functions Q N = tr N e −βĤ , where the trace is over the N -particle Hilbert space (see below). Computational framework.-To evaluate Q N , we implement a symmetric Suzuki-Trotter decomposition where we split β = τ N τ into N τ time steps. Thus, where the cyclic property of the trace was used. To proceed, we calculate the matrix elements of the factors inside the trace in coordinate space. For a single imaginary time step, those matrix elements define the factorized transfer matrix M ab , for a particles of spin-↑ and b particles of spin-↓. For example, in the 1 + 1 subspace, i.e. a = b = 1, we obtain where where 1 1 is a 3 × 3 unit matrix. While the above example does not involve identical particles, for the cases that do (e.g. the 2 + 1 subspace of the 3 particle Hilbert space), the (anti-)symmetrization can be carried out at the very end, i.e. after taking the N τ -th power of the distinguishable-particle transfer matrix. This property was already noted by Huang and Yang in 1959 [44,45] and is a consequence of the fact that the operators involved do not change the particles' statistics. Thus, there is no need to use (anti-)symmetrized intermediate states in the calculation, which greatly reduces the computational effort. In the Supplemental Materials we report on the generalization of the above result to M 21 , M 31 , M 22 , M 41 , and M 32 .
Armed with the above factorized transfer matrices M ab , we use automated algebra to symbolically expand [M ab ] Nτ for varying N τ . We then combine the results to obtain the relevant Q N and from them the b n , which are in turn extrapolated to the large-N τ limit. More explicitly, the interaction-induced change ∆b n , for n = 2, 3, 4, 5 is calculated as ∆b 2 = ∆b 11 , ∆b 3 = 2∆b 21 , ∆b 4 = 2∆b 31 + ∆b 22 , and ∆b 5 = 2∆b 41 + 2∆b 32 , where the subspace contributions are and the 4-and 5-particle subspaces are shown in the Supplemental Materials. Here, ∆X represents the change in X induced by the interactions and the Q ab are the canonical partition functions for a particles of spin-↑ and b particles of spin-↓. In the above expressions, the ∆b ab are intensive quantities, whereas the Q ab themselves scale as V a+b where V is the spatial volume. That property emphasizes the challenge in calculating ∆b ab numerically: the delicate cancellations must be resolved among the various terms involving different Q ab 's. It is for that reason that automated algebra methods are advocated here, where those cancellations can be resolved using arbitrary precision arithmetic, avoiding stochastic effects. Results: Approximate analytic expressions for ∆b n .-For N τ = 1, 2, we carry out calculations entirely analytically in which the coupling strength, the trapping frequency, and the spatial dimension appear as arbitrary variables (in principle, it is also possible to take this to even higher order, but the formulas become extremely long). The resulting formulas for ∆b 3 and ∆b 4 for N τ = 1, first shown in Ref. [35], qualitatively (and in some parameter regions quantitatively) capture the behavior of ∆b n . These formulas are also useful as checks for codes that implement higher values of N τ . Here, we provide results broken down by subspace for up to five particles, shown in full detail in the Supplemental Materials.
From those formulas we learn that (as shown in Fig. 1), increasing N τ does not immediately improve the quality of the final answer; rather, the results could move away from the N τ → ∞ limit before the asymptotic regime is reached, usually for N τ > 2. Simply put, as N τ is increased the results may worsen before they improve. Thus, it is important to investigate as large N τ as possible, even if low values are qualitatively correct. In our automated calculations, we explored up to N τ = 20 (for ∆b 21 ), 16 (∆b 31 ), 12 (∆b 22 ), 12 (∆b 41 ), and 8 (∆b 32 ), which we used to estimate the full ∆b 3 , ∆b 4 , and ∆b 5 , extrapolated to N τ → ∞.
Results: Virial coefficients in the unitary limit.-In our approach, we calculate ∆b 2 as a function of the bare coupling C, βω, and N τ , and renormalize by tuning C to the known result in the unitary limit [12], namely ∆b 2 = [4 cosh(βω/2)] −1 . Thus, the second-order VE is reproduced exactly by virtue of this renormalization condition, such that the line of constant physics is followed in the extrapolation to N τ → ∞, for each βω (see Supplemental Materials of Ref. [40]).
In Fig. 1 we show our results for ∆b 3 (top), ∆b 4 (center), and ∆b 5 (bottom) for a unitary Fermi gas in a harmonic trap as a function of βω. The error bars represent the uncertainty in the N τ → ∞ extrapolation, given by the difference between the maximum and minimum predictions of polynomial extrapolation schemes (degrees 2 to 5 for ∆b 3 and ∆b 4 , and degrees 2 and 3 for ∆b 5 , where the data is more limited; see Ref. [40]). Our results for ∆b 3 are in superb agreement with the quantum Monte Carlo data of Ref. [28] as well as with the homogeneous-limit answer of Ref. [40]. [The homogeneous limit is related to the results shown here by ∆b h n = n 3/2 ∆b n (βω → 0), see Refs. [12,22,23,46]]. The case of ∆b 4 is less clear cut: there is good agreement with Ref. [28] for βω ≥ 1, but a clear difference remains at low frequencies. We return to this issue below. Finally, we predict ∆b 5 as a function of βω, which to the best of our knowledge does not appear elsewhere in the literature. As the HO potential confines the system, it naturally increases its kinetic energy, effectively reducing the interaction effects. This suggests that, for a given interaction strength, the VE should enjoy better convergence properties when a trapping potential is turned on (as argued also in Ref. [12]). Indeed, although our results indicate that ∆b 4 ∆b 5 and, moreover, for 0.3 < βω < 1.4 we find ∆b 5 > |∆b 4 |, we also find that |∆b 2 | |∆b 3 | |∆b 4 |. To better understand the differences in ∆b 4 between our results and Ref. [28], we plot in Fig. 2 (top panel) the subspace contributions ∆b 31 and ∆b 22 . As pointed out in Ref. [28], these contributions partially cancel each other out, leading to the observed increased uncertainty in the final answer. Clearly, the largest differences arise [28] appears as red squares for ∆b3 and red circles for ∆b4, in both cases with error bars. The dashed-dotted line in the middle plot shows a high-temperature fit to the data of Ref. [28]. Black stars with error bars show the results by Hou and Drut from Ref. [40] calculated in the homogeneous gas limit. The dotted (dashed) line shows the Nτ = 1 (Nτ = 2) results given analytically in the Supplemental Materials. The latter show that, for ∆b3, increasing Nτ from 1 to 2 shows a dramatic improvement, whereas the case of ∆b4 is a cautionary tale: as Nτ goes from 1 to 2, the results move away from our final answer (blue crosses). In fact, it is not until Nτ = 5 that ∆b4 reaches the asymptotic regime one can use for extrapolation. Reference [29] presented a large-βω asymptotic formula for ∆bn, but its validity is well outside the 0 < βω < 3 region studied here.
in the determination of ∆b 22 , which is not unexpected as a contact interaction in that subspace is less susceptible to Pauli blocking than ∆b 31 . Figure 2 (bottom panel) shows our results for ∆b 41 and ∆b 32 , whose behavior parallels ∆b 31 and ∆b 22 in that they enter with different signs but similar magnitude, thus leading to increased uncertainty in the final result for ∆b 5   The black star and black pentagon show, respectively, the results for ∆b31 and −∆b22 in the βω → 0 limit, obtained in Ref. [40]. Bottom: −∆b41 (blue diamonds) and ∆b32 (green squares) as functions of βω.
The black star and black pentagon show, respectively, the results for −∆b41 and ∆b32 at βω = 0 from Ref. [40].
are able to resolve the fifth-order contribution, as shown already in the bottom panel of Fig. 1. Nevertheless, the size of the error bars of ∆b 32 is larger than that of ∆b 41 . This may come as a surprise given the results of Ref. [40], whose uncertainty at βω = 0 for ∆b 41 is larger than for ∆b 32 . Those results were calculated at the same N τ = 9 order for both coefficients, using an analytic cancellation of volume-dependent terms. In contrast, in the present work we achieved N τ = 12 for ∆b 41 but only N τ = 8 for ∆b 32 , due to the increasing computational cost of cancelling the volume-dependent terms, which is done numerically in the trapped case.
Results: Applications to thermodynamics.-Having obtained the precise form of ∆b 3 , ∆b 4 , and ∆b 5 as functions of βω for harmonically trapped fermions in the unitary limit, we apply those results to obtain thermodynamic information. As an example, we report here the compressibility and magnetic susceptibility, respectively χ n and χ s , defined as χ n,s = β −1 ∂ 2 ln Z/∂h 2 ± , where h ± = (µ ↑ ± µ ↓ )/2 and µ s is the chemical potential for spin-s particles. The interaction effects on χ s,n are where z s = e βµs is the fugacity for spin-s particles. Our results, shown in Fig. 3, indicate that the partial sums of the VE display large variations for ∆χ n,s as the VE order is increased, in particular for z ≥ 1. However, we also see that, using the high-order coefficients we calculated here, it is possible to carry out a Padé resummation [and related strategies (see e.g. [47])] to obtain sensible results for static response functions even as far as z = 3.
Conclusion and outlook.-In this work we have determined the frequency dependence of the virial coefficients b n of HO-trapped spin-1/2 fermions at unitarity. We used a discretization of the imaginary time direction and a Suzuki-Trotter factorization of the transfer matrix, together with automated algebra methods, to calculate canonical partition functions and from them the interaction induced change ∆b n , for n = 3, 4, 5, which we extrapolated to the continuous-time limit. To complement those numerical results, we provided analytic formulas for ∆b n in coarse lattices for arbitrary trap frequency and spatial dimension. Using our final ∆b n , we calculate the compressibility and susceptibility of the unitary Fermi gas and showed that the VE can be Padé-resummed to obtain sensible results even as far as z = 3.