Finite temperature auxiliary field-based quantum Monte Carlo methods, including determinant quantum Monte Carlo and Auxiliary Field Quantum Monte Carlo (AFQMC), have historically assumed pivotal roles in the investigation of the finite temperature phase diagrams of a wide variety of multidimensional lattice models and materials. Despite their utility, however, these techniques are typically formulated in the grand canonical ensemble, which makes them difficult to apply to condensates such as superfluids and difficult to benchmark against alternative methods that are formulated in the canonical ensemble. Working in the grand canonical ensemble is furthermore accompanied by the increased overhead associated with having to determine the chemical potentials that produce desired fillings. Given this backdrop, in this work, we present a new recursive approach for performing AFQMC simulations in the canonical ensemble that does not require knowledge of chemical potentials. To derive this approach, we exploit the convenient fact that AFQMC solves the many-body problem by decoupling many-body propagators into integrals over one-body problems to which non-interacting theories can be applied. We benchmark the accuracy of our technique on illustrative Bose and Fermi–Hubbard models and demonstrate that it can converge more quickly to the ground state than grand canonical AFQMC simulations. We believe that our novel use of HS-transformed operators to implement algorithms originally derived for non-interacting systems will motivate the development of a variety of other methods and anticipate that our technique will enable direct performance comparisons against other many-body approaches formulated in the canonical ensemble.

For decades, chemical physicists have focused on developing a constellation of electronic structure methods, from mean field1,2 to perturbation1,2 to coupled cluster theories,2,3 to describe the ground and low-lying excited states of molecules and solids. While it is undoubtedly true that many phenomena involve electrons that reside in their ground states, it is becoming increasingly apparent that there are a wealth of phenomena for which this assumption does not hold: atoms and molecules in the centers of large planets and stars can experience GPa of pressure and temperatures of over 10 000 K,4,5 and lasers can be used to heat and thereby steer chemical reactions along different mechanistic pathways to facilitate processes such as catalysis.6–8 Temperature is moreover one of the key parameters that can be used to tune the electronic properties of the materials that make up many of modern society’s most important technologies9–11 and is responsible for the pernicious loss of coherence in physical realizations of qubits.12 

Reflecting this growing list of applications, a growing number of methods have recently been developed to study them. One of the most straightforward ways of determining the finite temperature properties of quantum systems is by diagonalizing the system’s Hamiltonian to determine all of its many eigenvalues and weighting them to compute its partition function and related observables in a process termed exact diagonalization (ED).13 While this technique is numerically exact, as its name implies, its computational cost scales exponentially with system size. On the other end of the spectrum, finite temperature mean field theories, including finite temperature Hartree–Fock14,15 and Density Functional (DFT) theories,16,17 trade accuracy for computational expediency by approximating the many-body problem as a one-body problem in which an electron is coupled to an average or “mean” field of the other electrons. While such mean field theories continue to improve and hold great promise, they still struggle to achieve the accuracy often needed to correctly capture many-body physics. Between these two poles lie a variety of techniques that are generalizations of their ground state counterparts. The past few years, for example, have seen a resurgence of interest in finite temperature coupled cluster techniques18–22 and perturbation theories.23–28 Closer to the mean field end of the spectrum, finite temperature embedding theories that partition systems into correlated regions embedded within uncorrelated baths29 such as Dynamical Mean Field Theory (DMFT)30,31 and the finite temperature SEET and Second-order Green’s Function methods32–34 are distinctively capable of directly obtaining the full frequency-dependent spectra of the systems they model. Nevertheless, all of these methods struggle to balance computational cost with the need to account for the numerous electronic states that contribute to finite temperature expectation values.

Finite Temperature Quantum Monte Carlo (FT-QMC) methods are particularly advantageous in this regard because they have the exceptional ability to access numerous states without the exponential or high polynomial costs of other techniques by randomly sampling multidimensional state spaces for the most important states.35–38 One of the most successful of these techniques is Determinant Quantum Monte Carlo (DQMC)39–41 and its more recent extension, FT Auxiliary Field Quantum Monte Carlo (FT-AFQMC).42–44 DQMC has a long history of being employed to study a wide variety of lattice models. One of the key reasons why DQMC has been so widely adopted is because, for certain important classes of problems,39,45 it does not suffer from the infamous sign or phase problems46 that limit the practical utility of many other stochastic methods. AFQMC methods grew out of DQMC methods and were developed with the aim of being able to accurately model systems with clear sign problems by leveraging their more sophisticated sampling techniques and phaseless approximation.38,47–49 These methods have since shed light on such sign- and phaseful systems as the Hubbard model off of half-filling46,50,51 and many different ab initio molecular44,52–54 and solid state systems.55–58 Furthermore, these QMC techniques are versatile—they can be applied to virtually any second-quantized Hamiltonian—and, by sampling the overcomplete space of non-orthogonal determinants, they are able to sample large portions of Fock space with a high accuracy yet comparatively mild computational cost of O(N3)–O(N4), where N denotes the number of basis functions.38 

This said, one of the Janus-faced features of these methods is that they are formulated in the grand canonical ensemble, in which a system’s internal energy and particle number are allowed to fluctuate according to its fixed temperature and chemical potential.59 In many systems, it is more natural to fix intensive properties rather than extensive quantities (imagine the inherent difficulty involved with fixing the number of electrons in a solid), and auxiliary field-based methods formulated in the grand canonical ensemble may therefore be viewed as being better suited for modeling these systems. Nevertheless, it can be challenging to ascribe a meaningful chemical potential to systems such as condensates, superfluids, and superconductors that, by definition,60 can undergo large fluctuations in their particle numbers.61 Indeed, for this very reason, past attempts at using grand canonical FT-AFQMC to model bosons and Bose–Fermi mixtures were largely unsuccessful.62 However, even for the majority of systems for which the chemical potential is well-defined, the process of determining the correct chemical potential to achieve a desired filling or average particle number can be unwieldy. In current auxiliary field techniques that work in the grand canonical ensemble, practitioners must scan through tens to hundreds of possible chemical potentials before recovering the desired one for every Hamiltonian they study, which can accumulate into a substantial cost. Moreover, sampling in the grand canonical ensemble is not equivalent to sampling in the canonical ensemble, particularly at the particle numbers far from the thermodynamic limit used in many simulations. Sampling the grand canonical ensemble involves sampling a much larger Fock space of states, many of which minimally contribute to thermodynamic averages, than sampling the canonical ensemble. This makes comparing results from grand canonical ensemble approaches to results from canonical ensemble approaches challenging and can also introduce additional statistical noise that can make arriving at meaningful statistical averages more challenging than in the canonical ensemble. Because of the differences between these ensembles, properties measured in these ensembles may also converge to their limits in different manors, meaning that one ensemble may converge certain properties more rapidly than the other.

Given this context, in this work, we derive a new recursive formulation for performing finite temperature AFQMC simulations in the canonical ensemble. Unlike past canonical ensemble formalisms that relied upon Fourier extracting canonical results from grand canonical simulations and thus also depended upon a costly tuning of chemical potentials,63–67 our technique does not require knowledge of the chemical potential. To accomplish this, our method exploits one of the key yet often overlooked features of the AFQMC formalism: by decoupling two-body operators into an integral over one-body operators, the Hubbard–Stratonovich (HS) transformation38,68,69 used in AFQMC produces an ensemble of non-interacting systems to which theories developed for non-interacting systems can be applied. In particular, to derive our approach, we apply a recursive formalism for obtaining the partition functions of ideal gases to the one-body operators in AFQMC to ultimately yield a many-body recursive theory. In the following, we derive this formalism for systems of bosons, spinless fermions, and spinful fermions, and demonstrate its accuracy and practical advantages for several benchmark Hubbard models. Interestingly, we show that energies computed in the canonical ensemble converge more quickly to the ground state than energies computed in the grand canonical ensemble. Ultimately, we believe that this formalism will be most useful for studying condensates as well as systems of fermions approaching their ground states and will motivate the integration of other useful non-interacting theories into the AFQMC formalism.

We thus begin in Sec. II by deriving our formalism for bosons and fermions and contrasting it with the conventional grand canonical DQMC and AFQMC formalism. We then benchmark the accuracy and highlight interesting features of our method in Sec. III. Finally, in Sec. IV, we conclude with a discussion of the applications for which we expect our methodology to be of the greatest future use.

To contrast our canonical ensemble formalism with previous formalism, we begin by outlining the conventional grand canonical DQMC algorithm before describing our own formalism.

The central quantity in any grand canonical ensemble formalism is the grand partition function, Z, because it is from this quantity that all other thermodynamic quantities can be obtained. In the grand canonical ensemble, the internal energy and particle number are allowed to fluctuate around average values that can be tuned by the temperature and chemical potential, respectively.59 The grand partition function may thus be expressed as the trace over the Boltzmann factor containing both the temperature and the chemical potential.

