We investigate the viability of the phaseless finite-temperature auxiliary-field quantum Monte Carlo (ph-FT-AFQMC) method for *ab initio* systems using the uniform electron gas as a model. Through comparisons with exact results and FT coupled cluster theory, we find that ph-FT-AFQMC is sufficiently accurate at high to intermediate electronic densities. We show, both analytically and numerically, that the phaseless constraint at FT is fundamentally different from its zero-temperature counterpart (i.e., ph-ZT-AFQMC), and generally, one should not expect ph-FT-AFQMC to agree with ph-ZT-AFQMC in the low-temperature limit. With an efficient implementation, we are able to compare exchange-correlation energies to the existing results in the thermodynamic limit and find that the existing parameterizations are highly accurate. In particular, we found that ph-FT-AFQMC exchange-correlation energies are in better agreement with a known parameterization than is restricted path-integral MC in the regime of Θ ≤ 0.5 and *r*_{s} ≤ 2, which highlights the strength of ph-FT-AFQMC.

## I. INTRODUCTION

Temperature-dependent properties of interacting fermions are fundamentally important in both experiments and theory. Typical phenomena at finite temperature (FT) include the Bardeen Cooper Schrieffer (BCS)-Bose-Einstein Condensation (BEC) crossover in the attractive 2D Fermi model,^{1} the competition between stripe and superconducting orders in the 2D Hubbard model,^{2} and plasmonic catalysis.^{3} Furthermore, there is a growing interest in warm dense matter,^{4} an extreme state of matter found in planetary interiors^{5} that can be created with high intensity lasers.^{6,7} Understanding such phenomena at the theoretical level is challenging due to the delicate interplay between electron–electron, electron–ion, quantum mechanical, and thermal effects, all of which can be equally important and often cannot be treated perturbatively. Thus, accurate computational approaches are required that are capable of capturing these effects.

Density functional theory (DFT), as the workhorse of zero-temperature (ZT) electronic structure theory, is an ideal candidate: it is relatively cheap and accurate, it can be coupled with molecular dynamics to include ionic effects, and it can be rigorously formulated to incorporate thermal electronic effects.^{8} Indeed, DFT has proven itself effective in simulating warm dense matter.^{9–12} However, questions remain regarding the accuracy of using approximations made in thermal DFT, including the accuracy of the use of zero-temperature exchange-correlation functionals.^{13–16} Complementary approaches are therefore desired that can benchmark or supplement DFT results when necessary.

Quantum Monte Carlo (QMC) methods are a promising class of such computational methods. They are, in principle, exact in a finite supercell and can explicitly include many-body and thermal effects in an unbiased way. Of the many flavors of QMC, real space path-integral Monte Carlo (PIMC) is perhaps the best known and well established.^{17} PIMC, as a real space approach, works in the complete basis set limit, which is a considerable benefit at very high electronic temperatures. However, like all fermionic QMC methods, it suffers from the sign problem that can only be overcome using the uncontrolled restricted path approximation,^{18} which leads to the restricted PIMC (RPIMC) approach. This approach is similar in spirit to the fixed-node approximation in diffusion Monte Carlo (DMC),^{19} but now a constraint is enforced using a trial thermal density matrix. The quality of this constraint is *a priori* unknown; however, results with free fermion nodes in the uniform electron gas (UEG) suggest that it is unreliable at high densities and at lower temperatures.^{20–22} We note there have been promising developments in extending the scope of PIMC to higher densities and lower temperatures through algorithmic developments^{23} including the development of new algorithms such as the, in principle, exact real space permutation blocking PIMC^{24} (PB-PIMC) and second quantized configuration PIMC^{25} (CPIMC). In particular, PB-PIMC and CPIMC have been used as complementary approaches to simulate the warm dense UEG above half the Fermi temperature, with the fermion sign problem preventing simulations below this. We note that restricted CPIMC may be a promising route to access lower temperature, given its results on small UEG supercells and the 33-electron supercell spin-polarized case.^{26} Other interesting QMC approaches are the quantum chemistry inspired methods such as density matrix quantum Monte Carlo^{27–29} (DMQMC) and Krylov-projected full configuration interaction quantum Monte Carlo.^{30} Like CPIMC, these approaches work in a second quantized space and work well at high densities; however, they often struggle to reach the complete basis set limit. All of these methods offer unbiased exact thermal expectation values, in principle; however, they all ultimately scale exponentially with the number of electrons, in general. See Ref. 31 for a review of the parameter regimes accessible to these new QMC methods.

Recently, there has also been considerable interest in developing finite-temperature deterministic quantum chemistry methods that work in a finite basis set. These methods include second order perturbation theory,^{32–35} coupled cluster theory,^{36–38} and thermofield theory.^{39} These are promising and offer a systematic approach for including electronic temperature effects, but they often struggle to reach the continuum limit [also known as the complete basis set (CBS) limit] due to the steep computational scaling with respect to the number of basis functions, *M*. For example, FT-CCSD scales like $O(M6)$, which becomes prohibitively expensive to converge results to the CBS limit.^{37,38}

Finite-temperature auxiliary-field QMC^{40,41} (FT-AFQMC) is another promising QMC method. It works in the second quantized space and thereby suffers from the basis set incompleteness errors common to DMQMC and CPIMC. However, unlike DMQMC and CPIMC, it can be made to scale polynomially with system size at the cost of introducing an uncontrolled bias called the phaseless approximation.^{42,43} Moreover, zero-temperature phaseless AFQMC (ph-ZT-AFQMC) has proven itself as one of the most accurate and scalable post Hartree–Fock methods.^{44–48} In addition, the bias introduced by the phaseless approximation is typically much smaller than the fixed-node error in DMC,^{49} which, in principle, should serve as a rough upper bound to the bias in RPIMC. Thus, it is important to assess the quality of ph-FT-AFQMC for realistic systems as to date, the applications have largely focused on model systems^{50,51} or have not enforced the constraint,^{52–54} which is not a practical approach as the system size increases. Compared to FT-CCSD, ph-FT-AFQMC maintains favorable cubic scaling $O(M3)$ for each statistical sample,^{55,56} which makes it better-suited for large-scale warm dense matter simulations.

