## Third- and fourth-order virial coefficients of harmonically trapped fermions in a semiclassical approximation

K. J. Morrell, 1 C. E. Berger, 1 and J. E. Drut 1

1 Department of Physics and Astronomy, University of North Carolina, Chapel Hill, NC, 27599, USA (Dated: August 2, 2019)

Using a leading-order semiclassical approximation, we calculate the third- and fourth-order virial coefficients of nonrelativistic spin-1/2 fermions in a harmonic trapping potential in arbitrary spatial dimensions, and as functions of temperature, trapping frequency and coupling strength. Our simple, analytic results for the interaction-induced changes ∆ b 3 and ∆ b 4 agree qualitatively, and in some regimes quantitatively, with previous numerical calculations for the unitary limit of threedimensional Fermi gases.

## I. INTRODUCTION

where

The properties of fermions at finite temperature and density are in part governed by the dimensionless product βµ , where β = 1 ( / k B T ) is the inverse temperature and µ is the chemical potential. Typically, the region βµ glyph[similarequal] 0 displays a crossover between quantum and classical physics, while z = e βµ glyph[lessmuch] 1 indicates a dilute limit where the thermodynamics is given by the virial expansion, which expands a given physical quantity in powers of z . Since µ is coupled to the particle number N , the virial expansion at order N contains the physics of the N -body problem. In the simplest case, the coefficients b n of the virial expansion determine the pressure, density, and compressibility, as well as other elementary thermodynamic quantities such as energy and entropy. The change in b n due to interactions is usually denoted ∆ b n .

The previous work of Ref. [1] calculated the third- and fourth-order virial coefficients ∆ b 3 and ∆ b 4 , respectively, at leading-order (LO) in a semiclassical lattice approximation (SCLA), of homogeneous spin-1 / 2 fermions in arbitrary dimension. The follow-up work of Ref. [2] extended those results up to ∆ b 7 , while Ref. [3] carried out calculations up to next-to-next-to-leading order in the SCLA for up to ∆ b 5 . In this brief work we provide another piece of the puzzle by generalizing the calculations of Ref. [1] to systems in a harmonic trap of frequency ω . We present our derivations with intermediate steps in detail and give analytic formulas for ∆ b 3 and ∆ b 4 as functions of βω in arbitrary spatial dimension d . Our results, which will be given in terms of ∆ b 2 , are thus also functions of the coupling strength.

## II. HAMILTONIAN AND FORMALISM

As our focus is on systems with short-range interactions, such as dilute atomic gases or dilute neutron matter, the Hamiltonian reads

$$\hat { H } = \hat { H } _ { 0 } + \hat { V } _ { \text{int} }, \quad \ \ ( 1 ) \quad \stackrel { \cdot } { \ f u _ { l } }$$

$$\hat { H } _ { 0 } = \hat { T } + \hat { V } _ { \text{ext} },$$

and

$$\hat { T } = \sum _ { s = 1, 2 } \int d ^ { d } x \, \hat { \psi } _ { s } ^ { \dagger } ( { \mathbf x } ) \left ( - \frac { \hbar { ^ { 2 } } \nabla ^ { 2 } } { 2 m } \right ) \hat { \psi } _ { s } ( { \mathbf x } ), \quad ( 3 )$$

is the kinetic energy,

$$\hat { V } _ { \text{ext} } = \frac { 1 } { 2 } m \omega ^ { 2 } \int d ^ { d } x \, \mathbf x ^ { 2 } \left ( \hat { n } _ { 1 } ( \mathbf x ) + \hat { n } _ { 2 } ( \mathbf x ) \right ), \quad ( 4 )$$

is the spherically symmetric external trapping potential, and

$$\hat { V } _ { \text{int} } = - g _ { d D } \int d ^ { d } x \, \hat { n } _ { 1 } ( \mathbf x ) \hat { n } _ { 2 } ( \mathbf x ), \text{ \quad \ \ } ( 5 )$$

is the interaction.

In the above equations, the field operators ˆ ψ , ψ s ˆ † s correspond to particles of species s = 1 2, and ˆ ( , n s x ) are the coordinate-space densities. For the remainder of this work, we will set glyph[planckover2pi1] = k B = m = 1.

## A. Thermodynamics and the virial expansion

The equilibrium thermodynamics of our quantum many-body system can be captured by the grandcanonical partition function, namely

