Thermodynamics of spin-$1/2$ fermions on coarse temporal lattices
using automated algebra

We consider a non-relativistic, spin-$1/2$ fermionic system with Hamiltonian $\hat{H}=\hat{T}+\hat{V}$, where

is the kinetic energy in terms of spatial dimension $d$, particle mass $m$, and field operators $\hat{\psi}_{s}^{\dagger}$, $\hat{\psi}_{s}$ for spin $s$, and

is a contact interaction with coupling constant $g$ and particle density operators $\hat{n}_{s}=\hat{\psi}_{s}^{\dagger}\hat{\psi}_{s}$. To proceed, we let $\hbar=k_{B}=m=1$ and regularize the interaction by discretizing spacetime using a spatial lattice spacing $\ell$ such that the potential energy takes the form

where the operators are now dimensionless and $g_{lat}=g/\ell^{d}$.

To study thermodynamics, the central quantity of interest is the grand-canonical partition function of the unpolarized, interacting system,

where $\beta$ is the inverse temperature, $\mu$ is the chemical potential (common to both spin degrees of freedom, as we focus on unpolarized matter here), and $\hat{N}$ is the total particle number.

The central challenging point of quantum thermodynamics is that $\hat{T}$ and $\hat{V}$ do not commute. To address that problem, we follow a common route and introduce a temporal discretization $\beta=\tau N_{\tau}$ where $N_{\tau}$ is a positive integer and $\tau$ is the spacing, which allows us to separate the kinetic and potential energy exponential factors via a (symmetric) Trotter-Suzuki factorization, namely

where $\hat{K}_{0}=\hat{T}-\mu\hat{N}$. The net result on the partition function is the following approximation (valid up to $O(\tau^{2})$ once all $N_{\tau}$ time slices are accounted for):

where there are $N_{\tau}$ factors. Since $\tau=\beta/N_{\tau}$, it follows that increasing $N_{\tau}$ improves the accuracy of the final answer while also increasing the computational cost.

On a spatial lattice, we may use the property $\hat{n}^{2}=\hat{n}$, valid for dimensionless, lattice fermion density operators, to write

for $k\geq 0$ where we have introduced an arbitrary parameter $\lambda$ (see below) and

with $C=(e^{\tau g}-1)/\lambda$, $\hat{n}_{2}({\bf x})=\hat{n}_{\uparrow}({\bf x})\hat{n}_{\downarrow}({\bf x})$, and $\hat{V}_{0}=\openone$. When envisioning taking the large-$N_{\tau}$ limit at constant $\beta$, it is useful to set $\lambda=1/N_{\tau}$, such that, in that limit, one has $C\to\beta g$, i.e. $C$ approaches a well-defined limit in terms of the parameters of the theory. In practice, $C$ is set by renormalization (more on this below), such that the choice of $\lambda$ should be immaterial (up to, possibly, numerical artifacts). We will also see below that choosing $\lambda=1/N_{\tau}$ makes sense from the point of view of the scaling of the number of terms in the cumulant expansion presented below as $N_{\tau}$ is increased.

We make use of the interaction expansion of Eq. (7) by expanding the partition function into a quantum moment expansion:

where $\mathcal{Z}_{0}$ is the partition function of the non-interacting system and

are the quantum moments. Based on the above moment expansion, the following cumulant expansion results, which is one of the cornerstones of the QTCE, as described below:

where $\kappa_{0,0,\dots,0}=0$ for any value of $N_{\tau}$, because $W_{0,0,\dots,0}=1$ for all $N_{\tau}$. It is this expansion that gives us access to thermodynamic observables via the fundamental relation

where ($P_{0}$) $P$ is the pressure of the (non-)interacting system and $V$ is the $d$-dimensional volume.

The series form of Eq. (12) naturally allows one to study the effect of including successively higher-order terms, when the series is truncated at some order $N_{\kappa}$. Thus, our method has two parameters, namely $N_{\tau}$ and $N_{\kappa}$, which can be systematically increased. In this work, we favor increasing $N_{\kappa}$ while keeping $N_{\tau}$ relatively small (thus the “coarse temporal lattices” advertised in the title). Note that taking these limits in reverse order amounts to a perturbative approach, while increasing $N_{\kappa}$ first, at constant $N_{\tau}$ (as we propose), is formally closer to lattice Monte Carlo calculations.

Using

where $\lambda_{T}=\sqrt{2\pi\beta}$ is the thermal wavelength, and

we can reformulate Eq. (12) as

where

