Third-order perturbative lattice and complex Langevin analyses of the finite-temperature equation of state of non-relativistic fermions in one dimension

TL;DR

This paper develops and applies third-order lattice perturbation theory and complex Langevin methods to study the finite-temperature equation of state of one-dimensional non-relativistic fermions, achieving agreement with non-perturbative results and extending to strong couplings.

Contribution

It introduces an unconventional lattice perturbation approach using Hubbard-Stratonovich transformation and automates Wick's theorem, along with a dimension-independent technique for Matsubara sums, and applies these to fermions at finite temperature.

Findings

Excellent agreement between perturbative and non-perturbative results at weak couplings.

Predictions for strong couplings using complex Langevin methods.

Perturbative calculations of up to fifth-order virial coefficients.

Abstract

We analyze the pressure and density equations of state of unpolarized non-relativistic fermions at finite temperature in one spatial dimension. For attractively interacting regimes, we perform a third-order lattice perturbation theory calculation, assess its convergence properties by comparing with hybrid Monte Carlo results (there is no sign problem in this regime), and demonstrate agreement with real Langevin calculations. For repulsive interactions, we present lattice perturbation theory results as well as complex Langevin calculations, with a modified action to prevent uncontrolled excursions in the complex plane. Although perturbation theory is a common tool, our implementation of it is unconventional; we use a Hubbard-Stratonovich transformation to decouple the system and automate the application of Wick's theorem, thus generating the diagrammatic expansion, including symmetry factors, at any desired order. We also present an efficient technique to tackle nested Matsubara frequency sums without relying on contour integration, which is independent of dimension and applies to both relativistic and non-relativistic systems, as well as all energy-independent interactions. We find exceptional agreement between perturbative and non-perturbative results at weak couplings, and furnish predictions based on complex Langevin at strong couplings. We additionally present perturbative calculations of up to the fifth-order virial coefficient for repulsive and attractive couplings. Both the lattice perturbation theory and complex Langevin formalisms can easily be extended to a variety of situations including polarized systems, bosons, and higher dimension.

Figures

Tables

Table 1. Table 1: Detail of next-to-leading order (NLO), next-to-next-to-leading order (N2LO), and next-to-next-to-next-to-leading order (N3LO) contributions (respectively, order A 2 superscript 𝐴 2 A^{2} , A 4 superscript 𝐴 4 A^{4} , and A 6 superscript 𝐴 6 A^{6} ) to the grand-canonical partition function 𝒵 𝒵 \mathcal{Z} . The indicated diagram figure refers to the corresponding fully-connected Feynman diagram or product of such diagrams for each contribution. The volume V = N x ​ N τ 𝑉 subscript 𝑁 𝑥 subscript 𝑁 𝜏 V=N_{x}N_{\tau} that appears in all expressions refers to the spacetime volume of the lattice. The greatest number of loops that appear for each fully-connected diagram is also provided. It is implicit in the notation that a momentum-conserving Kronecker delta has been utilized to eliminate momentum sums: one in the case of S 4 subscript 𝑆 4 S_{4} and S 4 ( 1 ) superscript subscript 𝑆 4 1 S_{4}^{(1)} , but two in the case of S 6 subscript 𝑆 6 S_{6} .

Order	Diagram Figure	Prefactor	Diagram Frequency Sum	Loops
NLO	Fig. 3	$1/\left(2V\right)$	${\left[\sum _{{𝐩}_1}{S}_1\left(Q_1\right)\right]}^2$	1-loop
N2LO	${[\mathrm{Fig}.\text{<a><span>3</span></a>}]}^2$	$1/\left(8V^2\right)$	${\left[\sum _{{𝐩}_1}{S}_1\left(Q_1\right)\right]}^4$	1-loop
N2LO	Fig. 4(a)	$-1/\left(4V^2\right)$	$\left[\sum _{{𝐩}_1}{S}_1^{\left(1\right)}\left(Q_1\right)\right]{\left[\sum _{{𝐩}_1}{S}_1\left(Q_1\right)\right]}^2$	1-loop
N2LO	Fig. 4(b)	$1/\left(8V^2\right)$	$\sum _{{𝐩}_1,{𝐩}_2,{𝐩}_3}{S}_4(Q_1,Q_2,Q_3)$	3-loop
N3LO	${[\mathrm{Fig}.\text{<a><span>3</span></a>}]}^3$	$1/\left(48V^3\right)$	${\left[\sum _{{𝐩}_1}{S}_1\left(Q_1\right)\right]}^6$	1-loop
N3LO	$\mathrm{Fig}.\text{<a><span>4</span></a>}\left(a\right)\times \mathrm{Fig}.\text{<a><span>3</span></a>}$	$-1/\left(8V^3\right)$	$\left[\sum _{{𝐩}_1}{S}_1^{\left(1\right)}\left(Q_1\right)\right]{\left[\sum _{{𝐩}_1}{S}_1\left(Q_1\right)\right]}^4$	1-loop
N3LO	$\mathrm{Fig}.\text{<a><span>4</span></a>}\left(b\right)\times \mathrm{Fig}.\text{<a><span>3</span></a>}$	$1/\left(16V^3\right)$	$\left[\sum _{{𝐩}_1,{𝐩}_2,{𝐩}_3}{S}_4(Q_1,Q_2,Q_3)\right]{\left[\sum _{{𝐩}_1}{S}_1\left(Q_1\right)\right]}^2$	3-loop
N3LO	Fig. 5(a)	$1/\left(12V^3\right)$	$\left[\sum _{{𝐩}_1}{S}_1^{\left(2\right)}\left(Q_1\right)\right]{\left[\sum _{{𝐩}_1}{S}_1\left(Q_1\right)\right]}^3$	1-loop
N3LO	Fig. 5(b)	$1/\left(8V^3\right)$	${\left[\sum _{{𝐩}_1}{S}_1^{\left(1\right)}\left(Q_1\right)\right]}^2{\left[\sum _{{𝐩}_1}{S}_1\left(Q_1\right)\right]}^2$	1-loop
N3LO	Fig. 5(c)	$-1/\left(4V^3\right)$	$\left[\sum _{{𝐩}_1,{𝐩}_2,{𝐩}_3}{S}_4^{\left(1\right)}(Q_1,Q_2,Q_3)\right]\left[\sum _{{𝐩}_1}{S}_1\left(Q_1\right)\right]$	3-loop
N3LO	Fig. 5(d)	$1/\left(12V^3\right)$	$\sum _{{𝐩}_1,{𝐩}_2,{𝐩}_3,{𝐩}_4}{S}_6(Q_1,Q_2,Q_3,Q_4)$	4-loop

