Leading- and next-to-leading order semiclassical approximation to the first seven virial coefficients of spin-1/2 fermions across spatial dimensions

TL;DR

This paper develops a semiclassical approximation method to calculate virial coefficients of spin-1/2 fermions with contact interactions across different dimensions, providing analytic formulas and analyzing the virial expansion's convergence.

Contribution

It introduces leading- and next-to-leading order semiclassical formulas for virial coefficients up to seventh order, relating them to the second-order coefficient across dimensions.

Findings

Derived analytic expressions for virial coefficients up to seventh order.

Analyzed the virial expansion's convergence and its dependence on coupling strength.

Applied the results to calculate pressure and Tan's contact in 2D Fermi gases.

Abstract

Following up on recent calculations, we investigate the leading- and next-to-leading order semiclassical approximation to the virial coefficients of a two-species fermion system with a contact interaction. Using the analytic result for the second-order virial coefficient as a renormalization condition, we derive expressions for up to the seventh-order virial coefficient $\Delta b_7$. Our results at leading order, though approximate, furnish simple analytic formulas that relate $\Delta b_n$ to $\Delta b_2$ for arbitrary dimension, providing a glimpse into the behavior of the virial expansion across dimensions and coupling strengths. As an application, we calculate the pressure and Tan's contact of the 2D attractive Fermi gas and examine the radius of convergence of the virial expansion as a function of the coupling strength.

Figures

Tables

Table 1. Table 1: 𝒪 ​ ( C ) 𝒪 𝐶 \mathcal{O}(C) Terms for Δ ​ b 2 Δ subscript 𝑏 2 \Delta b_{2} to Δ ​ b 7 Δ subscript 𝑏 7 \Delta b_{7}

$\Delta b_n$	$𝒪\left(C\right)$
$\Delta b_2$	$\frac{Q\left(\beta \right)}{4}$
$\Delta b_3$	$-\frac{Q\left(2\beta \right)}{2}$
$\Delta b_4$	$\frac{Q\left(3\beta \right)}{2}+\frac{Q^2\left(2\beta \right)}{4Q\left(\beta \right)}$
$\Delta b_5$	$-\frac{Q\left(4\beta \right)}{2}-\frac{Q\left(2\beta \right)Q\left(3\beta \right)}{2Q\left(\beta \right)}$
$\Delta b_6$	$\frac{Q\left(5\beta \right)}{2}+\frac{Q\left(2\beta \right)Q\left(4\beta \right)}{2Q\left(\beta \right)}+\frac{Q{\left(3\beta \right)}^2}{4Q\left(\beta \right)}$
$\Delta b_7$	$-\frac{Q\left(6\beta \right)}{2}-\frac{Q\left(3\beta \right)Q\left(4\beta \right)+Q\left(2\beta \right)Q\left(5\beta \right)}{2Q\left(\beta \right)}$

Table 2. Table 2: 𝒪 ​ ( C 2 ) 𝒪 superscript 𝐶 2 \mathcal{O}(C^{2}) Terms for Δ ​ b 2 Δ subscript 𝑏 2 \Delta b_{2} to Δ ​ b 7 Δ subscript 𝑏 7 \Delta b_{7}

$\Delta b_n$	$𝒪\left(C^2\right)$
$\Delta b_4$	$\frac{G\left(M_4\right)}{2Q\left(\beta \right)}-\frac{Q\left(\beta \right)Q\left(2\beta \right)}{8}$
$\Delta b_5$	$-2\frac{G\left(M_5\right)}{Q\left(\beta \right)}+\frac{Q\left(\beta \right)Q\left(3\beta \right)}{4}+\frac{Q^2\left(2\beta \right)}{4}$
$\Delta b_6$	$\begin{array}{c}\frac{2G\left(M_{6a}\right)+2G\left(M_{6b}\right)+G\left(M_{6c}\right)}{Q\left(\beta \right)}\\ -\frac{3Q\left(4\beta \right)Q\left(\beta \right)}{8}-\frac{Q^3\left(2\beta \right)}{8Q\left(\beta \right)}-\frac{3Q\left(3\beta \right)Q\left(2\beta \right)}{4}\end{array}$
$\Delta b_7$	$\begin{array}{c}-2\frac{3G\left(M_{7a}^{\left(2\right)}\right)+G\left(M_{7b}^{\left(2\right)}\right)+G\left(M_{7c}^{\left(2\right)}\right)}{Q\left(\beta \right)}\\ +\frac{Q^2\left(3\beta \right)+Q\left(5\beta \right)Q\left(\beta \right)}{2}\\ +\frac{Q^2\left(2\beta \right)Q\left(3\beta \right)}{2Q\left(\beta \right)}+Q\left(4\beta \right)Q\left(2\beta \right)\end{array}$