and we define the “assembled cumulants” as

and so on.

Note that since $\bar{\kappa}_{1,0,\dots,0}=\bar{\kappa}_{0,1,0,\dots,0}=\dots=\bar{\kappa}_{0,0,\dots,1}$, we can write $\bar{K}_{1}=N_{\tau}\bar{\kappa}_{1,0,\dots,0}$, which explicitly shows that $\bar{K}_{1}$ is linear in $N_{\tau}$, and this dependence is then cancelled in the expansion by the factor $\lambda=1/N_{\tau}$ in Eq. (17). Similarly, it is easy to see that there are $\propto N^{2}_{\tau}$ terms in $\bar{K}_{2}$, whose number is balanced by the $\lambda^{2}=1/N^{2}_{\tau}$ dependence in Eq. (17). The same reasoning applies to all orders in the expansion of Eq. (17).

These $\bar{K}_{n}$ quantities are useful for analyzing the system and are the main output of QTCE, which returns results up to a given cumulant order $N_{\kappa}$, corresponding to the upper limit on the summation in Eq. (16). We expect the $\kappa_{n_{1},n_{2},\dots,n_{N_{\tau}}}$ to be extensive quantities proportional to $\left({L}/{\lambda_{T}}\right)^{d}$ and therefore all $\bar{\kappa}_{n_{1},n_{2},\dots,n_{N_{\tau}}}$ to be intensive.

The cumulants can generally be expressed in terms of the moments via recursive formulas, such that to calculate a cumulant at a given order, one needs all of the moments up to and including that order [12]. For example, for $N_{\tau}=1$,

where $C^{n}_{k}=n!/(k!(n-k)!)$ is the binomial coefficient.

Since Eq. (13) establishes that $\ln(\mathcal{Z}/\mathcal{Z}_{0})$ is proportional to the volume $V$, and $V$ and $\lambda$ are independent variables, it must be that within each cumulant $\kappa_{n}$ the only terms that contribute are those with the proper spatial volume scaling $\propto V$. Since this is true for all $n$, we see that in Eq. (21) the role of the terms in the sum is to cancel out the unnecessary terms in $W_{n}$. In practice, we therefore proceed simply by calculating $W_{n}$ and retaining only the terms with linear $V$ scaling. In terms of Feynman diagrams, this amounts to keeping only the connected diagrams, as indicated by the linked-cluster theorem (see e.g. Ref. [13]). Thus, we are able to eliminate a large number of terms early on in the calculation, thereby reducing the computational cost.

In Eq. (II.2) we saw that in order to calculate the moment $W_{n_{1},n_{2},\dots,n_{N_{\tau}}}$ we need

where we define

to represent the spatial integration that is part of the simplification and integration modules of QTCE.

Next, we turn to calculating the noninteracting expectation value

using a noninteracting generating functional and automated algebra.

The noninteracting generating functional of density correlation functions is defined as

where we couple our source $j({\bf x},t)$ to the density at each spacetime point via

which includes an imaginary-time coordinate $t$ to identify the various insertions of $e^{\hat{V}_{\text{ext}}}$ in $\mathcal{Z}_{0}[j]$. Furthermore, we have assumed that our system is placed on a lattice as a way to regularize a future interaction, such that all the quantities below, unless otherwise stated, are lattice quantities.

By taking functional derivatives of $\mathcal{Z}_{0}[j]$ with respect to $j$ at a given point in space and time, and then dividing by $\mathcal{Z}_{0}$, we can generate noninteracting expectation values of any combination of density operators. For example,

where the $1$ in $({\bf x},1)$ indicates the first time slice. By taking derivatives at different time slices, the above expression can be generalized to generate expectation values of the form in Eq. (24).

Because all the operators involved are now exponentials of one-body operators, the following essential simplification holds:

where $z=e^{\beta\mu}$ is the fugacity and $U=U_{1}U_{2}U_{3}\dots U_{N_{\tau}}$ is the product of the one-particle representation of the transfer matrices

where $|{\bf p}\rangle$ are single-particle states and $t$ corresponds to the time slice.

It is not difficult to show (in fact, this is a standard result in many-body physics; see e.g. Refs. [13, 1]) that there is another useful representation for the generating functional, namely

where

We will refer to Eqs. (28) and (30) as the $U$ and $M$ representations, respectively, but they are often also referred to as the spatial and spacetime representations.

To simplify our notation below, we define

which is a matrix as $T$ is the single-particle matrix representation of $\hat{T}$, i.e.

We also define