Table 2. Table 2: Results for the second, third, fourth, and fifth-order virial coefficients of the pressure P 𝑃 P [see Eq. ( 66 )] at NLO, N2LO, and N3LO for two repulsive and attractive couplings. All coefficients are computed for a spatial lattice size of N x = 100 subscript 𝑁 𝑥 100 N_{x}=100 , β = 8.0 𝛽 8.0 \beta=8.0 , and a temporal lattice spacing of τ = 0.05 𝜏 0.05 \tau=0.05 . The fifth-order virial expansion is displayed in Fig. 8 .

$\lambda$	Coefficient	NLO	N2LO	N3LO
-1.0	$b_2$	$-0.551$	$-0.473$	$-0.473$
	$b_3$	$0.472$	$0.373$	$0.342$
	$b_4$	$-0.452$	$-0.380$	$-0.280$
	$b_5$	$0.449$	$0.446$	$0.246$
-0.5	$b_2$	$-0.453$	$-0.433$	$-0.433$
	$b_3$	$0.333$	$0.308$	$0.304$
	$b_4$	$-0.289$	$-0.271$	$-0.259$
	$b_5$	$0.270$	$0.269$	$0.244$
0.5	$b_2$	$-0.253$	$-0.233$	$-0.233$
	$b_3$	$0.0508$	$0.0254$	$0.0295$
	$b_4$	$0.0408$	$0.0592$	$0.0461$
	$b_5$	$-0.0925$	$-0.0932$	$-0.0671$
1.0	$b_2$	$-0.152$	$-0.0712$	$-0.0712$
	$b_3$	$-0.0922$	$-0.194$	$-0.161$
	$b_4$	$0.208$	$0.282$	$0.176$
	$b_5$	$-0.276$	$-0.279$	$-0.0677$

Equations

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

$$\hat{H}=\hat{T}+\hat{V},$$

$$e^{-\tau \hat{H}}\simeq e^{-\tau \hat{T}/2}{e}^{-\tau \hat{V}}{e}^{-\tau \hat{T}/2},$$

$$e^{-\tau \hat{V}}=\int \mathcal{D}\sigma e^{-\tau \hat{W}[\sigma ]},$$

$$\mathcal{Z}=\int \mathcal{D}\sigma {\det }^{2}M[\sigma ],$$

$$M[\sigma ]\equiv \left(\begin{array}{cccc}1& 0& \cdots & {\mathcal{U}}_{N_{\tau }}[\sigma ]\\ -{\mathcal{U}}_1[\sigma ]& 1& \cdots & 0\\ 0& -{\mathcal{U}}_2[\sigma ]& \ddots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& \cdots & -{\mathcal{U}}_{N_{\tau }-1}[\sigma ]& 1\end{array}\right).$$

$${\left[{\mathcal{U}}_j[\sigma ]\right]}_{\mathbf{x},{\mathbf{x}}^{\prime }}\equiv {\left[e^{-\tau ({\mathbf{p}}^2/2m-\mu )}\right]}_{\mathbf{x},{\mathbf{x}}^{\prime }}\left[1+A\sin \sigma ({\mathbf{x}}^{\prime },t)\right]$$

$$\int \mathcal{D}\sigma =\prod _{\mathbf{x},\tau }{\int }_{-\pi }^{\pi }\frac{d\sigma (\mathbf{x},\tau )}{2\pi }.$$

$$\mathcal{Z}=\int \mathcal{D}\sigma \det M_{\uparrow }[\sigma ]\det M_{\downarrow }[\sigma ].$$

$$M=M_0+A\mathcal{T}\mathcal{S}[\sigma ]$$

$$\mathcal{T}\equiv \left(\begin{array}{cccc}0& 0& \cdots & T\\ -T& 0& \cdots & 0\\ 0& -T& \ddots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& \cdots & -T& 0\end{array}\right),$$

$${\left[\mathcal{S}[\sigma ]\right]}_{\mathbf{x},t;{\mathbf{x}}^{\prime },t^{\prime }}={\delta }_{t,t^{\prime }}{\delta }_{\mathbf{x},{\mathbf{x}}^{\prime }}\sin \sigma (\mathbf{x},t).$$

$$\det M=\det M_0\det (\text{openone}+AK\mathcal{S}[\sigma ]).$$

$$\det (\text{openone}+AK\mathcal{S}[\sigma ])$$

$$E_k=\frac{(-1{)}^{k+1}}{k}\mathrm{tr}\left[{\left(K\mathcal{S}[\sigma ]\right)}^k\right].$$

$$\exp \left[\sum _{k=1}^{\infty }{A}^k{E}_k\right]$$

$$=$$

$$\begin{gathered} \det (\text{openone}+AK\mathcal{S}[\sigma ])=1+A{E}_1+A^2\left(\frac{E_1^2}{2!}+E_2\right) \\ +A^3\left(\frac{E_1^3}{3!}+E_1{E}_2+E_3\right)+\cdots . \end{gathered}$$

$$\int \mathcal{D}\sigma E_2[\sigma ]$$

$$\begin{gathered} I_m(b_1,b_2,\dots ,b_m)= \\ \int \mathcal{D}\sigma \sin (\sigma [b_1])\sin (\sigma [b_2])\dots \sin (\sigma [b_m]), \end{gathered}$$

$$I_2(b_1,b_2)$$

$$I_4(b_1,b_2,b_3,b_4)=$$

$$+\frac{1}{4}$$

$${\int }_{-\pi }^{\pi }\frac{d\sigma }{2\pi }{\sin }^4\sigma =\frac{3}{8}.$$

$$[K{]}_{ab}=\sum _q{U}_{aq}^{†}[D_0{]}_q{U}_{qb}$$

$$[D_0{]}_q\equiv \frac{1}{1-\exp [i\omega +\tau ({\mathbf{k}}^2/2-\mu )]},$$

$$U_{aq}=\frac{1}{\sqrt{N_x^d{N}_{\tau }}}\exp [i(\omega t-\mathbf{k}\cdot \mathbf{x})].$$

$$S=\sum [K{]}_{a_1{a}_2}[K{]}_{a_3{a}_4}\dots [K{]}_{a_{2m-1}{a}_{2m}}{I}_{a_1\dots a_{2m}}$$

$$S=\sum [D_0{]}_{q_1}[D_0{]}_{q_2}\dots [D_0{]}_{q_m}{\tilde{I}}_{q_1,\dots ,q_m}$$

$$\sum _{q_1,q_2,\dots ,q_m}[D_0{]}_{q_1}[D_0{]}_{q_2}\dots [D_0{]}_{q_m}\delta (\{q_k\}),$$

$$S_1=\sum _{n=1}^{N_{\tau }}[D_0{]}_n=\sum _{n=1}^{N_{\tau }}\frac{1}{1-Q\exp (i{\omega }_n)},$$

Full text

I Introduction

