Toward an automated-algebra framework for high orders in the virial expansion of quantum matter

The virial expansion provides a non-perturbative view into the thermodynamics of quantum many-body systems in dilute regimes. While powerful, the expansion is challenging as calculating its coefficients at each order $n$ requires analyzing (if not solving) the quantum $n$-body problem. In this work, we present a comprehensive review of automated algebra methods, which we developed to calculate high-order virial coefficients. The methods are computational but non-stochastic, thus avoiding statistical effects; they are also for the most part analytic, not numerical, and amenable to massively parallel computer architectures. We show formalism and results for coefficients characterizing the thermodynamics (pressure, density, energy, static susceptibilities) of homogeneous and harmonically trapped systems, and explain how to generalize them to other observables such as the momentum distribution, Tan's contact, and the structure factor.


Introduction
Quantum many-body systems are notoriously difficult to compute. Whenever interactions play an important role, in atomic, condensed matter, or nuclear physics, most analytic approaches (if not all) are unable to give quantitatively reliable results and numerical methods often face challenges of their own as well, such as the infamous sign problem (notwithstanding remarkable progress over the last few years) [1].
The virial expansion (VE) (see e.g. Ref. [2] for an introduction to both the classical and quantum cases and Ref. [3] for a comprehensive review) aims to tackle the finite-temperature quantum many-body problem by breaking it down into contributions from subspaces of the full Fock space corresponding to fixed (and small) particle number. In this sense, the VE is effectively an expansion around a dilute limit in which the interparticle distance is much larger than every other scale in the system, in particular the thermal wavelength λ T = 2πβ, where β = 1/T is the inverse temperature and we have used units such thath = k B = m = 1. One might expect that, this being a high-temperature regime (as λ T is in the above sense small), interaction effects would play a quantitatively minor role; this, however, is not necessarily the case, as we will see. Another misconception is that quantum effects play a small role in such a regime, which is also not generally true. Both interaction and quantum effects are central in the calculation of the virial expansion and leave a clear imprint in the expansion coefficients. It is for that reason that such calculations are challenging and that specialized techniques are required to carry them out, as we explain in this review. But, we are getting ahead of ourselves; let us start from the beginning.
The origin of the VE is the classical cluster expansion of Mayer et al, developed in the 1930s [4], in which the classical grand-canonical partition function Z of a three-dimensional gas is expanded in powers of the fugacity z = exp(βµ), where µ is the chemical potential, namely where is the classical canonical partition function and V ij is the pairwise interaction potential between particle i and particle j. An expansion of both the pressure P and the number density n in powers of z is thus obtained, namely where the b k , originally called 'cluster coefficients', typically depend on temperature and the specific form of V ij . Substituting the density expansion into the pressure expansion, order by order in z, one obtains the old virial expansion of the imperfect gas equation of state in powers of the density, i.e.
where x = nλ 3 T and the a k are in older literature called 'virial coefficients' (more recently, that nomenclature has been used for the b k coefficients instead). Calculating the a k requires knowing all the b coefficients up to order k, which in turn requires calculating the canonical partition functions of up to k particles, as we will show in more detail below. The original work of Mayer et al proposed a diagrammatic technique for calculating the classical case, which is naturally simpler than its quantum counterpart simply because of the well known fact that kinetic and potential energies are commuting numbers in the classical case and noncommuting operators in the quantum case. The quantum virial expansion was first explored by Kahn and Uhlenbeck [5] and further developed by Lee and Yang [6][7][8][9][10]. As we review below, the calculation of b k for the quantum case is so challenging in practice that efforts to calculate the third-order coefficient and beyond were not successful until the 21st century, when computers became powerful enough to apply exact diagonalization techniques. From this point on, we focus entirely on the VE in the context of quantum systems.
Historically, applications of the VE have followed the above paradigm and thus centered on equations of state. However, since the VE is an expansion of the grand-canonical partition function, it is possible to apply it to any physical quantity. We will show techniques to calculate the VE for applications to pressure and density but also to the Tan contact [11][12][13] (relevant for systems with short-range interactions), the momentum distribution, and response functions such as the compressibility and the structure factor.
The recent developments of automated algebra, led by our group [14][15][16][17][18][19][20][21], have enabled the precise calculation of high-order coefficients (meaning beyond the third order, which can currently be addressed numerically). With such orders in hand, it becomes practical and meaningful to implement resummation techniques which, uncertainties notwithstanding, have been shown to substantially extend the domain of applicability of the virial expansion [18,19].
The remainder of this work is organized as follows. We begin in Sec. 2 by reviewing the formal elements of the virial expansion of the grand thermodynamic potential, as that is the simplest case, which also allows us to establish the basic identities and notation. In Sec. 3, we provide a brief review of conventional calculation methods for the virial coefficients. In Sec. 4, we present our method by example using a homogeneous system of Fermion gases with attractive interaction changing from non-interacting limit to the unitary limit. Encouraged by the agreements with and improvements over existing results, we further generalize the method to other systems in Sec. 5, including harmonically trapped systems in Sec. 5.1, neutron matter in Sec. 5.2, and the unitary Bose gas in Sec. 5.3. Finally, Sec. 6.2 demonstrates the applications to more complicated observables. Namely, one-body operators such as density or momentum distributions are discussed in Sec. 6.2, and two-body operators, which are of interest to quantities such as the structure factor or viscosity, are shown in Sec. 6.3.