In this paper, we investigate the viability of ph-FT-AFQMC both in terms of its accuracy and its ability to reach the complete basis set limit. We stress that good performance for *both* of these metrics is critical if the method is to be practically useful for finite-temperature *ab initio* systems. To address these problems, we use the uniform electron gas (UEG) model as a test bed for ph-FT-AFQMC. To the best of our knowledge, our work is the first to apply ph-FT-AFQMC at finite temperature beyond lattice models with short-range interactions. The main goal of our work is to assess the accuracy of ph-FT-AFQMC when applied to the UEG model.

Apart from being a foundational model in condensed matter physics,^{57–60} the UEG offers a number of useful features from a computational point of view. First, the model can be tuned from weak to strong correlation as a function of the density parameter, *r*_{s}, given that our previous work^{61} has shown that ph-ZT-AFQMC is highly accurate for *r*_{s} ≤ 3–4 at zero temperature and we have an idea of the magnitude of errors we might expect. Second, basis set convergence can be easily investigated in a plane wave basis set by increasing the energy cutoff. Finally, the UEG at warm dense matter conditions has been the subject of intense study over the past decade.^{62} In fact, many of the recent developments^{62} in finite-temperature fermionic QMC methods were spurred on by a discrepancy between RPIMC and CPIMC^{20,21} results for the warm dense UEG that has been incorporated into finite-temperature exchange-correlation functionals.^{63–65} Because of this effort, there is a considerable amount of essentially exact data for energetic,^{66,67} static,^{68–70} and dynamic properties^{71,72} of the model in a wide range of densities and temperatures. Nonetheless, despite this immense effort, there is a gap in accurate data below roughly half the Fermi temperature below, which no method could reach due to the Fermionic sign problem.^{73} In this paper, we use ph-FT-AFQMC to partially fill this gap.

This paper is organized as follows: In the first section, we outline the ph-FT-AFQMC method, paying careful attention to how to implement it efficiently in a numerically stable fashion. Next, we carefully benchmark the method at various temperatures against exact results in small basis sets and investigate the low-temperature limit in comparison to the ph-ZT-AFQMC results. We also compare the ph-FT-AFQMC results to FT-CCSD and PIMC where possible. Finally, we compare our results against other QMC methods and available exchange-correlation parameterization in the thermodynamic limit and finish by outlining our perspective for the future of the method.

## II. THEORY

We briefly summarize the phaseless approximation of AFQMC within the finite-temperature formalism (i.e., ph-FT-AFQMC).

### A. Finite-temperature AFQMC in the grand canonical ensemble

#### 1. General formalisms

The finite-temperature AFQMC algorithm aims to compute thermal expectation values based on the grand canonical partition function,

where $N^$ is a total number operator, *μ* is a chemical potential, and the Hamiltonian involves one-body ($\u01241$) and two-body ($\u01242$) terms,

The direct (deterministic) evaluation of trace in Eq. (1) scales exponentially since there are exponentially many states to consider. It is then natural to consider QMC algorithms to sample an instance of the terms in *Z*.

A particular flavor of QMC that we focus on here is AFQMC, where the two-body propagator is expressed by the Hubbard–Stratonovich transformation,^{74}

where $v^$ is related to the two-body Hamiltonian as a sum of squared operators,

where **x** is a vector of *N*_{α} auxiliary fields that are samples from the standard normal distribution, *p*(**x**). With the symmetric Trotter decomposition, the total propagator reads

where $B^$ is defined as

For a given number of imaginary time slices *n*, we sample the auxiliary fields via MC,