In DQMC39,41 and AFQMC48,62 methods that work in the grand canonical ensemble, the grand partition function may be re-expressed into a form amenable to sampling by first discretizing it into L imaginary time slices,

(1)

Here, Ĥ denotes the many-body Hamiltonian, μ denotes the chemical potential, N^=i,σĉiσĉiσ denotes the particle number operator, and Δτ = β/L. In order to facilitate its subsequent evaluation, the propagator, eΔτ(ĤμN^), is factored into short imaginary time kinetic and potential propagators via a Suzuki–Trotter factorization70,71 such that

(2)

As discussed in Sec. III A, in this work, we focus on models whose Hamiltonians can be written as Ĥ=K^+V^, where K^ denotes the collection of all one-body operators, including the chemical potential term, and V^ denotes that of all two-body operators. The exact grand partition function is recovered in the limit that Δτ → 0. This factorization enables us to treat the one-body and two-body operators separately. While one-body propagators, eΔτK^, may be neatly expressed as matrices in a given basis,40 two-body propagators, eΔτV^, may not be as easily expressed. Fortunately, two-body propagators of the form eΔτV^ can be re-expressed in terms of one-body propagators (i.e., linearized).68 In general, two-body operators such as V^ may be decomposed into combinations of squares of one-body operators,72 

(3)

where v^γ denote linear combinations of (and are therefore also linear) one-body operators and λγ weight each of those operators’ contributions to the potential operator. Using such decompositions, these propagators may then be re-expressed as one-body propagators using the continuous Hubbard–Stratonovich (HS) transformation69 

(4)

where ϕ is the so-called “auxiliary field.” As eϕ2/2/2π is a Gaussian, it may be viewed as a probability and the overall transformation may be viewed as re-expressing a two-body propagator into an integral over effective, one-body mean field propagators parameterized by different auxiliary fields. Note that, while this continuous transformation is more general and can be used to decouple a wide variety of fermion and boson Hamiltonians,38,69,73 it is often more efficient when treating Fermi–Hubbard models to employ the discrete version of the HS transformation.68,69 Nevertheless, because it is more general, we will proceed to base our remaining derivations off of the continuous HS transformation.

Substituting Eq. (4) into Eq. (1), the grand partition function may finally be written in a form that can be sampled,

(5)
(6)
(7)

as was our initial goal. In the above equation, ϕl{ϕlγ}, the vector of all auxiliary fields, p(ϕl), is the multivariate Gaussian distribution formed from the product of each field’s Gaussian distribution, and

(8)

One important advantage of working in the grand canonical ensemble is that, by taking the trace over fermions explicitly, the numerical complexities of evaluating Tr(lLB^(ϕl)) are avoided, and we obtain

(9)

where taking the trace over the product of B^(ϕl) leads to a product of determinants, one from each of the two spin sectors;40 Bσ(ϕl) is the matrix representation of B^σ(ϕl) in the single-particle basis; and I is the identity matrix. The grand partition function may thus be expressed as a multidimensional integral over auxiliary field space, which can be efficiently computed via Monte Carlo sampling. It is based on auxiliary field configurations sampled according to Z that observables such as energies and site/orbital occupancies may be determined. While this formalism has proven remarkably useful, it can be exceedingly costly to repeatedly determine the correct chemical potential at every temperature and system size one wants to model.

One can modify the grand canonical ensemble formalism described above to compute canonical ensemble quantities by replacing the Boltzmann factor containing the chemical potential term in Eq. (1) with eβĤ and taking the trace in Eq. (7) over all states in the canonical ensemble instead of in the grand canonical ensemble. Nevertheless, evaluating the trace over canonical ensemble states is far more complicated than evaluating the trace over grand canonical ensemble states as the particle number must be fixed.

This said, several approaches have been advanced over the years to transform sampling of the grand canonical ensemble into sampling of the canonical ensemble without explicitly taking the trace over canonical ensemble states. One technique that takes advantage of grand canonical treatments while still fixing particle numbers is the particle projection method.63–66 In this technique, the grand partition function, as written in Eq. (7), is multiplied by a constraint on the particle number in the form of a Kronecker delta, δN^,N. It turns out that, by re-expressing the Kronecker delta as a Fourier transform and integrating, one can arrive at a version of Eq. (9) in which the determinants are multiplied by a phase factor,63 

(10)

where ψm ≡ 2πm/Ns is the frequency of the Fourier transform. As the key grand canonical equations remain intact, one can sample the grand canonical ensemble as before but with a chemical potential that constrains the particle number to fluctuate around a predetermined value. Nevertheless, the same need to scan through chemical potentials to target particle numbers as in the grand canonical formalism makes this approach almost equally burdensome. A more direct canonical ensemble formalism is therefore crucial for boosting the computational efficiency, if not, feasibility, of many finite temperature auxiliary field simulations.

Given this backdrop, we have approached this problem from a very different angle: motivated by the simple recursive relations that can be used to construct the canonical partition functions of ideal gases,74,75 we have derived an analogous set of recursive relations to construct canonical partition functions in AFQMC. While such recursive relations are comparatively easy to construct for non-interacting systems because the energy spectra of such systems remain invariant under the addition/removal of particles, they are far harder to construct for interacting systems that possess many-body correlations. To surmount this fundamental problem, we have made use of the pivotal realization that the HS transformation of Eq. (4) essentially transforms our interacting problems into integrals over non-interacting problems. Therefore, so long as we can perform an HS transformation on a many-body Hamiltonian to decouple its interactions, we can then apply the same recursive relations that have been used for ideal gases to interacting systems.

We begin by deriving our recursive formalism for the canonical partition function and then one- and two-body observables for bosons and spinless fermions. As we will show in the following derivations, our canonical ensemble formulation is the same for both types of species except for the different signs that arise in their final expressions due to their disparate counting statistics.

1. Derivation of a recursive approach for obtaining the canonical partition function

Therefore, we denote â and â as particle creation and annihilation operators that can create/annihilate either fermions or bosons, respectively. The N-particle, canonical ensemble partition function may be expressed as

(11)

which we can factor and transform just as we did the grand partition function to obtain

(12)

In the above equation, we have added a subscript to the trace to differentiate it from that in Eq. (7) as the trace can only be taken over all N-particle quantum states in the canonical ensemble.

Using basic commutation/anticommutation rules, it can be proven40 that the product of B^(ϕl) exactly corresponds to an exponential of some one-body operator parameterized by the auxiliary field vector ϕ{ϕ1,,ϕL}, which we denote as Â(ϕ),

(13)

We can subsequently define a set of new boson/fermion coordinates in which Â(ϕ) is diagonal,40 

(14)

where we call the set of {ε̃γ(ϕ)} the effective single-particle spectrum because of the following relation:

(15)

and

(16)

Since eÂ(ϕ) is an independent-particle propagator that only depends on the auxiliary field vector, ϕ, the effective single-particle spectrum, {ε̃γ(ϕ)}, is independent of the particle number. For an N-particle, Ns-site system, taking the trace while constraining the particle number yields

(17)
(18)
(19)

Here, Γ is used to represent the set of N-particle states, and thus, Γn1++nNs=N and nγ denotes the number of particles in the γth eigenstate. For bosons, nγ = 0, 1, …, N, whereas for fermions, nγ = 0, 1. A more detailed derivation of these equations can be found in the  Appendix. The key implication of Eq. (19) is that, for a specified field ϕ, the single-particle spectrum can be decoupled from the particle number. Hence, the many-particle energy given such fields is simply the sum of all of the single-particle energies.

This key fact enables us to move beyond previous projection-based approaches and calculate Eq. (19) in a recursive fashion. To do so, we generalize the recursive approach to calculating canonical ensemble partition functions first developed for ideal gases.74,75 For an ideal gas of N particles whose total energy can be written as the sum of single-particle energies {εi}, we can express ZN in terms of its subsystem partition functions, ZNk,

(20)

with Z0 = 1 and zk, which can be identified as the partition function of the system at temperature , given by

(21)

Because we can write our decoupled many-body energy in terms of its single-particle energies, just as in the ideal case, we can then denote Trc(eβÂ(ϕ)) as ZN(ϕ) and re-express it using the same recursive formalism

(22)

with the same starting value Z0(ϕ)=1 and

(23)

In Eqs. (20) and (22), the plus sign should be employed whenever Bose statistics are imposed, while the minus sign should be employed whenever Fermi statistics are imposed. Otherwise, our formalism is identical for systems of bosons and spinless fermions since these systems involve only a single particle type and the total particle number can be constrained by a single value of N. Interestingly, this formalism bears a strong resemblance to a recursive formalism recently leveraged to incorporate particle statistics into path integral molecular dynamics simulations.76,77 However, whereas that formalism was derived with real space path integrals in mind, ours was designed with an eye toward Fock space.