Table 3. Table 3: 𝒪 ​ ( C 3 ) 𝒪 superscript 𝐶 3 \mathcal{O}(C^{3}) Terms for Δ ​ b 2 Δ subscript 𝑏 2 \Delta b_{2} to Δ ​ b 7 Δ subscript 𝑏 7 \Delta b_{7}

$\Delta b_n$	$𝒪\left(C^3\right)$
$\Delta b_6$	$\begin{array}{c}-G\left(M_5\right)+\frac{2F\left(N_1\right)}{3Q\left(\beta \right)}\\ +\frac{Q^2\left(2\beta \right)Q\left(\beta \right)}{16}+\frac{Q\left(3\beta \right)Q^2\left(\beta \right)}{24}\end{array}$
$\Delta b_7$	$\begin{array}{c}\frac{G\left(M_{7a}^{\left(3\right)}\right)Q\left(2\beta \right)}{Q\left(\beta \right)}-4\frac{F\left(N_2\right)}{Q\left(\beta \right)}\\ +3G\left(M_{7b}^{\left(3\right)}\right)+2G\left(M_{7c}^{\left(3\right)}\right)\\ -\frac{Q\left(4\beta \right)Q^2\left(\beta \right)+3Q\left(3\beta \right)Q\left(2\beta \right)Q\left(\beta \right)+Q^3\left(2\beta \right)}{8}\end{array}$

Table 4. Table 4: Full results in the LO-SCLA for Δ ​ b 3 Δ subscript 𝑏 3 \Delta b_{3} through Δ ​ b 7 Δ subscript 𝑏 7 \Delta b_{7} in powers of Δ ​ b 2 Δ subscript 𝑏 2 \Delta b_{2} in the continuum limit.

$\Delta b_n$	$𝒪\left(\Delta b_2\right)$	$𝒪{\left(\Delta b_2\right)}^2$	$𝒪{\left(\Delta b_2\right)}^3$
$\Delta b_3$	$-2^{1-d/2}$	–	–
$\Delta b_4$	$2\left(3^{-d/2}+2^{-d-1}\right)$	$2^{1-d/2}\left(2^{-d/2-1}-1\right)$	–
$\Delta b_5$	$-2\left(2^{-d}+6^{-d/2}\right)$	$4\left(2^{-d}+3^{-d/2}-7^{-d/2}\right)$	–
$\Delta b_6$	$2^{1-3d/2}+3^{-d}+2\cdot 5^{-d/2}$	$\begin{array}{c}-2^{1-3d/2}-3\cdot 2^{1-d}\left(1-3^{-d/2}\right)\\ -2^{2-d/2}\left(3^{1-d/2}-5^{-d/2}\right)\end{array}$	$\begin{array}{c}{2}^{2-d}-8\cdot 7^{-d/2}\\ +8\cdot 3^{-d/2-1}\left(1+2^{-d}\right)\end{array}$
$\Delta b_7$	$-2^{1-d/2}\left(6^{-d/2}+5^{-d/2}+3^{-d/2}\right)$	$\begin{array}{c}{2}^{4-3d/2}+8\cdot 5^{-d/2}\\ +8\cdot 3^{-d/2}\left(3^{-d/2}+2^{-d}\right)\\ -4\left(2^{-d}{5}^{-d/2}-{13}^{-d/2}-3\cdot {17}^{-d/2}\right)\end{array}$	$\begin{array}{c}-2^{3-3d/2}-2^{3-d}\left(1-3^{1-d/2}\right)\\ -8\cdot 2^{-d/2}\left(3^{1-d/2}-7^{-d/2}\right)\\ -16\cdot {10}^{-d/2}\left(2^{-d/2}-1\right)\end{array}$

Equations

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

$$\hat{V}=-g_d\int d^{d}x{\hat{n}}_{\uparrow }(\mathbf{x}){\hat{n}}_{\downarrow }(\mathbf{x}),$$

$$\hat{V}=-g_d\sum _{\mathbf{x}}{\hat{n}}_{\uparrow }(\mathbf{x}){\hat{n}}_{\downarrow }(\mathbf{x}),$$

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

$$\mathcal{Z}=\text{tr}\left[e^{-\beta (\hat{H}-\mu \hat{N})}\right]=\sum _{N=0}^{\infty }{z}^N{Q}_N,$$