Basic formalism
As mentioned above, the VE is an expansion in powers of the fugacity z where β is the inverse temperature and µ is the chemical potential. In the presence of spin or other internal degrees of freedom, there will naturally be a fugacity attached to each (conserved) particle number. In this section we present the formalism for a single flavor for simplicity, but in upcoming sections we will generalize it to spin-1/2 fermions.
The thermodynamics is encoded in the grand-canonical partition function, which for a quantum system is given by whereĤ is the Hamiltonian,N is the particle number operator, and µ is the chemical potential. Expanding Z in powers of the fugacity we obtain where is the canonical N-particle partition function. Then, the grand thermodynamic potential is expanded in powers of z using the above expressions to obtain the conventional expression where b k is the k-th order virial coefficient, which is an intensive, dimensionless quantity. Through Eq. (8) and Eq. (10), the b k are related to the canonical partition functions; for example, the first few b k are For a system of non-interacting, non-relativistic fermions in d dimensions, the virial coefficients are b whereas for bosons the result is simply k −(d+2)/2 . In the presence of interactions, one usually expands the ratio of Z to its noninteracting counterpart Z 0 , to obtain k captures the interaction effects and is the quantity most often reported. In practice, the ∆b k coefficients are determined by the interaction-induced change in Q j for 1 < j ≤ k, and therein lies the difficulty: those Q j must be determined with enough accuracy to cancel out all the volume dependence (which contains terms that scale with power up to j) and obtain a volume-independent ∆b k . As we will see below, some methods focus on extracting ∆b k directly from grand-canonical quantities (typically the density), where the volume cancellations have already happened, while others such as ours (and similarly exact diagonalization) propose to calculate the interaction effects on Q j and use those to calculate ∆b k .
The above is the VE as applied to the pressure; from it, the expansion for the density is easily derived. As mentioned above, the VE can in fact be applied to any observable such as Tan's contact, the momentum distribution, and the structure factor (see Sec. 6.3). We will return to those in a later section.

Second order
Besides the trivial one-body contribution, fully captured by Q 1 and factored out of the expansion, the leading contribution accounting for interaction effects appears at second order, i.e. the two-particle subspace. The interaction effects on the two-particle spectrum are captured by the scattering properties (binding energies and phase shifts) and in those terms the second order virial coefficient b 2 was first calculated analytically by Beth and Uhlenbeck in the 1930's [22,23]. Specifically, their result relates the interaction change ∆b 2 to the two-body scattering phase shift δ(E) such that, for a spin-1/2 Fermi gas in three spatial dimensions (3D), where the first summation is over all bound states and the second summation over all partial waves.
One may take the above expression in different dimensions and relate it to the parameters of the corresponding effective range expansion (i.e. scattering length, effective range, etc.) and obtain, for a zero-range interaction (see e.g. [24,25] for results in 1D, [26,27] for 2D, and [28] for 3D): where λ d is the physical coupling strength in d dimensions, defined as where a 0 is the s-wave scattering length and E B is the binding energy of the single two-body bound state of the 2D case. While the above results are sufficient for simple systems such as dilute gases of ultracold atoms, where the interaction can be modeled very precisely as being purely zero-range s-wave, a deeper analysis is needed to account for the complexities of nuclear systems (such as neutron and nuclear matter). References [29,30] presented such an extension of the 3D case to the richer scattering properties found in nuclear physics. There, one must account for not only finite range but also angular momentum channels beyond the simplest case of pure s-wave.
For pure neutron matter, Refs. [29,30] integrate by parts to rewrite the Beth-Uhlenbeck result as where δ tot neutrons (E) is the sum of all the scattering phase shifts at laboratory energy E, whose contributions from different angular momentum channels enter as δ tot where the partial wave terms δ2S+1 L J (E) are obtained from partial wave analyses of experimental data such as Nijmegen's [31]. For nuclear matter, on the other hand, one must account for the deuteron bound state, such that where E d is the binding energy of the deuteron and the −1 term comes from partial integration when accounting for the phase shift at zero energy being π times the number of bound states (see also Ref. [32]). The work of Refs. [29,30] also analyzed the contributions due to pure alphaparticle scattering and nucleon-alpha scattering, thus obtaining all possible contributions to the second-order virial expansion for nuclear matter composed of neutrons, protons, and alpha particles.