Combining Eqs. (12), (13), and (22), we thus arrive at our final recursive expression for the many-body partition function

(24)

which demonstrates that the many-body partition function may be recovered by taking the integral over auxiliary field space of all recursively constructed, field-dependent partition functions. As is customary, this integral may be evaluated through Monte Carlo sampling.

2. Derivation of a recursive approach for obtaining canonical ensemble observables

In the grand canonical ensemble, most observables of interest can be easily expressed in terms of elements of the one-body density matrix using thermal Wick’s theorem.78 Unfortunately, thermal Wick’s theorem cannot be applied to products of field operators in the canonical ensemble.79 In this section, we therefore present an alternative approach to obtaining the expectation values of such products in the canonical ensemble based on their occupation numbers, which is in the spirit of Schönhammer’s derivations for non-interacting fermions.80 

In the canonical ensemble, the (i, j)th element of the one-body density matrix may be expressed as

(25)

Since the denominator can be explicitly evaluated via Eq. (24), we must determine an approach for evaluating the numerator, the unnormalized version of the one-body density matrix element, which we denote as D̃ij. We can proceed to re-express the numerator as we did for the partition function by performing an imaginary-time breakup and HS transformation of the propagators, which leads to

(26)
(27)
(28)

As is done in Eq. (16), we can simplify this equation by applying a change of basis to âiâj to yield

(29)

where Uλμijλ|ij|μ is an overlap matrix. This formula allows the integrand to be expressed as

(30)
(31)
(32)

where the second equality comes from the orthogonality of the eigenstates. We also define the mean occupation number, nλ, of the λth eigenstate with the total particle number N as

(33)

which, as is illustrated in detail in the  Appendix, can also be calculated in a recursive fashion,

(34)

Here, the plus sign is again used to preserve Bose statistics, while the minus sign is used to preserve fermion statistics. As such, combining Eqs. (28), (32), and (34), we arrive at the final expression for the one-body matrix elements,

(35)

Following a similar scheme, the unnormalized (i, j, k, l)th element of the two-body density matrix can be expressed as

(36)

where

(37)
(38)
(39)

and the change of basis is performed in the following way:

(40)

Observe that the expression for n^λ2N only possesses a plus sign, which stems from the fact that two fermions cannot inhabit the same site. Again, the plus sign in the above equations corresponds to Bose statistics, while the minus sign corresponds to Fermi statistics. The elements of higher order density matrices can be evaluated using the same basic scheme, which essentially requires manipulating the order of operators using commutation/anti-commutation relations to create a sequence of occupation numbers that can be recursively calculated. With these one- and two-body matrix elements in hand, we can subsequently compute such important quantities as energies and correlation functions.

Developing a canonical ensemble treatment for spinful fermion systems is naturally more difficult than developing a treatment for spinless fermion systems due to the additional constraint(s) that must be imposed on the partition function to constrain the spin in each spin sector. Moreover, more complicated anticommutation relations make generalizing our formalism to arbitrary spinful Hamiltonians, such as ab initio Hamiltonians, far more challenging. Nevertheless, to simplify our discussion, here we focus on developing a formalism for treating the Fermi–Hubbard model given by Eq. (55) as a special but important case, given the model’s long history in condensed matter physics.

Since spin degrees of freedom may also be decoupled via an HS transformation for Hamiltonians of this type, Eq. (7) can be rewritten as

(41)

which means that spin up and spin down states are propagated separately, and we can accordingly define two one-body operators, Â(ϕ) and Â(ϕ), such that

(42)

and

(43)

As a result, the partition function under a specified auxiliary field can be decomposed into the product of two subpartition functions with N and N particles, respectively,

(44)

which implies that we can perform iterations over particle numbers for spin up and spin down partition functions separately. Thus, in line with our previous discussion, we end up with

(45)

Although expressions for the one- and two-body density matrix elements are more difficult to obtain for more general Hamiltonians, in the case of the Hubbard model, we can independently calculate matrix elements in different spin sectors, given a set of sampled auxiliary field vectors, ϕ. Thus, even though the field vectors that are ultimately sampled depend on both spin sectors, once these are obtained, the density matrix elements can be calculated separately. Hence, the unnormalized (i, j)th elements of the spin up/down one-body matrices take the same form as Eq. (32),

(46)

with the recursive relation being

(47)

Additionally, the unnormalized (i, j, k, l)th elements of the two-body density matrix may be obtained via

(48)

where the change of basis is applied in the following way:

(49)

Note that, for the Hubbard model in which repulsions occur locally, only the (i, i, i, i)th elements of Eq. (48) have nonzero contributions.

Although our formulas for the field-dependent partition functions are exact, as written, these formulas are not necessarily computationally expedient to evaluate. First, the numerical effort to evaluate Eq. (22) scales with the square of the particle number N. Moreover, zk(ϕ) grows exponentially with the inverse temperature β and the particle number N, which could cause severe numerical overflow issues at low temperatures.

To mitigate these issues, note that in all mean occupation number recursions, namely, Eqs. (34), (38), (39), and (47), the partition functions appear as ratios. Hence, we may circumvent these computational challenges by calculating the ratios of these partition functions instead of the individual partition functions themselves. To achieve this, we may sum Eq. (34) over λ and use the conservation of the particle number to obtain

(50)

Some reordering yields

(51)

Observables may be obtained for each sampled field at a dramatically reduced cost based on this equation as its cost only scales linearly with the particle number.

In order to calculate observables that are expressed in terms of combinations of particle operators, we need to average appropriate combinations of density matrices through a Monte Carlo technique. As an illustration, we use an arbitrary one-body operator T^=i,jtijâiâj and plug in Eq. (35) to yield

(52)

with n being the number of Monte Carlo samples and Φ being sampled from a Boltzmann-like distribution

(53)

This can be accomplished by using standard sampling algorithms such as the heat-bath or Metropolis algorithms.81 

In this paper, we elect to use the heat-bath algorithm to perform our integration over auxiliary fields because of its greater efficiency. If we define R as the ratio of the new field-dependent partition function ZN(Φ) to the old partition function ZN(Φ), the acceptance ratio of a random flip of discrete auxiliary fields or a random perturbation to continuous auxiliary fields may be written as

(54)

To benchmark this formalism against ED calculations and assess its performance relative to grand canonical ensemble algorithms, we implemented our own algorithms in Julia. Except where otherwise indicated, we set Δτ = 0.02 in our calculations to yield accurate results comparable to those produced by ED, which we also implemented ourselves for both canonical and grand canonical ensemble calculations. In general, smaller Δτs lead to more accurate results. However, the effective single-particle spectrum, {ε̃γ(ϕ)}, is Δτ-dependent and becomes extremely large when Δτ is too small, exacerbating the numerical sign problem illustrated in Sec. III B. Δτ = 0.02 thus is a balance point for us to investigate the models described in the following over physically meaningful parameter regimes without sign problems. All canonical ensemble AFQMC results are averaged over 50 walkers, which each sweep through the full lattice 2 × 105 times. The grand canonical ensemble results reported were obtained using a version of our previously described finite temperature ab initio AFQMC code modified to accommodate the Fermi–Hubbard model.44 All of the codes that accompany this work can be found in our project Github repository.82 

To test the accuracy and analyze the computational performance of our canonical ensemble AFQMC (CE-AFQMC) algorithm, we benchmarked it against ED results for both the Bose and Fermi–Hubbard models. These models have long been used to illustrate and study the effects of strong correlation in materials, free of complicating material-specific details. As its name suggests, the Fermi–Hubbard model (or, just the “Hubbard” model) exemplifies the two most basic behaviors of fermions, hopping and repulsion/attraction, on a lattice whose sites represent the locations of atoms within a material.83–85 The spinful version of this model is defined by the Hamiltonian

(55)

where the first, one-body term describes the hopping of fermions from site i to nearest-neighbor sites j (⟨ij⟩ denotes all pairs of nearest-neighbor sites ij), while the second term denotes the local, two-body fermion–fermion attraction/repulsion on the same site i. t denotes the hopping parameter, which quantifies the relative magnitude of the hopping contribution to the Hamiltonian; U is the attraction/repulsion constant that quantifies the relative magnitude of the two-body interaction to the Hamiltonian. If U > 0, the fermion–fermion interaction is repulsive and the fermions are discouraged from occupying the same site; if U < 0, the fermions are attracted to one another. and denote the up and down spins of the fermions.

The Bose–Hubbard model is a generalization of this model adapted to model bosons,

(56)

where b^i/b^i denote boson creation/annihilation operators at site i and n^i=b^ib^i denotes the boson number operator. As before, t and U modulate the relative contributions of the boson hopping and interaction terms. In the past, this model has been employed to describe the physics of bosons in optical lattices86 and superfluids.87 