where *p*(**x**_{1}, **x**_{2}, …, **x**_{n}) is the probability of sampling a specific path designated by auxiliary fields, **x**_{1}, **x**_{2}, …, **x**_{n}, and the time step, Δ*τ* = *β*/*n*, is determined by the temperature and the number of imaginary time slices. The evaluation of the trace in Eq. (7) is still difficult because it needs to consider all possible states in the Hilbert space despite the fact that every operator inside the trace is a one-body operator. One can show analytically that the trace can be written in terms of a determinant in the grand canonical ensemble,^{40,75,76}

where **B** is a matrix representation of $B^$ in the single-particle basis. For a later use, we define the product of **B** as

where we omit *μ* in the argument of **A** for simplicity. Note that for now, we defined **A** only at the end of the imaginary time propagation. Later, we will define **A** along the trajectory for an arbitrary imaginary time *τ*.

We are mostly interested in computing expectation values based on the partition function in Eq. (8), not the partition function itself,

The computation of expectation values can be easily achieved through an importance sampling procedure. To see this, we first write,

where **X** denotes the set of auxiliary fields along each imaginary path. The field configurations **X** are sampled by the MC algorithm where we write expectation values,

where the local expectation value is defined as

and walker weights *w*_{i} are updated via the importance sampling procedure based on the distribution, $tr(\xc2(X))/Z$.

In AFQMC, any expectation values are expressed as a function of the one-body Green’s function,

following the generalized Wick’s theorem.^{77} Therefore, computing **G** is sufficient to compute any expectation values. In terms of **B**(Δ*τ*, **x**, *μ*), it can be shown that a sample of the one-body Green’s function (i.e., the local quantity) is^{75,76}

We note that a sample of the one-body reduced density matrix (1RDM), **P**, is obtained by

Using Wick’s theorem, higher order reduced density matrices can be computed as products of the one-body Green’s function.

#### 2. Phaseless approximation

Since the value of determinants in Eq. (8) or equivalently $tr(\xc2(X))$ for a given **X** can be either positive or negative (in our case complex),^{78} the phase problem naturally arises in the finite-temperature algorithm. Therefore, it is important to impose the phaseless constraint^{43,51} to remove the phase problem similar to the case of the ph-ZT-AFQMC algorithm.^{42} This way, we can keep the overall scaling of the algorithm polynomial.

We introduce a trial density matrix, $B^T(\Delta \tau ,\mu T)$, and this is used to impose a phaseless constraint in the imaginary time evolution,^{43} which is defined as

where $\u0124T$ is some one-body operator and *μ*_{T} is the chemical potential for the trial density matrix. Throughout the paper, we will assume that *μ*_{T} is tuned so that the thermal one-body density matrix from $B^T(n\Delta \tau ,\mu T)$ has a trace of *N* with *N* being the desired number of particles.

We implement the imaginary time evolution of a path by the following algorithm. The central quantity in the constrained evolution is **A**. At an imaginary time *τ* = *k*Δ*τ*, we define

At each imaginary time step, we replace **B**_{T} in the middle with **B**, that is,

In general, the product of **B** and/or **B**_{T} requires a special care for numerical stability, especially when simulating low temperature. This can be achieved by the stabilization method described in Ref. 79, which we will describe further in Sec. II B.

During this propagation, the importance function is determined by the FT-AFQMC overlap ratio,

For a later use, we mention that we can equivalently write this ratio in terms of trace over all possible states in the grand canonical ensemble,

In practice, we employ the “optimal” force bias, which is a shift to the Gaussian distribution as well as the mean-field subtraction that enforces the normal ordering. The optimal force bias at *τ* = *k*Δ*τ* is

and the mean-field subtraction is

We can then define the importance function (in the hybrid form^{80}) as

The phaseless approximation ensures the reality and positivity of walkers using a modified importance function,

where

Using these, the *n*-th walker weight at *τ* = *k*Δ*τ* is updated via

This completes the description of the ph-AFQMC algorithm at finite temperature.

#### 3. The *T* → 0 limit

Converging ph-FT-AFQMC to the ph-ZT-AFQMC limit as decreasing *T* is highly desirable. This is due to the remarkable accuracy of ph-ZT-AFQMC in a variety of systems benchmarked to date.^{44–48} Reaching this zero-temperature limit would naturally suggest that the ph-FT-AFQMC algorithm is accurate at low temperature. Unfortunately, reaching the zero-temperature limit is, in fact, difficult due to the differences in the phaseless constraint between two algorithms as we shall see below.

We first define the zero-temperature overlap ratio at *τ* = *k*Δ*τ*,^{42}

where *N* is the number of particles and |*ψ*_{T}(*N*)⟩ is the trial wavefunction defined as

with |*ψ*_{0}(*N*)⟩ being some initial wavefunction and

We note that we set *μ* = *μ*_{T} = 0 because the ph-ZT-AFQMC algorithm typically works in a fixed particle number space. In other words, the number of particles in |*ψ*_{T}(*N*)⟩ and |*ψ*(*N*)⟩ is *N* as specified. Therefore, there is no need to have chemical potentials. This is then used to construct the phaseless importance function in Eq. (25), which defines the phaseless approximation at *T* = 0.

The correspondence between the *T* → 0 limit of ph-FT-AFQMC and ph-ZT-AFQMC can be understood by comparing Eqs. (21) and (28). It is then ultimately equivalent to showing $tr(\xc2(\tau ,{xi}i=1k\u22121))=\psi T(N)|\psi (\tau ,{xi}i=1k\u22121,N)$ where $\xc2$ is defined through Eq. (19). In the limit of *T* → 0, the number of time slices *n* becomes *∞*. First, we consider the propagation in the middle of a path. Namely, we can assume 1 ≪ *k* ≪ *n* in this case. One starts from

where the summation over *M* is done over states in different particle number sectors and the summation over *α* can be thought of as summing over all possible states in the basis of |*ψ*_{T}(*M*)⟩. Assuming that *μ*_{T} and *μ* are chosen so that it should pick out the *N*-particle sector, the action of $(B^T(\Delta \tau ,\mu T))n\u2212k+1$ to the left when *n* → *∞* (as well as 1 ≪ *k* ≪ *n*) yields only one term out of the summation, that is, using Eq. (29),

This then shows that the phaseless constraints are equivalent between ph-ZT-AFQMC and ph-FT-AFQMC algorithms when 1 ≪ *k* ≪ *n*.

Subtleties arise when 1 ≪ *k* ≪ *n* does not hold. In the ph-FT-AFQMC algorithm, such a case always happens toward the completion of a path. As an extreme case, let us consider the phaseless constraint at the last imaginary time step in the ph-FT-AFQMC algorithm. Namely, let us set *k* = *n*. With the same assumptions about *μ* and *μ*_{T}, we have

Not only does the summation over *α* not easily truncates but also this limit no longer corresponds to the zero-temperature limit unlike in Eq. (32). This is indeed why one should not expect the ph-FT-AFQMC energy to approach the ph-ZT-AFQMC energy, in general. Because of these issues, even the importance function may need to be modified, while the same importance function was advocated in Ref. 52. Further investigation on this in the future will be interesting.

Finally, there may be an additional complication when our assumptions about *μ* and *μ*_{T} do not hold exactly. In other words, it is possible to have fluctuation in the number of particles along an imaginary path. Such a fluctuation is not simply due to the fact that we are working in the grand canonical ensemble. Instead, it is due to the fact that the chemical potential used in $B^T$ (i.e., *μ*_{T}) is not necessarily the same as that of the many-body chemical potential used in $B^$. Therefore, on average, at a given imaginary time *τ*, **A** may not yield the same number of particles as our target *N*. However, this is generally only secondary compared to the issue discussed in Eq. (33) (i.e., the issue of not projecting to *α* = 0) unless the underlying system has a near-zero gap. When the system is metallic, the number of particle is very sensitive to both chemical potentials in $B^T$ and $B^$. Therefore, along an imaginary propagation, the number particle keeps changing even at very low temperature. Such a subtlety arises in other methods that work in the grand canonical ensemble such as low-order perturbation theory, which has been a subject of active research for some time.^{35} Nevertheless, when the gap is not so small, in the limit of *T* → 0, the number of particles changes very little as a function of *μ*. This makes the effect of particle fluctuation very small. In principle, one can work directly in the canonical ensemble that removes the need for the assumptions about *μ*_{T} and *μ*.^{54} Nonetheless, the illustration of Eq. (33) still holds and some modification to the constraint is necessary in the context as well.

These simple illustrations suggest that one may impose the phaseless constraint for *T* > 0 based on a modified overlap ratio (at *τ* = *k*Δ*τ*) (with the same assumptions about *μ*_{T} and *μ*),

where the trial density matrix is always multiplied up to the *n*-th power so that even in the last imaginary time step, one recovers the zero-temperature limit properly. However, the overall temperature in Eq. (34) is always lower than the physical temperature, which may break down in higher temperature regimes. Moreover, examining this constraint would be an interesting research topic in the future, but for the purpose of this work, we report numerical results based on the constraint in Eq. (20).

### B. Numerical stabilization

It is well known that the standard determinant QMC (DQMC) algorithm suffers from numerical instabilities resulting from the repeated multiplication of the **B** matrices.^{81} This issue can be overcome using the stratification method,^{79} which we will now briefly describe. We first write an **A** matrix via a column-pivoted QR decomposition (QRCP),

where **Q** is an orthogonal matrix, **R** is an upper-triangular matrix, and **Π** is a permutation matrix. We then define

where

The QRCP decomposition is more expensive than matrix multiplications, so we perform the decomposition once in a while, and this frequency is controlled by a “stack” size parameter *L*. The number of stacks is then specified by *n*_{stack} = *β*/*L*. The stack size *L* is set such that one can perform the product of **B***L*/Δ*τ* times without numerical instability.