By far most problems in quantum many-body physics, from quantum chemistry to condensed matter and high-energy areas, suffer from the so-called sign problem. This problem is an exponential deterioration of the signal-to-noise ratio that appears when estimating quantum expectation values stochastically (i.e. via some form of quantum Monte Carlo), and is present in nearly all non-relativistic models at finite polarization and relativistic ones at finite chemical potential (see e.g. SignProblemReview1 ; SignProblemReview2 ; SignProblemReview3 for reviews). While sufficient conditions for the avoidance of the sign problem have been under intense study recently ZhangEtAl ; LiEtAl , a general solution has been shown to be a formidable challenge TroyerWiese . Motivated in part by the above, we set out to examine a classic many-body problem normally plagued by a sign problem using both perturbative and non-perturbative methods, in all cases regulated by a spatial lattice. Specifically, we focus on a system of non-relativistic fermions in one spatial dimension (1D) with a contact interaction, for which a sign problem is present when interactions are repulsive.

On one hand, we calculate analytically the first few orders of the lattice perturbative expansion of the pressure, aiming to assess their reliability towards characterizing thermodynamic equations of state. The results of this calculation provide a posteriori motivation for our present work: First, the perturbative finite-temperature lattice results and formalism shown here are unconventional, as perturbative calculations are typically done in the continuum; second, we have found methodological aspects that also appear to have been previously overlooked, which greatly simplify the calculation on the lattice, and which apply to a very wide range of systems (including relativistic ones); third, we consider 1D as an initial benchmarking step towards more demanding, higher-dimensional calculations; finally, semi-analytic tools like perturbation theory can address the shortcomings of conventional quantum Monte Carlo calculations where the sign problem is present, at least in some regimes.

On the other hand, we extend the hybrid Monte Carlo calculations of Ref. EoS1D and study the repulsive case by applying complex Langevin (CL) techniques. The latter have garnered considerable attention in the last few years in the area of finite-density lattice QCD FiniteDQCD0 ; FiniteDQCD1 ; FiniteDQCD2 ; FiniteDQCD3 ; FiniteDQCD4 ; FiniteDQCD5 ; FiniteDQCD6 ; FiniteDQCD7 ; FiniteDQCD8 , but their application to non-relativistic matter is almost nonexistent (see however CLYamamoto1 ; CLYamamoto2 ). We use the CL method with a modified action to prevent uncontrolled excursions into the complex plane, which would otherwise lead the method to converge to an incorrect result. While we focus here on a 1D problem, all the methods used extend easily to arbitrary dimensions, as well as to polarized and multi-component systems. However, a comment is in order that is specific to 1D: As pointed out in Ref. BA , in 1D one may apply the thermodynamic Bethe ansatz in certain finite-temperature regimes, but it is generally an uncontrolled approximation as it involves solving an infinite tower of non-linear integral equations (see however Ref. KakashviliBolech ).

Our motivation for studying 1D fermions is in part computational convenience, as mentioned above, and in part due to the fact that the general 1D quantum many-body problem remains an area of great interest for condensed matter, atomic, and high-energy physics. Substantial progress was made in the last few decades in particular for the 1D Hubbard model (see 1DHubbardBook ) using numerical methods such as quantum Monte Carlo (see e.g. 1DMC1 ; 1DMC2 ; 1DMC3 ), density matrix renormalization group (see e.g. DMRG1 ; DMRG2 ; DMRG3 ; DMRG4 ), and the Bethe ansatz (see BA and BatchelorEtAl ; BatchelorEtAlGaudinYang1 ; BatchelorEtAlGaudinYang2 ), as well as with analytic approaches (such as bosonization Bosonization and beyond-mean-field approaches RBD ). Remarkably, there are aspects of this problem that remain elusive in the sense that they lie beyond the Luttinger liquid paradigm, as pointed out in Ref. RMPGlazman , which justifies more detailed studies. On the experimental side, ultracold atoms in optical lattices continue to provide an unparalleled realization of clean, malleable, fermionic and bosonic systems, and thus these systems continue to shed light on multiple aspects of strongly coupled quantum dynamics (see ExpReviewLattices ; 1DExp ; 1DSUNExp ).

Below we present the lattice perturbation theory formalism (Sec. II), which is unconventional in that it uses a Hubbard-Stratonovich (HS) transformation HS1 ; HS2 to facilitate the automation of generating a diagrammatic expansion. With the diagrams in hand, the evaluation of physical observables proceeds in two steps: first, the Matsubara-frequency sums are done with a technique that allows us to carry them out simultaneously (rather than iteratively, as is the case when using countour integrals, see e.g. Ref. Nieto ); second, the remaining sum over spatial momenta is calculated numerically. Performing the frequency sums analytically is primarily avantageous in that it yields a vastly improved scaling of the computation time, as the total number of nested sums is greatly reduced. In Sec. III we present a brief overview of the CL method, contrast it with conventional hybrid Monte Carlo, and introduce a modified action to stabilize the calculation. In Sec. IV we show our results for the pressure and density equations of state up to next-to-next-to-next-to (N3LO) leading order in perturbation theory for both attractive and repulsive couplings, and compare with hybrid Monte Carlo for the case of attractive interactions where no sign problem is present. We also demonstrate results for a modified-action form of real and complex Langevin calculations for both attractive and repulsive couplings, and compare with our perturbative results. Using particle-number projection, we also provide the corresponding perturbative expansion for the first few virial coefficients. In Sec. V we summarize and present our conclusions. The appendices contain further details on the Matsubara sums involved and how to calculate them.

II Perturbation Theory Formalism

We first present the formalism for our method of computing the perturbative expansion for the grand-canonical partition function of interacting quantum gases on the lattice. In a nutshell, the procedure involves analytically evaluating expressions for the diagrammatic expansion and the resulting path integral on a computer using an object-oriented approach. The final expressions, which contain sums over the complete frequency-momentum basis, are evaluated numerically to determine the pressure equation of state, or other physical observable, such as the density equation of state or projected virial coefficients.

II.1 Path integral form of the partition function

The starting point of our lattice perturbative expansion of the equation of state is the grand-canonical partition function,

[TABLE]

where, as usual, $\hat{H}$ is the Hamiltonian, $\hat{N}$ is the particle number operator, $\beta$ the inverse temperature, and $\mu$ the chemical potential; the trace is over all the multiparticle states in the Fock space. We will assume that

[TABLE]

where $\hat{T}$ is the kinetic energy operator, and $\hat{V}$ is the potential energy operator. Rather than expanding directly in powers of $\hat{V}$, as is conventional in perturbation theory, we will introduce a HS transformation. We first discretize the imaginary time direction as $\beta=\tau N_{\tau}$, and apply a first-order Suzuki-Trotter factorization such that

[TABLE]

and next perform the HS transformation,

[TABLE]

where $\hat{W}[\sigma]$ is a one-body operator representing an external potential set by the auxiliary field $\sigma$, and the field integral is over all possible configurations of that field. The specific form of $\hat{W}[\sigma]$ depends on the choice of the HS transformation; in the case of fermions it can be either discrete or continuous, and in the latter case it can be compact or non-compact MCReview . For the purposes of this work, the specific form of the HS transform is immaterial, as we will undo the transformation order-by-order in the perturbative expansion.

