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

TL;DR

This paper derives analytical expressions for third- and fourth-order virial coefficients of trapped fermions using a semiclassical approximation, providing insights into their temperature, trapping, and interaction dependence.

Contribution

It introduces a simple semiclassical method to calculate higher-order virial coefficients for trapped fermions across various dimensions and interaction strengths.

Findings

Analytic formulas for virial coefficients agree with numerical results in certain regimes.

Interaction effects on virial coefficients are characterized as functions of temperature and coupling.

Results are applicable to understanding thermodynamics of trapped Fermi gases.

Abstract

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.

Equations

$$\hat{H}={\hat{H}}_0+{\hat{V}}_{\text{int}},$$

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

$$\hat{T}=\sum _{s=1,2}\int d^{d}x{\hat{\psi }}_s^{†}(\mathbf{x})\left(-\frac{{\hslash }^2{\nabla }^2}{2m}\right){\hat{\psi }}_s(\mathbf{x}),$$

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

$${\hat{V}}_{\text{int}}=-g_{dD}\int d^{d}x{\hat{n}}_1(\mathbf{x}){\hat{n}}_2(\mathbf{x}),$$

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

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

$$\mathcal{Z}=\sum _{N=0}^{\infty }{z}^N{Q}_N,$$

$$b_1$$

$$b_2$$

$$b_3$$

$$b_4$$

$$Q_1$$

$$Q_2$$

$$Q_3$$

$$Q_4$$

$$\Delta b_2$$

$$\Delta b_3$$

$$\Delta b_4$$

$$\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).$$

$$Q_1\sum _{k=1}^{\infty }{b}_k^0{z}^k$$

$$Q_1$$

$$=$$

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

$$b_n^0\to (-1{)}^{n+1}{\left(\frac{1}{n}\right)}^{d+1},$$

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

$$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 ,$$

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

$$Q_{1,1}$$

$$M_{{\mathbf{x}}_1,{\mathbf{x}}_2}=1+C{\delta }_{{\mathbf{x}}_1,{\mathbf{x}}_2},$$

$$\Delta Q_{1,1}=C\sum _{{\mathbf{k}}_1,{\mathbf{k}}_2}{e}^{-\beta (E_{{\mathbf{k}}_1}+E_{{\mathbf{k}}_2})}|{\varphi }_{{\mathbf{k}}_1}(\mathbf{x}){|}^2|{\varphi }_{{\mathbf{k}}_2}(\mathbf{x}){|}^2,$$

$$\phi_{n}(x)=\frac{1}{\sqrt{2^{n}n!}}\bigg{(}\frac{\omega}{\pi}\bigg{)}^{1/4}e^{-\omega x^{2}/2}H_{n}(\sqrt{\omega}x),$$

$$\sum _{n=0}^{\infty }{e}^{-\beta E_n}|{\varphi }_n(x){|}^2=\sqrt{\frac{\omega }{\pi }}{e}^{-\beta \omega /2-\omega x^2}G(\beta \omega ,\sqrt{\omega },x),$$

$$G(k,b,x)\equiv \sum _{n=0}^{\infty }\frac{e^{-kn}}{2^{n}n!}[H_n(bx){]}^2$$

$$\Delta Q_{11}$$

Full text

I Introduction

The properties of fermions at finite temperature and density are in part governed by the dimensionless product $\beta\mu$, where $\beta=1/(k_{B}T)$ is the inverse temperature and $\mu$ is the chemical potential. Typically, the region $\beta\mu\simeq 0$ displays a crossover between quantum and classical physics, while $z=e^{\beta\mu}\ll 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 $\mu$ 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 $\Delta b_{n}$.

The previous work of Ref. ShillDrut calculated the third- and fourth-order virial coefficients $\Delta b_{3}$ and $\Delta 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. HouEtAl extended those results up to $\Delta b_{7}$, while Ref. HouDrut carried out calculations up to next-to-next-to-leading order in the SCLA for up to $\Delta b_{5}$. In this brief work we provide another piece of the puzzle by generalizing the calculations of Ref. ShillDrut to systems in a harmonic trap of frequency $\omega$. We present our derivations with intermediate steps in detail and give analytic formulas for $\Delta b_{3}$ and $\Delta b_{4}$ as functions of $\beta\omega$ in arbitrary spatial dimension $d$. Our results, which will be given in terms of $\Delta 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

[TABLE]

where

[TABLE]

and

[TABLE]

is the kinetic energy,

[TABLE]

is the spherically symmetric external trapping potential, and