$$Q_1{b}_2$$

$$Q_1{b}_3$$

$$Q_1{b}_4$$

$$Q_1=2\sum _{\mathbf{p}}{e}^{-\beta \frac{{\mathbf{p}}^2}{2m}},$$

$$Q_1\Delta b_2$$

$$Q_1\Delta b_3$$

$$Q_1\Delta b_4$$

$$\Delta Q_2$$

$$\Delta Q_3$$

$$\Delta Q_4$$

$$\Delta Q_5$$

$$\Delta Q_6$$

$$\Delta Q_7$$

$$e^{-\beta (\hat{T}+\hat{V})}=e^{-\beta \hat{T}}{e}^{-\beta \hat{V}},$$

$$e^{-\beta (\hat{T}+\hat{V})}=\underset{n\to \infty }{\lim }{\left(e^{-\beta \hat{T}/n}{e}^{-\beta \hat{V}/n}\right)}^n.$$

$$Q_{11}$$

$$e^{-\beta \hat{V}}|{\mathbf{x}}_1{\mathbf{x}}_2⟩$$

$$\Delta Q_{11}=C\sum _{{\mathbf{p}}_1{\mathbf{p}}_2,{\mathbf{x}}_1{\mathbf{x}}_2}{e}^{-\beta (p_1^2+p_2^2)/2m}{\delta }_{{\mathbf{x}}_1,{\mathbf{x}}_2}|⟨{\mathbf{x}}_1{\mathbf{x}}_2|{\mathbf{p}}_1{\mathbf{p}}_2⟩{|}^2,$$

$$\Delta Q_{11}=C\frac{Q_{10}^2}{V},$$

$$Q_{10}=\sum _{{\mathbf{p}}_1}{e}^{-\beta p_1^2/2m}.$$

$$Q_{22}=\sum _{\bar{\mathbf{P}}}{e}^{-\beta {\bar{\mathbf{P}}}^2/2m}⟨\bar{\mathbf{P}}|e^{-\beta \hat{V}}|\bar{\mathbf{P}}⟩.$$

$$e^{-\beta \hat{V}}|\bar{\mathbf{X}}⟩$$

$$f_1(\bar{\mathbf{X}})={\delta }_{{\mathbf{x}}_1,{\mathbf{x}}_3}+{\delta }_{{\mathbf{x}}_1,{\mathbf{x}}_4}+{\delta }_{{\mathbf{x}}_2,{\mathbf{x}}_3}+{\delta }_{{\mathbf{x}}_2,{\mathbf{x}}_4},$$

$$f_2(\bar{\mathbf{X}})=2[{\delta }_{{\mathbf{x}}_1,{\mathbf{x}}_3}{\delta }_{{\mathbf{x}}_2,{\mathbf{x}}_4}+{\delta }_{{\mathbf{x}}_1,{\mathbf{x}}_4}{\delta }_{{\mathbf{x}}_2,{\mathbf{x}}_3}].$$

$$\Delta Q_{22}=\sum _{\bar{\mathbf{P}},\bar{\mathbf{X}}}{e}^{-\beta {\bar{\mathbf{P}}}^2/2m}|⟨\bar{\mathbf{X}}|\bar{\mathbf{P}}⟩{|}^2(C{f}_1(\bar{\mathbf{X}})+C^2{f}_2(\bar{\mathbf{X}})).$$

$$|⟨{\mathbf{y}}_1{\mathbf{y}}_2|{\mathbf{q}}_1{\mathbf{q}}_2⟩{|}^2=\frac{1}{V^2}(1-\cos (({\mathbf{y}}_1-{\mathbf{y}}_2)({\mathbf{q}}_1-{\mathbf{q}}_2))).$$

$$\Delta Q_{MN}=\sum _{\bar{\mathbf{P}},\bar{\mathbf{X}}}{e}^{-\beta {\bar{\mathbf{P}}}^2/2m}|⟨\bar{\mathbf{X}}|\bar{\mathbf{P}}⟩{|}^2(C{f}_a(\bar{\mathbf{X}})+C^2{f}_b(\bar{\mathbf{X}})+\dots ),$$

$$\Delta Q_{M1}=MC\frac{Q_{10}}{V}\sum _{\bar{\mathbf{P}},\bar{\mathbf{X}}}{e}^{-\beta {\bar{\mathbf{P}}}^2/2m}|⟨\bar{\mathbf{X}}|\bar{\mathbf{P}}⟩{|}^2,$$

$$\mathcal{Z}$$