Naturally, the use of the HS transform is much more common in Monte Carlo approaches than in perturbative ones, and it is the natural path to mean-field theory. Here, however, we will use it as a device to recover the results of Wick’s theorem but bypassing the operator algebra completely, as we show below. In addition, this approach has the advantage of facilitating the use of many-body forces, as it is relatively simple to write down HS transforms for them and the steps after the transform are quite mechanical (which is the main reason for our choice).

Collecting all the field integrals and performing the trace over the resulting product of exponentials of one-body operators yields a path-integral form of the grand-canonical partition function $\mathcal{Z}$,

[TABLE]

where in the above form we have assumed that the system contains two identical species and the matrix $M$ of dimension $N_{\tau}\times V$ by $N_{\tau}\times V$ (where $N_{\tau}$ is the temporal volume and $V$ is the spatial volume) encodes all relevant system dynamics, such that

[TABLE]

For the case of contact interactions, such as in the Gaudin-Yang model, the matrix ${\mathcal{U}}[\sigma]$ can be written as

[TABLE]

where ${\bf x}$ and ${\bf x}^{\prime}$ are spatial indices, $t$ is a temporal index, and where $A\equiv\sqrt{2(e^{\tau g}-1)}$, $g$ is the zero-range coupling, and $\tau$ is the temporal lattice spacing (see e.g. EoS1D ). We have chosen a realization of the HS transformation in which the continuous field $\sigma$ takes values between $-\pi$ and $\pi$ at each point in spacetime, such that

[TABLE]

This particular kind of HS transformation was first explored in Refs. HSLee1 ; HSLee2 .

Note that in the case of non-identical species, such as in polarized or mass imbalanced systems of two flavors, two different determinants are present in the partition function, i.e.

[TABLE]

At this point, we have eliminated from the problem all the quantum operators, and this is one of the main advantages of the method. Moreover, we have accomplished this by “integrating out” the fermionic degrees of freedom, which is an unorthodox route to perturbation theory. We will see, however, that this is a useful way to proceed in the sense that it mechanically generates the correct answers, and therefore it is amenable to automation.

II.2 Expanding the fermion determinant

To obtain the perturbative expansion of the grand-canonical partition function on the lattice, we expand the determinant of the matrix $M$ in powers of the dimensionless parameter $A$ defined above. To this end, first note that

[TABLE]

where $M_{0}$ is the non-interacting counterpart, and

[TABLE]

with $\left[T\right]_{{\bf x},{\bf x}^{\prime}}\equiv\left[e^{-\tau(\mathbf{p}^{2}/2-\mu)}\right]_{{\bf x},{\bf x}^{\prime}}$, (where we have taken $\hbar=k_{B}=m=1$), and

[TABLE]

Thus, we further obtain that $M=M_{0}(\openone+AK\mathcal{S}[\sigma])$, where $K\equiv M_{0}^{-1}\mathcal{T}$ will play the role of the free propagator, as we explain below, and of course

[TABLE]

At this point, we set aside the non-interacting factor and make use of the identity $\det=\exp\mathrm{tr}\ln$ to formally expand the logarithm in powers of $A$, such that

[TABLE]

where

[TABLE]

Note that each power of the fermion determinant translates simply as a prefactor in the definition of $E_{k}$. Thus, the expansion coefficients for an unpolarized system of $N_{f}$ species are simply $N_{f}E_{k}$. We further write the determinant as

[TABLE]

Eq. (15) begins to expose the powers of $A$ in the perturbative expansion. The task ahead is to isolate those powers and carry out the path integral exactly. Using these expressions, by fully expanding and combining terms of identical powers of $A$, it is straightforward to see that

[TABLE]

The above expansion corresponds to one of the spin degrees of freedom and therefore it is to be combined with a corresponding expansion for the other spin, which could be done for polarized (different $\mu$ for each species) or unpolarized (same $\mu$) cases. Here we focus on the latter with a two-component system, which simply means that we are expanding the square of the determinant. As mentioned above, this can be accomplished easily by setting $E_{k}\to N_{f}E_{k}$ for the general case of $N_{f}$ flavors.

II.3 Recovering Wick’s theorem by calculating the path integral exactly at each order

Once the determinants are expanded to the desired power of $A$, which is easy to automate, one obtains products of powers of the various $E_{k}$ defined in Eq. (15). The path integral of each of these products must be evaluated to obtain the expansion of the partition function. For instance, at second order, one of the terms is

[TABLE]

where all the (collective, spacetime) indices $i,j,k,m$ on the right-hand side are assumed to be summed over.

Writing out the indices explicitly at each order, we see that the main problem is computing path integrals of the form

[TABLE]

where each $b_{k}$ represents a spacetime coordinate. These integrals vanish if $m$ is odd, but otherwise the result is generally finite and positive. For instance, for $m=2$,

[TABLE]

and for $m=4$,

[TABLE]

where $3/8$ is the result for the case of all four coordinates coinciding, i.e.

[TABLE]

It is the tensor expressions like Eq. (II.3) that automatically implement Wick’s theorem when contracting with the various $K$ propagators. As with the Taylor expansion of the determinant, the calculation of the $I_{m}$ needed for each of the terms at a given order was also automated, as was the subsequent index contraction with the propagators. The result of that process is that not only the diagrams themselves, but also all the symmetry factors are generated automatically, minimizing the amount of manual bookkeeping. Naturally, the topology of the diagrams enters through the $I_{m}$ tensors, which encode the multiple ways in which the corresponding path integral can give a non-vanishing result.

The complexity of the expression for $I_{m}$ grows with each order in the perturbative expansion, and causes the number of terms in the expansion of $\mathcal{Z}$ (i.e. after contracting with the $K$’s) to grow very quickly. Although the next-to-leading (NLO) contribution can easily be verified by hand, the next-to-next-to-next-to-leading (N3LO) (i.e. $A^{6}$) order produces (naively) on the order of $10^{3}$ terms, all of which must be simplified to obtain the final results. Moreover, in the case of polarized systems, multiple products of the determinant must be expanded. This scaling underscores why the method we proposed here is well suited for automation but it is otherwise not ideal for manual calculation, especially if a high-order, finite-temperature expansion is the goal.

II.4 Transforming to frequency-momentum space on the lattice

Naturally, the matrix $K$ can be diagonalized in the momentum basis with a discrete Fourier transform, such that

[TABLE]