At an imaginary time *τ*, we are interested in computing

The product of stack {**A**_{i}} needs to be performed following the stratification algorithm:

Compute the QRCP A_{1} = Q_{1}D_{1}T_{1} |

for ( do 2 ≤ i ≤ n_{stack}) |

Compute C_{i} = (A_{i}Q_{i−1})D_{i−1} |

Compute the QRCP C_{i} = Q_{i}R_{i}Π_{i} |

Set D_{i} = diag(R_{i}) |

Compute $Ti=(Di\u22121Ri)(\Pi iTi\u22121)$ |

end for |

Compute the QRCP A_{1} = Q_{1}D_{1}T_{1} |

for ( do 2 ≤ i ≤ n_{stack}) |

Compute C_{i} = (A_{i}Q_{i−1})D_{i−1} |

Compute the QRCP C_{i} = Q_{i}R_{i}Π_{i} |

Set D_{i} = diag(R_{i}) |

Compute $Ti=(Di\u22121Ri)(\Pi iTi\u22121)$ |

end for |

In the end of this algorithm, we achieve

Finally, the Green’s function is also computed via a numerically stable form,

where

and

This can be derived from

### C. Exploiting the stack structure

It is possible to reuse the QDT factorization of all stacks but one which is the central stack we are propagating. We write **A**(*τ*) as a product of two matrices (left and right blocks),

We note that the QDT factorizations of the left and right blocks are already given from the previous propagation. Namely, we already have

Then, the numerically stable formation of **A**(*τ*) requires only two QRCP calls (as opposed to *n*_{stack} calls described previously) and also far less matrix–matrix multiplications to do. The algorithm for propagating one time step is as follows:

Compute $TL=TLBT\u22121$ |

Compute C_{CR} = (BQ_{R})D_{R} |

Compute the QRCP C_{CR} = Q_{CR}R_{CR}Π_{CR} |

Set D_{CR} = diag(R_{CR}) |

Compute $TCR=(DCR\u22121RCR)(\Pi CRTR)$ |

Compute C_{LCR} = Q_{L}D_{L}T_{L}Q_{CR}D_{CR} |

Compute the QRCP C_{LCR} = Q_{LCR}R_{LCR}Π_{LCR} |

Set D_{LCR} = diag(R_{LCR}) |

Compute $TLCR=(DLCR\u22121RLCR)(\Pi LCRTCR)$ |

Compute $TL=TLBT\u22121$ |

Compute C_{CR} = (BQ_{R})D_{R} |

Compute the QRCP C_{CR} = Q_{CR}R_{CR}Π_{CR} |

Set D_{CR} = diag(R_{CR}) |

Compute $TCR=(DCR\u22121RCR)(\Pi CRTR)$ |

Compute C_{LCR} = Q_{L}D_{L}T_{L}Q_{CR}D_{CR} |

Compute the QRCP C_{LCR} = Q_{LCR}R_{LCR}Π_{LCR} |

Set D_{LCR} = diag(R_{LCR}) |

Compute $TLCR=(DLCR\u22121RLCR)(\Pi LCRTCR)$ |

We, therefore, achieve

as in Eq. (40). The Green’s function can then be computed as before.

### D. Exploiting the low-rank structure

He and co-workers^{82} found that at low temperature, both **A**_{L} and **A**_{R} are low-rank, which can enable significant savings ($O(M/N)$). Such low-rank structures are reflected on **D**_{L} and **D**_{R} where with a certain threshold, we can approximate them as

where

We denote the rank of **d**_{L/R} to be *m*_{L/R} and show how scaling reduction can be achieved in terms of these ranks. From now on, we will only work in the reduced dimension provided by **d**_{L/R}. **d**_{L/R} is an *m*_{L/R} × *m*_{L/R} matrix that is much smaller than the original *M* × *M* matrix.

We write

where **q**_{L/R} is a matrix of dimension *M* × *m*_{L/R} and **t**_{L/R} is a matrix of dimension *m*_{L/R} × *M*. To maximize cost saving, one needs to modify the stratification algorithm further. The most efficient propagation with stratification can be done as follows (assuming that we sampled **B**):

Compute $tL=tLBT\u22121$ |

Compute c_{CR} = Bq_{R}d_{R} |

Compute c_{LCR} = d_{L}t_{L}c_{CR} |

Compute the QRCP c_{LCR} = q_{LCR}r_{LCR}π_{LCR} (q_{LCR} |

is m_{L} × m_{T} and r_{LCR} is m_{T} × m_{R} where m_{T} = min(m_{R}, m_{L})) |

Compute q_{LCR} = q_{L}q_{LCR} |

Set d_{LCR} = diag(r_{LCR}) |

Compute $tLCR=(dLCR\u22121rLCR)(\pi LCRtR)$ |

Compute $tL=tLBT\u22121$ |

Compute c_{CR} = Bq_{R}d_{R} |

Compute c_{LCR} = d_{L}t_{L}c_{CR} |

Compute the QRCP c_{LCR} = q_{LCR}r_{LCR}π_{LCR} (q_{LCR} |

is m_{L} × m_{T} and r_{LCR} is m_{T} × m_{R} where m_{T} = min(m_{R}, m_{L})) |

Compute q_{LCR} = q_{L}q_{LCR} |

Set d_{LCR} = diag(r_{LCR}) |

Compute $tLCR=(dLCR\u22121rLCR)(\pi LCRtR)$ |

We then achieve