In Fig. 1, we compare our canonical ensemble AFQMC results against canonical ensemble results obtained using ED for a 1D, three-site Bose–Hubbard model containing three bosons, while in Fig. 2, we compare our CE-AFQMC and ED results for a 1D, six-site Hubbard model at half-filling. In order to obtain the total energies of these systems, we must compute their kinetic energies, Ek, which are one-body quantities, and potential energies, Ep, which are two-body quantities, which serve as direct tests of the formalism we presented for computing many-body density matrices in Secs. II C and II D. For both systems, we make these comparisons for several different positive values of U from high temperatures down to sufficiently low temperatures that their energies begin to plateau to their ground state values. It is apparent from Figs. 1 and 2 that our algorithm yields exact results, within statistical error bars, for both of these models across all of the U values and temperatures surveyed. Direct numerical comparisons are presented in the tables in the supplementary material. Interestingly, whereas the Bose–Hubbard results monotonically converge to their final plateaus with decreasing temperature, the Fermi–Hubbard results manifest kinks at larger U values in their potential energy curves between 0.75 and 2 1/kBT, illustrating the competition that occurs between the model’s kinetic and potential contributions. These results verify that our formalism can achieve exact results using the heat bath algorithm described above over physically meaningful parameter regimes.

FIG. 1.

Total (Etot), kinetic (Ek), and potential (Ep) energies of a three-site Bose–Hubbard model simulated for several values of U at t = 1 and N = 3. Energies are given in units of t.

FIG. 1.

Total (Etot), kinetic (Ek), and potential (Ep) energies of a three-site Bose–Hubbard model simulated for several values of U at t = 1 and N = 3. Energies are given in units of t.

Close modal
FIG. 2.

Total (Etot), kinetic (Ek), and potential (Ep) energies of a six-site Hubbard model simulated for several values of U at t = 1 and N = N = 3. Energies are given in units of t.

FIG. 2.

Total (Etot), kinetic (Ek), and potential (Ep) energies of a six-site Hubbard model simulated for several values of U at t = 1 and N = N = 3. Energies are given in units of t.

Close modal

The results presented for the Fermi–Hubbard model above naturally beg the question how does the sign problem arise in this formalism. After all, the sign problem should arise in this QMC approach, just as it arises in all previous QMC approaches.88 We observe two origins of the sign problem in our formalism—one “physical” and one “numerical.” The unavoidable “physical” origin of the sign problem in our approach stems from the fact that the expectation value of the propagator for the states of certain systems in Eq. (19) can become negative, resulting in a mixture of positive and negative terms that can cancel one another out during sampling. This physical sign problem corresponds to that which emerges in grand canonical ensemble algorithms when the trace over grand canonical states yields products of determinants that can also be negative.46 

While we did not observe such a sign problem in the Fermi–Hubbard simulations presented above, we did observe the second numerical form of the sign problem. This numerical sign problem arises from rounding errors that can accumulate during the calculation of Eq. (34) and has been reported in several previous papers describing recursive treatments of ideal Fermi gases.79,80,89 More specifically, when a large eβελ(ϕ) appears in the recursive calculation, its corresponding mean occupation number, nλN, tends to approach 1, yielding 1nλN values as small as 10−15. As illustrated in Table IV of the supplementary material, when this occurs, the corresponding occupation number cannot be correctly represented by finite precision arithmetic and corrupts the recursion. As a result, there exists a critical particle number, Nc, at which nλNc begins to exceed 1 for the λ that corresponds to the largest eigenvalue of eβA(ϕ), which violates the Pauli exclusion principle and causes the ratio of partition functions in Eq. (51) to become negative in turn. This is illustrated for a 16-site Hubbard model in Fig. 3, from which we can clearly see a rapid collapse of the average sign as the number of particles surpasses Nc at β = 2 and ξ = (meaning that no approximations were invoked). In contrast with the physical sign problem described above, if left untreated, this numerical sign problem renders the simulation either completely sign-problem-free or completely noisy as a function of filling.

FIG. 3.

The average sign of ZN(ϕ) as a function of filling for a 16-site, U = 4 Hubbard model for various β’s and ξ cutoff values. ξ = means that no cutoff approximation is employed.

FIG. 3.

The average sign of ZN(ϕ) as a function of filling for a 16-site, U = 4 Hubbard model for various β’s and ξ cutoff values. ξ = means that no cutoff approximation is employed.

Close modal

To curb this problem, we have developed an approximation similar in spirit to the selection of active spaces in other electronic structure techniques that strikes a balance between alleviating the numerical sign problem and biasing the simulation. This approximation is motivated by the observation that the numerical instabilities our simulations incur stem from the fact that, at low temperatures, eβελ(ϕ)’s become extremely large for small λ’s and extremely small for large λ’s. Thus, to prevent our results from being corrupted by round-off errors, the influence of these large and small eβελ(ϕ)’s on the other eβελ(ϕ)’s must be reduced. To accomplish this, note that εN(ϕ) is the energy level closest to the Fermi level and, at low temperatures, the occupancies of energy levels far below the Fermi level tend to be “frozen” around 1. Based on this observation, a coarse-grained scheme can be devised in which the occupancies of the energy levels ελl(ϕ) (where l denotes lower) are set to 1 (nλlN=1) if eβ(εN(ϕ)ελl(ϕ))>ξ, where ξ is a cutoff value. A similar approximation can then be made by setting nλuN=0 for energy levels ελu(ϕ) (where u denotes upper) well above the Fermi level such that eβ(ελu(ϕ)εN(ϕ))>ξ. The smaller the ξ, the smaller the numerical sign problem, but the larger the bias. Conversely, the larger the ξ, the larger the numerical sign problem, but the smaller the bias.

That said, we benchmarked the performance of this approximation against ED for the Hubbard model for different ξ values at different temperatures and U’s, as reported in Table I. As expected, a finer approximation, ξ = 104, yields more accurate results, and with this cutoff value, less than 1% of the samples are corrupted by the numerical sign problem. In contrast, a coarser approximation with ξ = 103 completely eliminates the numerical sign problem for all parameter regimes we investigated, but its results are more biased. To better analyze the applicability of this scheme, we performed several simulations on the 16-site Hubbard model with U = 4 and plotted the average sign of the partition function ZN(ϕ) as a function of filling, n=N+NNs. We chose ξ = 103 to get rid of any effects from the numerical sign problem, in order to better investigate the emergence of the physical sign problem in our algorithm. As illustrated in Fig. 3, we find that when the cutoff approximation is applied, the overall sign problem is greatly improved at β = 2, and the physical sign problem emerges at ⟨n⟩ = 1, which is a consequence of the fact that as one iterates over the particle number, some extreme auxiliary field configurations are more likely to occur. For the β = 4 curve, one finds that the physical sign problem emerges at ⟨n⟩ = 0.75 and decreases as the filling increases. These examples demonstrate that a reasonable choice of ξ can mitigate the numerical sign problem without significantly biasing the algorithm’s results, highlighting the algorithm’s overall scalability.

TABLE I.

Comparison of QMC calculations with different cutoff thresholds (ξ) against Exact Diagonalization (ED) results for the U = 4 Hubbard model at β = 5. The λlth eigenstate is discarded in the recursive calculation if expβ(εN(ϕ)ελl(ϕ))>ξ holds, where N=12nNs.

EDQMC, ξ = 103QMC, ξ = 104
Ns = 6, ⟨n⟩ = 1 −3.6429 −3.61(3) −3.67(6) 
Ns = 8, ⟨n⟩ = 0.75 −6.6298 −6.61(1) −6.64(5) 
Ns = 8, ⟨n⟩ = 1 −4.5002 −4.44(3) −4.51(5) 
EDQMC, ξ = 103QMC, ξ = 104
Ns = 6, ⟨n⟩ = 1 −3.6429 −3.61(3) −3.67(6) 
Ns = 8, ⟨n⟩ = 0.75 −6.6298 −6.61(1) −6.64(5) 
Ns = 8, ⟨n⟩ = 1 −4.5002 −4.44(3) −4.51(5) 

One of the advantages that accompanies simulating in a different ensemble is that properties computed in that ensemble may converge to the ground state as a function of temperature in a different manner than in the original ensemble. Any differences in convergence to the ground state between the canonical and grand canonical ensembles would be particularly important to note in the context of AFQMC because the higher the temperature at which properties approximate ground state properties, the more likely that those properties can be reliably modeled without a significant sign problem. To draw this comparison between ensembles, we therefore simulated the same six-site Fermi–Hubbard model described in Sec. III A using both our CE and GCE44,62 algorithms from high to low temperatures at which the model’s properties approach their ground state values. Note that, although CE and GCE properties may differ at finite temperatures, they should converge to the same ground state values. This is because while states with particle numbers that differ from the average may be readily accessed in grand canonical simulations at finite temperatures, the step-like nature of the particle number vs chemical potential curve that arises at low temperatures ensures that only states with a fixed N are sampled in the grand canonical ensemble at temperatures approaching the ground state. This consequently ensures that the same ground state energy will be attained in both ensembles. As illustrated in Fig. 4 for U = 2, we find that our CE-AFQMC algorithm converges significantly more rapidly than its GCE counterpart—by about a factor of a few 1/kBT units. Although this may seem small, a few 1/kBT units can make the difference between being able to secure coveted insights into low-temperature physics and not. The intuitive reason why the canonical ensemble converges more rapidly is because, as alluded to above, it accesses fewer higher-energy states at any given temperature than the grand canonical ensemble. We expect this convergence difference to grow with increasing U, as evidenced in the supplementary material. This finding suggests that our canonical algorithm may assume an essential role in future research focused on understanding how finite temperature algorithms and properties converge to the ground state—especially given that ground state AFQMC algorithms are performed at fixed particle numbers.