where $a,b$ are collective spacetime indices of the form $(t,x)$, and $q=(\omega,k)$ is a collective frequency-momentum index, with $\omega=(2n_{\omega}+1)/N_{\tau}$ being a fermionic Matsubara frequency ($n_{\omega}=1,\dots,N_{\tau}$), and $\mathbf{k}=2\mathbf{n_{k}}\pi/N_{x}$ ($\mathbf{n}_{\mathbf{k},i}=1,\dots,N_{x}$) a $d$-dimensional spatial wavevector. The free propagator $D_{0}$ is then

[TABLE]

while the Fourier transform matrices take the form

[TABLE]

By computing the path integral over the auxiliary field $\sigma$ for all the terms in the expansion of the determinant at a particular order, we are left with an object of the generic schematic form

[TABLE]

where the sum is taken over all free spacetime indices, and the object $I_{a_{1}\dots a_{2m}}$ results from the path integral over the interaction terms $\sin(\sigma_{i})$.

Upon inserting the Fourier representation of the propagator [c.f. Eq. (24)], we obtain a form that is described by a collection of indices in the frequency-momentum basis:

[TABLE]

where $\tilde{I}_{q_{1},\dots,q_{m}}$ results from appropriately contracting the various Fourier tensors $U$ and $U^{\dagger}$ with $I_{a_{1}\dots a_{2m}}$. The $\tilde{I}$ tensors represent momentum conservation laws for each specific term.

The advantage of going to frequency-momentum space is that $S$ can be obtained by performing $m$ frequency-momentum sums instead of the $2m$ spacetime sums. This optimization, however, is not enough; it is crucial to carry out the frequency sums analytically in order to have a numerically manageable expression at the end. We turn to those sums next.

II.5 Computing finite Matsubara frequency sums analytically: two tricks

In the evaluation of expressions at a given order $n$, we are faced with nested frequency sums that schematically look like the expression

[TABLE]

where we have left out sums over momenta, which are to be carried out afterwards, and the delta factor represents an energy-momentum conservation law that is derived from each diagram’s topology. Note that, here and below, we will use the delta notation to represent the discrete Kronecker symbol rather than the Dirac symbol (we have no need for Dirac deltas here, as all our expressions are discrete). In this section, we show how we performed the frequency sums in a way that does not use complex contour integration and, moreover, allows us to treat all the sums simultaneously.

We begin with an example. The simplest case is that of a single factor, which may seem trivial but is nevertheless instructive:

[TABLE]

where we have encoded all the momentum and chemical potential dependence in the factor $Q$, and where $\omega_{n}=(2n+1)\pi/N_{\tau}$. In this case we expand the expression inside the sum as a geometric series:

[TABLE]

where we used the fact that $k\geq 0$ and

[TABLE]

The $S_{1}$ sum is useful per se, but also because high-order calculations need the more general sums, such as

[TABLE]

such that $S_{1}^{(0)}=S_{1}$. For general $k$, this is easy to compute by introducing a new parameter $\lambda$ via

[TABLE]

and differentiate with respect to $\lambda$ as needed,

[TABLE]

Clearly, $S_{1}(\lambda)$ is a generating function for frequency sums of higher order. In particular, for instance,

[TABLE]

The method of expanding the numerator as a power series applies as well when multiple sums are present, as we show in Appendix A with an example that corresponds to the “beach-ball” diagram of Fig. 4(b).

The general procedure of the method is as follows. First, take a step backward and use the Fourier-sum representation of all the frequency Kronecker deltas, i.e.

[TABLE]

in combination with the geometric series expression of the denominators. This first step is similar to that of techniques used in the continuum, where an integral representation of a Dirac delta function is utilized. Beyond this point, however, the calculations differ considerably.

The second step is to sum over each fermionic frequency $\omega_{n}$, which yields as many delta functions as denominators, with their corresponding “boundary” sums and factors of $(-1)^{m}$, i.e. use

[TABLE]

Here it is important to keep the terms for all $m$, as $y$ will generally vary over a semi-infinite range due to the geometric series expansion.

Third, sum over the geometric sum index to saturate the delta functions generated in the previous step. To this end, it is useful to extend the geometric sum index to $-\infty$, such that Heaviside step functions should be inserted accordingly.

Fourth, implement the constraints over the boundary sums and evaluate such sums, which lead to the expected $(1+Q_{k}^{N_{\tau}})^{-1}$ factors. Finally, sum over the Fourier index of the energy-conserving delta functions. As is clear from the above, the operations to be performed are rather elementary and do not involve complex analysis, only a small amount of bookkeeping. It is the order of the operations that is crucial for the simplicity of the method.

III Complex Langevin Approach

The second method we present in this work involves the use of CL dynamics in the context of the lattice quantum Monte Carlo technique, which is, of course, a nonperturbative approach. Recent developments have shed light on several aspects of the CL method, which had previously hindered progress but which are now starting to be better understood and in some cases even resolved. In this section we discuss our implementation of CL dynamics, including a kind of dynamical stabilization that is reminiscent of the one put forward in Refs. DynamicStabilisation1 ; DynamicStabilisation2 .

The starting point is once again the partition function

[TABLE]

where one normally identifies $P[\sigma]={\det}^{2}M[\sigma]$ as the (unnormalized) probability measure to be used in a Metropolis-based Monte Carlo calculation. One may define an effective action $S$ via

[TABLE]

Observables, at least simple ones, are then shown to take the form

[TABLE]

such that the expectation value can be determined by sampling the auxiliary field $\sigma$ according to $P[\sigma]$. Two well known ways of carrying out that sampling, without the presence of a sign problem, are the hybrid Monte Carlo algorithm (HMC) (see Refs. HMC1 ; HMC2 ) and the real Langevin method (RL) (also known simply as stochastic quantization; see Refs. RL1 ; RL2 ; RL3 ).

In HMC, one defines an auxiliary field variable $\pi$ conjugate to the HS field $\sigma$ along with global (molecular-dynamics type) equations of motion in a fictitious phase-space time $t$ such that

[TABLE]

where the right-hand side of the last equation naturally represents a molecular dynamics force. These equations are usually integrated via the leapfrog method (or more sophisticated variants) to ensure reversibility, which in turn is essential for maintaining detailed balance. The global updates that this method allows are essential for lattice QCD calculations. The new field obtained at the end of the trajectory is accepted or rejected using the criteria of a Metropolis step.

In RL, on the other hand, there is no Metropolis step and the equations of motion are different:

[TABLE]

where we note that there is no auxiliary momentum field $\pi$ but a $t$-dependent noise field $\eta$ appears instead. The latter satisfies $\langle\eta(x,\tau)\rangle=0$ and $\langle\eta(x,\tau)\eta(x^{\prime},\tau^{\prime})\rangle=2\delta_{x,x^{\prime}}\delta_{\tau,\tau^{\prime}}$.