$$⟨{\hat{n}}_{\uparrow }{\hat{n}}_{\downarrow }{⟩}_0=\frac{1}{V^2}\sum _{\mathbf{p},\mathbf{q}}\frac{z_{\uparrow }{e}^{-\beta p^2/2m}}{1+z_{\uparrow }{e}^{-\beta p^2/2m}}\frac{z_{\downarrow }{e}^{-\beta q^2/2m}}{1+z_{\downarrow }{e}^{-\beta q^2/2m}}.$$

Full text

I Introduction

In a recent paper ShillDrut , two of us presented results for virial coefficients in a semiclassical lattice approximation (SCLA), at leading order (LO), for spin-$1/2$ fermions with a contact two-body interaction. We found that, in spite of the crudeness of the approximation, the results for $\Delta b_{3}$ and $\Delta b_{4}$ were surprisingly good when written in terms of the exact $\Delta b_{2}$, which amounted to using the latter as a renormalized coupling. Specifically, quantitative or at least qualitative agreement was found between the LO-SCLA and diagrammatic and Monte Carlo results for 1D and 2D Fermi gases with attractive interactions.

In this work, we explore the LO-SCLA further by carrying out the evaluation of virial coefficients up to $b_{7}$, and furthermore extending our previous analysis to next-to-leading order (NLO). Our (approximate) analytic answers, obtained partially by algebra automation, provide insight into the behavior of the virial expansion as a function of the spatial dimension of the problem. As an application, we compare with Monte Carlo results for the density equation of state of an attractive 2D Fermi gas. As an additional example, we calculate the many-body contribution to Tan’s contact for the same system.

II Hamiltonian and virial expansion

We assume a non-relativistic kinetic energy and a two-body contact interaction, such that the Hamiltonian for two flavors $\uparrow,\downarrow$ is $\hat{H}=\hat{T}+\hat{V}$, where

[TABLE]

and

[TABLE]

where the field operators $\hat{\psi}_{s},\hat{\psi}^{\dagger}_{s}$ are fermionic fields for particles of spin $s=\uparrow,\downarrow$ (summed over $s$ above), and $\hat{n}_{s}({\bf x})$ are the coordinate-space densities. In the remainder of this work, we will take $\hbar=k_{\text{B}}=m=1$ and discretize spacetime using the spatial lattice spacing $\ell$ to set the scale for all quantities. In particular, we will define the lattice kinetic energy exactly as above by using a momentum-space representation with periodic boundary conditions (rather than a local three-point formula for the second derivative in coordinate space), and the lattice potential energy will take the form

[TABLE]

where now all of the operators and constants on the right-hand side represent dimensionless lattice quantities, and we have omitted an overall prefactor $\ell^{-2}$ that gives $\hat{V}$ its physical units.

One way to characterize the thermodynamics of this system is through the virial expansion VirialReview , 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. The corresponding coefficients accompanying the powers of $z$ in the expansion of the grand-canonical potential $\Omega$ are the virial coeffiecients; specifically,

[TABLE]

where

[TABLE]

is the grand-canonical partition function and $Q_{N}$ is the $N$-body partition function. By definition, $b_{1}=1$ and the higher-order coefficients require solving the corresponding few-body problems:

[TABLE]

and so forth (see Appendix). For completeness and future reference, we note here the values of the noninteracting virial coefficients for nonrelativistic fermions in $d$ spatial dimensions: $b^{(0)}_{n}=(-1)^{n+1}n^{-(d+2)/2}$.

For our system, in arbitrary spatial dimensions,

[TABLE]

which in the continuum limit becomes $Q_{1}=2{V}/{\lambda_{T}^{d}}$, where $V$ is the $d$-dimensional spatial volume and $\lambda_{T}=\sqrt{2\pi\beta}$ is the de Broglie thermal wavelength. Since $Q_{1}\propto V$, the above expressions for $b_{n}$ display precisely how the volume dependence should cancel out to yield volume-independent coefficients. In particular, the highest power of $Q_{1}$ does not involve the interaction and therefore always disappears in the interaction-induced change $\Delta b_{n}$:

[TABLE]

and so forth. In the Appendix we show the corresponding expressions up to $b_{7}$, but the pattern is repeated: $\Delta b_{n}$ involves a contribution from $\Delta Q_{n}$ and several contributions involving the previous virial coefficients $b_{m}$, $m<n$ and powers of $Q_{1}$; the latter always cancel against specific terms within $\Delta Q_{n}$ to yield a volume-independent $b_{n}$. Those cancellations are a challenging feature for stochastic approaches, but they become a useful check for our calculations.