where **q** is an *M* × *m*_{T} matrix, **d** is an *m*_{T} × *m*_{T} diagonal matrix, and **t** is an *m*_{T} × *M* matrix. Using these reduced dimension matrices and the Woodbury identity, we can compute the Green’s function at a reduced cost,

where $(d\u22121+tq)\u22121$ occurs in the dimension of *m*_{T} × *m*_{T} and other matrix multiplications are done at the cost of $O(M2mT)$. Similarly, using the matrix determinant lemma,

The determinant evaluation occurs in the dimension of *m*_{T} × *m*_{T} as well.

Furthermore, we note that the evaluation of $ImT+tqd$ needs to be done by the usual stabilization algorithm,

From this, we evaluate the determinant as well as the Green’s function.

For measurements, we use 1RDM, **P**, which is now expressed as

This **P** is also of low-rank, and this structure can be exploited to accelerate the local energy evaluation and other measurements.

### E. Uniform electron gas

The uniform electron gas (UEG) model is usually defined and solved in the plane wave basis set. We will follow this convention in this work as well. The kinetic energy operator is defined as

where **K** here is a plane wave vector. The electron–electron interaction operator is (in a spin-orbital basis)

where Ω is the volume of the unit cell. Finally, the Madelung energy *E*_{M} should be included to account for the self-interaction of the Ewald sum under periodic boundary conditions and^{21}

where *N* is the number of electrons in the unit cell and *r*_{s} is the Wigner-Seitz radius. We define the UEG Hamiltonian as a sum of these three terms,

The two-body Hamiltonian $V^ee$ needs to be rewritten as a sum of squares to employ the AFQMC algorithm. It was shown in Ref. 83 that

where

and

with a momentum transfer operator $\rho ^$ defined as

where Θ is the Heaviside step function and *E*_{cut} is the kinetic energy cutoff in the simulation. The HS operators $v^$ are now $\xc2(Q)$ and $B^(Q)$.

The local energy *ϵ*_{n}(*τ*) for the UEG then reads

where the one-body energy, *E*_{1}, is

and the two-body energy, *E*_{2}, is

with the Coulomb two-body density matrix Γ_{Q},

and the exchange two-body density matrix Λ_{Q},

The formation of Γ_{Q} costs $O(M2)$, whereas the formation of Λ_{Q} takes $O(M3)$ amount of work. Therefore, the evaluation of the exchange contribution is the bottleneck in the local energy evaluation. As noted in Ref. 84, the evaluation of the energy (and propagation) can be accelerated using fast Fourier transforms, and we found this to be slower than an algorithm based on sparse linear algebra for system sizes considered in this work. Therefore, our implementation is exclusively based on sparsity.

Unlike the ph-ZT-AFQMC algorithm, finite-temperature estimators are all “pure” as opposed to “mixed,” which does not require any special treatments to measure expectation values of operators that do not commute with $\u0124$. For instance, the one-body and two-body energies can be straightforwardly read off from the local energy evaluation.

## III. COMPUTATIONAL DETAILS

All simulations were performed with the open source PAUXY code.^{85} HANDE was used to perform FCI calculations.^{86,87} Q-Chem was used to crosscheck some of the mean-field finite-temperature calculations.^{88}

We used 640 walkers along with the “comb” population control method.^{89,90} For Θ ≤ 0.25, we found that the comb algorithm introduced significant population control biases, so we switched here to the pair-branch algorithm.^{91} Multiple time steps were employed depending on *r*_{s}, which are Δ*τ* = 0.05 for *r*_{s} ≤ 1, Δ*τ* = 0.025 for *r*_{s} = 2, and Δ*τ* = 0.005 for *r*_{s} > 2 all in reduced units. This choice was made to maintain the time step to be roughly the same throughout all *r*_{s} values in the atomic unit. Given our time step error and population control bias study in the 54-electron super cell at zero temperature,^{61} we expect that the time step and population control biases will not affect the conclusions of our study.