The (conventional) mathematical underpinnings of HMC and RL depend on $P[\sigma]$ being positive semi-definite. As is well known, this property generally fails to hold e.g. for repulsive interactions, polarized systems, etc. as mentioned above, and the calculation is then said to have a sign problem (or more generally a complex phase problem). In the case of HMC, this means that the Metropolis step is simply no longer well-defined and thus the algorithm is no longer available. For RL, on the other hand, a generalization is possible, namely CL, as first noted in Ref. Parisi1983 . In CL, one complexifies the HS field $\sigma$ into

[TABLE]

where $\sigma_{R}$ and $\sigma_{I}$ are both real fields. The CL equations of motion are given by

[TABLE]

where now $S[\sigma]$ is to be understood as a complex function of the complex variable $\sigma$. Note that, when the action is real, the imaginary part of the force is zero and then CL reduces to RL.

Under certain conditions, which have lately received much attention (see CLSufficientConditions00 ; CLSufficientConditions01 ; CLSufficientConditions02 ; CLSufficientConditions03 ), the CL method can be guaranteed to converge to the right answer. In those cases, expectation values $\langle\mathcal{O}\rangle$ are correctly obtained by averaging over the real part of $\mathcal{O}[\sigma]$ with complex fields $\sigma$ sampled throughout the CL dynamics evolution.

In the process of making CL a viable solution to the sign problem in lattice QCD, two crucial challenges were identified: the appearance of numerical instabilities in the CL evolution and the uncontrolled excursions of $\sigma$ into the complex plane due to singularities in $P[\sigma]$. The former was largely resolved by implementing adaptive time-step solvers CLAdaptiveSolvers ; the latter, on the other hand, is currently under investigation and a few approaches have been proposed (see e.g. RCL1 ; RCL2 ). In the calculations presented here, those excursions are highly problematic because the dependence of $P[\sigma]$ on $\sigma_{I}$ is through hyperbolic functions; indeed, the HS transformation we use depends on $\sin\sigma$, and

[TABLE]

Thus, a growing (positive or negative) $\sigma_{I}$ effectively increases the coupling at an exponential rate, which can completely stall the calculation or result in a converged but wrong answer (as we have observed in our tests). This exponential growth is similar to the problem found in gauge theories, as the complexified link variables representing the gauge field become unbounded in the same fashion in those theories.

III.1 Modified action

To overcome the problem of large excursions in the complex plane in the CL algorithm, we modified the action in a way reminiscent of the dynamical stabilization approach of Ref. DynamicStabilisation1 ; DynamicStabilisation2 . In the latter, a new term was added to the CL dynamics which vanishes in the continuum limit and renders the calculation stable. We propose here to modify the CL equations by adding a regulating term controlled by a parameter $\xi$, such that the new equations, in their discretized form, are

[TABLE]

This change amounts to modifying the action by

[TABLE]

The rationale for adding such a regulating term was originally based on our understanding of ${\delta S[\sigma]}/{\delta\sigma}$ as a molecular dynamics force in the context of the HMC algorithm (although it is typically referred to as the “drift” in the context of the CL method). In the HMC sense, the new term in the action can be understood as a harmonic oscillator trapping potential, i.e. a restoring force that prevents the field from running away. However, HMC does not apply when $\sigma$ is complex. A more appropriate interpretation is obtained by keeping only those new terms in the CL equations and neglecting the rest, which results in the decoupled form

[TABLE]

whose solution is a decaying exponential (assuming $\xi>0$) for both $\sigma_{R}$ and $\sigma_{I}$. Thus, this new ad hoc term represents a damping force.

Naturally, the proposed modification introduces a systematic effect that needs to be studied for each quantity of interest, i.e. it is crucial to understand the $\xi$ dependence of the output. The results presented below correspond to $\xi=0.1$. To gain specific insight into the variations of the density with $\xi$, we show in Fig. 1 a plot of the running average of the density as a function of the Langevin time $t$ for several values of $\xi$ in the neighborhood of [math]. As evident in that figure, there is a sizable window of small values of $\xi$ where CL converges. Additionally, we show in Fig. 2 a plot of the phase and magnitude of $e^{-S}$, where $S$ is the complex action, for different values of $\xi$. The effect of the $\xi$ term on the CL dynamics can be clearly seen in that figure (note the change of scale for the $\xi<0$ plot, where the results are expected to diverge).

In the next section we present our results. Whenever RL or CL results are shown, we have kept $\xi=0.1$. The agreement of RL with HMC (where the regulating term should not be needed) and the weak dependence of the density on $\xi$ (at least within a window close to $\xi=0$) support the idea that the damping term has negligible impact on the results within the precision studied here.

IV Results

In this section we show the results of implementing the above two formalisms and methods to the case of unpolarized, spin-$1/2$ fermions in one spatial dimension. The specific Hamiltonian we explore is the Gaudin-Yang form GaudinYang1 ; GaudinYang2 where $\hat{H}=\hat{T}+\hat{V}$ with

[TABLE]

and

[TABLE]

where $\hat{\psi}_{s}^{\dagger},\hat{\psi}_{s}$ are the creation and annihilation operators in coordinate space for particles of spin $s$, and $\hat{n}_{s}=\hat{\psi}_{s}^{\dagger}\hat{\psi}_{s}$ are the corresponding density operators. There are only two dimensionless parameters characterizing this problem at finite temperature: the fugacity $z=e^{\beta\mu}$ and the coupling $\lambda=\sqrt{\beta}g$.

In all of our results below, perturbative as well as non-perturbative, the system was placed on a lattices of spatial size $N_{x}=60,80,100$ and temporal size $N_{\tau}=160$. The spatial lattice spacing was set to $1$, thus setting the length scale for the problem, and the temporal lattice spacing was chosen as $\tau=0.05$ (in units of the spatial spacing). Previous studies show that, with those parameters (particularly with $N_{x}\geq 60$), the equations of state of this system are within $\leq 5\%$ of the continuum limit result for the values of $\lambda$ studied here.

IV.1 Analytic expressions for the perturbative expansion

Although the results shown here correspond to unpolarized fermions with two internal degrees of freedom (e.g., spin-1/2 particles), the perturbative and complex Langevin formalisms can easily be extended to greater degrees of freedom (or, “flavors”) by introducing additional determinants in the integrand of the path integral. Additionally, although all observables computed here correspond to the one-dimensional system, the perturbative expansion itself is independent of the spatial dimension; one must simply modify the definition of the wavevector $\mathbf{k}$ to extend the system to higher dimension.

Among the main results of this work are the final expressions for the perturbative contributions to the partition function $\mathcal{Z}$ at each order $n$ up to next-to-next-to-next-to-leading order (N3LO), which are shown in Table 1 (explained in detail below). Each order $n$ is to be accompanied by a factor of $A^{2n}$, but all the information required to reproduce the numerical results in this work are otherwise presented in that table and in the appendices.