which is also a matrix, where $x_{k}=({\bf x}_{k},t_{k})$, and similarly

It is very easy to see from the above form of $M$ that any derivative of $M$ beyond the first one vanishes when evaluated at unequal times, i.e.

for $t_{1}\neq t_{2}$ and arbitrary $j$.

Similarly, it is not difficult to see that, for $t_{1}=t_{2}$,

for ${\bf x}_{1}\neq{\bf x}_{2}$ and arbitrary $j$, because both $\hat{V}_{\text{ext},t}$ and its exponential are diagonal matrices in coordinate space.

As usual when working with generating functionals, we set $j=0$ at the end of the calculation. With these definitions one finds, for example,

and

where the last term vanishes if $t_{1}\neq t_{2}$. In the $U$ representation, the above becomes

and so forth, where the traces are computed in momentum space for the homogeneous system. [In the above expressions, the time dependence appearing in the operators $\hat{n}$ simply indicates the time slice where they are to be inserted, in expressions such as Eq. (24).]

The case of equal times $t_{1}=t_{2}$ requires some care for general interactions, as it may generate terms that diverge in the continuum limit [namely the third term in the right-hand side of Eqs. (39) and (40)]. For a contact interaction, as is our interest here, those terms are simply not present, by virtue of Eq. (9), i.e. the required equal-time correlations are only needed at unequal spatial separation.

Thus, we see that to calculate the moments (and thus access the cumulants), one can begin by generating trace expressions for density correlations that simply involve the propagator of Eq. (32) and the first derivative of $U$. This is one of the first steps of the automated algebra code implemented by the QTCE. Note that at this stage, we have focused on two-body local interactions (i.e. of the density-density form); generalizations of the method are discussed in Section VI.

The QTCE uses the above formalism and computational approach to evaluate the quantum cumulants of Eq. (12) following the steps outlined below.

1. 1.

   Expectation value identification: Determine which expectation values need to be calculated. These correspond to Eq. (24) but at this stage are simply identified by the set of indices $\{n_{i}\}$.

2. 2.

   Term generation: Generate expressions for the required non-interacting expectation values using functional derivatives. At the end of this step, we have expressions as in Eq. (40).

3. 3.

   Simplification: Carry out coordinate sums of Eq. (23) as well as other simplifications before assembling the cumulants. The outputs of this step are multidimensional momentum integrals, each corresponding to a diagram. They have the general form

   $$\int{\mathcal{D}}{\bf P}\prod_{k=1}^{N}G[{\bf P},{\bf Q}_{k},z],$$

   (41)

   where (for $N_{\tau}=1)$

   $$G[{\bf P},{\bf Q}_{k},z]=\frac{ze^{-\frac{1}{2}|{\bf Q}_{k}^{T}\cdot{\bf P}|^{2}}}{1+ze^{-\frac{1}{2}|{\bf Q}_{k}^{T}\cdot{\bf P}|^{2}}},$$

   (42)

   ${\bf P}^{T}=({\bf p}_{1},{\bf p}_{2},...,{\bf p}_{K})$ is a vector made out of $d$-dimensional vectors, ${\bf Q}_{k}$ is a $K$-dimensional vector, and

   $${\mathcal{D}}{\bf P}\equiv d{\bf p}_{1}d{\bf p}_{2}...d{\bf p}_{K}.$$

   (43)

4. 4.

   Evaluation: Evaluate the diagrams, i.e. the multidimensional integrals obtained in the previous step. This can be done with a variety of methods including a pure numerical evaluation, such as the VEGAS algorithm [14], or a semi-analytic approach, such as an expansion in powers of $z$, i.e. the virial expansion. Another semi-analytic approach which we have found to be useful involves an expansion in powers of $\frac{z}{1+z}$ (which is always less than 1 for all $z>0$), as further discussed below.

5. 5.

   Assembly: Once the diagrams are evaluated, their contribution to each assembled cumulant $\bar{K}_{n}$ (and therefore to the pressure) are finally assembled to obtain physical results.

When considering integrals of products of $G(z)$ functions, the main object of interest is the function

where $x^{2}$ represents the dispersion relation of noninteracting, nonrelativistic matter.

In the approach we advocate here, we aim to generate Gaussian integrals only. To that end, one may consider expanding $F(x,z)$ in powers of $z$, which effectively yields the virial expansion:

As the $e^{-\frac{k}{2}x^{2}}$ factor is always positive and less than or equal to $1$, this expansion is restricted to $z<1$, as expected and is well-known. Since our interest is in arbitrary $z$, in particular larger than $1$, we put aside this expansion and proceed in a different direction.