$$\mathcal { Z } = \text{tr} \left [ e ^ { - \beta ( \hat { H } - \mu \hat { N } ) } \right ] = e ^ { - \beta$$

where β is the inverse temperature, Ω is the grand thermodynamic potential, ˆ is the total particle number opN erator, and µ is the overall chemical potential (we will not consider polarized systems in this work).

As the calculation of Z is a formidable problem in the presence of interactions, we resort to approximations or numerical evaluations in order to access the thermodynamics. To that end, in this work we will use the virial expansion, which is an expansion around the dilute limit z → 0, where z = e βµ is the fugacity, i.e. it is a lowfugacity expansion (see Ref. [4] for a review on recent

applications of the virial expansion to ultracold atoms). The coefficients accompanying the powers of z in the expansion Ω are the virial coefficients b n :

$$- \beta \Omega = \ln \mathcal { Z } = Q _ { 1 } \sum _ { n = 1 } ^ { \infty } b _ { n } z ^ { n }, \$$

where Q 1 is the one-body partition function. Using the fact that Z is itself a sum over canonical partition functions Q N of all possible particle numbers N , namely

$$\mathcal { Z } = \sum _ { N = 0 } ^ { \infty } z ^ { N } Q _ { N }, \quad \quad ( 8 ) \quad \quad \quad$$

we obtain expressions for the virial coefficients

$$b _ { 1 } & = 1, & ( 9 ) \quad Q _ { 1 } \, \mathfrak { m }$$

b

Q

2

Q

1

2

=

-

,

(10)

1

Q

2!

$$b _ { 3 } = \frac { \overset { \cdot } { Q } _ { 3 } } { Q _ { 1 } } - b _ { 2 } Q _ { 1 } - \frac {$$

$$b _ { 4 } = \frac { \stackrel { \cdot \cdot } { Q _ { 4 } } } { Q _ { 1 } } - \left ( b _ { 3 } + \frac { b _ { 2 } ^ { 2 } } { 2 } \right ) Q _ { 1 } - b _ { 2 } \frac { Q _ { 1 } ^ { 2 } } { 2! } - \frac { Q _ { 1 } ^ { 3 } } { 4! }, \quad ( 1 2 ) \quad \text{ need}$$

and so on. In this work we will not pursue the virial expansion beyond b 4 . The Q N can themselves be written in terms of the partition functions Q a,b for a particles of type 1 and b particles of type 2:

$$Q _ { 1 } & = 2 Q _ { 1, 0 }, & ( 1 3 ) \quad \text{Th}$$

$$Q _ { 2 } = 2 Q _ { 2, 0 } + Q _ { 1, 1 }, \quad \quad \ \ ( 1 4 ) \quad \text{spi}$$

$$Q _ { 3 } = 2 Q _ { 3, 0 } + 2 Q _ { 2, 1 },$$

$$Q _ { 4 } = 2 Q _ { 4, 0 } + 2 Q _ { 3, 1 } + Q _ { 2, 2 }, \text{ \quad \ \ } ( 1 6 )$$

and so on for higher orders. In the absence of intraspecies interactions, only Q 1 1 , , Q 2 1 , , Q 3 1 , , and Q 2 2 , are affected, such that the change in b 2 , b 3 , and b 4 due to interactions is entirely given by

$$\Delta b _ { 2 } = & \frac { \Delta Q _ { 1, 1 } } { Q _ { 1 } }, & ( 1 7 ) \quad \text{which} \\ \quad \text{$\text{$\text{$\text{$Ea.\,T}$} }$}$$

$$\Delta b _ { 3 } = \frac { 2 \dot { \Delta } Q _ { 2, 1 } } { Q _ { 1 } } - \Delta b _ { 2 } Q _ { 1 }, \quad \ \ ( 1 8 ) \quad \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
 \quad \quad \quad$$

$$\Delta b _ { 4 } = \frac { \stackrel { \mu \mu } { 2 \Delta Q _ { 3, 1 } } + \Delta Q _ { 2, 2 } } { Q _ { 1 } } - \Delta \left ( b _ { 3 } + \frac { b _ { 2 } ^ { 2 } } { 2 } \right ) Q _ { 1 } - \frac { \Delta b _ { 2 } } { 2 } Q _ { 1 } ^ { 2 } ( 1 9 ) \quad \stackrel { \infty \infty } { \text{non} }$$

To calculate ∆ Q m,n , we implement a semiclassical approximation, as described in the next section. Once we obtain the virial coefficients, one may rebuild the grandcanonical potential Ω to access the thermodynamics of the system as a function of the various parameters.

In order to connect to the physical parameters of the systems at hand, we will use the value of ∆ b 2 as a renormalization condition by relying on the exact answers as functions of βω and the physical coupling λ . These exact answers are not always known analytically, but they can easily be obtained numerically by solving the two-body problem of interest.

Although in this work we will focus on systems in a harmonic trap, thus far the identities presented in this section are more general. As a reference for the trapped case, we present here the calculation of the noninteracting virial coefficients for arbitrary βω . [We note that such a calculation, while simple, does not appear in the literature.] Starting from the logarithm of the noninteracting partition function in d spatial dimensions, we have, for two fermion species,

$$\ln \mathcal { Z } = 2 \sum _ { \mathbf n } \ln \left ( 1 + z e ^ { - \beta d / 2 } \prod _ { i = 1 } ^ { d } e ^ { - \beta \omega n _ { i } } \right ). \quad ( 2 0 )$$

Expanding in powers of z on both sides, and switching the order of the sums, we obtain

$$& \text{$one order on use sums, we odom} \\ ( 9 ) & & Q _ { 1 } \sum _ { k = 1 } ^ { \infty } b _ { k } ^ { 0 } z ^ { k } = 2 \sum _ { k = 1 } ^ { \infty } z ^ { k } \frac { ( - 1 ) ^ { k + 1 } } { k } e ^ { - \beta \omega d k / 2 } \left ( \sum _ { n = 0 } ^ { \infty } e ^ { - \beta \omega k n } \right ) ^ { d } \\ & & = 2 \sum _ { k = 1 } ^ { \infty } z ^ { k } \frac { ( - 1 ) ^ { k + 1 } } { 2 ^ { d } k } \left ( \frac { 1 } { \sinh ( \beta \omega k / 2 ) } \right ) ^ { d }. \quad ( 2 1 ) \\ & & \text{To identify the noninteraction $\text{$variial coefficient} \text{ $h^{0}$} \text{ } we}$$

To identify the noninteracting virial coefficients b 0 n , we need Q 1 :

$$Q _ { 1 } = 2 \sum _ { \mathbf n } e ^ { - \beta E _ { n } } = 2 e ^ { - \beta \omega d / 2 } \left ( \frac { 1 } { 1 - e ^ { - \beta \omega } } \right ) ^ { d } \ \ ( 2 2 )$$

$$= 2 \left ( \frac { 1 } { 2 \sinh ( \beta \omega / 2 ) } \right ) ^ { d }.$$

Thus, the virial coefficients of a trapped noninteracting spin-1/2 Fermi gas in d dimensions are

$$b _ { n } ^ { 0 } = \frac { ( - 1 ) ^ { n + 1 } } { n } \left ( \frac { \sinh ( \beta \omega / 2 ) } { \sinh ( \beta \omega n / 2 ) } \right ) ^ { d }. \quad \ ( 2 4 )$$

Notably, in the limit βω glyph[lessmuch] 1, we obtain

$$b _ { n } ^ { 0 } \rightarrow ( - 1 ) ^ { n + 1 } \left ( \frac { 1 } { n } \right ) ^ { d + 1 }, \text{ \quad \ \ } ( 2 5 )$$

which agrees in d = 3 with the local density approximation result quoted in Ref. [4]. The simple result of Eq. (24) should be a textbook calculation, but it does not appear elsewhere, to the best of our knowledge. Note that for the homogeneous (i.e. untrapped) system, the noninteracting virial coefficients in d dimensions are

$$b _ { n } ^ { 0, \hom } = ( - 1 ) ^ { n + 1 } \left ( \frac { 1 } { n } \right ) ^ { \frac { d } { 2 } + 1 }, \text{ \quad \ \ } ( 2 6 )$$

such that b 0 n = b 0 hom , n n -d 2 for βω glyph[lessmuch] 1.

## B. Semiclassical lattice approximation

To calculate the interaction-induced change ∆ Q m,n , we implement an approximation which consists in keeping the leading term in the commutator expansion:

$$\text{oay} \quad e ^ { - \beta ( \hat { H } _ { 0 } + \hat { V } _ { \text{int} } ) } = e ^ { - \beta \hat { H } _ { 0 } } e ^ { - \beta \hat { V } _ { \text{int} } } \times e ^ { - \frac { \beta ^ { 2 } } { 2 } [ \hat { H } _ { 0 }, \hat { V } _ { \text{int} } ] } \times \dots, \ ( 2 7 )$$

where the higher orders involve exponentials of nested commutators of ˆ H 0 with ˆ V int . Thus, the leading order in this expansion consists in setting [ ˆ H ,V 0 ˆ int ] = 0, which corresponds to a semiclassical approximation. Another way to see this approximation is in terms of a TrotterSuzuki factorization, i.e.

$$e ^ { - \beta ( \hat { H } _ { 0 } + \hat { V } _ { \text{int} } ) } = \lim _ { n \to \infty } \left ( e ^ { - \beta \hat { H } _ { 0 } / n } e ^ { - \beta \hat { V } _ { \text{int} } / n } \right ) ^ { n }, \quad ( 2 8 )$$

where the n = 1 case is the leading order we pursue in this work and higher orders can be defined by increasing n .

## C. Example: Calculation of ∆ Q 1 1 , and ∆ b 2

In the approximation proposed above, the two-particle problem is analyzed as follows:

$$& \text{prod} \text{em is anayze} \text{ as follows:} \\ & Q _ { 1, 1 } = \text{Tr} \left [ e ^ { - \beta \hat { H } _ { 0 } } e ^ { - \beta \hat { V } _ { \text{int} } } \right ] & \text{Using} \quad & \text{Using} \\ & \quad = \sum _ { k _ { 1 }, k _ { 2 }, x _ { 1 }, x _ { 2 } } \langle k _ { 1 } k _ { 2 } | e ^ { - \beta \hat { H } _ { 0 } } | x _ { 1 } x _ { 2 } \rangle \langle x _ { 1 } x _ { 2 } | e ^ { - \beta \hat { V } _ { \text{int} } } | k _ { 1 } k _ { 2 } \rangle & \quad \Delta b _ { 2 } \\ & \quad = \sum _ { k _ { 1 }, k _ { 2 }, x _ { 1 }, x _ { 2 } } e ^ { - \beta ( E _ { k _ { 1 } } + E _ { k _ { 2 } } ) } M _ { x _ { 1 }, x _ { 2 } } | \langle k _ { 1 } k _ { 2 } | x _ { 1 } x _ { 2 } \rangle | ^ { 2 } ( 2 9 ) & \quad \text{where} \quad & \text{to} \quad & \text{wi} \quad & \Lambda \circ \,.$$

where we have inserted complete sets of states in coordinate space {| x x 1 2 〉} and in the basis | k k 1 2 〉 of eigenstates of ˆ H 0 , whose single-particle eigenstates | k 〉 have eigenvalues E k . We have also made use of the fact that ˆ V int is diagonal in coordinate space, such that

$$M _ { \mathbf x _ { 1 }, \mathbf x _ { 2 } } = 1 + C \delta _ { \mathbf x _ { 1 }, \mathbf x _ { 2 } }, \text{ \quad \ \ } ( 3 0 )$$

where C = ( e βg dD -1) glyph[lscript] d and glyph[lscript] is an ultraviolet length scale to be defined by our renormalization condition (see below).

Thus,

$$\Delta Q _ { 1, 1 } = C \sum _ { k _ { 1 }, k _ { 2 } } e ^ { - \beta ( E _ { k _ { 1 } } + E _ { k _ { 2 } } ) } | \phi _ { k _ { 1 } } ( { \mathbf x } ) | ^ { 2 } | \phi _ { k _ { 2 } } ( { \mathbf x } ) | ^ { 2 }, \ \ ( 3 1 ) \quad \ U \text{ Using}$$

and we will use normalized single-particle wavefunctions in cartesian coordinates which in 1D take the form

$$\phi _ { n } ( x ) = \frac { 1 } { \sqrt { 2 ^ { n } n! } } \left ( \frac { \omega } { \pi } \right ) ^ { 1 / 4 } e ^ { - \omega x ^ { 2 } / 2 } H _ { n } ( \sqrt { \omega } x ), \quad ( 3 2 ) \quad \text{As w} \\ \text{scrip} \quad.$$

where the H n are Hermite polynomials. The sums over k 1 and k 2 in Eq. (31) are independent and identical and take the form

$$\sum _ { n = 0 } ^ { \infty } e ^ { - \beta E _ { n } } | \phi _ { n } ( x ) | ^ { 2 } = \sqrt { \frac { \omega } { \pi } } e ^ { - \beta \omega / 2 - \omega x ^ { 2 } } G ( \beta \omega, \sqrt { \omega }, x ), \\ \text{for each certain dimension where the function } C \, \text{$C$ can}$$

for each cartesian dimension, where the function G can be calculated as a special case of Mehler's formula [5]:

$$G ( k, b, x ) \equiv \sum _ { n = 0 } ^ { \infty } \frac { e ^ { - k n } } { 2 ^ { n } n! } [ H _ { n } ( b x ) ] ^ { 2 } = \frac { \exp \left [ 2 b ^ { 2 } x ^ { 2 } / ( 1 + e ^ { k } ) \right ] } { \sqrt { 1 - e ^ { - 2 k } } }. \quad \text{the u} \quad \text{effect} \quad \text{Gaus}$$

This formula encodes the finite-temperature, singleparticle density matrix of a noninteracting, nonrelativistic system in a harmonic trapping potential, and therefore its use is essential in the calculations that follow.

Squaring the result, we obtain, in 1 spatial dimension,

$$( 2 8 ) & \quad \text{ouquing the re-un, we you can, in $1$ spavail unmemsion}, \\ \quad \Delta Q _ { 1 1 } & = 2 C \int _ { 0 } ^ { \infty } d x \frac { \omega } { 2 \pi \sinh ( \beta \omega ) } \exp \left [ - 2 \omega x ^ { 2 } \tanh ( \beta \omega / 2 ) \right ] \\ \text{use in} \\ \text{ easing} & \quad, \quad \text{$\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, \quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $ \quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad,$\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad, $\quad. \quad, $\quad. \quad, $\quad. \quad, $\quad. \quad, $\quad. \quad, $\quad. \quad, $\quad. \quad, $\quad. \quad, $\quad. \quad, $\quad. \quad, $\quad. \quad, $\quad. \quad, $\quad. \quad. \quad, $\quad. \quad. \quad, $\quad. \quad. \quad, $\quad. \quad. \quad. \quad, $\quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \ \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \ { \text{$3}$}$$

where we have performed the last Gaussian integral along with some hyperbolic function simplifications. Generalizing to d spatial dimensions is very simple in this case:

$$\Delta Q _ { 1 1 } = C \left [ \frac { \sqrt { \omega } } { 4 } \frac { 1 } { \sqrt { \pi \sinh ( \beta \omega ) } \sinh ( \beta \omega / 2 ) } \right ] ^ { d }.$$

Using Eq. (22), we find that Q 1 cancels exactly in the final expression, as expected, such that

$$\mathbf k _ { 2 } \rangle \quad & \Delta b _ { 2 } \equiv \frac { C } { 2 } \left [ \frac { \sqrt { \omega } } { 2 \sqrt { \pi } \sinh ( \beta \omega ) } \right ] ^ { d } = \frac { C } { \lambda _ { T } ^ { d } } \, \frac { 1 } { 2 } \left [ \frac { \beta \omega } { 2 \sinh ( \beta \omega ) } \right ] ^ { \frac { d } { 2 } }, \\ \quad & \text{where we have used the thermal wavelength } \lambda _ { T } = \sqrt { 2 \pi \beta }$$

where we have used the thermal wavelength λ T = 2 πβ to write the result in dimensionless form.

As mentioned above, we will use this result to connect to the physical coupling λ of a given system, as a renormalization condition at a given value of βω . To that end, we first solve for C/λ d T :

$$\frac { C } { \lambda _ { T } ^ { d } } = 2 \Delta b _ { 2 } \left [ \frac { 2 \sinh ( \beta \omega ) } { \beta \omega } \right ] ^ { \frac { d } { 2 } }. \quad \quad ( 3 7 )$$

In the unitary limit of the 3D Fermi gas [6], for instance, the exact answer for ∆ b 2 is known [4]:

$$\Delta b _ { 2 } = \frac { 1 } { 2 } \left ( \frac { e ^ { - \beta \omega / 2 } } { 1 + e ^ { - \beta \omega } } \right ) = \frac { 1 } { 4 \cosh ( \beta \omega / 2 ) }. \quad ( 3 8 )$$

Using that result in Eq. (37) yields

$$\frac { C } { \lambda _ { T } ^ { 3 } } = \frac { e ^ { - \beta \omega / 2 } } { 1 + e ^ { - \beta \omega } } \left [ \frac { 2 \sinh ( \beta \omega ) } { \beta \omega } \right ] ^ { 3 / 2 }. \quad \ ( 3 9 )$$

As we will show below, this type of renormalization prescription is very practical as our results for ∆ b 3 and ∆ b 4 are simple quadratic functions of C/λ d T with βω -dependent coefficients.

## III. RESULTS

Following the steps outlined above in the example calculation of ∆ b 2 , we have calculated the various contributions to ∆ b 3 and ∆ b 4 , which we present in this section. In all cases, the central component of the calculation is the use of the analytic form of Mehler's kernel, which effectively reduces the calculation to a small number of Gaussian integrals.

## A. Result for ∆ Q 2 1 , and ∆ b 3

With small modifications to the example for ∆ b 2 , it is straightforward to show that

$$\Delta Q _ { 2, 1 } & = \frac { C } { \lambda _ { T } ^ { d } } \left [ \frac { \beta \omega } { 4 \sinh ^ { 3 } ( \beta \omega ) } \right ] ^ { \frac { d } { 2 } } \times \\ & \left [ \frac { 1 } { 2 ^ { \frac { d } { 2 } } \tanh ^ { d } \left ( \frac { \beta \omega } { 2 } \right ) } - \left ( \frac { 2 \cosh ^ { 2 } ( \beta \omega / 2 ) } { 4 \cosh ^ { 2 } ( \beta \omega / 2 ) - 1 } \right ) ^ { d / 2 } \right ] ( 4 0 ) \quad \text{We} \\ & \text{case,}$$

and using the results of the previous section, it is easy to assemble the final answer for ∆ b 3 in our approximation using

$$\Delta b _ { 3 } = \frac { 2 \Delta Q _ { 2, 1 } } { Q _ { 1 } } - \Delta b _ { 2 } Q _ { 1 }. \quad \quad ( 4 1 ) \quad \text{exp}$$

Indeed, after simplifications, we obtain

$$\text{it is} \quad \Delta b _ { 3 } = - \frac { C } { \lambda _ { T } ^ { d } } \left [ \frac { \beta \omega } { 2 \sinh ( \beta \omega ) } \right ] ^ { \frac { d } { 2 } } \left ( \frac { 1 } { 4 \cosh ^ { 2 } ( \beta \omega / 2 ) - 1 } \right ) ^ { d / 2 }, \\ \text{with} \quad \text{with} \text{is manipulate} \text{$two$ finite in $the$ $\mathcal{A}$} \rightharpoonup 0 limit $10$ that}$$

which is manifestly finite in the βω → 0 limit. In that limit,

$$\Delta b _ { 3 } \rightarrow - \frac { C } { \lambda _ { T } ^ { d } } \frac { 1 } { 6 ^ { 4 } }.$$

We recall the result of Ref. [1] for the homogeneous case, namely

$$\Delta b _ { 3 } ^ { \hom } = - \frac { C } { \lambda _ { T } ^ { d } } \frac { 1 } { 2 ^ { \frac { d } { 2 } } },$$

which shows that the relationship between the homogeneous and trapped cases, pointed out in the introduction, is also satisfied once interactions are turned on, as expected.

## B. Result for ∆ Q 3 1 , , ∆ Q 2 2 , , and ∆ b 4 .

Note that Q 1 diverges in the βω → 0 limit, but that divergence will cancel out in the final expression for ∆ b 3 .

Again following the steps outlined above, we obtain

$$\frac { 2 \Delta Q _ { 3, 1 } } { Q _ { 1 } } = \frac { C } { 2 \lambda _ { T } ^ { d } } \left [ \frac { \sqrt { 2 \beta \omega \sinh ( \beta \omega / 2 ) } } { 2 \sinh ^ { 2 } ( \beta \omega ) } \right ] ^ { d } \times \left \{ \left ( \frac { 1 } { 2 \tanh ^ { 3 } ( \beta \omega / 2 ) } \right ) ^ { \frac { 4 } { 4 } } + 2 \tanh ^ { \frac { 4 } { 2 } } ( \beta \omega ) \left ( \frac { \cosh ( \beta \omega ) + 1 } { 2 \cosh ( \beta \omega ) + 1 } \right ) ^ { \frac { 4 } { 4 } } \\ - 2 \coth ^ { \frac { 4 } { 4 } } ( \beta \omega / 2 ) \left ( \frac { \cosh ( \beta \omega ) + 1 } { 2 \cosh ( \beta \omega ) + 1 } \right ) ^ { \frac { 4 } { 4 } } - \left ( \frac { 1 } { 2 \tanh ( \beta \omega / 2 ) } \right ) ^ { \frac { 4 } { 4 } } \right \},$$

$$\frac { \Delta Q _ { 2, 2 } } { Q _ { 1 } } = \frac { C } { 2 \lambda _ { T } ^ { d } } \left [ \frac { \sqrt { 2 \beta \omega } \sinh ( \beta \omega / 2 ) } { 2 \sinh ^ { 2 } ( \beta \omega ) } \right ] ^ { d } \times \left \{ \left ( \frac { 1 } { 2 \tanh ^ { 3 } ( \beta \omega / 2 ) } \right ) ^ { \frac { 4 } { 2 } } - 2 \coth ^ { \frac { 4 } { 2 } } ( \beta \omega / 2 ) \left ( \frac { \cosh ( \beta \omega ) + 1 } { 2 \cosh ( \beta \omega ) + 1 } \right ) ^ { \frac { 4 } { 2 } } + \left ( \frac { \tanh ( \beta \omega ) } { 2 } \right ) ^ { \frac { 4 } { 2 } } \\ + \left ( \frac { C } { 2 \lambda _ { T } ^ { d } } \right ) \frac { ( 2 \beta \omega ) ^ { \frac { 4 } { 2 } } } { 2 ^ { d } } \left [ \frac { 1 } { \tanh ^ { d } ( \beta \omega / 2 ) } + 1 - 2 \left ( \sech ( \beta \omega ) + 1 \right ) ^ { \frac { 4 } { 2 } } \right ] \right \}.$$

Combining these with results from the previous sections, the final answer for ∆ b 4 can be assembled using while the final result for ∆ b 4 is finite), we obtain

$$\Delta b _ { 4 } = \frac { 2 \Delta Q _ { 3, 1 } + \Delta Q _ { 2, 2 } } { Q _ { 1 } } - \Delta \left ( b _ { 3 } + \frac { b _ { 2 } ^ { 2 } } { 2 } \right ) Q _ { 1 } - \frac { \Delta b _ { 2 } } { 2 } Q _ { 1 } ^ { 2 } ( 4 7 ) \quad \Delta b _ { 4 } =$$

After several simplifications and cancellations (which can be tracked by their degree of divergence as βω → 0,

$$\hat { \Omega } _ { 1 } ^ { 2 } ( 4 7 ) \quad & \Delta b _ { 4 } = \frac { C } { 2 \lambda _ { T } ^ { d } } \left [ \frac { \sqrt { 2 \beta \omega } \sinh ( \beta \omega / 2 ) } { 2 \sinh ^ { 2 } ( \beta \omega ) } \right ] ^ { d } \times \\ & \quad \left \{ 2 \tanh ^ { \frac { d } { 2 } } ( \beta \omega ) \left ( \frac { \cosh ( \beta \omega ) + 1 } { 2 \cosh ( \beta \omega ) + 1 } \right ) ^ { \frac { d } { 2 } } + \left ( \frac { \tanh ( \beta \omega ) } { 2 } \right ) ^ { \frac { d } { 2 } } \\ \text{which} \quad & \quad + \left ( \frac { C } { 2 \lambda _ { T } ^ { d } } \right ) \frac { ( 2 \beta \omega ) ^ { \frac { d } { 2 } } } { 2 ^ { d } } \left [ 1 - 2 \left ( \text{sech} ( \beta \omega ) + 1 \right ) ^ { \frac { d } { 2 } } \right ] \right \}. \quad ( 4 8 )$$

In this case, the limit βω → 0 yields

$$\Delta b _ { 4 } & \rightarrow \frac { C } { \lambda _ { T } ^ { d } } 4 ^ { - \frac { d } { 2 } } ( 3 ^ { - \frac { d } { 2 } } + 2 ^ { - d - 1 } ) \\ & \quad + \left ( \frac { C } { \lambda _ { T } ^ { d } } \right ) ^ { 2 } 4 ^ { - d - 1 } \left ( 1 - 2 ^ { \frac { d } { 2 } + 1 } \right ). \quad \left ( 4 9 \right ) \quad \left ( \frac { \tau } { \tau } \right ) ^ { 2 }.$$

Once again, we recall the homogeneous result:

$$\Delta b _ { 4 } ^ { \hom } & = \frac { C } { \lambda _ { T } ^ { d } } ( 3 ^ { - \frac { c } { 2 } } + 2 ^ { - d - 1 } ) &. \\ & \quad + \left ( \frac { C } { \lambda _ { T } ^ { d } } \right ) ^ { 2 } 4 ^ { - \frac { c } { 2 } - 1 } ( 1 - 2 ^ { \frac { c } { 2 } + 1 } ), \quad ( 5 0 )$$

and find that ∆ b 4 = ∆ b hom 4 4 -d 2 , as expected.

## C. Results in terms of ∆ b 2 .

Finally, collecting our results for ∆ b 3 and ∆ b 4 and expressing them in terms of ∆ b 2 [via Eq. (37)], we obtain

$$\Delta b _ { 3 } = - 2 \left ( \frac { 1 } { 4 \cosh ^ { 2 } ( \beta \omega / 2 ) - 1 } \right ) ^ { \frac { d } { 2 } } \Delta b _ { 2 }, \quad ( 5 1 )$$

and

$$\Delta b _ { 4 } & = \left ( \frac { 1 } { 2 \cosh ( \beta \omega ) } \right ) ^ { \frac { d } { 2 } } f _ { 1 } ( \beta \omega, d ) \Delta b _ { 2 } + f _ { 2 } ( \beta \omega, d ) ( \Delta b _ { 2 } ) ^ { 2 }, \quad \begin{matrix} \text{FIG. 1} \\ \text{of the} \\ \Delta b _ { 2 } \, \text{co} \end{matrix}$$

where

$$& f _ { 1 } ( \beta \omega, d ) = \left ( \frac { 1 } { 2 \cosh ( \beta \omega / 2 ) } \right ) ^ { d } + 2 \left ( \frac { 1 } { 4 \cosh ^ { 2 } ( \beta \omega / 2 ) - 1 } \right ) ^ { \frac { d } { 2 } }, \\ & \text{and}$$

and

$$f _ { 2 } ( \beta \omega, d ) = \left ( \frac { 1 } { 2 \cosh ( \beta \omega / 2 ) } \right ) ^ { d } - 2 \left ( \frac { 1 } { 4 \cosh ^ { 2 } ( \beta \omega / 2 ) - 2 } \right ) ^ { \frac { d } { 2 } } \text{.} \quad \text{sure} \\ \text{$\mathfrak{T}$} \colon \text{$\mathfrak{T}$} \dots \text{$\mathfrak{T}$} \dots \text{$\mathfrak{T}$} \dots \text{$\mathfrak{T}$} \dots \text{$\mathfrak{T}$} \dots \text{$\mathfrak{T}$} \dots \text{$\mathfrak{T}$} \dots \text{$\mathfrak{T}$} \dots \text{$\mathfrak{T}$} \dots \text{$\mathfrak{T}$} \dots \text{$\mathfrak{T}$} \dots \text{$\mathfrak{T}$} \dots \text{$\mathfrak{T}$} \dots \text{$\mathf
fixed}.$$

The above formulas for ∆ b 3 and ∆ b 4 are the main result of this work. In the following we explore their behavior as a function of d and βω , focusing in particular on the unitary limit of the 3D Fermi gas. While numerical results exist for these quantities in some cases, in particular in 2D [9-12] (see also [13-15]) and in 3D at unitarity [16-18, 20], most of those correspond to homogeneous systems and do not feature explicit, analytic dependence on the dimension nor on βω , as shown here. Our results are therefore useful in that they are able to provide analytic insight into the behavior of virial coefficients across dimensions, and as a function of the temperature (or trapping frequency) as well as the coupling strength. Below, we evaluate our formulas and discuss the resulting answers.

.

FIG. 1. ∆ b 3 / ∆ b 2 (top) and ∆ b 4 / ∆ b 2 (bottom) as functions of the spatial dimension d at fixed βω = 0 1, 1 0, and 5 0, for . . . ∆ b 2 corresponding to the unitary limit in d = 3, Eq. (38).

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00070v1-figure-a31dd8f9-c97c-41a5-a142-8fa7a7f5016c.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00070v1-figure-35cc3149-145d-4a4f-99c4-f3644b708ad5.jpg)

## D. Qualitative behavior.

To illustrate our analytic results, in Fig. 1 we show ∆ b 3 / ∆ b 2 and ∆ b 4 / ∆ b 2 as a function of the spatial dimension d , at various βω , fixing ∆ b 2 to its value in the unitary limit (as a reference point). We find that, as d increases, the magnitude of the interaction-induced change ∆ b n decreases. This suggests that, using ∆ b 2 as the fixed, dimension-independent coupling, the radius of convergence of the virial expansion increases with d . This is consistent with the idea that, in higher dimensions, the kinetic energy dominates over the interactions and mean-field type of approaches capture the behavior of the system correctly.

As a comparison with previous calculations, we show in Fig. 2 our results in 3D at unitarity as a function of βω , superimposed with the data from Ref. [7]. While we do not expect, a priori, good quantitative agreement in this strong coupling regime, we find at least qualitative agreement for both ∆ b 3 and ∆ b 4 , and surprisingly good agreement at the quantitative level for ∆ b 4 . Clearly, the LO-SCLA is able to capture more than just the shape of the βω dependence of the virial coefficients.

FIG. 2. Comparison of our LO-SCLA results for ∆ b 3 (top) and ∆ b 4 (bottom) as a function of βω . The high-temperature fits and PIMC results are from Ref. [7]; the sum-over-states results are from Ref. [8].

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00070v1-figure-3ebfdd52-9fdc-4079-9211-ed0dec80c771.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00070v1-figure-430672a4-4efe-4724-9e04-db29b229fd3e.jpg)