FIG. 4.

Six-site Fermi–Hubbard model simulated within the canonical ensemble with N = N = 3 and the grand canonical ensemble with ⟨N⟩ = 6 (half-filling) at U = 2. Δτ is set to 0.05 at β = 8 to avoid numerical sign problems. The horizontal dashed line at E = −5.4095 is the exact ground state energy.

FIG. 4.

Six-site Fermi–Hubbard model simulated within the canonical ensemble with N = N = 3 and the grand canonical ensemble with ⟨N⟩ = 6 (half-filling) at U = 2. Δτ is set to 0.05 at β = 8 to avoid numerical sign problems. The horizontal dashed line at E = −5.4095 is the exact ground state energy.

Close modal

Given our algorithm’s encouraging ground state convergence properties, one may ask how it scales with sites and filling relative to the grand canonical algorithm. Because it relies upon the same imaginary time propagation of matrices as the DQMC and GCE-AFQMC algorithms before it,42,90 it requires O(LNs3) operations to construct its full propagator. If we group all L imaginary time steps into m sets, the algorithm moreover requires O(L2m2Ns3) operations for numerical stabilizations, as have previous algorithms.90 The distinguishing feature of our algorithm is its use of recursions. To initialize recursions, we must diagonalize the matrix eβA(ϕ), and to account for N particles, O(N) recursions are required. These two operations thus involve an additional computational cost that scales as O(Ns3+N). Therefore, unless NNs, the cost of performing recursions is negligible compared to the usual cost of propagating matrices. However, previous DQMC algorithms can rely on the Sherman–Morrison formula, which only requires O(Ns2) operations to update the determinant of a matrix with a rank-one perturbation, to quickly update their Green’s functions.37,39 Because no analogous formula exists for updating the eigenvalues of a matrix with a rank-one perturbation, as is needed in our algorithm, the numerical cost of updating quantities after each Monte Carlo step in our approach scales as O(Ns3), making it the most expensive part of our algorithm.

In this work, we have presented a new finite temperature AFQMC formalism for computing observables within the canonical ensemble. Unlike previous canonical ensemble algorithms, which relied upon projecting canonical ensemble properties out from grand canonical simulations, our formalism is based on recursions among purely canonical ensemble quantities, eliminating the costly need to scan through chemical potentials to fix average particle numbers in most DMC and AFQMC simulations. To verify the accuracy of our approach while also demonstrating its broad applicability, we used our algorithm to calculate thermal properties of the Bose and Fermi–Hubbard models, quintessential models of strong correlation, and benchmarked it against ED results. We moreover demonstrated that our canonical algorithm converges to ground states more rapidly with decreasing temperature than conventional grand canonical AFQMC algorithms, which illustrates its promise for studying how low-temperature phenomena ultimately give rise to ground state phenomena.

Although we presented a rather limited range of illustrations of our algorithm here, we foresee our algorithm having a wide variety of practical applications. First and foremost, because our canonical ensemble approach definitively fixes the particle number, we believe that it will eliminate the “rogue eigenvalue problem” that curtails the direct application of grand canonical AFQMC algorithms to boson and fermion condensates whose average particle numbers become challenging to fix at low temperatures.62 This will enable the study of low temperature Bose and Bose–Fermi condensates62 as well as superconductors91 with the high accuracy typical of AFQMC techniques. In addition, our algorithm will enable apples-to-apples, canonical ensemble comparisons between FT-AFQMC and other canonical ensemble finite temperature electronic structure techniques, such as density matrix QMC92 and finite temperature coupled cluster approaches,18–21 for the first time. Some of these comparisons will necessitate generalizing our formalism to arbitrary ab initio Hamiltonians, which is an ongoing effort. This will provide the community with the critical ability to compare and refine finite temperature methods, much as has been recently done for ground state electronic structure techniques.93,94

Supplementary data including additional benchmarks and an illustration of the numerical sign problem may be found in the supplementary material.

The authors thank Richard Stratt, Hongxia Hao, James Shepherd, and members of the Center for Predictive Simulation of Functional Materials for stimulating conversations. T.S. received partial support from the Brown Open Graduate Education program, while Y.L. was graciously supported by the Brown Presidential Fellows program. B.M.R. acknowledges funding from the U.S. Department of Energy, Office of Science, Basic Energy Sciences, Materials Sciences and Engineering Division, as part of the Computational Materials Sciences Program and Center for Predictive Simulation of Functional Materials. Y.Y. contributed to this work while a part of the International Exchange Program for Excellent Students of the University of Science and Technology of China. All calculations were conducted using computational resources and services at the Brown University Center for Computation and Visualization.

The data that support the findings of this study are also openly available in the Zenodo repository at https://doi.org/10.5281/zenodo.3991899.95 

In this appendix, we provide detailed derivations of several recursive relations that are essential to obtaining the partition functions and density matrices for a specific auxiliary field. Some of these equations have been derived in other contexts before.62,96–98

First, we derive Eqs. (17)–(19), which undergird our expressions for the canonical partition function. Suppose we have an arbitrary N × N matrix A expressed in the single-particle basis and form two Fock-space operators that operate on boson and fermion Fock states, respectively,

(A1)
(A2)

We then consider the operations Âbb^ and Âfĉ. Using the elementary commutation/anticommutation relations, it is straightforward to show that

(A3)
(A4)

which gives

(A5)
(A6)

where I is the N × N unit matrix. Note that the elements of the matrices A and I are c numbers in Fock space and  is a scalar relative to the N × N matrices. As a result, I and A commute as matrices. Repeatedly applying these two equations yields

(A7)
(A8)

for any positive integer n. Hence,

(A9)
(A10)

can be shown inductively, where in the last step of both equations, the exponential is allowed to be broken up because the two terms in the exponent commute, as noted above. These results allow eβÂ(ϕ) in Eq. (17) to “walk through” the string of creation operators that defines the N-particle state |Γ⟩. More specifically, by expanding |Γ⟩ in terms of particle creation operators, we have

(A11)

where, in the last step, we have used the fact eÂ(ϕ)|0=|0. Thus, sandwiching eÂ(ϕ) between ⟨Γ| and |Γ⟩ yields

(A12)

which enables us to proceed from Eqs. (18) and (19). Note that Eqs. (A9) and (A10) have exactly the same form. As such, Eq. (A12) holds true for both boson and fermion propagators as they propagate the corresponding states in the same manner.

Next, we show how to iteratively obtain the mean occupation numbers of the eigenstates, as well as the products of those occupancies, i.e., Eqs. (34), (38), and (39). These are needed to evaluate a system’s energy and correlation functions. Since bosons and fermions possess different particle statistics, we provide their derivations separately, starting with those for bosons.

To arrive at Eq. (34), without loss of generality, we consider the λ = 1 case of Eq. (33). By definition,

(A13)

where we have introduced the shorthand {nγ}1N to denote all states that conserve their particle number such that γ=1Nsnγ=N and have omitted all field dependence for clarity. Since n1 can take values from 0 up to N for bosons, the above equation can be re-expressed as

(A14)

where in the first equality, we can perform the sum over n1 first and leave the rest constrained by the particle number conservation constraint {nγ}2Nn1γ=2Nsnγ=Nn1, while in the second equality, we have used the fact that the enumerations of {n2,,nNs} are independent of {n1 = 0}. This allows us to count n1 from n1 = 1, instead of n1 = 0. Replacing n^1 with n^12 in the expectation values above gives

(A15)

The procedure leading to Eq. (A14) can be easily extended to calculate expectation values of products of different occupation numbers, i.e., ninjN with ij. Replacing n^1 with n^1n^2 and observing that the enumerations of {n3,,nNs} are independent from {n1 = 0, n2 = 0}, one obtains

(A16)

which is equivalent to Eq. (39) as the indices in that equation are dummy variables, and hence, the equation can be rewritten into a symmetric form.

Since nγ can only be 0 and 1 for fermions, we arrive at the following two identities:

(A17)
(A18)