To capture the general $z$ behavior, we note that regardless of the value of $z$, at large $x$ the leading term is

while at small $x$,

To capture these limits, we add and subtract $z$ in the denominator to obtain

which we then expand as a series in powers of $z/(1+z)$, such that

As desired, this is a collection of Gaussians, as the virial expansion proposed earlier, but in a different form. Moreover, the above expression is within the radius of convergence of the geometric series for all values of $z$ and $x$ of interest. Since the binomial theorem dictates that

we may write the final form

When carrying out integrals involving products of $G$ functions, one may expand such products in powers of $\frac{z}{1+z}$ using the above expression, such that the resulting integrals are always sums of Gaussian integrals. We will refer to such an approach to integration as the large-$z$ expansion (LZE), although it is not an expansion around $1/z=0$. One of the crucial advantages of the LZE over purely numerical integration techniques is that the latter are usually based on stochastic methods and therefore incur a statistical uncertainty. Furthermore, such methods are usually confined to integer spatial dimensions. The LZE, on the other hand, is systematic and can be used to evaluate integrals in arbitrary dimensions, even complex. We make use of this property below. To further display the appealing properties of the LZE, we show in Fig. 1 its behavior at second and fourth orders when compared with the virial expansion.

As a first application of the QTCE, we show in Fig. 2 the pressure equation of state of spin-$1/2$ fermions at unitarity, compared with prior data from experiments as well as selected numerical approaches (including self-consistent diagrammatic [15] and Monte Carlo methods [16]). Specifically, we show in this plot the variation in our results when renormalizing using the second-order virial expansion (VE, shown in dashed gold line) or its fifth-order counterpart (solid gold line); the corresponding QTCE results are shown as light- and dark-shaded blue bands, respectively. The limits of the QTCE bands are set by renormalizing at $\beta\mu=-3.0$ and $\beta\mu=-1.0$; we find that the fifth-order VE renormalization scheme yields a much more constrained band than the second-order scheme. Remarkably, our results match the MIT [17] experimental results up to at least $\beta\mu=1.0$ (with a slight undershooting deviation around $\beta\mu=0$). We emphasize that the QTCE answer is given by a plain (no resummation) partial sum of the cumulant expansion of Eq. (12) up to $N_{\kappa}=6$ with the most rudimentary approximation for the imaginary-time direction, namely $N_{\tau}=1$.

Extending our results beyond unitarity, in Fig. 2 we show the pressure as a function of $\beta\mu>0$, across the BCS-BEC crossover, for varying couplings parametrized by the ratio $\Delta b_{2}/\Delta b_{2}^{\text{UFG}}$ (shown with different colors). Here, $\Delta b_{2}$ is the interaction-induced change in the second-order virial coefficient, which is in one-to-one correspondence with the interaction strength as measured usually by the scattering length; $\Delta b_{2}^{\text{UFG}}=1/\sqrt{2}$ is the value in the unitary limit. In addition, the plot shows partial sums of the cumulant expansion up to order $N_{\kappa}$ (shown with different shades of the same color), indicating the highest power of $\lambda$ summed in the series of Eq. (12). All of the results shown in this plot correspond to $N_{\tau}=1$, i.e. the simplest possible Suzuki-Trotter factorization (i.e. the most coarse temporal lattice). For $\beta\mu<0$, there is no noticeable difference from varying $N_{\kappa}$, which means that the first cumulant, properly renormalized, is enough to capture the behavior of the curve.

The weakest coupling we considered is $\Delta b_{2}/\Delta b_{2}^{\text{UFG}}=0.01$, which is indistinguishable from the noninteracting limit in the scale shown in the plot. As the coupling is increased, the partial sums display oscillatory behavior at sufficiently large $\beta\mu$, indicating that the series must be resummed using methods such as Pade approximants. In particular, we note that the splitting of the partial sums coincides roughly (and not surprisingly) with the location of the superfluid phase transition (gray vertical band) around $\beta\mu=2.4$ in the unitary limit $\Delta b_{2}/\Delta b_{2}^{\text{UFG}}=1$ (blue lines).

To characterize the thermodynamics of a system in a more detailed fashion, here we take a closer look at derivatives of the pressure, which yield the density equation of state and the elementary response functions given by the compressibility and heat capacity.

The density equation of state is accessed from Eq. (13) using the thermodynamic relation

Taking further derivatives of the pressure, one may obtain the isothermal compressibility via the thermodynamic relation

