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.

## I. INTRODUCTION

For decades, chemical physicists have focused on developing a constellation of electronic structure methods, from mean field^{1,2} to perturbation^{1,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 technologies^{9–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–Fock^{14,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 techniques^{18–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 baths^{29} such as Dynamical Mean Field Theory (DMFT)^{30,31} and the finite temperature SEET and Second-order Green’s Function methods^{32–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 problems^{46} 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-filling^{46,50,51} and many different *ab initio* molecular^{44,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*(*N*^{3})–*O*(*N*^{4}), 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) transformation^{38,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.

## II. METHODS

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

### A. Review of determinant quantum Monte Carlo (DQMC) in the grand canonical ensemble

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 DQMC^{39,41} and AFQMC^{48,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,

Here, $\u0124$ denotes the many-body Hamiltonian, *μ* denotes the chemical potential, $N^=\u2211i,\sigma \u0109i\sigma \u2020\u0109i\sigma $ denotes the particle number operator, and Δ*τ* = *β*/*L*. In order to facilitate its subsequent evaluation, the propagator, $e\u2212\Delta \tau (\u0124\u2212\mu N^)$, is factored into short imaginary time kinetic and potential propagators via a Suzuki–Trotter factorization^{70,71} such that

As discussed in Sec. III A, in this work, we focus on models whose Hamiltonians can be written as $\u0124=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\u2212\Delta \tau K^$, may be neatly expressed as matrices in a given basis,^{40} two-body propagators, $e\u2212\Delta \tau V^$, may not be as easily expressed. Fortunately, two-body propagators of the form $e\u2212\Delta \tau 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}

where $v^\gamma $ 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) transformation^{69}

where *ϕ* is the so-called “auxiliary field.” As $e\u2212\varphi 2/2/2\pi $ 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,

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

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(\u220flLB^(\varphi \u2192l))$ are avoided, and we obtain

where taking the trace over the product of $B^(\varphi \u2192l)$ leads to a product of determinants, one from each of the two spin sectors;^{40} $B\sigma (\varphi \u2192l)$ is the matrix representation of $B^\sigma (\varphi \u2192l)$ 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.

### B. Obtaining canonical ensemble properties from the grand canonical formalism

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\u2212\beta \u0124$ 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, $\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}

where *ψ*_{m} ≡ 2*πm*/*N*_{s} 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.

### C. Canonical ensemble AFQMC algorithm for spinless bosons and fermions

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 $\xe2\u2020$ and $\xe2$ 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

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

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 proven^{40} that the product of $B^(\varphi \u2192l)$ exactly corresponds to an exponential of some one-body operator parameterized by the auxiliary field vector $\varphi \u2192\u2261{\varphi \u21921,\u2026,\varphi \u2192L}$, which we denote as $\xc2(\varphi \u2192)$,

We can subsequently define a set of new boson/fermion coordinates in which $\xc2(\varphi \u2192)$ is diagonal,^{40}

where we call the set of ${\epsilon \u0303\gamma (\varphi \u2192)}$ the effective single-particle spectrum because of the following relation:

and

Since $e\u2212\xc2(\varphi \u2192)$ is an independent-particle propagator that only depends on the auxiliary field vector, $\varphi \u2192$, the effective single-particle spectrum, ${\epsilon \u0303\gamma (\varphi \u2192)}$, is independent of the particle number. For an *N*-particle, *N*_{s}-site system, taking the trace while constraining the particle number yields

Here, Γ is used to represent the set of *N*-particle states, and thus, $\u2211\Gamma \u2261\u2211n1+\cdots +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 $\varphi \u2192$, 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 *Z*_{N} in terms of its subsystem partition functions, *Z*_{N−k},

with *Z*_{0} = 1 and *z*_{k}, which can be identified as the partition function of the system at temperature *kβ*, given by

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\u2212\beta \xc2(\varphi \u2192))$ as $ZN(\varphi \u2192)$ and re-express it using the same recursive formalism

with the same starting value $Z0(\varphi \u2192)=1$ and

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

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

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\u0303ij$. 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

As is done in Eq. (16), we can simplify this equation by applying a change of basis to $\xe2i\u2020\xe2j\u2020$ to yield

where $U\lambda \mu ij\u2261\lambda |i\u2009j|\mu $ is an overlap matrix. This formula allows the integrand to be expressed as

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

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

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,

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

where

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

Observe that the expression for $\u27e8n^\lambda 2\u27e9N$ 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.

### D. Canonical ensemble AFQMC algorithm for spinful fermions

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

which means that spin up and spin down states are propagated separately, and we can accordingly define two one-body operators, $\xc2\u2191(\varphi \u2192)$ and $\xc2\u2193(\varphi \u2192)$, such that

and

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,

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

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, $\varphi \u2192$. 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),

with the recursive relation being

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

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

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.

### E. Computational details

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(\varphi \u2192)$ 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

Some reordering yields

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^=\u2212\u2211i,jtij\xe2i\u2020\xe2j$ and plug in Eq. (35) to yield

with *n* being the number of Monte Carlo samples and $\Phi \u2192$ being sampled from a Boltzmann-like distribution

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(\Phi \u2192\u2032)$ to the old partition function $ZN(\Phi \u2192)$, the acceptance ratio of a random flip of discrete auxiliary fields or a random perturbation to continuous auxiliary fields may be written as

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, ${\epsilon \u0303\gamma (\varphi \u2192)}$, 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 × 10^{5} 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}

## III. RESULTS AND DISCUSSION

### A. Benchmarks against exact diagonalization results for the Bose and Fermi–Hubbard models

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

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,

where $b^i\u2020$/$b^i$ denote boson creation/annihilation operators at site *i* and $n^i=b^i\u2020b^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 lattices^{86} 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, *E*_{k}, which are one-body quantities, and potential energies, *E*_{p}, 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/*k*_{B}*T*, 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.

### B. Emergence of the sign problem

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\u2212\beta \epsilon \lambda (\varphi \u2192)$ appears in the recursive calculation, its corresponding mean occupation number, $\u27e8n\lambda \u27e9N$, tends to approach 1, yielding $\u27e81\u2212n\lambda \u27e9N$ 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, *N*_{c}, at which $\u27e8n\lambda \u27e9Nc$ begins to exceed 1 for the *λ* that corresponds to the largest eigenvalue of $e\u2212\beta A(\varphi \u2192)$, 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 *N*_{c} 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.

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\u2212\beta \epsilon \lambda (\varphi \u2192)$’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\u2212\beta \epsilon \lambda (\varphi \u2192)$’s on the other $e\u2212\beta \epsilon \lambda (\varphi \u2192)$’s must be reduced. To accomplish this, note that $\epsilon N(\varphi \u2192)$ 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 $\epsilon \lambda l(\varphi \u2192)$ (where *l* denotes lower) are set to 1 ($\u27e8n\lambda l\u27e9N=1$) if $e\u2212\beta (\epsilon N(\varphi \u2192)\u2212\epsilon \lambda l(\varphi \u2192))>\xi $, where *ξ* is a cutoff value. A similar approximation can then be made by setting $\u27e8n\lambda u\u27e9N=0$ for energy levels $\epsilon \lambda u(\varphi \u2192)$ (where *u* denotes upper) well above the Fermi level such that $e\u2212\beta (\epsilon \lambda u(\varphi \u2192)\u2212\epsilon N(\varphi \u2192))>\xi $. 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, *ξ* = 10^{4}, 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 *ξ* = 10^{3} 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(\varphi \u2192)$ as a function of filling, $\u27e8n\u27e9=N\u2191+N\u2193Ns$. We chose *ξ* = 10^{3} 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.

. | ED . | QMC, ξ = 10^{3}
. | QMC, ξ = 10^{4}
. |
---|---|---|---|

N_{s} = 6, ⟨n⟩ = 1 | −3.6429 | −3.61(3) | −3.67(6) |

N_{s} = 8, ⟨n⟩ = 0.75 | −6.6298 | −6.61(1) | −6.64(5) |

N_{s} = 8, ⟨n⟩ = 1 | −4.5002 | −4.44(3) | −4.51(5) |

. | ED . | QMC, ξ = 10^{3}
. | QMC, ξ = 10^{4}
. |
---|---|---|---|

N_{s} = 6, ⟨n⟩ = 1 | −3.6429 | −3.61(3) | −3.67(6) |

N_{s} = 8, ⟨n⟩ = 0.75 | −6.6298 | −6.61(1) | −6.64(5) |

N_{s} = 8, ⟨n⟩ = 1 | −4.5002 | −4.44(3) | −4.51(5) |

### C. Convergence to the ground state

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 GCE^{44,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/*k*_{B}*T* units. Although this may seem small, a few 1/*k*_{B}*T* 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.

### D. Scaling relative to the grand canonical AFQMC algorithm

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\u2212\beta A(\varphi \u2192)$, 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 *N* ≫ *N*_{s}, 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.

## IV. CONCLUSIONS

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 condensates^{62} as well as superconductors^{91} 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 QMC^{92} 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 MATERIAL

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

## ACKNOWLEDGMENTS

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.

## DATA AVAILABILITY

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}

### APPENDIX: DETAILED DERIVATION OF PARTITION FUNCTIONS AND UNNORMALIZED DENSITY MATRICES

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,

We then consider the operations $\xc2bb^\u2020$ and $\xc2f\u0109\u2020$. Using the elementary commutation/anticommutation relations, it is straightforward to show that

which gives

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 $\xc2$ is a scalar relative to the *N* × *N* matrices. As a result, $I\xc2$ and *A* commute as matrices. Repeatedly applying these two equations yields

for any positive integer *n*. Hence,

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\u2212\beta \xc2(\varphi \u2192)$ 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

where, in the last step, we have used the fact $e\u2212\xc2(\varphi \u2192)|0=|0$. Thus, sandwiching $e\u2212\xc2(\varphi \u2192)$ between ⟨Γ| and |Γ⟩ yields

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,

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

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

The procedure leading to Eq. (A14) can be easily extended to calculate expectation values of products of different occupation numbers, i.e., $\u27e8ninj\u27e9N$ with *i* ≠ *j*. Replacing $n^1$ with $n^1n^2$ and observing that the enumerations of ${n3,\u2026,nNs}$ are independent from {*n*_{1} = 0, *n*_{2} = 0}, one obtains

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:

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

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

These are precisely the equations referenced in the text.

## REFERENCES

_{2}on Au