In terms of the partition functions $Q_{MN}$ of $M$ particles of one type and $N$ of the other type, we have

[TABLE]

We thus see that the number of non-trivial contributions to each virial coefficient is actually small. The challenge is in determining each of these terms and for that purpose we implement the semiclassical approximation advertised above, which we describe in detail next.

III The semiclassical approximation at leading order

III.1 Basic formalism

In a wide range of many-body methods, the grand-canonical partition function $\mathcal{Z}$ is expressed as a path integral over an auxiliary Hubbard-Stratonovich field. Here we use a different route, but with the same first step: we introduce a Trotter-Suzuki (TS) factorization of the Boltzmann weight. At the lowest non-trivial order in such a factorization,

[TABLE]

where the higher orders involve exponentials of nested commutators of $\hat{T}$ with $\hat{V}$. Thus, the LO in this expansion consists in setting $[\hat{T},\hat{V}]=0$, which becomes exact in the limit where either $\hat{T}$ or $\hat{V}$ can be ignored (i.e. respectively the strong- and weak-coupling limits). Orders beyond LO can be reached using a factorization based on the Trotter identity

[TABLE]

Indeed, the leading order can simply be viewed as the most coarse possible TS factorization, i.e. with time step $\tau=\beta$. Higher orders $n>1$ will be defined by using progressively finer discretizations $\tau=\beta/n$. We leave such explorations to future work.

III.2 A simple example

As the simplest example of the LO-SCLA, we calculate $Q_{11}$:

[TABLE]

The kinetic energy operator piece is thus trivially evaluated. The central step is to insert a coordinate-space completeness relation to evaluate the potential energy piece, which we do using the following identity:

[TABLE]

where $C=(e^{\beta g_{d}}-1)\ell^{d}$ and we used the fermionic relation $\hat{n}^{2}_{s}=\hat{n}_{s}$. The $C$-independent term yields the noninteracting result, such that we may write

[TABLE]

which simplifies dramatically in this particular case when using a plane wave basis, because $|\langle{\bf x}_{1}{\bf x}_{2}|{\bf p}_{1}{\bf p}_{2}\rangle|^{2}=1/V^{2}$. We then find

[TABLE]

where

[TABLE]

Thus, $\Delta b_{2}=C{Q_{10}^{2}}/{(VQ_{1})}=C{Q_{1}}/{(4V)}$, where we used $Q_{1}=2Q_{10}$. Following essentially the same steps, it is not difficult to see that $\Delta b_{3}=-CQ_{1}(2\beta)/V$, where $Q_{1}(2\beta)$ is $Q_{1}$ evaluated at $\beta\to 2\beta$.

III.3 A more difficult example

To display the complexity of the calculation in a less trivial case, we show $Q_{22}$ as another example. Using the notation ${\bar{\bf P}}=({\bf p}_{1},{\bf p}_{2},{\bf p}_{3},{\bf p}_{4})$, where $1,2$ refer to spin-up particles and $3,4$ to spin-down particles, we have

[TABLE]

As before, we must insert a complete set of coordinate eigenstates to evaluate the remaining matrix element. To that end, we note that, using the notation ${\bar{\bf X}}=({\bf x}_{1},{\bf x}_{2},{\bf x}_{3},{\bf x}_{4})$,

[TABLE]

where

[TABLE]

and

[TABLE]

Thus,

[TABLE]

Note that $\langle{\bar{\bf X}}|{\bar{\bf P}}\rangle$ factorizes across spins where, in this case, each of the factors involved takes the form

[TABLE]

Based on the above examples, it is easy to glean that the general form of the change $\Delta Q_{M,N}$ in the partition function for $M$ spin-up particles and $N$ spin-down particles, with a contact interaction, is given by

[TABLE]

where ${\bar{\bf P}},{\bar{\bf X}}$ represent all momenta and positions of the $M+N$ particles, and the functions $f_{a}$, $f_{b}$, $\dots$, which encode the matrix element of $e^{-\beta\hat{V}}$, depend on the specific case being considered. In particular, the case of $\Delta Q_{M1}$ is particularly simple and reduces to

[TABLE]

where ${\bar{\bf P}},{\bar{\bf X}}$ go over the momenta and positions of the $M$ identical particles.