Our results for the density in the unitary limit are shown in Fig. 4. It is easy to see that, as for the pressure, the QTCE results, even at $N_{\tau}=1$, are a dramatic improvement over the fifth-order VE. As in Fig. 3, it is easy to see that as the superfluid phase transition is approached from low fugacity, the cumulant expansion displays wild oscillations that must be resummed.

Finally, the heat capacity per particle at unitarity is given by

where the first equality is the definition and the second equality is valid only at unitarity, where one may use the pressure-energy relation $E=3PV/2$; furthermore, we used $P_{0,T=0}={2nE_{F}}/{5}$, where $E_{F}=(3\pi^{2}n)^{2/3}/2$ is the Fermi energy of a noninteracting gas at density $n$.

Our results for the compressibility and heat capacity are shown in Fig. 5. Both of these response functions show the correct qualitative features for temperatures above the superfluid transition. Moreover, the agreement with other approaches is arguably not merely qualitative. However, as with other quantities discussed above, dramatic oscillations appear in the cumulant expansion in low-temperature phase. This indicates that future analyses will require resummation approaches to interpret QTCE data when approaching a phase transition, which is not surprising.

One of the distinguishing aspects of QTCE is that it encodes the spatial dimension $d$ as an analytic variable that can be tuned continuously. To exemplify that capability, we show in Fig. 6 the results obtained for the pressure ratio $P/P_{0}$ as a function of $d$, for various chemical potentials, tuned to the unitary limit by renormalizing using the condition of matching at low fugacities with the fifth-order VE, as calculated in Ref. [8], across dimensions. We emphasize that we do not imply that this condition defines a unitary point in different dimensions; rather, it simply provides a way to formally renormalize the problem as $d$ is varied.

As seen in Fig. 6, the interaction effects on $P$ disappear as $d$ is decreased. Put differently, weaker interactions are more easily able to reproduce the condition $\Delta b_{2}=1/\sqrt{2}$ in low dimensions than in high dimensions. In a real system, the physical picture behind this property is that a two-body bound state appears in 1D and 2D with vanishingly small attraction, which can dramatically increase $\Delta b_{2}$ already at weak coupling, as $\Delta b_{2}$ contains a term that increases exponentially with the binding energy.

The above observation suggests that, when renormalizing using the VE, one may expect the cumulant expansion (and more generally the QTCE results in coarse temporal lattices) to present better convergence properties in low dimensions. To explore that intuition, we show in Fig. 7 (top) the pressure as a function of $\beta\mu$ at unitarity, obtained with QTCE by dimensional extrapolation using data from $d=1.125$ to $d=2.5$, compared with QTCE exactly at $d=3.0$, as well as with the experimental results from MIT [17]. The results obtained with QTCE via extrapolation from low dimensions are clearly in much better quantitative agreement with the experimental numbers than those obtained by calculating directly in 3D.

In the bottom panel of Fig. 7, we provide a more detailed look at the dimension dependence of the pressure, at fixed $\beta\mu=2.5$, for $N_{\tau}=1$ and $N_{\kappa}=3$. There is a clear quantitative difference in the final results between taking the QTCE results at face value at $d=3$ (regardless of the resummation strategy or integrator used, as shown) and extrapolating from low dimensions (in this case between $d=1.125$ and $d=2.25$.

For completeness, we show in Fig. 8 a comparison of QTCE results at $N_{\tau}=1$ with a third-order perturbation theory calculation (in 1D; Ref. [24]) and auxiliary field quantum Monte Carlo calculations (in 2D; Ref. [25]) for the pressure as a function of $\beta\mu$. In 1D, the perturbation theory results are used to renormalize the QTCE calculation at $\beta\mu=-3.0$. In 2D, on the other hand, the second order virial expansion is used to renormalize at $\beta\mu=-3.0$. Also in both cases, it is clear that the results of QTCE at $N_{\tau}=1$ are superior to those of the virial expansion when increasing $\beta\mu$, even in the simplest possible case of $N_{\kappa}=1$. At sufficiently large $\beta\mu$, however, the partial sums of the QTCE cumulant expansion break down; in particular, they do so earlier in $\beta\mu$ for stronger couplings, as measured by the parameter $\lambda=2\sqrt{\beta}/a_{0}$ in 1D (unrelated to the QTCE parameter $\lambda$) and $\beta\epsilon_{B}$ in 2D, where $a_{0}$ is the 1D scattering length and $\epsilon_{B}$ is the binding energy of the two-body system.