## IV. SUMMARY AND CONCLUSIONS

In this work we have implemented a semiclassical approximation, at leading order, to calculate the virial coefficients ∆ b 3 and ∆ b 4 of harmonically trapped Fermi gases. Our calculations yield analytic answers as functions of βω and, by a renormalization prescription that matches ∆ b 2 to the known exact result, we also obtain

[1] C. R. Shill, J. E. Drut, Virial coefficients of 1D and 2D Fermi gases by stochastic methods and a semiclassical lattice approximation , Phys. Rev. A 98 , 053615 (2018).

[2] Y. Hou, A. Czejdo, J. DeChant, C. R. Shill, J. E. Drut, Leading-order semiclassical approximation to the first seven virial coefficients of spin-1/2 fermions across spatial dimensions , arXiv:1907.10120.

[3] Y. Hou, J. E. Drut, Semiclassical approximation to virial coefficients beyond the leading order , arXiv:1908.xxxx.

[4] X.-J. Liu, Virial expansion for a strongly correlated Fermi system and its application to ultracold atomic Fermi gases , Phys. Rep. 524 , 37 (2013).

the dependence on the physical coupling strength. Notably, our results are also analytic functions of the spatial dimension d , allowing us to study the behavior of the virial expansion across dimensions. We find that, at fixed ∆ b 2 , the magnitude of ∆ b n decreases as d increases, for all βω .