[TABLE]

is the interaction.

In the above equations, the field operators $\hat{\psi}_{s},\hat{\psi}^{\dagger}_{s}$ correspond to particles of species $s=1,2$, and $\hat{n}_{s}({\bf x})$ are the coordinate-space densities. For the remainder of this work, we will set $\hbar=k_{\text{B}}=m=1$.

II.1 Thermodynamics and the virial expansion

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

[TABLE]

where $\beta$ is the inverse temperature, $\Omega$ is the grand thermodynamic potential, $\hat{N}$ is the total particle number operator, and $\mu$ is the overall chemical potential (we will not consider polarized systems in this work).

As the calculation of $\mathcal{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\to 0$, where $z=e^{\beta\mu}$ is the fugacity, i.e. it is a low-fugacity expansion (see Ref. VirialReview for a review on recent applications of the virial expansion to ultracold atoms). The coefficients accompanying the powers of $z$ in the expansion $\Omega$ are the virial coefficients $b_{n}$:

[TABLE]

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

[TABLE]

we obtain expressions for the virial coefficients

[TABLE]

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:

[TABLE]

and so on for higher orders. In the absence of intra-species 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

[TABLE]

To calculate $\Delta Q_{m,n}$, we implement a semiclassical approximation, as described in the next section. Once we obtain the virial coefficients, one may rebuild the grand-canonical potential $\Omega$ 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 $\Delta b_{2}$ as a renormalization condition by relying on the exact answers as functions of $\beta\omega$ and the physical coupling $\lambda$. 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 $\beta\omega$. [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,

[TABLE]

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

[TABLE]

To identify the noninteracting virial coefficients $b^{0}_{n}$, we need $Q_{1}$:

[TABLE]

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

[TABLE]

Notably, in the limit $\beta\omega\ll 1$, we obtain

[TABLE]

which agrees in $d=3$ with the local density approximation result quoted in Ref. VirialReview . 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

[TABLE]

such that $b_{n}^{0}=b_{n}^{0,\text{hom}}n^{-\frac{d}{2}}$ for $\beta\omega\ll 1$.

II.2 Semiclassical lattice approximation

To calculate the interaction-induced change $\Delta Q_{m,n}$, we implement an approximation which consists in keeping the leading term in the commutator expansion:

[TABLE]

where the higher orders involve exponentials of nested commutators of $\hat{H}_{0}$ with $\hat{V}_{\text{int}}$. Thus, the leading order in this expansion consists in setting $[\hat{H}_{0},\hat{V}_{\text{int}}]=0$, which corresponds to a semiclassical approximation. Another way to see this approximation is in terms of a Trotter-Suzuki factorization, i.e.

[TABLE]

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

II.3 Example: Calculation of $\Delta Q_{1,1}$ and $\Delta b_{2}$

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

[TABLE]

where we have inserted complete sets of states in coordinate space $\{|{\bf x}_{1}{\bf x}_{2}\rangle\}$ and in the basis $|{\bf k}_{1}{\bf k}_{2}\rangle$ of eigenstates of $\hat{H}_{0}$, whose single-particle eigenstates $|{\bf k}\rangle$ have eigenvalues $E_{\bf k}$. We have also made use of the fact that $\hat{V}_{\text{int}}$ is diagonal in coordinate space, such that

[TABLE]

where $C=(e^{\beta g_{dD}}-1)\ell^{d}$ and $\ell$ is an ultraviolet length scale to be defined by our renormalization condition (see below).

Thus,

[TABLE]

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

[TABLE]

where the $H_{n}$ are Hermite polynomials. The sums over ${\bf k}_{1}$ and ${\bf k}_{2}$ in Eq. (31) are independent and identical and take the form

[TABLE]

for each cartesian dimension, where the function $G$ can be calculated as a special case of Mehler’s formula MehlerFormula :

[TABLE]

This formula encodes the finite-temperature, single-particle 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,

[TABLE]

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:

[TABLE]

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

[TABLE]

where we have used the thermal wavelength $\lambda_{T}=\sqrt{2\pi\beta}$ to write the result in dimensionless form.

As mentioned above, we will use this result to connect to the physical coupling $\lambda$ of a given system, as a renormalization condition at a given value of $\beta\omega$. To that end, we first solve for $C/\lambda_{T}^{d}$:

[TABLE]

In the unitary limit of the 3D Fermi gas ZwergerBook , for instance, the exact answer for $\Delta b_{2}$ is known VirialReview :

[TABLE]

Using that result in Eq. (37) yields

[TABLE]

As we will show below, this type of renormalization prescription is very practical as our results for $\Delta b_{3}$ and $\Delta b_{4}$ are simple quadratic functions of ${C}/{\lambda_{T}^{d}}$ with $\beta\omega$-dependent coefficients.

III Results

Following the steps outlined above in the example calculation of $\Delta b_{2}$, we have calculated the various contributions to $\Delta b_{3}$ and $\Delta 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.

III.1 Result for $\Delta Q_{2,1}$ and $\Delta b_{3}$

With small modifications to the example for $\Delta b_{2}$, it is straightforward to show that

[TABLE]

and using the results of the previous section, it is easy to assemble the final answer for $\Delta b_{3}$ in our approximation using

[TABLE]

Note that $Q_{1}$ diverges in the $\beta\omega\to 0$ limit, but that divergence will cancel out in the final expression for $\Delta b_{3}$. Indeed, after simplifications, we obtain

[TABLE]

which is manifestly finite in the $\beta\omega\to 0$ limit. In that limit,

[TABLE]

We recall the result of Ref. ShillDrut for the homogeneous case, namely

[TABLE]

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.

III.2 Result for $\Delta Q_{3,1}$, $\Delta Q_{2,2}$, and $\Delta b_{4}$.

Again following the steps outlined above, we obtain

[TABLE]

[TABLE]

Combining these with results from the previous sections, the final answer for $\Delta b_{4}$ can be assembled using

[TABLE]

After several simplifications and cancellations (which can be tracked by their degree of divergence as $\beta\omega\to 0$, while the final result for $\Delta b_{4}$ is finite), we obtain

[TABLE]

In this case, the limit $\beta\omega\to 0$ yields

[TABLE]

Once again, we recall the homogeneous result:

[TABLE]

and find that $\Delta b_{4}=\Delta b_{4}^{\text{hom}}4^{-\frac{d}{2}}$, as expected.

III.3 Results in terms of $\Delta b_{2}$.

Finally, collecting our results for $\Delta b_{3}$ and $\Delta b_{4}$ and expressing them in terms of $\Delta b_{2}$ [via Eq. (37)], we obtain

[TABLE]

and

[TABLE]

where

[TABLE]

and

[TABLE]

The above formulas for $\Delta b_{3}$ and $\Delta b_{4}$ are the main result of this work. In the following we explore their behavior as a function of $d$ and $\beta\omega$, 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 DrummondVirial2D ; virial2D2 ; PhysRevA.89.013614 ; Ngampruetikorn (see also virial2D ; Daza2D ; Ordo ) and in 3D at unitarity LeeSchaeferPRC1 ; LiuHuDrummond ; LiuHuDrummond2 ; DBK , most of those correspond to homogeneous systems and do not feature explicit, analytic dependence on the dimension nor on $\beta\omega$, 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.

III.4 Qualitative behavior.

To illustrate our analytic results, in Fig. 1 we show $\Delta b_{3}/\Delta b_{2}$ and $\Delta b_{4}/\Delta b_{2}$ as a function of the spatial dimension $d$, at various $\beta\omega$, fixing $\Delta 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 $\Delta b_{n}$ decreases. This suggests that, using $\Delta 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 $\beta\omega$, superimposed with the data from Ref. Doerte . While we do not expect, a priori, good quantitative agreement in this strong coupling regime, we find at least qualitative agreement for both $\Delta b_{3}$ and $\Delta b_{4}$, and surprisingly good agreement at the quantitative level for $\Delta b_{4}$. Clearly, the LO-SCLA is able to capture more than just the shape of the $\beta\omega$ dependence of the virial coefficients.

IV Summary and Conclusions

In this work we have implemented a semiclassical approximation, at leading order, to calculate the virial coefficients $\Delta b_{3}$ and $\Delta b_{4}$ of harmonically trapped Fermi gases. Our calculations yield analytic answers as functions of $\beta\omega$ and, by a renormalization prescription that matches $\Delta b_{2}$ to the known exact result, we also obtain 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 $\Delta b_{2}$, the magnitude of $\Delta b_{n}$ decreases as $d$ increases, for all $\beta\omega$.

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 $\Delta b_{3}$ and $\Delta 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 $\Delta b_{2}$ does not by itself eliminate all the lattice artifacts. Future studies should explore the use of improved actions (see e.g. DrutNicholson ; Drut ), potentially making use of prior knowledge of $\Delta b_{3}$ were available, to enhance the quality of the expansion.

Acknowledgements.

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-FG02-97ER25308.