Each contribution to the perturbative expansion of the partition function corresponds to a fully connected Feynman diagram of $n$ vertices (or a product of two or more such diagrams), and can be described by the product of some scalar prefactor with a $1/(N_{x}N_{\tau})^{n}$ dependence, and one or more sums over the complete momentum-frequency basis. The Feynman diagrams that appear at NLO, N2LO, and N3LO are provided in Figs. 3, 4, 5, respectively. The analytic expressions for each diagram are provided in Table 1, where the Matsubara-frequency sums have already been carried out using the technique outlined in Sec. II.5; therefore all sums and indices that remain in these expressions are over spatial momenta. The frequency sums always contain $Q_{n}$ as parameters that encode the spatial momentum dependence:

[TABLE]

A variety of frequency sums for sums over various products of propagators and momentum conservation conditions appear in these expressions. The expressions for $S_{1}$ and $S_{1}^{(1)}$ are given earlier in the text by Eqs. (33) and (38), respectively, $S_{4}$ is derived in Appendix A, and all others are given in Appendix B.

The number of loops, or the number of nested momentum sums that appear for each expression, is also provided explicitly in Table 1. Since the momentum sums must be computed in full, this number provides an estimate for the computational scaling that is required to compute the value of each diagram; the number of terms that must be computed scales as $N_{x}^{\ell d}$, where $\ell$ is the number of loops and $d$ is the spatial dimension.

IV.2 Pressure equation of state via perturbation theory

To determine the pressure $P$, we use the perturbative expansion of the interacting partition function $\mathcal{Z}$, which we expand to order $2n$ in the parameter $A$ and write as

[TABLE]

where $\Delta_{n}$ is the N$n$LO contribution, as given by the expressions in Table 1, and where $\mathcal{Z}_{0}$ is the non-interacting grand-canonical partition function,

[TABLE]

with $\epsilon_{k}=(2\pi k/N_{x})^{2}/2$.

Thus, the interacting pressure $P$, relative to the non-interacting result $P_{0}$, becomes

[TABLE]

To keep the computed value of $\ln\mathcal{Z}$ consistent with the highest order of $A$ in expansion of $\mathcal{Z}$, we expand the numerator $\ln(1+\sum_{n}A^{2n}\Delta_{n})$ in a Taylor series about $A=0$, such that the expanded form up to N3LO is

[TABLE]

where

[TABLE]

Note that while this procedure offers consistency in the order of the coupling at all stages of the calculation, in performing this expansion it is important to be mindful of the validity of choosing to expand about the non-interacting limit. Since the partition function is an extensive quantity, this expansion may yield unphysical results in cases where $\sum_{n}A^{2n}\Delta_{n}\gg 1$, which may occur in cases where stronger effective interactions are present, as is the case when $N_{f}>2$. However, in all cases shown here, the observables demonstrate physical behavior and are in agreement with alternative methods where available.

Using Eq. (62), we computed the pressure at NLO, N2LO, and N3LO, as a function of the dimensionless parameter $\beta\mu=\ln z$. Our results for this quantity, for the case of attractive interactions, are shown for a variety of interaction strengths in Fig. 6 (top). Remarkably, we see evidence of convergence for $\beta\mu\geq 2$, even for $\lambda=2$. For $\beta\mu\leq 2$, on the other hand, our results are qualitatively correct but fail to match the Monte Carlo answers by roughly $10\%$ in the worst case of $\lambda=2$.

A possible explanation of this behavior is that perturbation theory simply fails to capture the effect of the two-body bound state on the virial coefficients as $\lambda$ is increased, which in the case of $b_{2}$ involves a dominant inverse gaussian contribution $\propto\exp(\lambda^{2}/4)$. This effect could be absorbed into the definition of $\lambda$ by way of a new renormalization scheme defined by matching to an exact calculation of $b_{2}$ (lattice or continuum), instead of identifying the coupling $g$ with the inverse scattering length. Clearly, such a renormalization scheme would improve the agreement with the Monte Carlo results in the semiclassical region of negative $\beta\mu$, but by the same token it would spoil it for $\beta\mu\geq 2$. Nevertheless, as the perturbative expansion is extended beyond third-order, the convergence properties of the regime is expected to improve.

In all cases, the NLO results are quantitatively disappointing, but N2LO and N3LO display substantial improvement both in convergence and in approaching the Monte Carlo results. Naturally, this improvement is not free: the computational effort increases dramatically for N2LO and N3LO relative to NLO. Similar behavior is seen for the repulsive case in the bottom of Fig. 6 but with an important difference: the results oscillate as the perturbative order is increased, whereas in the attractive case convergence appears to be monotonic.

IV.3 Density equation of state via perturbation theory and complex Langevin

Based on our knowledge of the pressure, differentiation with respect to $\beta\mu$ gives access to the density, for which results from perturbation theory can be compared more directly with quantum Monte Carlo results (the pressure is not typically a quantity computed directly in Monte Carlo calculations, but is rather obtained by integrating the density; see Ref. EoS1D for details). We show such a comparison in Fig. 7 (top) for the case of attractive interactions where Monte Carlo calculations are possible without a sign problem. For repulsive interactions, our calculations yield the results shown in Fig. 7 (bottom). In both cases, we have numerically differentiated with a point spacing of $d\mu=0.01$.

Both figures also display density results obtained using CL techniques. Remarkably, this method demonstrates excellent agreement with both hybrid Monte Carlo and perturbation theory for attractive interactions (where CL becomes RL), and perturbative results for repulsive interactions. The auxiliary field $\sigma$ was evolved with CL dynamics for $10^{5}$ iterations with an adaptive timestep $\delta t$, adjusted by monitoring the magnitude of the CL drift.

IV.4 Perturbative virial coefficients

Based on the results shown above, we implement a method of particle-number projection, which is well known in the area of nuclear physics, to calculate virial coefficients perturbatively. To explain the method, we recall that the virial expansion for the pressure is given by

[TABLE]

where $\mathcal{Q}_{1}=N_{f}L/\lambda_{T}$ is the one-particle canonical partition function (here $\lambda_{T}=\sqrt{2\pi\beta}$ is the thermal de Broglie wavelength), and $b_{n}$ are the virial coefficients we want to extract.

Using our (semi-)analytic expressions for $\ln\mathcal{Z}(z)$, it is possible to extract the coefficients $b_{n}$ via a numerical Fourier projection,

[TABLE]

where the partition function is evaluated at $z=\alpha\exp(i\phi)$ and $\alpha$ is an arbitrary real parameter that is independent of the result for $b_{n}$, but is used to ensure that a well-behaved integrand is used. A value of $\alpha=0.01$ is well-suited for the evaluations performed here.

The perturbative expansion of virial coefficients $b_{n}$ displays signs of convergence at weak coupling $|\lambda|\leq 1$ and for low enough virial order $n$. We find that for $\lambda=\pm 1$ the predictions for $b_{5}$ cease to converge, at least at N3LO. For comparison, an exact lattice calculation of the second-order virial coefficient at $N_{x}=100$ and $\beta=8$ under attractive interactions (see Ref. EoS1D ) shows that for $\lambda=0.5$, $b_{2}=-0.230$, and for $\lambda=1.0$, $b_{2}=-0.0351$. This discrepancy compared with the perturbative value for $b_{2}$, particularly for larger values of the coupling, underscores the difficulty in quantitatively capturing the behavior in the virial region using perturbation theory.