Although there have been many (and very precise) determinations of virial coefficients in the literature, they are mostly numerical and focus on specific dimensions or couplings (and most of them are for homogeneous systems). Our approach and results are, in that sense, complementary: we do not expect high precision from the LO-SCLA, but through it we are able to study, explicitly, the variations with the parameters of the problem, which yield qualitative analytic insight into the properties of the virial expansion. We have demonstrated the quality of our leading-order (!) results for ∆ b 3 and ∆ b 4 in the unitary limit by showing that they qualitatively follow the expected answers, which is an encouraging sign to proceed to next-to-leading order in future work. Furthermore, the approximate agreement with prior results at unitarity suggests that, between the noninteracting regime and the unitary point, that agreement should be even better than shown here.

Finally, it should be pointed out that the renormalization prescription based on ∆ b 2 does not by itself eliminate all the lattice artifacts. Future studies should explore the use of improved actions (see e.g. [21, 22]), potentially making use of prior knowledge of ∆ b 3 were available, to enhance the quality of the expansion.

## ACKNOWLEDGMENTS

This material is based upon work supported by the National Science Foundation under Grant No. PHY1452635 (Computational Physics Program). C.E.B. acknowledges support from the United States Department of Energy through the Computational Science Graduate Fellowship (DOE CSGF) under grant number DE-FG0297ER25308.