The wavefunction $\langle{\bar{\bf X}}|{\bar{\bf P}}\rangle$ is a product of two Slater determinants which, if using a plane-wave single-particle basis, leads to simple Gaussian integrals over the momenta $\bar{\bf P}$. The only challenge is in carrying out the sum over $\bar{\bf X}$ before the sum over $\bar{\bf P}$, as the Slater determinants will naively lead to a large number of terms; to that end, it is crucial to use the interaction matrix elements to simplify the determinants before carrying out any momentum sums or integrals. In all cases, the integrals involved will be multidimensional Gaussian integrals.

III.4 Generating $\mathcal{O}(C)$ results for all $n$

As a check for our calculations and a useful result in itself, we show here how to calculate the contributions at $\mathcal{O}(C)$ in the LO-SCLA for all $n$. We begin with the grand-canonical partition function at LO, generalized to arbitrary chemical potentials $\mu_{\uparrow}$ and $\mu_{\downarrow}$ (sum over $s=\uparrow,\downarrow$ implied below):

[TABLE]

where $\mathcal{Z}_{0}=\tr\left[e^{-\beta(\hat{T}-\mu_{s}\hat{N}_{s})}\right]$ is the noninteracting partition function and $\langle\cdot\rangle_{0}$ denotes a noninteracting thermal expectation value. This is, of course, nothing other than leading-order perturbation theory. In the general, asymmetric case we have

[TABLE]

By differentiation of the above expression with respect to $z_{\uparrow}=e^{\beta\mu_{\uparrow}}$ and $z_{\downarrow}=e^{\beta\mu_{\downarrow}}$, it is straightforward to derive the $\mathcal{O}(C)$ change in $\Delta Q_{M,N}$ for arbitrary $M,N$.

IV Results

IV.1 Virial coefficients on the lattice

Using the above formalism, we obtained expressions for the virial coefficients on the lattice as shown in Tables 1, 2, and 3, in the LO-SCLA. The results shown in those tables, up to $\Delta b_{5}$, were obtained on paper; the remaining answers were obtained using an automated computer algebra code of our own design.

As mentioned in the Introduction, the explicit calculation of $\Delta b_{n}$ by way of $\Delta Q_{MN}$ involves delicate cancellations of various volume-dependent contributions which scale as $V$, $V^{2}$, $\dots,V^{n-1}$, which yields a volume-independent result for $\Delta b_{n}$. As a check for our automated calculations, we have verified that those cancellations take place exactly as expected.

The functions $Q$, $F$, and $G$ appearing in 1, 2, and 3 are defined by

[TABLE]

[TABLE]

[TABLE]

where ${\bf P}^{T}=({\bf p}_{1},{\bf p}_{2},{\bf p}_{3},\cdots)$ is a vector collecting all the $d$-dimensional momentum variables appearing in the sums, $\bf M$ and $\bf N$ are block matrices of the form

[TABLE]

for $d=3$, where the explicit form of the block entries $X_{i}$ = ${M}$ and ${N}$, respectively of size $3\times 3$ and $4\times 4$, are shown in the Appendix. In the cases studied here, the matrices $X_{i}$ will be the same for all cartesian components.

In the continuum limit, where the above sums turn into integrals, we have

[TABLE]

[TABLE]

and

[TABLE]

Using these formulas, we present our continuum results in the next section.

Note that, as a feature of the LO-SCLA, the expressions for $\Delta b_{2}$ and $\Delta b_{3}$ stop at $\mathcal{O}(C)$; the results for $\Delta b_{4}$ and $\Delta b_{5}$ terminate at $\mathcal{O}(C^{2})$; and finally $\Delta b_{6}$ and $\Delta b_{7}$ go only up to $\mathcal{O}(C^{3})$. We emphasize that those are full results within the LO-SCLA rather than approximations in powers of $C$.

IV.2 Virial coefficients in the continuum limit

In Ref. ShillDrut , it was shown that the LO-SCLA gives

[TABLE]

for a fermionic two-species system with a contact interaction, in $d$ spatial dimensions. [Note that we have corrected the coefficient of $(\Delta b_{2})^{2}$ relative to Ref. ShillDrut .]

The main result of this work is the extension of the above formulas to up to seventh order in the virial expansion and up to NLO in the SCLA. We collect the LO results in Table 4 and provide the NLO answers as Supplemental Material SuppMatt . In the LO-SCLA, $\Delta b_{2}\propto C$, and so $\Delta b_{2}$ tracks the order of $C$ appearing in each virial coefficient. At the same level in the approximation, $\Delta b_{3}$ only displays contributions up to $\mathcal{O}(C)$, while $\Delta b_{4}$ and $\Delta b_{5}$ stop at $\mathcal{O}(C^{2})$; $\Delta b_{6}$ and $\Delta b_{7}$ contain terms up to $\mathcal{O}(C^{3})$.