V Summary and Conclusions

In summary, we have determined, using both up to third order in lattice perturbation theory and CL techniques, the pressure and density equations of state of non-relativistic spin-$1/2$ fermions with contact interactions in both attractive and repulsive regimes. We have used those results in combination with particle-number projection techniques to obtain a perturbative expansion for virial coefficients.

Our results for the perturbative pressure and density display some remarkable features: even for couplings as high as $\lambda=2$ there is a region of fugacity values (namely for $\beta\mu\geq 1$) that displays clear signs of convergence toward excellent agreement with extant non-perturbative Monte Carlo results (for the case of attractive interactions).

We performed our perturbative lattice calculations by following an unconventional route based on first decoupling the interaction via a Hubbard-Stratonovich transformation, and then undoing that transformation after expanding the non-interacting fermion determinant under the path integral. The approach is interesting because it is amenable to automation in a way that seems easier to handle than conventional approaches that use operators throughout.

Additionally, we have presented a way to carry out the Matsubara-frequency sums that is both simple and efficient, in particular in that it does not involve contour integration in the complex plane and in that it treats all the nested sums simultaneously. While techniques that are similar in spirit exist in the continuum, we are not aware of this type of approach for systems on the lattice. Carrying out those frequency sums is more than a matter of convenience: it is essential in order to obtain numerically manageable expressions, namely expressions that only contain spatial momentum sums. Furthermore, we have shown how a common trick, namely differentiation with respect to an auxiliary parameter, can generate expressions needed for certain diagrams using expressions for lower-order diagrams. The frequency sums presented above and in the appendices are universal in that they apply to systems regardless of the dispersion relation, external potential, and interaction potential, as long as the latter is energy-independent.

To access repulsive interactions in a non-perturbative fashion (a well-known case suffering from the sign problem), we implemented the CL method and proposed a modification to the action that prevents uncontrolled excursions into the complex plane. This modified action yields a damping term in the CL dynamics which, in turn, introduces systematic effects of a priori unknown size. Our investigations of those effects, both on the attractive and repulsive sides, support the idea that the impact on the density equation of state is within the uncertainties of our calculations; other observables would require independent investigations. The effectiveness of this approach in other situations, like polarized matter and higher dimensions, remains to be studied; however, nothing seems to prevent us from using the same idea in those cases. For the 1D system studied here, the CL results agree with the perturbative ones in regions where perturbation theory appears to converge (within the orders studied here); everywhere else, the CL answers represent non-perturbative predictions.

Acknowledgements.

We thank G. Aarts, J. Braun, and A. Vuorinen for comments on the manuscript. This material is based upon work supported by the National Science Foundation under Grants No. DGE1144081 (Graduate Research Fellowship Program), PHY1306520 (Nuclear Theory Program), and PHY1452635 (Computational Physics Program).

Appendix A Matsubara frequency sum derivation example of $S_{4}$

To further illustrate our method of analytically evaluating Matsubara frequency sums on the lattice, we will show in detail how to calculate the case of $S_{4}$. Note that higher-order frequency sums, such as $S_{4}^{(1)}$, can also be determined from the expressions derived here by inserting a parameter $\lambda$ such that, for example, $Q_{1}\rightarrow Q_{1}/\lambda$, and proceed to calculate derivatives with respect to $\lambda$, evaluated at $\lambda=1$. This is a generalization of the procedure used for the case of $S_{1}^{(1)}$, as illustrated in Eq. (36).

The frequency sum $S_{4}$ is defined as the sum over four free propagators, where the indices of the first three propagators are free, and the fourth is constrained by the momentum conservation condition $\omega_{4}=\omega_{1}-\omega_{2}+\omega_{3}$. Therefore, we write the expression for $S_{4}$ with a Kronecker delta such that

[TABLE]

We first represent the Kronecker delta with a sum over the complete frequency-space basis and expand the denominators using a geometric series,

[TABLE]

Next, we reorder the sums in preparation to carry out all the frequency sums, which result in delta functions. Note that in the process, the appearance of powers of $(-1)$, which results from the fermionic Matsubara frequencies.

[TABLE]

To saturate the delta functions when summing over $\{k\}$, we must account for the fact that all $k$ are positive. To this end, we extend the sums to negative values of $k$ and insert Heaviside functions accordingly (defined to saturate for non-negative values of the argument), which results in

[TABLE]

Because of the range of $p$, the $\theta$ functions imply the following conditions: $m_{1}\geq 0$ if $p=0$, but otherwise $m_{1}>0$; an identical condition for $m_{3}$; $m_{2}\geq 0$; and $m_{4}\geq 0$. Using the last two conditions, we obtain a first simplification:

[TABLE]

Applying the remaining two conditions amounts to (note we separate the $p=0$ and $p\neq 0$ cases in the second line)

[TABLE]

where $X={Q_{2}Q_{4}}/({Q_{1}Q_{3}})$. Thus, the final result is simply

[TABLE]

Note that, as anticipated, we obtained this result without resorting to contour integration, and moreover we addressed all the frequency sums simultaneously.

Appendix B Listing of Higher-Order Matsubara Frequency Sums

In this section, we will provide a listing of all Matsubara frequency sums that appear in the perturbative contributions to the grand-canonical partition function at third-order. Note that the first- and second-order sums $S_{1}$, $S_{1}^{(1)}$, and $S_{4}$ are provided elsewhere in the text. In all cases, our expressions were checked by comparing with a direct evaluation of the defining sums in a small spacetime volume.

The frequency sum $S_{1}^{(2)}$ is defined as the sum over the product of three identical propagators, and the derived expression is given by

[TABLE]

The frequency sum $S_{4}^{(1)}$ can be obtained, as explained in the main text, by using a differentiation trick:

[TABLE]

where $S_{4}(\lambda)\equiv\lambda^{-1}\left.S_{4}\right|_{Q_{1}\to Q_{1}/\lambda}$. Thus,

[TABLE]

where $X\equiv Q_{2}Q_{4}/(Q_{1}Q_{3})$.

Finally, the frequency sum $S_{6}$ is defined as

[TABLE]

such that

[TABLE]

[TABLE]

where we also define

[TABLE]

as well as $X\equiv Q_{1}Q_{6}/(Q_{2}Q_{5})$ and $Y\equiv Q_{3}Q_{5}/(Q_{4}Q_{6})$. As evident from the above expressions, it is important to consider limiting cases near poles (e.g. $X=1$) when numerically evaluating these sums; such limits are straightforward to compute and implement.