[5] H. Bateman, A. Erd´ elyi Higher transcendental functions. II , (McGraw Hill, New York, 1953)

- [6] The BCS-BEC Crossover and the Unitary Fermi Gas , edited by W. Zwerger (Springer-Verlag, Berlin, 2012).

[7] Y. Yan, D. Blume, Path integral Monte Carlo determination of the fourth-order virial coefficient for unitary twocomponent Fermi gas with zero-range interactions , Phys. Rev. Lett. 116 , 230401 (2016).

[8] D. Rakshit, K. M. Daily, and D. Blume, Natural and unnatural parity states of small trapped equal-mass twocomponent Fermi gases at unitarity and fourth-order virial coefficient , Phys. Rev. A 85 , 033634 (2012).

- [9] X.-J. Liu, H. Hu, and P. D. Drummond, Exact few-body results for strongly correlated quantum gases in two dimensions , Phys. Rev. B 82 , 054524 (2010)
- [10] V. Ngampruetikorn, J. Levinsen, and M. M. Parish, Pair correlations in the two-dimensional Fermi gas , Phys. Rev. Lett. 111 , 265301 (2013).
- [11] M. Barth and J. Hofmann, Pairing effects in the nondegenerate limit of the two-dimensional Fermi gas , Phys. Rev. A. 89 , 013614 (2014).
- [12] V. Ngampruetikorn, M. M. Parish, and J. Levinsen, High-temperature limit of the resonant Fermi gas , Phys. Rev. A 91 , 013606 (2015).
- [13] C. Chaffin and T. Sch¨fer, a Scale breaking and fluid dynamics in a dilute two-dimensional Fermi gas , Phys. Rev. A 88 , 043636 (2013).
- [14] W. S. Daza, J. E. Drut, C. L. Lin, and C. R. Ord´˜ez, on Virial expansion for the Tan contact and Beth-Uhlenbeck formula from two-dimensional SO(2,1) anomalies Phys. Rev. A 97 , 033630 (2018).
- [15] C. R. Ordo˜ nez, Path-integral Fujikawa approach to anomalous virial theorems and equations of state for sys-
- tems with SO(2,1) symmetry , Physica, 446 , 64 (2016).
- [16] D. Lee and T. Sch¨fer, a Cold dilute neutron matter on the lattice. I. Lattice virial coefficients and large scattering lengths Phys. Rev. C 73 , 015201 (2006).
- [17] X.-J. Liu, H. Hu, and P. D. Drummond, Virial expansion for a strongly correlated Fermi gas , Phys. Rev. Lett. 102 , 160401 (2009).
- [18] X.-J. Liu, H. Hu, and P. D. Drummond, Three attractively interacting fermions in a harmonic trap: Exact solution, ferromagnetism, and high-temperature thermodynamics , Phys. Rev. A 82 , 023619 (2010).
- [19] X. Leyronas, Virial expansion with Feynman diagrams , Phys. Rev. A 84 , 053633 (2011).
- [20] D. B. Kaplan, S. Sun, A new field theoretic method for the virial expansion , Phys. Rev. Lett. 107 , 030601 (2011).
- [21] J. E. Drut and A. N. Nicholson, Lattice methods for strongly interacting many-body systems , J. Phys. G: Nucl. Part. Phys. 40 , 043101 (2013).
- [22] J. E. Drut Improved lattice operators for non-relativistic fermions , Phys. Rev. A 86 , 013604 (2012).
K. J. Morrell C. E. Berger J. E. Drut
Third- and fourth-order virial coefficients of harmonically trapped fermions in a semiclassical approximation
Using a leading-order semiclassical approximation, we calculate the third-
and fourth-order virial coefficients of nonrelativistic spin-1/2 fermions in a
harmonic trapping potential in arbitrary spatial dimensions, and as functions
of temperature, trapping frequency and coupling strength. Our simple, analytic
results for the interaction-induced changes $\Delta b_3$ and $\Delta b_4$ agree
qualitatively, and in some regimes quantitatively, with previous numerical
calculations for the unitary limit of three-dimensional Fermi gases.