A phaseless auxiliary-field quantum Monte Carlo perspective on the uniform electron gas at finite temperatures: Issues, observations, and benchmark study

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,¹ the competition between stripe and superconducting orders in the 2D Hubbard model,² and plasmonic catalysis.³ Furthermore, there is a growing interest in warm dense matter,⁴ an extreme state of matter found in planetary interiors⁵ 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.⁸ 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.¹⁷ 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,¹⁸ which leads to the restricted PIMC (RPIMC) approach. This approach is similar in spirit to the fixed-node approximation in diffusion Monte Carlo (DMC),¹⁹ 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²³ including the development of new algorithms such as the, in principle, exact real space permutation blocking PIMC²⁴ (PB-PIMC) and second quantized configuration PIMC²⁵ (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.²⁶ 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.³⁰ 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.³⁹ 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,⁴⁹ 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⁶¹ 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.⁶² In fact, many of the recent developments⁶² 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.⁷³ 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.

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

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 (⁠$Ĥ1$⁠) and two-body (⁠$Ĥ2$⁠) 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,⁷⁴ 

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₁, x₂, …, x_n) is the probability of sampling a specific path designated by auxiliary fields, x₁, x₂, …, 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(Â(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.⁷⁷ 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.

Since the value of determinants in Eq. (8) or equivalently $tr(Â(X))$ for a given X can be either positive or negative (in our case complex),⁷⁸ 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.⁴² This way, we can keep the overall scaling of the algorithm polynomial.

We introduce a trial density matrix, $B^T(Δτ,μT)$⁠, and this is used to impose a phaseless constraint in the imaginary time evolution,⁴³ which is defined as

where $ĤT$ 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Δτ,μ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⁸⁰) 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.

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Δτ,⁴² 

where N is the number of particles and |ψ_T(N)⟩ is the trial wavefunction defined as

with |ψ₀(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(Â(τ,{xi}i=1k−1))=ψT(N)|ψ(τ,{xi}i=1k−1,N)$ where $Â$ 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(Δτ,μT))n−k+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.³⁵ 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 μ.⁵⁴ 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).

It is well known that the standard determinant QMC (DQMC) algorithm suffers from numerical instabilities resulting from the repeated multiplication of the B matrices.⁸¹ This issue can be overcome using the stratification method,⁷⁹ 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 BL/Δτ 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:

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

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:

We, therefore, achieve

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

He and co-workers⁸² 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):

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−1+tq)−1$ 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.

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²¹ 

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 $ρ^$ 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 $Â(Q)$ and $B^(Q)$⁠.

The local energy ϵ_n(τ) for the UEG then reads

where the one-body energy, E₁, is

and the two-body energy, E₂, 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 $Ĥ$⁠. For instance, the one-body and two-body energies can be straightforwardly read off from the local energy evaluation.

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

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.⁹¹ 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,⁶¹ 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 $⟨N^⟩T=N¯$⁠. The interacting chemical potential was then tuned to also match $N¯$⁠. 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 $μ̃$⁠. Following this, we perform five poorly converged simulations (320 walkers averaged over 10 simulations) in the interval of $[0.975μ̃,1.025μ̃]$⁠. 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πN¯1/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_BT)] in atomic unit to increase as $O(rs2)$⁠, which makes the higher r_s values computationally costly. A low-rank threshold of 10⁻⁶ in Eq. (54) was used in every calculation.

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.

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.⁹³ 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¯=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,⁶¹ 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,⁵⁰ we mention a total of four possible reasons for this disagreement:

1. 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.

2. 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.

3. 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;⁵⁰ 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;⁹⁴ however, we see from Fig. 4 that this choice makes essentially no practical difference.

4. Time-reversal symmetry breaking: The final possibility is that as the phaseless constraint breaks imaginary time symmetry in the estimators⁵⁰ 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,⁵⁰ 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).

We now study a larger UEG supercell considering $N¯=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³⁷ 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¯=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¯=2$ and we observe consistently accurate results at a larger supercell (⁠$N¯=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,⁹⁵ 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.⁹⁵ 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,⁶¹ this result is not surprising.

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¯≃14$⁠, 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.