The results in Table 4 correspond to taking the continuum limit of the lattice expressions and using $\Delta b_{2}$ as the renormalized coupling to replace the bare lattice coupling $C$. Using those results, we show in Fig. 1 a plot of $\Delta b_{3},\dots,\Delta b_{7}$ as functions of the spatial dimension $d$ at $\Delta b_{2}=1/\sqrt{2}$, which is the value corresponding to the 3D unitary Fermi gas. We compare those answers with the $\Delta b_{3}$ results of Ref. ShillDrut in 1D, the diagrammatic results of Ref. Ngampruetikorn in 2D (see also Ref. DrummondVirial2D ), and the known results in 3D (Ref. Leyronas calculated the exact answer, while the work of Ref. LiuHuDrummond calculated it numerically, and Ref. DBK semi-analytically). We also compare our results for $\Delta b_{4}$ at unitarity with Ref. YanBlume (see also Refs. Rakshit and Ngampruetikorn2 ). We note that, while the LO-SCLA is quite rudimentary, the answers it provides are qualitatively correct as a function of $d$ but can be far from the expected numbers (in the sense and scale of Fig. 1). At NLO, on the other hand, the agreement improves considerably for $\Delta b_{3}$ but again deteriorates for $\Delta b_{4}$ (when comparing, in the latter, case, with the only data point available, which is at $d=3$).

While the change in the progression from LO to NLO in Fig. 1 is substantial, more weakly coupled regimes than unitarity feature much improved behavior. As an example, we show a plot of $\Delta b_{3},\dots,\Delta b_{7}$ as functions of the spatial dimension $d$ at $\Delta b_{2}=1/(5\sqrt{2})$ in Fig. 2, which corresponds to the BCS side of the 3D resonant Fermi gas. In all cases in that figure, the changes when going from LO to NLO are small on the overall scale of the plot, across all dimensions, whereas the $\Delta b_{n}$ are of the same order as the noninteracting values $b^{(0)}_{n}=(-1)^{n+1}n^{-(d+2)/2}$. This shows that the SCLA is able to capture the behavior of virial coefficients even in regimes where the interaction effects are not small.

To give a more precise sense of the behavior described above, we show in Fig. 3 a plot of the LO-NLO change in $\Delta b_{3}$ and $\Delta b_{4}$, namely $|\Delta b_{3}^{\text{NLO}}$-$\Delta b_{3}^{\text{LO}}|$ and $|\Delta b_{4}^{\text{NLO}}$-$\Delta b_{4}^{\text{LO}}|$, as percentages relative to the LO result, versus $\Delta b_{2}/\Delta b_{2}^{\mathrm{UFG}}$, where $\Delta b_{2}^{\mathrm{UFG}}=1/\sqrt{2}$ corresponds to the unitary limit of the 3D Fermi gas. Based purely on this LO-NLO analysis, we conclude that, in all cases, the convergence properties of the SCLA deteriorate both as a function of the coupling and as a function of the virial order. In other words, achieving a desired convergence level across a set of $\Delta b_{n}$ would require higher SCLA orders for higher $n$. This is not unexpected; in fact, it would be surprising to find the opposite behavior (i.e. improved convergence for higher $n$). Interestingly, lower dimensions display better convergence properties than the higher dimensional counterparts across all couplings, which is unexpected given that interaction effects are typically enhanced in low dimensions.

IV.3 Application: pressure and Tan’s contact in 2D

In this section we apply our estimates of the virial coefficients to two simple thermodynamic observables: the pressure and Tan’s contact. For concreteness, we focus on the 2D attractive Fermi gas, but we emphasize that our results are explicit analytic functions of the dimension and can therefore be evaluated for arbitrary $d$.

To access the pressure, we combine the calculated virial coefficients according to

[TABLE]

where $\Delta\Omega=-\Delta PV$ is the change in the pressure $P$ due to interaction effects. From this equation, it is straightforward to determine the density change $\Delta n$ and the compressibility change $\Delta\chi$ by differentiation with respect to $z$. In Fig. 4 we show the pressure, in units of its noninteracting counterpart $P_{0}$, as a function of $\beta\mu=\ln z$, for the 2D Fermi gas with attractive interactions, for the LO- (top) and NLO-SCLA (bottom). The results at NLO show somewhat improved agreement with previous results at all couplings studied when considering the full expressions that include $\Delta b_{7}$.