Third order and beyond
The complexity of the quantum many-body problem for three particles and beyond forces one to switch to a combination of analytic and numerical approaches to calculate virial coefficients beyond second order. In cases of great interest such as spin-1/2 fermions, one has to further break up the problem into the subspaces of fixed particle number; for example, calculating b 4 requires solving the problem of 3 + 1 particles (i.e. 3 spin up and 1 spin down, and viceversa) as well as 2 + 2 particles. The number of such subspaces naturally proliferates with higher total particle number, and furthermore each subspace may present its own difficulties.
To calculate ∆b 3 there is only one distinct subspace that matters, namely 2 + 1 particles (assuming that both particles have the same mass) and impressive exact analytic progress was made by several authors, notably the work of Leyronas [33], Kaplan and Sun [34], and Castin and colleagues [35][36][37][38], as well as the large effective range expansion of Ngampruetikorn et al [39] (see also the early work of Ref. [40] focusing on the unitary limit).
Leyronas organizes the calculation of ∆b k around the VE of the density equation of state, itself expressed as an integral over all momenta of the equal-time single-particle Green's function in momentum space (i.e. the momentum distribution). Diagrams of various types are then identified at each order in z (up to order z 3 ), where contributions from the 2-and 3-body T matrix appear (the latter describing the atom-dimer scattering). The resulting time integrals are converted into energy integrals, which are then evaluated analytically where possible and otherwise numerically. The resulting approach is thus for the most part analytic and in principle exact and is in remarkable agreement in the unitary limit with prior purely numerical results for b 3 [41].
Other diagrammatic approaches also made interesting contributions. The work of Kaplan and Sun [34], which preceded Leyronas, starts from the density equation of state written as a momentum integral over the single-particle Green's function (as Leyronas does), but rather than carrying out the Matsubara sum from the outset, it uses a Poisson summation to express the propagator directly as a power series in z. The latter is then interpreted as a sum over winding number of worldlines around the compact imaginary time direction. The diagrams associated with each term in that expansion are referred to as 'chronographs'. Adding the contributions from such chronographs and accounting for systematic effects by extrapolation, very good agreement with prior numerical results [41] for b 3 was obtained in the unitary limit.
Ngampruetikorn et al [39], used an expansion around large effective range R * (compared to the thermal wavelength λ T ), which allows them to examine up to the four-particle subspace diagrammatically, thus obtaining numerical estimates for up to ∆b 4 . They focused on the unitary Fermi gas by interpolating between R * λ T and λ T /|a| 1, where a is the scattering length, and applied their method to the pressure, density, entropy, and spectral functions. Their interpolation results for ∆b 3 and ∆b 4 at unitarity agree with those obtained by other groups, including those presented here. In Ref. [42], Ngampruetikorn et al also studied the pairing correlations of the 2D Fermi gas up to third order in the virial expansion, additionally obtaining Tan's contact (we return to the expansion of this and other quantities below).
In Ref. [43] Werner and Castin analyzed the (2+1)-body problem of harmonically trapped spin-1/2 fermions at unitarity, obtaining their exact spectrum and eigenstates. Generalizing that work to the problem of 3+1 and 2+2 particles, Endo and Castin [36,37] (see also [38]) calculated the value of ∆b 4 as a function of the trapping frequency βω. As we will show in Sec. 5.1, our non-perturbative determination turned out to be in remarkable agreement with their result, which they considered to be only a conjecture.
On the numerical side, some of the early works used exact diagonalization in hyperspherical coordinates [35,41], whereby a large number of eigenstates can be calculated and their energies summed over to calculate canonical partition functions, thus providing access to b 3 (and to some extent b 4 , and with low accuracy b 5 ) for systems of cold atoms with short-range interactions.
In an outstanding numerical feat, Yan and Blume [44,45] designed an ad hoc Monte Carlo method to tackle the calculation of ∆b 4 for fermions at unitarity, resulting in the first determination of this quantity with stochastic methods. Their calculation featured a harmonic trapping potential, which induces a temperature dependence in ∆b 4 , which is expected to be temperature-independent in the unitary limit. (Below we will show a comparison between Yan and Blume's results and ours, when our method is generalized to include a harmonic trap.) The only important drawback of this work was the large uncertainty in the final result, induced by the increased stochastic noise as the trapping potential is removed. We will return to a discussion of this result below.
In spite of all of the above remarkable progress, it is evident that, due to the complexity of the n-particle quantum mechanical problem, a different kind of approach is needed if one is to determine high-order virial coefficients with well controlled systematic error. In particular, stochastic methods tend to have too large uncertainties to yield the accurate estimates needed to implement resummation techniques (we return to these below). On the other hand, direct numerical methods such as exact diagonalization can be very powerful in providing detailed information (furnishing not only energies but also the associated eigenstates), but have not yet succeeded in accurately determining virial coefficients beyond third order. One of the main objectives of this paper is to present our work in developing and applying a non-perturbative, semi-analytic, computational approach that is free of stochastic effects, beginning in the next section.

Homogeneous Fermi gases with a zero-range interaction
In this section we explain our method in detail and review its application to the simplest case, namely that of Fermi gases in homogeneous space, focusing on the VE for the pressure.

Factorizing the transfer matrix
The cornerstone of the approach is the Suzuki-Trotter factorization of the transfer matrix (i.e. the quantum version of the Boltzmann weight). To that end, the Hamiltonian is split into kinetic and potential energy terms, i.e.Ĥ =T +V, (20) such that the simplest symmetric Suzuki-Trotter factorization is where τ is in principle an arbitrary parameter, but we will define it such that for some integer N τ , one has β = τN τ , i.e. τ defines the imaginary time discretization. Note that, since our interest is in taking the trace of powers of e −τĤ , the same accuracy is obtained for the symmetric decomposition as for its asymmetric counterpart (Note: When calculating expectation values of operators, not mere traces, using the symmetric decomposition does make a difference.) The above factorization step is always needed in our method, regardless of whether the target system is in homogeneous space or in a trapping potential; in this section we focus on the former, returning to the latter in a later section.
Using the above factorization, the objective is to calculate Q N for the desired particle content at progressively larger values of N τ . The resulting Q N are then used to calculate the b k , and the limit of large N τ is taken at the end by extrapolation. The latter is an extrapolation to the continuous imaginary-time limit.
The fundamental building blocks in the calculation of Q N are the matrix elements of the factorized transfer matrix. For homogeneous non-relativistic spin-1/2 fermions with a zero-range interaction one hasT wheren σ (p) is the number density operator for particles of spin σ and momentum p, and we use m ≡ 1 for simplicity; moreover,V wheren σ (r) is the number density operator for particles of spin σ at position r, and we have regularized the problem by putting it on a spatial lattice of spacing . In practice, we renormalize this interaction by tuning it to reproduce the exact value of b 2 , as set by the Beth-Uhlenbeck formula reviewed above.
To calculate Q 2 , for example, the desired matrix elements are given by for the (1 + 1)-particle problem, where K(p) is the non-interacting, factorized Boltzmann weight C is the coupling strength and we use the notation |P = |p 1 p 2 · · · p a p a+1 · · · p a+b , for the (a + b)-particle system, where p 1 to p a are for spin-↑ particles and p a+1 to p a+b for spin-↓ particles. Similarly, for the (2 + 1)-particle problem, which enters in calculating Q 3 , one obtains Naturally, the complexity of the above matrix elements results mainly from the interaction elements P|e −τV |Q and rises rapidly with the particle number. Another example is the (2 + 2)-particle problem, which pertains to Q 4 , for which we obtain

From transfer matrices to canonical partition functions
In the above expressions for M 21 and M 22 we have not implemented any symmetrization or antisymmetrization, which is needed to account for quantum statistics. Without approximation or loss of generality, that operation can be carried out at the end of the calculation, i.e. upon taking the N τ -th power of the desired transfer matrix M ab , because the operators involved preserve the particle statistics. For example, in the fermionic case, the antisymmetrization yields whereas in the bosonic case one would have a symmetric form, namely To take the N τ -th power as shown above, we use automated tensor algebra (further details on this automation are presented below). Although the resulting number of terms is very large (on the order of 10 9 for the cases we have explored), it is manageable. Furthermore, each term is given by a multidimensional Gaussian integral (once the continuum and largevolume limits are taken) with an easily identifiable quadratic form. Since those integrals are easily evaluated using determinants, the entire process of algebraic manipulation and evaluation can be farmed out to massively distributed computing architectures as a large set of independent processes. In practice, we have been able to explore subspaces with up to 9 particles with several days of calculations on the Open Science Grid [46,47]. The number of particles one can analyze is of course limited by N τ ; for instance with 10 6 CPU hours, it is possible to study up to N τ = 23 for 3 particles, but only up to N τ = 4 for 9 particles.
Combining the results for the Q ab , one obtains the desired b k . It should be noted that, in practice, one does not use these quantities directly but rather the changes induced by the interaction as given by ∆Q ab and ∆b k .
The main advantage of this method is that it is not a stochastic approach; in fact, it is closer to an analytic approach, as it amounts to an automated, direct evaluation of a lattice field theory calculation. The automated algebra allows us to resolve the volume cancellations that plague the evaluation of virial coefficients, which stem from combining Q ab for varying a and b. As the latter scale as V a+b plus sub-leading terms, whereas the b k are volume independent, resolving those cancellations is both crucial and very difficult to achieve with stochastic approaches.

Computational details of automated algebra
In this section we present a more detailed technical discussion of our automated algebra method to capture the general idea represented in our code. The ultimate goal of the method is to evaluate the canonical partition functions as shown in Eq. (31) and (32), which involves three steps:

1.
Term generation: Expand the product M N τ mj symbolically, which will yield a large number of terms as N τ is increased.

2.
Delta crunch: Contract indices to saturate all Kronecker deltas, thus simplifying each term into a product of Gaussian functions, namely the propagator K(p), by integrating out a subset of variables. This is the most computationally expensive step.

3.
Gaussian integration: For each term, take the summation over the rest of the variables and take the continuum limit, ultimately turning each term into a multidimensional Gaussian integral whose results are analytically available as the well-known formula where n is the dimension of vector x.
We now proceed to elaborate on the above steps.
Step 1 As shown for instance in Eq. (30), the kinetic energy appears in our transfer matrix expressions as a prefactor, fully factorized across particles. The complexity of the problem lies in the interaction operators, of course. In this first step, we expand the product that results from such operators when N τ factors are present: where To evaluate the matrix sums with different subscripts [e.g. as in Eq. (31)], thus implementing the pertinent quantum statistics, we implement different boundary conditions on P (1) and P (N τ ) . For example, in Eq. (31), the first term on the right-hand side is a normal trace, which we obtain by imposing a periodic boundary condition p On the other hand, the second term is a kind of "shifted" trace, obtained by setting the boundary condition p We have considered so far the complete expansion of the product. However, thanks to the cyclic property of the trace, only a subset of the complete expansion needs to be evaluated. For example, for N τ = 2, there are two terms f 1 (P (1) , P (2) ) f 2 (P (2) , P (3) ) and f 2 (P (1) , P (2) ) f 1 (P (2) , P (3) ) at O(C 3 ), which are equivalent under the cyclic variable substitution P (1) → P (2) , P (2) → P (3) , P (3) → P (1) . Such a property defines a mathematical object called "combinatorial necklace" and the goal of this first step is then to generate all possible necklaces of functions f i , for which we used the algorithm developed in Ref. [48].
Step 2 Once we have expanded the product of interaction operators, each term is in the form of a product of propagators and δ's as K(P (1) )K(P (2) )K(P (3) ) · · · × ∆(P (1) , P (2) , P (3) , · · · ) (36) where K(P (i) ) is a shorthand for the kinetic-energy product K(p M+N ), and ∆(P (1) , P (2) , P (3) , · · · ) is a product of Kronecker δ's from the combination of f i (P (i) , P (i+1) ) functions. Equations (25), (29) or (30) exemplify what one would obtain in the simplest case of N τ = 1. Thus, the output of this step is effectively the tensor contraction (for all internal indices; not the trace indices) of such N τ = 1 results for N τ as large as needed.
The second step of the method is to carry out the sums on such a term over all momentum variables from P (2) to P (N τ −1) , i.e. all intermediate complete sets inserted. This operation is called "Delta crunch" as it will reduce the ∆ function by substituting all available momentum variables, i.e. "crunching" the δ's into the K factors. To this end, we loop through the δ's in the ∆ function one at a time, and perform variable substitution in both K and ∆. To provide efficiency in variable substitutions, the K and the δ's are represented as a hashmap, which offers a O(1) time in lookup and modification.
Step 3 Once the Delta crunch step is complete, the resulting expression will be a product of kinetic energy factors with corresponding variable substitutions, e.g.
Recalling that K(p) is a Gaussian function, the evaluation of the summation over the remaining momentum variables is easiest carried out in the continuum limit. To that end, we convert the above product into a quadratic form where p contains all the momentum variables. The matrix A is symmetric and positivedefinite, such that one can use Cholesky decomposition to evaluate the determinant, which is computationally more efficient and numerically more stable than LU decomposition.

Parallelization
Before concluding this section, we want to add one more technical note on the parallelization of our method. Compared to more conventional methods like QMC, ours is much easier to parallelize. All the terms generated in the first step are independent from each other, which means they can be evaluated in fully parallel fashion with little or no communication overhead among processes. Moreover, the evaluation of each term is cheap as it does not involve complicated linear algebra operations, and so it is suitable to run on any number of CPU cores. These features make our method ideal to run on a distributed, heterogeneous computing cluster, such as the Open Science Grid or the Folding@home project, where the computational power is unevenly distributed across nodes, in contrast with traditional supercomputers, while the number of available cores can be much higher.

Selected results
Using the method presented above, which resulted from a sequence of prior studies [14][15][16][17], we tackled the problem of calculating, with as high precision as possible, the virial coefficients up to ∆b 5 for homogeneous spin-1/2 Fermi gases with an attractive contact interaction, for varying coupling strengths. For illustrative purposes, we focus here on the 3D case [18], but extensions to other dimensions were also studied [19].
In Fig. 1, we present the coefficients ∆b k for k = 3, 4, 5 (left panel) and the corresponding subspace contribution ∆b ij (right panel). The inset on the left panel shows a comparison with experiments [49,50] and theory [36,39,45]. Besides the excellent agreement with the ∆b 3 results of Ref. [33], the above provide a precise non-perturbative, non-stochastic determination of ∆b 4 and ∆b 5 . The results show that, due to the competing subspace contributions (right panel), both ∆b 4 and ∆b 5 are nonmonotonic as a function of the coupling strength. In addition, when approaching the unitary limit it is evident that ∆b 5 becomes comparable in magnitude to ∆b 4 . It is this property that complicates the experimental determination of ∆b 4 at strong coupling (dashed error bars in inset of left panel), as one is forced to assume that ∆b 5 is comparatively small, which is not the case. The agreement for ∆b 4 with the estimates of Refs [36,39,45], which use entirely different methods (among them and with the present work), further supports the accuracy and precision of the method. Note, however, the wide disparities in results for ∆b 22 compared to ∆b 31 , which indicate that the (2 + 2)-particle subspace is a considerably more difficult problem to tackle than the polarized (3 + 1) subspace.

Harmonic traps
In the case of harmonically trapped particles, it is convenient to split the Hamiltonian into a noninteracting harmonic oscillator pieceĤ 0 , and the interactionV in order to perform the Suzuki-Trotter factorization, such that whereĤ 0 =T +V HO andV is the harmonic trapping potential, and therefore The main advantage of this choice of splitting ofĤ is that, as we will see below, the resulting factorization of the transfer matrices can be easily written in terms of the so-called Mehler kernel [15]. It should be pointed out that the other possible factorization where one combineŝ V HO withV rather than withT is an interesting possibility that one could explore entirely in momentum space (as in the homogeneous cases of the previous section). However, such an approach effectively results in a more complicated interaction made out of a one-body piece and a two-body piece with independent coupling constants, which must be kept track of throughout the calculation, with the concomitant increase in computational complexity. We will return to a similar situation below, namely the case of attractive Bose gases, which require two-and three-body forces for stability reasons. As anticipated, the transfer matrices can be written conveniently in terms of the Mehler kernel. For example, when examining the (1 + 1)-particle subspace, we obtain M 11 = x 1 , x 2 |e −τ(T+V HO ) e −τV |y 1 , y 2 = ρ(x 1 , y 1 )ρ(x 2 , y 2 )[1 + Cδ(y 1 − y 2 )], (42) where ρ(x, y) is the Mehler kernel which encodes the coordinate representation of the transfer matrix of a single-particle in the harmonic potentialV HO . Here, Z T = (x T /λ T , y T /λ T ), and where 1 is a 3 × 3 unit matrix. Similarly, for the (2 + 1)-particle subspace, we find Note that, as in a previous section, we do not include any symmetrization or antisymmetrization at this stage, as quantum statistics can be accounted for after the N τ -th power is taken. Further transfer matrices for up to 5 particles can be found in Ref. [20].
As in the homogeneous case, we rely on known, analytic results for ∆b 2 in order to renormalize the interaction. In a harmonic trapping potential ∆b T 2 is given at unitarity by [3] ∆b We note that, for arbitrary order, the trapped virial coefficients ∆b T n in the high-temperature (or low-frequency) limit βω → 0 is connected with their homogeneous counterparts via which is useful to know when comparing homogeneous results with low-frequency extrapolations from the trapped case.
Using our framework at leading order (N τ = 1) in d spatial dimensions, it is not difficult to obtain analytic formulas such as and Following the same approach outlined in previous sections, we worked out the next-to-leading order (N τ = 2) formulas, which the reader can find in Ref. [20]. With automated algebra, these can be pushed to much higher N τ to obtain estimates for ∆b k which are then extrapolated to large N τ . In Fig. 2, we show our results at N τ = 1, 2 for the full-space contribution ∆b k , and those extrapolated to N τ → ∞ limit for both ∆b k and the corresponding subspace contribution ∆b ij . We find very good agreement between the extrapolated results and the Monte Carlo calculations of Ref. [45] for large trapping frequency βω. As βω approaches 0, our results are smoother compared to the Monte Carlo results, likely due to the volume cancellations that we are able to resolve analytically but which induce noise in the Monte Carlo method. Although the leading-and next-to-leading-order results deviate from the expected curve, they capture most qualitative features and may shed light into higher-order virial coefficients where the computational costs are too large for a full calculation at large N τ . Lastly, note that for ∆b 4 , the leading-order result is closer to the extrapolated curve compared to the next-leading-order result. This indicates that the discretization error from the Suzuki-Trotter decomposition is not monotonic in N τ . In other words, as N τ increases, we may expect worse results before the asymptotic regime is approached. The onset of that regime is of course coupling dependent. Therefore, at a given interacting strength it is essential to investigate N τ as high as possible in order to obtain quantitatively correct results.   [36], and the open circle is to emphasize the results in the βω → 0 limit. The PIMC results from Ref. [45] are shown as red circles.

Neutron Matter
The outermost crust of a neutron star is a low-density regime for which the unitary Fermi gas results presented above can be regarded as an approximation. To go beyond the unitary limit, a realistic model for neutron matter must include interaction range effects. In this section we elaborate on a possible implementation of a finite-range interaction in the context of our method. To that end, we anticipate that it will be essential to expand the matrix elements of e −τV as a sum of Gaussians, as that would take advantage of the computational framework developed in the simpler cases discussed above.
As a reminder, for a contact interaction in the two-particle subspace one has q 1 q 2 |e −τV |p 1 p 2 = δ q 1 ,p 1 δ q 2 ,p 2 + C δ q 1 +q 2 ,p 1 +p 2 , where C includes the lattice spacings and the bare coupling. Naturally, the first term represents a noninteracting piece and the second term encodes the interaction effects. A contact interaction in coordinate space has constant matrix elements in momentum space, but in the presence of a finite range there will be nontrivial momentum dependence and therefore we generalize the above via where q r = q 1 − q 2 and p r = p 1 − p 2 . In the limit λ 1 , λ 2 → 0, we reproduce the zero-range model. This type of expansion can be further generalized by including powers of q 2 r , p 2 r and q r · p r as a prefactor rather than in the exponent, which can be used to account for more features of the interaction. However, such a generalization should be used sparingly or not at all if possible, as it would dramatically (specifically, factorially) increase the number of terms that result in the expansion, by virtue of Wick's theorem.

Unitary Bose gas
The case of Bose gases with attractive interactions, in particular close to the universal regime of the unitary limit, requires special care. Given the absence of Pauli exclusion, attractively interacting bosons undergo Thomas collapse [51], i.e. they are unstable. Moreover, close to unitarity they display the Efimov effect [52,53]. To properly renormalize such a system, a repulsive three-body interaction is required [54]. The latter introduces a new dimensionful parameter that is sensitive to the ultraviolet cutoff and must be fixed by using some known physical quantity.
To account for the above, we consider an interaction of the form where we have used the normal-ordered form andâ † r ,â r are the bosonic creation and annihilation operators at point r.
In previous sections we used ∆b 2 to renormalize a two-body contact interaction. It is therefore natural in the present case to use ∆b 2 and ∆b 3 for the same purpose. Since ∆b 2 involves only the two-particle subspace, we use it to renormalize the two-body interaction in the usual way before proceeding to ∆b 3 , which depends on both the two-and three-particle subspaces and thus determines the three-body coupling. In the bosonic unitary limit, the threebody parameter induces a temperature dependence on ∆b 3 , which was calculated exactly in Ref. [55].
Calculating ∆b 2 using the above interaction in d spatial dimensions and N τ = 1 we obtain where V = L d is the d-dimensional volume, such that in the (spatial) continuum limit where we used that in that limit Q 1 /V → λ −d T . On the other hand, a calculation of ∆b 3 , also in d spatial dimensions and N τ = 1, yields which in the (spatial) continuum limit becomes We thus see that, as anticipated, ∆b 2 and ∆b 3 determine g 2 and g 3 . At the N τ = 1 level of this calculation, the temperature dependence of ∆b 2 and ∆b 3 will typically induce a temperature dependence in g 2 and g 3 . More explicitly, we may invert the above results to obtain the dimensionless form Armed with these answers, the renormalization program is complete and we may proceed to calculate and predict higher-order virial coefficients. Our leading order (N τ = 1) for the fourth-order bosonic virial coefficient ∆b 4 is neatly expressed in terms of a relationship with ∆b 2 and ∆b 3 : In Fig. 3 we show the above result applied to d = 3 at unitarity, where ∆b 2 = √ 2 and ∆b 3 is as determined in Ref. [55]. Our N τ = 1 approximation is likely a rather crude one, but it should be asymptotically correct for small βE T .

Pressure, density, and generalization to other observables
So far, we limited our scope to the expansion for pressure, but the virial expansion is much more general and can be applied to a wide range of observables. In this section, we present the generalization of our approach to the Tan contact, the momentum distribution, and the structure factor.

From pressure to density and Tan's contact
By way of the grand thermodynamic potential, the b n grant access to all thermodynamic quantities, at least in principle. One of the most basic quantities besides the pressure is the density, which is given by where n 0 is the density at absence of interaction, and similarly in the presence of finite polarization (i.e. different fugacities z ↑ and z ↓ for each spin), where s =↑, ↓ for each spin.
Beyond global thermodynamic quantities, the Tan contact can also be investigated through the virial expansion. For systems with short-range two-body forces, Tan's contact represents the probability of finding two particles at the same spatial location. For that reason, it captures the short-distance behavior of all correlation functions. In particular, the momentum distribution in such systems decays at large momentum k as where I is the contact. Thus, one way to measure or calculate the contact is to determine the large-momentum tail of the momentum distribution. In practice, however, it has proven more efficient to use the so-called adiabatic relation, whereby I is obtained from the variation of the grand thermodynamic potential with the coupling strength. Below we focus on the 3D case, but similar equations can be derived in one-and two-spatial dimensions and readers are referred to Ref. [19] for more details.
In three spatial dimensions, the VE for the contact is obtained via where are the virial coefficients of I. To evaluate the partial derivative, the most straightforward method is to apply the chain rule as where the first derivative can be calculated numerically and the second is analytically given by the Beth-Uhlenbeck formula as Thanks to the analytical nature of our method, one can improve the accuracy of the numerical derivative without repeating calculations. Alternatively, one can take advantage of the analytic form of ∆b k from our method, which is a polynomial in terms of coupling strength C where l max = min(M, N) · N τ is the degree of the polynomial. The derivative can then be treated as a polynomial in C of degree l max − 1 as where D = ∂C ∂l and B l (D) = l A l D. Now, we treat C and D as two independent variables and tune their values in two passes: firstly the value of C is renormalized to reproduce the ∆b 2 ; and then we plug the resulting C in Eq. (70) and tune D to produce the expected second-order result as given in Eq. (68).
In Fig. 4, we plot the density (left panel) and Tan's contact (right panel) in the homogeneous case. On the left panel, the solid lines (with uncertainty bands) are the results of using the VE at face value, i.e. truncated at a given order: blue for third order, red for fourth order, and green for fifth order. In each case, a corresponding dashed line of the same color shows the results of Padé resummation at that order. Finally, the second-order results are shown as a black dashed line for reference, while the black dots come from the MIT experiment of Ref. [50]. For both the truncated VE and the resummed results, we found improved agreement as higher-order contributions were included, the effect being most notable when the resummation is applied. This inspires an open question under investigation on the effect of the resummation when using even higher orders. This is particularly intriguing as access to higher-order coefficients would allow more freedom in exploring a range of resummation techniques.
In the right panel of Fig. 4 we show the dimensionless contact in the form as a function of the reduced temperature T/T F , which is related to the density via Therefore, the resulting curve is an interplay between two independent series for the density and the contact. In particular, the curve will be more sensitive to the density estimate as both quantities implicitly depend on it. In other words, in the region where we observed deviation in the left panel between our results and experimental determinations, we expect errors in both horizontal and vertical directions. This approximately corresponds to z 0.75, or T/T F 1.4, for the truncated VE results. At higher temperature, where our results coincide with the experimental measurement, the main source of error is the expansion series for I. To better show I as given in Eq. (65), and isolate the influence on the density, we use the experimental values [56] to calculate the reduced temperature T/T F (thus providing an essentially exact density) and evaluate the series ∑ ∞ k=2 c k z k using Padé resummation at different orders (dashed lines). Overall, we observed the same trend as for the density where the results are improved as higher-order contributions are included. Although the resummed fifth-order VE (red dashed line) seems to give very good results below the critical temperature, one should take such an agreement with a grain of salt, as the error in the dimensionless contact is greatly reduced by a factor of (T/T F ) 2 . Furthermore, this result is free from errors in the density because we have used the experimental data for the latter.
To represent the more practical situation where we have no access to accurate experimental measurement of the density, we use the best density estimate (green dashed line in the left panel) and the final result is plotted as the thick red solid line. Even though the approximation begins to fail when approaching the critical temperature, we find impressive agreement until T/T F ≈ 0.4, roughly corresponding to the diverging point between the Padé and experimental results at z = 5.0 for the density.    [50]. Right: Dimensionless contact I/(Nk F ) as a function of reduced temperature T/T F . The blue, orange and green solid lines are the results using truncated VE for both the contact and the density, which is connected to the reduced temperature T/T F . The same color code is used. Colorful dotted lines are the results using the experimental density measurement and VE contact. The red curves are the results using resummation: the dotted is with experimental density and the solid with Padé resummed density. The light gray and brown points are experimental determinations from Ref. [56] and Ref. [57] respectively.

One-body operator example: the momentum distribution
Besides the density and pressure equations of state, the next simplest quantity one can calculate is the expectation value of a one-body operator, for instance the momentum distributionn σ (p) of particles with spin σ. Of course, here the only added difficulty is that the operator singles out a specific momentum p, but it warrants special attention.
The starting point is the thermal expectation value Having available the virial expansion for Z already, we focus on the numerator Thus, after expanding the denominator Z, where and so forth. (Note that, naturally, the content of m σ,1 (k) is entirely noninteracting, as it corresponds to a single-particle subspace.) The new "virial coefficients" m σ,k (p) requires the calculation of the Hilbert trace tr N [e −βĤn σ (p)], denoted as W N [n σ (p)] or just W N,σ for shorthand. In our method, its evaluation shares the same formalism as the usual partition function Q N as in Eqs. (31) and (32), i.e. the trace W N,σ encodes the particle statistics. Taking the (2, 1)-system for example, where N σ,p is the matrix elements of momentum density operator N 21,σ (p) = p 1 p 2 p 3 |n σ (p)|q 1 q 2 q 3 = δ pp 1 + δ pp 2 p 1 p 2 p 3 |q 1 q 2 q 3 = δ pp 1 + δ pp 2 δ p 1 q 1 δ p 2 q 2 δ p 3 q 3 . (81) At leading-order up to the third order in fugacity, we obtain the change of momentum distribution with respect to that in the non-interacting case, wherek = λ T k is the dimensionless momentum. As in previous cases, we emphasize that this is merely an approximation at the lowest non-trivial order in C, but it provides a qualitative guide and a reasonable answer at weak coupling.

Two-body operators and real-time evolution
An advantage of our automated algebra method is its versatility. The formalism and implementation can be adapted to more complicated observables without introducing significant new complexity. In this section, we present the first steps for a few such examples, which show the way for future work.
As a natural extension of the case of one-body operators, Eq. (74) is generalized to two-body operators: From this point on, the evaluation of the N-particle Hilbert space trace W(Ô 1 ,Ô 2 ) = tr N [e −βĤÔ 1Ô2 ] follows the same methodology explained in the previous section.
An interesting two-body quantity is the density-density correlation n(r)n(0) , which is the central ingredient of the static structure factor: S(q) = d 3 re −iq·r δn(r)δn(0) where δX =X − X . The latter approximately describes the in-medium scattering rate in certain contexts [58]. Recent studies [59,60] examined this quantity within the VE up to second order, and it is interesting to examine the effect of higher-order contributions. Yet another example that has been intensely investigated recently is the viscosity. Specifically, in the work of Refs. [61][62][63], the VE was applied to the calculation of both bulk and shear viscosity, which requires the commutator of the contact operator.
Our formalism can also be generalized beyond equilibrium systems by including realtime evolution. The recent work of Ref. [64] investigated the effects of an interaction quench in a bosonic system, in an attempt to explain the results of the experiment of Ref. [65]. Here, the system starts as noninteracting and in thermal equilibrium, and then an interparticle interaction is suddenly switched on. The process was investigated in Ref. [64] using the VE of the momentum distribution n(k) up to second order and varying behaviors at low and high momenta k were found. Although this is an interesting result, the effect of higher-order contributions remain an open problem and may not be small, which further motivates the extension of our framework for more general dynamic processes.
The quench process mentioned above is described by the Hamiltonian whereĤ 0 is the free Hamiltonian andV is the interaction, which is turned on at t = 0. To study the evolution of a given operatorÔ one starts with expressions similar to Eq. (83), with the exception that the Hilbert space trace now becomes time-dependent W N (Ô, t) = tr N (e −βĤ 0 e itĤÔ e −itĤ ).
Work in investigating these kinds of expressions with automated algebra, specifically for the interaction quench, is underway.

Summary and conclusions
The VE is an expansion of the quantum many-body thermodynamics in powers of the fugacity z that is capable of non-perturbatively characterizing such many-body systems in dilute, high-temperature regimes. The challenge of the VE is that calculating its coefficients b k at n-th order in the expansion requires analyzing the quantum n-body problem in some detail.
In this brief review, we have outlined some of the main methods used to calculate b k and advocated an approach based on discretizing imaginary time, thus obtaining a factorized transfer matrix, and automating the resulting algebra to calculate the trace of the N τ -th power of the transfer matrix. We have shown in detail how such a method is structured and how it is able to provide reliable results for up to b 5 for strongly coupled nonrelativistic fermions, focusing on the unitary limit. The automated algebra algorithm is fully parallelizable and may thus be extended beyond b 5 with increased computational power.
We have also shown how to generalize the approach to account for harmonic trapping potentials, neutron matter and attractive Bose gases (which require three-body forces). Furthermore, the method is also straightforwardly generalized to observables other than the pressure, such as the Tan contact, the momentum distribution, and the structure factor.