We typically averaged results over between 50 and 100 independent simulations, which typically were run with 160 CPUs for 10 hours. We used a free-electron trial density matrix and tuned the chemical potential *μ*_{T} such that $\u27e8N^\u27e9T=N\xaf$. The interacting chemical potential was then tuned to also match $N\xaf$. To determine the interacting chemical potential, we first use root bracketing based on a single step of the ph-FT-AFQMC algorithm (i.e., without considering the error bars) until the chemical potential was determined to about an accuracy of 1%, which we call $\mu \u0303$. Following this, we perform five poorly converged simulations (320 walkers averaged over 10 simulations) in the interval of $[0.975\mu \u0303,1.025\mu \u0303]$. We next perform a weighted least squares fit to these five data points to determine the optimal chemical potential *μ*^{*}. Finally, we ran at this *μ*^{*} using a larger walker number and average over a large number of independent simulations to get the final results reported here. We find that the procedure works most of the time at the lower temperatures considered here, although some care needs to be taken at higher temperatures where the electron number varies more significantly with *μ*. The simulation input and output are available in Ref. 92. We will primarily focus on the regime of 0.5 ≤ *r*_{s} ≤ 4 and Θ ≤ 1, i.e., the warm dense regime, where Θ = *T*/*T*_{F} and *T*_{F} is the Fermi temperature of an unpolarized UEG. Box sizes were fixed via $L=43\pi N\xaf1/3rs$ to enable comparison to results in the canonical ensemble. We note that the reduced unit system depends on the value of *r*_{s} because the Fermi temperature depends on *r*_{s} as

where *k*_{B} is the Boltzmann constant. Due to this, for a given reduced temperature Θ, we observe the total path length [i.e., 1/(*k*_{B}*T*)] in atomic unit to increase as $O(rs2)$, which makes the higher *r*_{s} values computationally costly. A low-rank threshold of 10^{−6} in Eq. (54) was used in every calculation.

## IV. RESULTS

### A. Assessing the phaseless constraint

In ph-FT-AFQMC, the main source of error is the bias introduced by the phaseless constraints. The impact of this bias is heavily dependent on the quality of trial density matrices. Here, we will employ a simple non-interacting one-body trial density matrix based on the kinetic energy operator. Namely,

which we call a free-electron (FE) trial density matrix with a trial chemical potential, *μ*_{T}.

#### 1. A study of $N\xaf=2$ with *M* = 7

To assess the quality of ph-FT-AFQMC with the FE density matrix, we can compare to exact full configuration interaction results (FT-FCI), which is possible in a very small finite basis set.^{93} We study the UEG model with only seven plane waves (*M* = 7) and tune the chemical potential to reach an on-average two-electron system.

The results of this comparison are plotted in Fig. 1 where we compare to a small system size ($N\xaf=2$) as a function of *r*_{s} and Θ. At first look, the figure looks reasonable: ph-FT-AFQMC agrees with FT-FCI very well for low *r*_{s} (high densities), which is consistent with results at zero temperature,^{61} and the results disagree more as *r*_{s} increases and as the temperature decreases. Since the *T* = 0 ph-ZT-AFQMC results are essentially exact for this system at all *r*_{s} values considered here, these results unfortunately are a manifestation of the discrepancy between ph-FT-AFQMC and ph-ZT-AFQMC at the zero-temperature limit. For example, at *r*_{s} = 3 and *T* = 0, the ph-ZT-AFQMC energy is −0.239 68(3) *E*_{h}/*e*, which is statistically identical to the FCI result of −0.239 68 *E*_{h}/*e*. This suggests that the phaseless constraint at finite temperature differs from that at zero temperature and that it is potentially considerably larger.

In Sec. II A 3, we analytically showed that two factors (constraint and chemical potential) can make the zero-temperature limit unreachable using the ph-FT-AFQMC algorithm. Incorporating other observations made in prior works,^{50} we mention a total of four possible reasons for this disagreement:

*The issue of chemical potential mismatch*: Since we work in different ensembles (FT with the grand canonical ensemble and ZT with the canonical ensemble), it is important to tune*μ*and*μ*_{T}properly so that the number of electrons does not fluctuate so much along an imaginary time path as assumed throughout the discussion in Sec. II A 3. Based on Fig. 2, we suspect that this effect is quite minor since despite the fact that the average particle number deviates from the desired value along the path, the effect of this deviation in other physical observables is negligible.*The fundamental differences between the two constraints*: A second possibility is that the constraints are fundamentally different even if*μ*and*μ*_{T}are chosen properly. This is then the direct consequence of the illustration suggested by Eq. (33) in Sec. II A 3. Indeed, as can be seen in Fig. 3, we see that this appears to be the case. Toward the end of the path, we see that the ph-FT-AFQMC results begin to deviate substantially from the FCI result, with the magnitude of the deviation increasing with*r*_{s}. This is consistent with the fact that larger*r*_{s}requires longer imaginary time to project out the ground state. This longer projection time results in a larger portion of the path using a determinant ratio that does not resemble the zero-temperature overlap ratio. In the middle of the path, results are close to exact and the algorithm more closely resembles ph-ZT-AFQMC.*The difference between the trial density matrix and trial wavefunction*: A third possibility is that the FE trial density matrix is not appropriate to reproduce the correct zero-temperature trial wavefunction used in ph-AFQMC. This is an important concern, in general;^{50}however, at the high densities and small system sizes considered here, we do not expect there to be any unrestricted Hartree–Fock solutions. Furthermore, for the UEG, RHF is equivalent to a free-electron trial in a closed shell system. Nevertheless, to verify this, we tested the thermal Hartree–Fock density matrix;^{94}however, we see from Fig. 4 that this choice makes essentially no practical difference.*Time-reversal symmetry breaking*: The final possibility is that as the phaseless constraint breaks imaginary time symmetry in the estimators^{50}and we should instead average across time slices, as suggested in Ref. 50. Interestingly, we find that at*r*_{s}= 3, this averaging procedure bias results to be above the exact value (see Fig. 5) and again does not resolve this discrepancy with the zero-temperature algorithm. This bias can also be understood in the context of Fig. 3, where the value of the energy along the imaginary time path is not symmetric with respect to*τ*. While the time-averaged estimators were found to be more accurate in the Hubbard model,^{50}given our results, one should not expect this to be a general solution to this problem.

Despite these concerns, it is also clear from Fig. 1 that ph-FT-AFQMC performs quite well with only small deviations seen from the exact result in the regimes we are interested in (i.e., Θ < 1 and *r*_{s} < 3).

#### 2. A study of $N\xaf=66$ with *M* = 57

We now study a larger UEG supercell considering $N\xaf=66$ with *M* = 57. Obviously, such a parameter set-up would have a very large basis set incompleteness error and, therefore, should not be used to draw any physically meaningful results. Nevertheless, for this basis set size, FT-CCSD results are available for several Θ and *r*_{s} values^{37} to which we compare ph-FT-AFQMC results. We use the FE trial density matrix as before since we do not have any UHF solutions at zero temperature below *r*_{s} = 3, as shown in the Appendix.

The result of this comparison is shown in Fig. 6. We see that ph-FT-AFQMC agrees well with FT-CCSD for low *r*_{s} in all temperatures considered here. As we observed in the study of $N\xaf=2$, the performance of ph-FT-AFQMC at low temperature may not be as good as the ph-ZT-AFQMC in the zero-temperature limit. Nonetheless, we found accurate results at low *r*_{s} in the case of $N\xaf=2$ and we observe consistently accurate results at a larger supercell ($N\xaf=66$) when compared against FT-CCSD. CCSD was shown to be accurate for low *r*_{s} such as *r*_{s} = 0.5 in the zero-temperature benchmark,^{95} so this helps us build confidence on the performance of ph-FT-AFQMC at low *r*_{s}.

However, the quality of FT-CCSD is expected to gradually degrade as *r*_{s} increases and Θ decreases, as indicated by its zero-temperature benchmark study.^{95} Consequently, at *r*_{s} = 4 FT-CCSD, exchange-correlation energies show erratic behavior changing the energy trend as a function of *r*_{s} completely. ph-FT-AFQMC, on the other hand, shows a smooth monotonic behavior as a function of *r*_{s} in all temperatures. Given the superior performance of ph-ZT-AFQMC compared to CCSD in the UEG model,^{61} this result is not surprising.

### B. Efficacy of the low-rank truncation

In Fig. 7, we show the practical effectiveness of the low-rank truncation discussed in Sec. II D. We evaluate the compression percentage based on the ratio between the rank of the left (*m*_{L}), right (*m*_{R}), and total (*m*_{T}) and the number of basis functions (*M*). Namely,

where *i* ∈ {*L*, *R*, *T*}. The set of parameters that we chose to look at this is $N\xaf\u224314$, *M* = 1189, *r*_{s} = 0.5, 4.0, and Θ = 0.5, 0.125. For the purpose of demonstration, tuning chemical potential is not important, so we set *μ* = *μ*_{T}, which ensures a correct average number of particles in the trial density matrix.

Comparing Figs. 7(a) and 7(b), we see the effect of the temperature change at *r*_{s} = 0.5. As expected, the lower the temperature, we observe the higher compression ratio. In fact, we are in the limit where $M\u226bN\xaf$, so the low-rank compression becomes very effective. As we move along an imaginary path, we see that the compression efficiency decreases for the left (trial density matrix blocks) and increases for the right (sampled blocks). This is the feature of our propagation algorithm, which replaces $B^T$ by $B^$ each time step. We observe nearly constant compression efficiency for the total block. We obtain about 80% compression in *m*_{T} at Θ = 0.5 and nearly 100% compression in *m*_{T} at Θ = 0.125. A similar conclusion can be drawn from Figs. 7(c) and 7(d). We see only small changes in the compression efficiency when changing from *r*_{s} = 0.5 to *r*_{s} = 4.0.

While this technique is useful in speeding up the calculation, in general, one needs a small $N\xaf/M$ ratio and low temperature for a large saving. In our work, among larger calculations, the biggest saving was made in the case of $N\xaf=66$ and *M* = 485 at Θ = 0.125, which will be presented below. For this particular example, we did not find the cost saving to be substantial due to its sizable $N\xaf/M$, but the compression algorithm was still used for the numerical stability.

### C. Comparison to other approaches

In Sec. IV A, we focused on assessing the accuracy of ph-FT-AFQMC on a very small supercell for a given finite basis set. The results from Sec. IV A suggest that ph-FT-AFQMC is quite accurate for *r*_{s} ≤ 3 and it is now important to assess the utility of this method in more realistic calculations specifically toward the CBS limit.

In Fig. 8, we investigate the magnitude of the basis set error as a function of *r*_{s} at Θ = 0.5 where there exist previous exact CPIMC results as well as RPIMC data. We see that the ph-FT-AFQMC results are indistinguishable from the CPIMC results at *r*_{s} = 0.5 when for *M* ≥ 257 on the plotted scale. We also see that *M* = 485 is sufficient to converge results up to *r*_{s} = 2 on this scale. As seen in previous studies, we find the RPIMC results to be too low; however, the magnitude of this bias seems to decrease with increasing *r*_{s}. We can also ascribe the deviation of the FT-CCSD result from the PIMC data to the basis set size. We see that FT-CCSD and ph-FT-AFQMC agree well for *M* = 123 (the largest basis set considered in Ref. 38) up to *r*_{s} = 1 with the FT-CCSD results deviating from the expected smooth trend with *r*_{s} beyond this. Unfortunately, we find obtaining ph-FT-AFQMC data beyond *r*_{s} = 2 to be too expensive beyond *M* = 257 due to the smaller time steps required [equivalently the larger values of *β* required as *r*_{s} increases as in Eq. (76)]. Nevertheless, the results look promising and it is possible that reliable results could be obtained for even larger *r*_{s} values if more computational resources were expended.

With the confidence that the basis set error and phaseless error are under control in this parameter regime, we next look toward the thermodynamic limit. To reach the thermodynamic limit, we use finite size corrections calculated from the finite-temperature random phase approximation.^{57,96} The subject of finite size corrections is treated exhaustively elsewhere,^{62,66,97} and we will not discuss them in any detail. In essence, we add a correction Δ*u*_{xc}(*r*_{s}, Θ, *N*) to our QMC results where

where $uxcRPA$ can be computed numerically through derivatives of the exchange-correlation free energy. Detailed equations and the code necessary for this are available in Refs. 97 and 98. We also verified the accuracy of these corrections across *r*_{s} and Θ ≥ 1 independently in our supplementary material using the existing CPIMC and PB-PIMC data.^{67}

In Fig. 9, we compare our finite size corrected ph-FT-AFQMC results with *M* = 485 to the fit of Ref. 65 denoted here as GDSMFB and the RPIMC data of Ref. 20. For *r*_{s} < 2.0, we find an excellent agreement between our data and the GDSMFB fit, which is reassuring as the GDSMFB fit was generated using QMC data for the interaction energy, not the internal energy, and also relied exclusively on finite-temperature PIMC data above Θ = 0.5, with a correction being applied to bridge the gap to the known zero-temperature result.^{58,99} We also compare to the KSDT fit from Karasiev *et al.*, which was originally fitted to the RPIMC data and subsequently reparameterized^{100} (corrKSDT) to correct the zero-temperature limit and incorporate the more accurate CPIMC and PB-PIMC data from Ref. 73. Similarly to previous work,^{65,100} we find an excellent agreement between our QMC data and GDSMFB and corrKSDT, while we see some slight deviations at Θ = 0.25 from the KSDT fit. Note that we applied the same size corrections to the ph-FT-AFQMC and RPIMC data points and not the original size corrections from Ref. 20, which were subsequently found to be incorrect.^{66} Thus, the RPIMC data do not generally agree with the parameterizations. Finally, we note that we observe a visible deviation of the ph-FT-AFQMC energy from the fits at Θ = 0.5 and *r*_{s} = 3, which can be attributed to phaseless constraint bias and basis set incompleteness error.

The agreement between ph-FT-AFQMC and the GDSMFB fit gives added confidence that this procedure was well founded and the fit of GDSMFB is highly accurate, particularly in the warm dense matter regime (Θ ≤ 0.5 and *r*_{s} ≤ 2). Similar to other studies,^{21,22,67} we find that the RPIMC data are biased and generally too low; however, this bias becomes small past *r*_{s} = 4. We emphasize that in this parameter regime, previously no other methods than RPIMC could run and no reliable verification of the GDSMFB fit was available. This gives us confidence that the ph-FT-AFQMC method will be useful for *ab initio* simulations of warm dense matters in the future.

Finally, to assess the reliability of ph-FT-AFQMC for the calculation of properties other than the total energy, we investigated the accuracy of the static structure factor

## V. CONCLUSIONS

In this work, we have examined the accuracy of phaseless finite-temperature auxiliary-field quantum Monte Carlo (ph-FT-AFQMC) on the uniform electron gas (UEG) model over various Wigner-Seitz radii *r*_{s} and temperatures Θ. The ultimate goal of this work was to access the regime where Θ ≤ 0.5 and *r*_{s} ≤ 2.0. This is the regime where other commonly used QMC methods such as configuration path-integral MC (CPIMC) and permutation blocking PIMC (PB-PIMC) cannot be easily run. Furthermore, another popular flavor of QMC, restricted PIMC (RPIMC), typically exhibits a large bias for *r*_{s} ≤ 2.0, which is currently the only many-body calculations available in this regime.

We summarize the findings and conclusions of our work as follows:

*Difficulties in reaching the accuracy of the zero-temperature ph-AFQMC (ph-ZT-AFQMC) algorithm*: We showed both analytically and numerically that one should not expect the ph-FT-AFQMC energy to match the ph-ZT-AFQMC energy in the low-temperature limit even when the number of particles is tuned properly.*Utility of low-rank truncation*: We demonstrated that the low-rank truncation method discussed in Ref. 82 is effective in the UEG Hamiltonian as well, especially when Θ < 0.5.*Accuracy of ph-FT-AFQMC energies*: We were able to use the ph-FT-AFQMC reliably to investigate Θ ≤ 0.5 and*r*_{s}≤ 2.0. Given the benchmark results on small basis sets that compare favorably to finite-temperature coupled cluster with singles and doubles (FT-CCSD), ph-FT-AFQMC energies are expected to be accurate for a given basis set. Furthermore, we were able to run a large enough basis set so that the basis set incompleteness error is insignificant for the purpose of our study. Our work suggests that the bias in RPIMC is not negligible for*r*_{s}< 2, and therefore, one must be cautious when using RPIMC for dense electron gas simulations as previously suggested by others.^{21,22,67}*Validation of the GDSMFB fit*^{65}*of the exchange-correlation energies at difficult parameter regimes*: In the regime of*r*_{s}≤ 2.0 and Θ ≤ 0.5, no many-body methods have been able to verify the GDSMFB fit whose parameterization relies on the unbiased PIMC data for Θ ≥ 0.5. We found an excellent agreement between the GDSMFB fit and our ph-FT-AFQMC results.*Accuracy of ph-FT-AFQMC static structure factors*: One benefit of working directly in the second quantized space is to be able to compute properties straightforwardly. As an example, we computed the static structure factor of a super cell of $N\xaf=34$ and compared our results to the numerical exact PIMC results. We found a nearly perfect agreement between ph-FT-AFQMC and PIMC.

Given our results, we are cautiously optimistic that ph-FT-AFQMC is a useful tool that is more scalable than other many-body methods based on the second quantization and can provide accurate results for the regimes where other methods either cannot run at all or cannot perform well.

However, the remaining issues in ph-FT-AFQMC should not be ignored. Most notably, the inability to access the same zero-temperature limit as ph-ZT-AFQMC may become a more serious issue in the future since ph-ZT-AFQMC has been shown to be accurate in many benchmark systems. In light of Eq. (33) and Ref. 54, it will be interesting to investigate different types of constraints in the canonical ensemble, which can guarantee the same zero-temperature limit as ph-ZT-AFQMC.

Furthermore, reaching the basis set limit for higher temperature than Θ = 1.0 is excruciating, even though the imaginary time propagation is short. Since there is no more obvious low-rank structure in the propagator, there seems not much that can be done to speed-up. Nevertheless, given the exceptional speed-up shown by the ph-ZT-AFQMC implementation for solids using graphics processing units (GPUs),^{101} we expect that the analogous ph-FT-AFQMC implementation using GPUs will help greatly to ameliorate this situation.

Looking to the future, one obvious extension will be to explore the dynamical structure factor that can be computed from analytically continued imaginary time displaced correlation functions in ph-FT-AFQMC.^{102–104} These results would build upon recent promising advances in the analytic continuation of PIMC data^{72,105,106} in the warm dense regime, where again we could potentially bridge the gap to lower temperatures. Another interesting avenue would be to explore the ability to compute free energy differences much like is done in the interaction picture DMQMC.^{22} While the extension of ph-FT-AFQMC to Fermi-Bose mixtures when the number of bosons is conserved was presented,^{51} it will be interesting to see further development for the cases with a non-conserving boson number.^{107}

## SUPPLEMENTARY MATERIAL

The supplementary material of this work is available, which contains the test of the RPA size correction used in this correction to the existing PIMC data.

## ACKNOWLEDGMENTS

The work of J.L. was in part supported by the CCMS summer internship at the Lawrence Livermore National Lab in 2018. J.L. thanks Martin Head-Gordon and David Reichman for encouragement and support. We thank Hao Shi, Yuan-Yao He, and Shiwei Zhang for useful conversations on the FT-AFQMC algorithm. This work was performed under the auspices of the U.S. Department of Energy (DOE) by LLNL under Contract No. DE-AC52-07NA27344. The work of F.D.M and M.A.M. was supported by the U.S. DOE, 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 (CPSFM). Computing support for this work came from the LLNL Institutional Computing Grand Challenge program.

## DATA AVAILABILITY

The data that support the findings of this study are available within the article. Raw data are published through a Zenodo repository.^{92}

### APPENDIX: SYMMETRY BREAKING AT ZERO TEMPERATURE

We studied the instability of RHF to UHF for the 66-electron supercell at zero temperature using the algorithm presented in Ref. 61. In Fig. 11, we present the spin expectation value (⟨*S*^{2}⟩) of UHF as a function of *r*_{s}. We found that there is no spin polarization occuring for *r*_{s} < 4.0 at zero temperature. As the focus of this study was the 66-electron supercell model for *r*_{s} ≤ 4.0, we did not consider UHF trial density matrices. Since spin polarization is only small at *r*_{s} = 4.0 and no spin polarization occurs for *r*_{s} < 4.0 at zero temperature, it is expected that spin polarization would not occur at *T* > 0, justifying our choice of the free-electron trial density matrix in this study.