Within the context of the LO- and NLO-SCLA results for $\Delta b_{n}$, Fig. 4 provides an indication of the range of validity of the virial expansion. In those plots, the large oscillations observed as $\beta\mu$ is increased show that the radius of convergence $r$ of the virial expansion, as a function of $z$, is notably reduced as the coupling is increased. For $\beta\epsilon_{B}=1$, we have $r\simeq\exp(-1.25)$, whereas for $\beta\epsilon_{B}=3$, we have $r\simeq\exp(-2.75)$. Such an effect is expected but, to understand to what extent that reduction is due to the LO approximation, higher orders in the approximation must be investigated.

To obtain Tan’s contact Tan , we differentiate with respect to the coupling $\lambda$ (i.e. we use the so-called adiabatic relation). In this case, since we use $\Delta b_{2}$ as our physical dimensionless coupling (or as a proxy to scattering parameters or binding energies), we simply differentiate $\Delta\Omega$ with respect to that parameter and obtain the dimensionless form

[TABLE]

To connect the above expression to the conventional form of Tan’s contact we only need the overall factor $\partial b_{2}/\partial\lambda$, where $\lambda$ is the coupling. Such a factor contains only two-body physics and can therefore be calculated explicitly using the well-known Beth-Uhlenbeck formula BU . In Fig. 5 we show our results for ${\Delta\mathcal{C}}/{Q_{1}}$ as a function of $\beta\mu=\ln z$, for the 2D Fermi gas with attractive interactions. Although somewhat more difficult to visualize, we glean from this figure that the virial expansion breaks down at lower values of $\beta\mu$ at stronger couplings.

V Summary and Conclusions

In this work, we calculated the virial coefficients of spin-$1/2$ fermions in the LO-SCLA, up to $\Delta b_{7}$. We have presented analytic results on the lattice and, much more succinctly, in the continuum limit, where they feature an explicit analytic dependence on the number of spatial dimensions $d$.

As a renormalization prescription, we fixed the bare constant $C$ by using the fact that $\Delta b_{2}$ is in many cases known analytically through the Beth-Uhlenbeck formula BU (see e.g. Refs. EoS1D ; virial2D ; virial2D2 ; PhysRevA.89.013614 ; Daza2D ; LeeSchaeferPRC1 ). That choice allowed us to express our results in powers of $\Delta b_{2}$ and to perform cross-dimensional comparisons by varying $d$ at fixed $\Delta b_{2}$. In turn, that comparison shows that the SCLA behaves qualitatively as expected. Notably, the agreement for $\Delta b_{3}$ improves across dimensions when going from LO to NLO, but it appears to deteriorate for $\Delta b_{4}$ at unitarity. To better understand that feature, we showed results at weaker couplings, which show much better convergence (i.e. much smaller changes) when going from LO to NLO for all the $\Delta b_{n}$ studied. These results are encouraging towards exploring higher orders in the SCLA, which will be carried out elsewhere HouDrut .

The fact that we have used the continuum limit of momentum sums in our final answers, together with the continuum $\Delta b_{2}$ result to fix the coupling, does not eliminate all the lattice artifacts. Indeed, implicit in our derivations is the spatial lattice spacing $\ell$, with which the coupling runs and which induces finite-range effects. To avoid such effects, future studies should use improved actions (see e.g. DrutNicholson ; Drut ).

Finally, it should be stressed that the applicability of the SCLA goes beyond the approximation of virial coefficients. In the form implemented here, it can be used to access the thermodynamics of finite systems. For those, the SCLA simply represents the non-perturbative analytic evaluation of a finite-temperature lattice calculation which would not be possible due to the sign problem, and which is carried out in a coarse temporal lattice.

Acknowledgements.

This material is based upon work supported by the National Science Foundation under Grant No. PHY1452635 (Computational Physics Program).

Appendix A Matrices

In this section we provide the detailed form of the matrices that appear in the evaluation of the functions in Tables 1, 2, and 3. To that end, it is useful to define the following notation:

[TABLE]

Using the above, we have

[TABLE]

[TABLE]

[TABLE]

Appendix B High-order virial expansion formulas

For completeness, we provide here some of the formulas which we omitted in the main text for the sake of brevity and clarity. These are model independent, except as noted below. The complete expressions for $b_{5}$, $b_{6}$, and $b_{7}$ in terms of the corresponding canonical partition functions and prior virial coefficients can be written as

[TABLE]

whereas the change in the above due to interactions (assuming here two-body interactions) are given by

[TABLE]

To use these, it is useful to have the following:

[TABLE]

[TABLE]