Substituting these identities into the recursive relation for the mean occupation number [Eq. (A13)] yields

(A19)

Proceeding along similar lines yields the fermion version of Eq. (A16),

(A20)

These are precisely the equations referenced in the text.

1.
A.
Szabo
and
N.
Ostlund
,
Modern Quantum Chemistry: Introduction to Advanced Electronic Structure Theory
, Dover Books on Chemistry (
Dover Publications
,
1996
).
2.
I.
Shavitt
and
R.
Bartlett
,
Many-Body Methods in Chemistry and Physics: MBPT and Coupled-Cluster Theory
, Cambridge Molecular Science (
Cambridge University Press
,
2009
).
3.
R. J.
Bartlett
and
M.
Musiał
, “
Coupled-cluster theory in quantum chemistry
,”
Rev. Mod. Phys.
79
,
291
352
(
2007
).
4.
W.
Lorenzen
,
A.
Becker
, and
R.
Redmer
, “
Progress in warm dense matter and planetary physics
,” in
Frontiers and Challenges in Warm Dense Matter
(
Springer International Publishing
,
2014
), pp.
203
234
.
5.
J. M.
McMahon
,
M. A.
Morales
,
C.
Pierleoni
, and
D. M.
Ceperley
, “
The properties of hydrogen and helium under extreme conditions
,”
Rev. Mod. Phys.
84
(
4
),
1607
(
2012
).
6.
S.
Mukherjee
,
F.
Libisch
,
N.
Large
,
O.
Neumann
,
L. V.
Brown
,
J.
Cheng
,
J. B.
Lassiter
,
E. A.
Carter
,
P.
Nordlander
, and
N. J.
Halas
, “
Hot electrons do the impossible: Plasmon-induced dissociation of H2 on Au
,”
Nano Lett.
13
,
240
247
(
2013
).
7.
L.
Zhou
,
C.
Zhang
,
M. J.
McClain
,
A.
Manjavacas
,
C. M.
Krauter
,
S.
Tian
,
F.
Berg
,
H. O.
Everitt
,
E. A.
Carter
,
P.
Nordlander
, and
N. J.
Halas
, “
Aluminum nanocrystals as a plasmonic photocatalyst for hydrogen dissociation
,”
Nano Lett.
16
,
1478
1484
(
2016
).
8.
H.
Yong
,
N.
Zotev
,
J. M.
Ruddock
,
B.
Stankus
,
M.
Simmermacher
,
A. M.
Carrascosa
,
W.
Du
,
N.
Goff
,
Y.
Chang
,
D.
Bellshaw
,
M.
Liang
,
S.
Carbajo
,
J. E.
Koglin
,
J. S.
Robinson
,
S.
Boutet
,
M. P.
Minitti
,
A.
Kirrander
, and
P. M.
Weber
, “
Observation of the molecular response to light upon photoexcitation
,”
Nat. Commun.
11
,
2157
(
2020
).
9.
M.
Imada
,
A.
Fujimori
, and
Y.
Tokura
, “
Metal–insulator transitions
,”
Rev. Mod. Phys.
70
,
1039
1263
(
1998
).
10.
J.
Orenstein
and
A. J.
Millis
, “
Advances in the physics of high-temperature superconductivity
,”
Science
288
,
468
474
(
2000
).
11.
Q.
Si
,
R.
Yu
, and
E.
Abrahams
, “
High-temperature superconductivity in iron pnictides and chalcogenides
,”
Nat. Rev. Mater.
1
,
16017
(
2016
).
12.
T. D.
Ladd
,
F.
Jelezko
,
R.
Laflamme
,
Y.
Nakamura
,
C.
Monroe
, and
J. L.
O’Brien
, “
Quantum computers
,”
Nature
464
,
45
53
(
2010
).
13.
S. A.
Jafari
, “
Introduction to Hubbard model and exact diagonalization
,”
Iranian J. Phys. Res.
8
,
113
(
2008
); arXiv:0807.4878 [cond-mat.str-el].
14.
N. D.
Mermin
, “
Stability of the thermal Hartree–Fock approximation
,”
Ann. Phys.
21
,
99
121
(
1963
).
15.
G.
Sanyal
,
S. H.
Mandal
, and
D.
Mukherjee
, “
On optimal mean-field descriptions in finite-temperature many-body theories: Use of thermal Brillouin and Bruckner conditions
,”
Proc. Indiana Acad. Sci. - Chem. Sci.
106
,
407
(
1994
).
16.
V. V.
Karasiev
,
T.
Sjostrom
,
J.
Dufty
, and
S. B.
Trickey
, “
Accurate homogeneous electron gas exchange–correlation free energy for local spin-density calculations
,”
Phys. Rev. Lett.
112
,
076403
(
2014
).
17.
A.
Pribram-Jones
,
S.
Pittalis
,
E. K. U.
Gross
, and
K.
Burke
, “
Thermal density functional theory in context
,” in
Frontiers and Challenges in Warm Dense Matter
, edited by
F.
Graziani
,
M. P.
Desjarlais
,
R.
Redmer
, and
S. B.
Trickey
(
Springer International Publishing
,
Cham
,
2014
), pp.
25
60
.
18.
G.
Harsha
,
T. M.
Henderson
, and
G. E.
Scuseria
, “
Thermofield theory for finite-temperature quantum chemistry
,”
J. Chem. Phys.
150
,
154109
(
2019
).
19.
G.
Harsha
,
T. M.
Henderson
, and
G. E.
Scuseria
, “
Thermofield theory for finite-temperature coupled cluster
,”
J. Chem. Theory Comput.
15
,
6127
6136
(
2019
).
20.
A. F.
White
and
G. K.-L.
Chan
, “
A time-dependent formulation of coupled-cluster theory for many-fermion systems at finite temperature
,”
J. Chem. Theory Comput.
14
,
5690
5700
(
2018
).
21.
A. F.
White
and
G.
Kin-Lic Chan
, “
Finite-temperature coupled cluster: Efficient implementation and application to prototypical systems
,”
J. Chem. Phys.
152
,
224104
(
2020
).
22.
P.
Shushkov
and
T. F.
Miller
, “
Real-time density-matrix coupled-cluster approach for closed and open systems at finite temperature
,”
J. Chem. Phys.
151
,
134107
(
2019
).
23.
P. K.
Jha
and
S.
Hirata
, “
Numerical evidence invalidating finite-temperature many-body perturbation theory
,” in
Annu. Rep. Comput. Chem.
(
Elsevier
,
2019
), Chap. 1, pp.
3
15
.
24.
P. K.
Jha
and
S.
Hirata
, “
Finite-temperature many-body perturbation theory in the canonical ensemble
,”
Phys. Rev. E
101
,
022106
(
2020
).
25.
A. F.
White
and
G.
Kin-Lic Chan
, “
Comment on “Numerical evidence falsifying finite-temperature many-body perturbation theory”
,” arXiv:1810.03653 (
2018
).
26.
S.
Hirata
and
P. K.
Jha
, “
Converging finite-temperature many-body perturbation theory in the grand canonicals ensemble that conserves the average number of electrons
,” in
Annu. Rep. Comput. Chem.
(
Elsevier
,
2019
), Chap. 2, pp.
17
37
.
27.
S.
Hirata
and
P. K.
Jha
, “
Finite-temperature many-body perturbation theory in the grand canonical ensemble
,”
J. Chem. Phys.
153
,
014103
(
2020
).
28.
R.
Santra
and
J.
Schirmer
, “
Finite-temperature second-order many-body perturbation theory revisited
,”
Chem. Phys.
482
,
355
361
(
2017
), part of Special Issue: Electrons and nuclei in motion - correlation and dynamics in molecules (on the occasion of the 70th birthday of Lorenz S. Cederbaum).
29.
D.
Zgid
and
E.
Gull
, “
Finite temperature quantum embedding theories for correlated systems
,”
New J. Phys.
19
,
023047
(
2017
).
30.
A.
Georges
,
G.
Kotliar
,
W.
Krauth
, and
M. J.
Rozenberg
, “
Dynamical mean-field theory of strongly correlated fermion systems and the limit of infinite dimensions
,”
Rev. Mod. Phys.
68
,
13
125
(
1996
).
31.
G.
Kotliar
,
S. Y.
Savrasov
,
K.
Haule
,
V. S.
Oudovenko
,
O.
Parcollet
, and
C. A.
Marianetti
, “
Electronic structure calculations with dynamical mean-field theory
,”
Rev. Mod. Phys.
78
,
865
951
(
2006
).
32.
A. R.
Welden
,
A. A.
Rusakov
, and
D.
Zgid
, “
Exploring connections between statistical mechanics and Green’s functions for realistic systems: Temperature dependent electronic entropy and internal energy from a self-consistent second-order Green’s function
,”
J. Chem. Phys.
145
,
204106
(
2016
).
33.
A. A.
Kananenka
,
A. R.
Welden
,
T. N.
Lan
,
E.
Gull
, and
D.
Zgid
, “
Efficient temperature-dependent Green’s function methods for realistic systems: Using cubic spline interpolation to approximate Matsubara Green’s functions
,”
J. Chem. Theory Comput.
12
,
2250
2259
(
2016
).
34.
A. A.
Kananenka
,
J. J.
Phillips
, and
D.
Zgid
, “
Efficient temperature-dependent Green’s functions methods for realistic systems: Compact grids for orthogonal polynomial transforms
,”
J. Chem. Theory Comput.
12
,
564
571
(
2016
).
35.
W. M. C.
Foulkes
,
L.
Mitas
,
R. J.
Needs
, and
G.
Rajagopal
, “
Quantum Monte Carlo simulations of solids
,”
Rev. Mod. Phys.
73
,
33
83
(
2001
).
36.
J.
Toulouse
,
R.
Assaraf
, and
C. J.
Umrigar
, “
Introduction to the variational and diffusion Monte Carlo methods
,” in
Electron Correlation in Molecules - Ab Initio Beyond Gaussian Quantum Chemistry
, Advances in Quantum Chemistry, edited by
P. E.
Hoggan
and
T.
Ozdogan
(
Academic Press
,
2016
), Chap. 15, Vol. 73, pp.
285
314
.
37.
J.
Gubernatis
,
N.
Kawashima
, and
P.
Werner
,
Quantum Monte Carlo Methods: Algorithms for Lattice Models
(
Cambridge University Press
,
2016
).
38.
M.
Motta
and
S.
Zhang
, “
Ab initio computations of molecular systems by the auxiliary-field quantum Monte Carlo method
,”
Wiley Interdiscip. Rev.: Comput. Mol. Sci.
8
,
e1364
(
2018
).
39.
S. R.
White
,
D. J.
Scalapino
,
R. L.
Sugar
,
E. Y.
Loh
,
J. E.
Gubernatis
, and
R. T.
Scalettar
, “
Numerical study of the two-dimensional Hubbard model
,”
Phys. Rev. B
40
,
506
516
(
1989
).
40.
J. E.
Hirsch
, “
Two-dimensional Hubbard model: Numerical simulation study
,”
Phys. Rev. B
31
,
4403
4419
(
1985
).
41.
Z.
Bai
,
W.
Chen
,
R.
Scalettar
, and
I.
Yamazaki
, “
Numerical methods for quantum Monte Carlo simulations of the Hubbard model
,” in
Multi-Scale Phenomena in Complex Fluids
(
World Scientific
,
2009
), pp.
1
110
.
42.
S.
Zhang
, “
Finite-temperature Monte Carlo calculations for systems with fermions
,”
Phys. Rev. Lett.
83
,
2777
2780
(
1999
).
43.
Y.-Y.
He
,
H.
Shi
, and
S.
Zhang
, “
Reaching the continuum limit in finite-temperature ab initio field-theory computations in many-fermion systems
,”
Phys. Rev. Lett.
123
,
136402
(
2019
).
44.
Y.
Liu
,
M.
Cho
, and
B.
Rubenstein
, “
Ab initio finite temperature auxiliary field quantum Monte Carlo
,”
J. Chem. Theory Comput.
14
,
4722
4732
(
2018
).
45.
Z.-X.
Li
and
H.
Yao
, “
Sign-problem-free fermionic quantum Monte Carlo: Developments and applications
,”
Annu. Rev. Condens. Matter Phys.
10
,
337
356
(
2019
).
46.
E. Y.
Loh
,
J. E.
Gubernatis
,
R. T.
Scalettar
,
S. R.
White
,
D. J.
Scalapino
, and
R. L.
Sugar
, “
Sign problem in the numerical simulation of many-electron systems
,”
Phys. Rev. B
41
,
9301
9307
(
1990
).
47.
S.
Zhang
,
J.
Carlson
, and
J. E.
Gubernatis
, “
Constrained path Monte Carlo method for fermion ground states
,”
Phys. Rev. B
55
,
7464
7477
(
1997
).
48.
S.
Zhang
and
H.
Krakauer
, “
Quantum Monte Carlo method using phase-free random walks with Slater determinants
,”
Phys. Rev. Lett.
90
,
136401
(
2003
).
49.
S.
Zhang
,
J.
Carlson
, and
J. E.
Gubernatis
, “
Constrained path quantum Monte Carlo method for fermion ground states
,”
Phys. Rev. Lett.
74
,
3652
3655
(
1995
).
50.
B.-X.
Zheng
,
C.-M.
Chung
,
P.
Corboz
,
G.
Ehlers
,
M.-P.
Qin
,
R. M.
Noack
,
H.
Shi
,
S. R.
White
,
S.
Zhang
, and
G. K.-L.
Chan
, “
Stripe order in the underdoped region of the two-dimensional Hubbard model
,”
Science
358
,
1155
1160
(
2017
).
51.
J. P. F.
LeBlanc
,
A. E.
Antipov
,
F.
Becca
,
I. W.
Bulik
,
G. K.-L.
Chan
,
C.-M.
Chung
,
Y.
Deng
,
M.
Ferrero
,
T. M.
Henderson
,
C. A.
Jiménez-Hoyos
,
E.
Kozik
,
X.-W.
Liu
,
A. J.
Millis
,
N. V.
Prokof’ev
,
M.
Qin
,
G. E.
Scuseria
,
H.
Shi
,
B. V.
Svistunov
,
L. F.
Tocchio
,
I. S.
Tupitsyn
,
S. R.
White
,
S.
Zhang
,
B.-X.
Zheng
,
Z.
Zhu
, and
E.
Gull (Simons Collaboration on the Many-Electron Problem)
, “
Solutions of the two-dimensional Hubbard model: Benchmarks and results from a wide range of numerical algorithms
,”
Phys. Rev. X
5
,
041041
(
2015
).
52.
W.
Purwanto
,
S.
Zhang
, and
H.
Krakauer
, “
An auxiliary-field quantum Monte Carlo study of the chromium dimer
,”
J. Chem. Phys.
142
,
064302
(
2015
).
53.
M.
Suewattana
,
W.
Purwanto
,
S.
Zhang
,
H.
Krakauer
, and
E. J.
Walter
, “
Phaseless auxiliary-field quantum Monte Carlo calculations with plane waves and pseudopotentials: Applications to atoms and molecules
,”
Phys. Rev. B
75
,
245123
(
2007
).
54.
J.
Shee
,
S.
Zhang
,
D. R.
Reichman
, and
R. A.
Friesner
, “
Chemical transformations approaching chemical accuracy via correlated sampling in auxiliary-field quantum Monte Carlo
,”
J. Chem. Theory Comput.
13
,
2667
2680
(
2017
).
55.
Y.
Liu
,
T.
Shen
,
H.
Zhang
, and
B.
Rubenstein
, “
Unveiling the finite temperature physics of hydrogen chains via auxiliary field quantum Monte Carlo
,”
J. Chem. Theory Comput.
16
,
4298
4314
(
2020
).
56.
F.
Ma
,
W.
Purwanto
,
S.
Zhang
, and
H.
Krakauer
, “
Quantum Monte Carlo calculations in solids with downfolded Hamiltonians
,”
Phys. Rev. Lett.
114
,
226401
(
2015
).
57.
W.
Purwanto
,
S.
Zhang
, and
H.
Krakauer
, “
Frozen-orbital and downfolding calculations with auxiliary-field quantum Monte Carlo
,”
J. Chem. Theory Comput.
9
,
4825
4833
(
2013
).
58.
S.
Zhang
,
F. D.
Malone
, and
M. A.
Morales
, “
Auxiliary-field quantum Monte Carlo calculations of the structural properties of nickel oxide
,”
J. Chem. Phys.
149
,
164102
(
2018
).
59.
D.
McQuarrie
,
Statistical Mechanics
(
University Science Books
,
2000
).
60.
J. F.
Annett
,
Superconductivity, Superfluids and Condensates
, Oxford Master Series in Physics (
Oxford University Press
,
Oxford
,
2004
).
61.
W. J.
Mullin
and
J. P.
Fernández
, “
Bose–Einstein condensation, fluctuations, and recurrence relations in statistical mechanics
,”
Am. J. Phys.
71
,
661
669
(
2003
).
62.
B. M.
Rubenstein
,
S.
Zhang
, and
D. R.
Reichman
, “
Finite-temperature auxiliary-field quantum Monte Carlo technique for Bose–Fermi mixtures
,”
Phys. Rev. A
86
,
053606
(
2012
).
63.
W. E.
Ormand
,
D. J.
Dean
,
C. W.
Johnson
,
G. H.
Lang
, and
S. E.
Koonin
, “
Demonstration of the auxiliary-field Monte Carlo approach for sd–shell nuclei
,”
Phys. Rev. C
49
,
1422
1427
(
1994
).
64.
R. D.
Sedgewick
,
D. J.
Scalapino
,
R. L.
Sugar
, and
L.
Capriotti
, “
Canonical and grand canonical ensemble expectation values from quantum Monte Carlo simulations
,”
Phys. Rev. B
68
,
045120
(
2003
).
65.
C. N.
Gilbreth
,
S.
Jensen
, and
Y.
Alhassid
, “
Reducing the complexity of finite-temperature auxiliary-field quantum Monte Carlo
,” arXiv:1907.10596 [physics.comp-ph] (
2019
).
66.
S.
Rombouts
and
K.
Heyde
, “
An accurate and efficient algorithm for the computation of the characteristic polynomial of a general square matrix
,”
J. Comput. Phys.
140
,
453
458
(
1998
).
67.
Z.
Wang
,
F. F.
Assaad
, and
F.
Parisen Toldin
, “
Finite-size effects in canonical and grand-canonical quantum Monte Carlo simulations for fermions
,”
Phys. Rev. E
96
,
042131
(
2017
).
68.
J. E.
Hirsch
, “
Discrete Hubbard–Stratonovich transformation for fermion lattice models
,”
Phys. Rev. B
28
,
4059
4061
(
1983
).
69.
G. M.
Buendia
, “
Comparative study of the discrete and the continuous Hubbard–Stratonovich transformation for a one-dimensional spinless fermion model
,”
Phys. Rev. B
33
,
3519
3521
(
1986
).
70.
H. F.
Trotter
, “
On the product of semi-groups of operators
,”
Proc. Am. Math. Soc.
10
,
545
551
(
1959
).
71.
M.
Suzuki
, “
Relationship between d-dimensional quantal spin systems and (d + 1)-dimensional ising systems: Equivalence, critical exponents and systematic approximants of the partition function and spin correlations
,”
Prog. Theor. Phys.
56
,
1454
1469
(
1976
).
72.
S.
Zhang
, “
Ab initio electronic structure calculations by auxiliary-field quantum Monte Carlo
,” arXiv:1807.06688 [cond-mat.str-el] (
2018
).
73.
H.
Shi
and
S.
Zhang
, “
Symmetry in auxiliary-field quantum Monte Carlo calculations
,”
Phys. Rev. B
88
,
125132
(
2013
).
74.
P.
Borrmann
and
G.
Franke
, “
Recursion formulas for quantum statistical partition functions
,”
J. Chem. Phys.
98
,
2484
(
1993
).
75.
P.
Borrmann
,
J.
Harting
,
O.
Mülken
, and
E. R.
Hilf
, “
Calculation of thermodynamic properties of finite Bose–Einstein systems
,”
Phys. Rev. A
60
,
1519
(
1999
).
76.
B.
Hirshberg
,
V.
Rizzi
, and
M.
Parrinello
, “
Path integral molecular dynamics for bosons
,”
Proc. Natl. Acad. Sci. U. S. A.
116
,
21445
21449
(
2019
).
77.
B.
Hirshberg
,
M.
Invernizzi
, and
M.
Parrinello
, “
Path integral molecular dynamics for fermions: Alleviating the sign problem with the Bogoliubov inequality
,”
J. Chem. Phys.
152
,
171102
(
2020
).
78.
A. L.
Fetter
and
J. D.
Walecka
,
Quantum Theory of Many-Particle Systems
(
McGraw-Hill
,
Boston
,
1971
).
79.
K.
Schönhammer
, “
Deviations from Wick’s theorem in the canonical ensemble
,”
Phys. Rev. A
96
,
012102
(
2017
).
80.
K.
Schönhammer
, “
Thermodynamics and occupation numbers of a Fermi gas in the canonical ensemble
,”
Am. J. Phys.
68
,
1032
1037
(
2000
).
81.
D.
Landau
and
K.
Binder
,
A Guide to Monte Carlo Simulations in Statistical Physics
(
Cambridge University Press
,
USA
,
2005
).
82.
T.
Shen
,“
Canonical Ensemble AFQMC for Bose/Fermi Hubbard Models
,” Zendodo, v1.0, Dataset
83.
J.
Hubbard
and
B. H.
Flowers
, “
Electron correlations in narrow energy bands
,”
Proc. R. Soc. London, Ser. A
276
,
238
257
(
1963
).
84.
J.
Hubbard
and
B. H.
Flowers
, “
Electron correlations in narrow energy bands. II. The degenerate band case
,”
Proc. R. Soc. London, Ser. A
277
,
237
259
(
1964
).
85.
T.
Esslinger
, “
Fermi-Hubbard physics with atoms in an optical lattice
,”
Annu. Rev. Condens. Matter Phys.
1
,
129
152
(
2010
).
86.
M.
Greiner
,
O.
Mandel
,
T.
Esslinger
,
T. W.
Hänsch
, and
I.
Bloch
, “
Quantum phase transition from a superfluid to a Mott insulator in a gas of ultracold atoms
,”
Nature
415
,
39
44
(
2002
).
87.
M. P. A.
Fisher
,
P. B.
Weichman
,
G.
Grinstein
, and
D. S.
Fisher
, “
Boson localization and the superfluid–insulator transition
,”
Phys. Rev. B
40
,
546
570
(
1989
).
88.
M.
Troyer
and
U.-J.
Wiese
, “
Computational complexity and fundamental limitations to fermionic quantum Monte Carlo simulations
,”
Phys. Rev. Lett.
94
,
170201
(
2005
).
89.
W.
Magnus
and
F.
Brosens
, “
Occupation numbers in a quantum canonical ensemble: A projection operator approach
,”
Physica A
518
,
253
264
(
2019
).
90.
R.
Blankenbecler
,
D. J.
Scalapino
, and
R. L.
Sugar
, “
Monte Carlo calculations of coupled boson-fermion systems. I
,”
Phys. Rev. D
24
,
2278
2286
(
1981
).
91.
C. N.
Gilbreth
and
Y.
Alhassid
, “
Pair condensation in a finite trapped Fermi gas
,”
Phys. Rev. A
88
,
063643
(
2013
).
92.
H. R.
Petras
,
S. K.
Ramadugu
,
F. D.
Malone
, and
J. J.
Shepherd
, “
Using density matrix quantum Monte Carlo for calculating exact-on-average energies for ab initio Hamiltonians in a finite basis set
,”
J. Chem. Theory Comput.
16
,
1029
1038
(
2020
).
93.
K. T.
Williams
,
Y.
Yao
,
J.
Li
,
L.
Chen
,
H.
Shi
,
M.
Motta
,
C.
Niu
,
U.
Ray
,
S.
Guo
,
R. J.
Anderson
,
J.
Li
,
L. N.
Tran
,
C.-N.
Yeh
,
B.
Mussard
,
S.
Sharma
,
F.
Bruneval
,
M.
van Schilfgaarde
,
G. H.
Booth
,
G. K.-L.
Chan
,
S.
Zhang
,
E.
Gull
,
D.
Zgid
,
A.
Millis
,
C. J.
Umrigar
, and
L. K.
Wagner (Simons Collaboration on the Many-Electron Problem)
, “
Direct comparison of many-body methods for realistic electronic Hamiltonians
,”
Phys. Rev. X
10
,
011041
(
2020
).
94.
M.
Motta
,
D. M.
Ceperley
,
G. K.-L.
Chan
,
J. A.
Gomez
,
E.
Gull
,
S.
Guo
,
C. A.
Jiménez-Hoyos
,
T. N.
Lan
,
J.
Li
,
F.
Ma
,
A. J.
Millis
,
N. V.
Prokof’ev
,
U.
Ray
,
G. E.
Scuseria
,
S.
Sorella
,
E. M.
Stoudenmire
,
Q.
Sun
,
I. S.
Tupitsyn
,
S. R.
White
,
D.
Zgid
, and
S.
Zhang (Simons Collaboration on the Many-Electron Problem)
, “
Towards the solution of the many-electron problem in real materials: Equation of state of the hydrogen chain with state-of-the-art many-body methods
,”
Phys. Rev. X
7
,
031059
(
2017
).
95.
T.
Shen
(
2020
). “
Canonical ensemble AFQMC for Bose/Fermi Hubbard models data
,” Zenodo.
96.
D. R.
Hamann
and
S. B.
Fahy
, “
Energy measurement in auxiliary-field many-electron calculations
,”
Phys. Rev. B
41
,
11352
(
1990
).
97.
H.
Schmidt
, “
Eine einfache herleitung der verteilungsfunktionen für bose-und fermi-statistik
,”
Z. Phys.
134
,
430
431
(
1953
).
98.
H.
Schmidt
, “
A simple derivation of distribution functions for Bose and Fermi statistics
,”
Am. J. Phys.
57
,
1150
1151
(
1989
).

Supplementary Material