A Pseudo-BCS Wavefunction from Density Matrix Decomposition:Application in Auxiliary-Field Quantum Monte Carlo

We present a method to construct pseudo-BCS wave functions from the one-body density matrix. The resulting many-body wave function, which can be produced for any fermion systems, including those with purely repulsive interactions, has the form of a number-projected BCS form, or antisymmetrized germinal power (AGP). Such wave functions provide a better ansatz for correlated fermion systems than a single Slater determinant, and often better than a linear combination of Slater determinants (for example from a truncated active space calculation). We describe a procedure to build such a wave function conveniently from a given reduced density matrix of the system, rather than from a mean-field solution (which gives a Slater determinant for repulsive interactions). The pseudo-BCS wave function thus obtained reproduces the density matrix or minimizes the difference between the input and resulting density matrices. One application of the pseudo-BCS wave function is in auxiliary-field quantum Monte Carlo (AFQMC) calculations as the trial wave function to control the sign/phase problem. AFQMC is often among the most accurate general methods for correlated fermion systems. We show that the pseudo-BCS form further reduces the constraint bias and leads to improved accuracy compared to the usual Slater determinant trial wave functions, using the two-dimensional Hubbard model as an example. Furthermore, the pseudo-BCS trial wave function allows a new systematically improvable self-consistent approach, with pseudo-BCS trial wave function iteratively generated by AFQMC via the one-body density matrix.


I. INTRODUCTION
The study of strongly-correlated quantum many-body systems is highly challenging. A general approach does not yet exist to compress the complexity of the manybody wave functions that is widely applicable and yields systematic accuracy across different ranges of many-body models and materials 1-3 . Methodological developments thus have a key role in the study of interacting quantum systems, which spans several subfields in physics, including condensed matter physics, nuclear physics, cold atoms physics, as well as in quantum chemistry and materials science. Recent work and collaborations on method developments have lead to significant progress with computational approaches.
The simplest approaches to many-fermion systems are based on the independent-particle framework. The basic entity in this framework is the Slater determinant. The Slater determinant can be the wave function ansatz itself, as in a Hartree-Fock (HF) mean-field calculation. Alternatively and more commonly, it is used as a vehicle either to capture some property of the system, for example the electronic density and gradients in a densityfunctional theory (DFT) calculations 4,5 , or as the starting point or reference state for many-body calculations such as perturbation theory (e.g. MP2) or coupled cluster [e.g. CCSD(T)] calculations in quantum chemistry 6 , or to impose fermion antisymmetry in quantum Monte Carlo (QMC) calculations [7][8][9] . The Slater determinant can be manipulated with low-polynomial computational cost while fully accounting for permutation antisymmetry, which is a great advantage in treating fermion systems. The Slater determinant has a major shortcoming in these roles, however. It is fundamentally a Fermi liquid picture that contains only occupied orbitals and discards any information on the virtual orbitals.
The BCS wave function 10,11 is a simple ansatz that allows one to overcome this shortcoming of the Slater determinant, and describe a non-trivial modification to the topology of the independent-electron Fermi surface and momentum distribution. A number-projected BCS wave function can be viewed as a linear combination of Slater determinants. It maintains fermion antisymmetry and requires only marginally more computational cost to manipulate. This form can therefore be advantageous in various contexts in the study of interacting fermion systems. For example, the BCS form is clearly more desirable in treating systems with superconductivity. Even in more conventional systems, the BCS wave function (usually in its particle-number projected form) has been found to lead to better trial wave functions [12][13][14][15] in QMC calculations. One can also imagine that a BCS wave function might serve as a better reference state for many-body (e.g., coupled-cluster) calculations.
The BCS wave function can be generated via a meanfield calculation, as the solution to a Hamiltonian which contains attractive interactions. Conceptually, for a system with attractive interactions, the BCS wave function can be thought of as a natural way to go beyond the Fermi liquid picture, by letting pairs of spin-up and spindown particles occupy those orbitals which are virtual in the so-called restricted Hartree-Fock (RHF) solution, thereby creating a better variational wave function with lower energy. (For simplicity, we consider a spin-balanced system with singlet pairing only.) For electronic systems with purely repulsive Coulomb interactions (such as in most electronic models or atoms, molecules, and most real materials), this picture breaks down. If we think of the occupied and virtual orbitals defined by an RHF calculation, the BCS ansatz includes a linear combination of electronic configurations consisting of the RHF and all possible paired excitations from it. The pair occupancy of the higher virtuals raises the oneparticle energy without lowering the interaction energy. The lowest energy mean-field solution is thus a Slater determinant, not BCS.
We consider in this paper a generalized form of the BCS wave function designed for use in electronic systems. The pseudo-BCS wave function, as we shall refer to it by, is based on the concept of the one-body density matrix instead of the variational energy. The wave function is constructed in the structure of a particle-number-projected BCS wave function, with the goal of reproducing the desired one-body density matrix of the many-body system. The pseudo-BCS wave function can be thought of as a general pairing form involving the occupied and virtual orbitals of, for example, an unrestricted Hartree-Fock (UHF) solution (as opposed to RHF in standard BCS). As is well known, the UHF determinant, by allowing spin-symmetry-breaking, often provides a better reference description than RHF for so-called correlated systems (e.g., in Hubbard-like models 16 , or bond breaking in molecules 17 ). In our pseudo-BCS form, we allow the coefficients for orbital pairing occupancy to be complex numbers. The complex phase provides an additional variational degree of freedom to lower the interaction energy, whose physical origin will be discussed below. We introduce a low-cost algorithm to construct the pseudo-BCS wave function by a decomposition of the density matrix.
Our approach allows a straightforward way to incorporate the pseudo-BCS form into a computational framework beyond mean-field. For example, by coupling it with the auxiliary field quantum Monte Carlo (AFQMC) method 8,9 , we obtain a self-consistent procedure in which the pseudo-BCS wave function serves as a trial wave function for the constraint to control the sign/phase problem. The AFQMC calculation computes a one-body density matrix, to which we apply our density matrix decomposition procedure to produce a new pseudo-BCS wave function, and the process is iterated until convergence. We describe the approach below, and illustrate its im-plementation and performance in the Hubbard model. We show that this self-consistent AFQMC can generate a systematically improved trial wave function during the simulation. With the converged pseudo-BCS trial wave function, the AFQMC calculation leads to smaller projection time for reaching the ground state, and smaller systematic bias from the constraint compared with Slater determinant trial wave functions. In the Hubbard model with next-nearest-neighbor hopping, the improvement with self-consistent pseudo-BCS trial wave function allows AFQMC to accurately distinguish the subtle spin orders as the hopping amplitudes are varied.
The rest of this paper is organized as follows. In Sec. II, we introduce the formalism of the pseudo-BCS wave function, and show how it can be obtained from the onebody density matrix of the many-body state. In Sec. III, we briefly introduce AFQMC and present details to couple the AFQMC to pseudo-BCS trial wave functions to achieve self-consistency. In Sec. IV, we show illustrative results, and demonstrate the improvement of the pseudo-BCS wave function and the self-consistent AFQMC with it. Then in Sec. V we discuss several additional points and summarize.

II. PSEUDO-BCS WAVEFUNCTION
The approach discussed in this paper applies generally to any electronic Hamiltonian. To help make the discussion more concrete, we use the Hubbard model to introduce the pseudo-BCS wave function and the density matrix decomposition algorithm. Also the results in Sec. IV will all be for this model.
The Hubbard Hamiltonian is where c † i,σ (c i,σ ) creates (annihilates) an electron with spin σ (σ =↑, ↓) at lattice site i defined in two dimensions on a rectangular lattice of N s ≡ L x × L y sites, and n i,σ ≡ c † i,σ c i,σ is the density operator. Periodic or open boundary condition will be applied. The hopping matrix elements {t i,j } will contain near-neighbor terms with amplitude t (which is used to set the units of energy) and next-near-neighbor terms with amplitude t . The parameters t /t and U/t and the boundary conditions will be specified explicitly later for each system. In this paper, we use N σ to denote the number of electrons with spin σ. As mentioned, we focus in this work on spin-balanced systems, with N ↑ = N ↓ , and no spin-flip terms in the Hamiltonian. We comment briefly on generalizations to other cases in Sec. V.

A. Formalism
In the following, we introduce the formalism and notations, before discussing how to obtain a pseudo-BCS wave function from a given one-body density matrix of the targeted many-body state. A general Slater determinant, for example the UHF solution, can be written as where |0 represents the vacuum state and the operator creates an electron of spin-σ in an orbital described by the n-th column vector (Φ σ ) i,n of the N s × N σ matrix Φ σ , with σ =↑ or ↓. The one-body density matrix of the Slater determinant wave function is given by The structure of a particle-number-projected BCS wave function, also known as anti-symmetrized germinal power (AGP) 18,19 , is where ψ † is a pair creation operator, defined as The N s × N s matrix F can be written, for example by a singular value decomposition (SVD), in the form F = U DV † , with where D is a diagonal matrix whose elements are given by {d n }, and the N s × N s matrices U and V contain single-particle orbitals. Our pseudo-BCS wave function will have the same form, and will allow the diagonal elements to be complex: The pseudo-BCS wave function will couple orbitals described by U and V in a way that represents different principles of pairing from standard BCS, as we discuss in the next subsection.

B. Pseudo-BCS wave function from density matrix
For the reverse process of Eq. (4), namely to find the best Slater determinant wave function given the density matrix G σ , one can construct so-called natural orbitals. The density matrix can be diagonalized to obtain the natural orbitals (eigenvectors) and occupancies (eigenvalues). The matrix Φ σ for the Slater determinant is then be formed by the natural orbitals with the N σ highest occupancies. Clearly the resulting density matrix from Eq. (4) is an approximation and does not reproduce G σ exactly, with information of the higher (virtual) orbitals lost. Now we consider this process for pseudo-BCS wave functions: If we have a best estimate of the one-body density matrix G σ , where |Ψ 0 representing the the many-body ground-state wave function. how to obtain a pseudo-BCS wave function which best reproduces the density matrix. This is accomplished using a decomposition of the density matrix with no truncation of the occupancy, as we illustrate next.
If G ↑ and G ↓ share the same eigenvalues, there exists a pseudo-BCS wave function which can reproduce the density matrix exactly. To demonstrate the procedure, we work with the natural orbitals and occupancy for each spin sector: where Λ is a diagonal matrix with eigenvalues {λ 1 , λ 2 , · · · , λ Ns }, the common occupancy numbers. Using the natural orbitals: we can write down a pseudo-BCS wave function: where {n} = {n 1 , n 2 , · · · , n Nσ } is a set of N σ nonrepetitive indices from {1, 2, · · · , n, · · · , N s }.
In the basis of these natural orbitals, the density matrix of |Ψ is It is straightforward to verify thatG σ is diagonal. To match it with the occupancy of the target density matrix Λ, we have the following condition: where in the numerator the sum is restricted to sets of {n} which contain i. This gives N s equations which can be solved to determine the N s elements of the diagonal matrix D. The resulting pseudo-BCS wave function |Ψ reproduces the target density matrix: The pseudo-BCS wave function |Ψ in Eq. (12) has the characteristic matrix Comparing it to Eq. (7), we see that U = P and V = Q . It is useful to consider the special case of AGP wave functions. If the natural orbitals are real, as in RHF-type orbitals in quantum chemistry, we have P = Q. Otherwise, if the orbitals are complex, P should be equal to Q for BCS pairing. That is, we must organize the degenerate natural orbitals in P and Q such that ↑-spin and ↓-spin orbitals with the same occupancy are arranged in complex conjugate pairs. For example, for the attractive Hubbard model (U/t < 0 in Eq. (1)), the column-vectors P i,n and Q j,n should give a plane-wave with momentum k and −k, respectively, so as to have |k for ↑-spin and |−k for ↓-spin be coupled, and the corresponding F matrix in the form P DP † . In the reverse direction, if we obtain F via a BCS mean-field calculation, which yields diagonal matrix elements in momentum space given by v k /u k 15 , the SVD of F will result in a form U DU † , consistent with the construction from the one-body density matrix discussed above.
In general, when a degeneracy is present in the density matrix eigenvalues in Eq. (10), there are extra degrees of freedom in how the eigenvectors (i.e., natural orbitals) among the degenerate set can be paired in constructing the pseudo-BCS wave function. In the AGP case, this is determined by BCS theory as discussed above. With pseudo-BCS, we have adopted the approach of minimizing the variational energy.
We next discuss the solution of Eq. (14) for D. For wave functions in BCS form without number projection, that is, within grand canonical ensemble, the density matrix decomposition procedure in Eq. (14) has the following exact solution 20 : For our pseudo-BCS wave functions, this is an approximate solution which should become more accurate as the system becomes larger. In our applications below, we use Eq. (17) to determine the amplitude of d i .
Only the amplitude of d i is determined by Eq. (14) or (17), however. Any complex phase factor can be assigned to d i without affecting the density matrixG. We use this degree of freedom to minimize the variational energy of the resulting pseudo-BCS wave function. For the Hubbard Hamiltonian, this turns out to be equivalent to minimizing the double occupancy i n i,↑ n i,↓ , which can be done efficiently. This is further discussed in Sec. V. In general, one can also replace the solution of Eq. (14) with a minimization procedure, to find the set of complex diagonal elements {d n } which minimize the difference between the density matrix and the target, as well as the variational energy, of the resulting pseudo-BCS wave function.
In our applications below, we sometimes deliberately break the spin symmetry in the Hubbard model, by applying a spin-dependent pinning field to induce antiferromagnetic correlations 21 . Similar to the situation in UHF, this preserves the condition that G ↑ and G ↓ share the same eigenvalues but have different eigenvectors (orbitals). For systems that do not satisfy this condition, one could still enforce the condition and adopt the approach of seeking a pseudo-BCS wave function whose density matrices are closest to the target, but that would of course be an additional level of approximation.

III. SELF-CONSISTENT AFQMC WITH PSEUDO-BCS WAVE FUNCTION
As mentioned earlier, the pseudo-BCS wave function can provide a better ansatz not only for systems with superconducting correlations, but for systems with purely repulsive electron-electron interactions. By generalizing the concept of AGP, the pseudo-BCS wave function allows paired states in the sense of UHF orbitals. For many systems with antiferromagnetic correlations, the UHF provides a better description (at the cost of symmetry breaking). In such systems, the pseudo-BCS with UHF orbitals is analogous to AGP with RHF-like orbitals in BCS systems.
In this section, we illustrate one application of the pseudo-BCS form by coupling it to the AFQMC framework. We input into the AFQMC a trial wave function chosen in the form of a pseudo-BCS wave function, and then self-consistently improve it with the density matrix from AFQMC, using the decomposition discussed in the previous section. Below we first briefly review the concept and algorithms of AFQMC, before introducing the procedure for realizing the self-consistency with the pseudo-BCS wave function.

A. AFQMC and self-consistency
The ground-state AFQMC method relies on imaginary time evolution from an initial Slater determinant |Φ I (or any linear combination of Slater determinants): which will project to the ground state of H if the overlap Ψ 0 |Φ I is non-zero. To realize this imaginary time evolution, we first use Suzuki-Trotter decomposition to break up the imaginary-time evolution operator where δτ = τ /N . We then apply the Hubbard-Stratonovich transformation to the two-body term; whereĥ(x) is a one-body operator which dependents on the auxiliary-field vector x, p(x) is a probability density function, and eĥ (x) propagate a Slater determinant |Φ to another Slater determinant |Φ . Putting these together 8 , we have The integrals on the right-hand side are in many dimensions and will require Monte Carlo. In a sense, the imaginary-time evolution in AFQMC can be viewed as an ensemble of random walks in a space of Slater determinants. This population of Slater determinants in the random walks are each orthonormal, but not orthogonal to each other, as the branching random walks occur in an over-complete determinant space 22 .
The anti-symmetry of fermions will cause an arbitrary sign or phase to develop in each of the Slater determinants during the stochastic propagation. If a walker propagates to become perpendicular to the ground state, this walker will effectively cease to contribute under further projection. The number of such walkers will in general grow exponentially with projection time, which will result in large statistical fluctuations. Eliminating them is an exact condition which removes the sign or phase problem 8,9 ; in practice the condition is implemented using a trial wave function |Ψ T , which is applied to constrain the random walk paths. In Hubbard-like models, this constraint is Ψ T |Φ > 0 for ground-state calculations, but for a more general form of interaction a condition involving the phase is needed 9 . If the trial wave function is the exact ground state, then this constraint is unbiased. Many studies have shown that even a freeelectron or HF wave function typically gives high accuracy in a variety of correlated systems 1,2,23 .
To reduce the dependence on the trial wave function in the constrained path or phaseless approximation, a self-consistent procedure was introduced based on Slaterdeterminant-type of wave functions 21 . Here we introduce a new self-consistent AFQMC based on the pseudo-BCS decomposition. A trial wave function in the form of a pseudo-BCS is used as a constraint to carry out the AFQMC calculation, from which a density matrix is computed. We then apply the density matrix decomposition to the result to obtain a new pseudo-BCS wave function, which is used as the trial wave function in the next iteration of AFQMC calculation. The procedure is repeated until the density matrix or other physical observables computed from AFQMC are converged. The approach proposed here achieves self-consistency via the one-body density matrix without needing to involve a fictitious mean-field calculation 21,23 .

B. Additional details: application of pseudo-BCS trial wave function in AFQMC
We provide some of the formalism and details 15,20 necessary to apply a pseudo-BCS trial wave function in AFQMC and to realize the self-consistent procedure described above. We define the overlap matrix between the pseudo-BCS wave function of Eq. (5) and a Slater determinant of Eq. (2) as: The overlap between |Ψ and |Φ is then (The global sign and coefficient above will have no affect in the calculations in this paper.) The one-body term mixed-estimate 22 is and In order to evaluate the pure estimator (as opposed to the mixed estimator) of the density matrix or other observables, back-propagation is needed 8,24 . The backpropagation of a pseudo-BCS wave function is more subtle, and a scheme to ensure numerical stability has recently been proposed 14 , which we adopt here.

IV. RESULTS
In this section, we use the Hubbard model as an example to illustrate the method described above and show the improvement of the self-consistent AFQMC with a pseudo-BCS form of the trial wave function. In Sec. IV A we show results in the pure Hubbard model with t = 0. Here extensive results exist from previous studies which have shown that AFQMC with the usual Slater determinant trial wave functions is quite accurate, and we use this case as a benchmark. Then in Sec. IV B we apply the method to the case with t = 0 where it is shown that our pseudo-BCS leads to an improvement in parameter regimes where the single determinant trial wave function is less accurate.

A. Illustration in the pure Hubbard model
We first study the Hubbard model with only nearestneighbor hopping, i.e., with t i,j = t for near-neighbors i, j and t i,j = 0 otherwise. We will work with a 4 × 4 lattice under periodic boundary condition as a test case, in which exact diagonalization (ED) can be performed straightforwardly. We choose U = 4t, and doping h ≡ 1 − (N ↑ + N ↓ )/N s = 1/8. This doping represents the parameter regime where the sign problem is the most severe, and the constraint error is the largest because of both open-shell and interaction effects. The ground state of this system has three fold degeneracy, distinguished by their overall momentum, (0, 0) and (0, π) or (π, 0) 25 . The (0, 0) case has a degeneracy in the occupancy of the natural orbitals. We focus on the latter situation below in order to compare with single-determinant (SD) cases with no ambiguity. In Fig. 1, we illustrate the density matrix decomposition procedure discussed in Sec. II B. The pseudo-BCS wave function is generated from the exact one-body density matrix from ED. The resulting eigenvalues of the pseudo-BCS wave function from the decomposition is compared with the exact results. The result of a single Slater determinant wave function constructed from natural orbitals (which is equivalent to the free-electron wave function in this case) is also shown. For this test, we obtain the exact ground state density matrix from ED, and apply the density matrix decomposition to obtain a pseudo-BCS wave function. This wave function is then applied in AFQMC as the trial wave function. In Fig. 2, the density matrix computed from AFQMC using this trial wave function is compared with that from ED. For reference, the density matrix com- puted from AFQMC using a SD trial wave function which is formed with the natural orbitals with the N σ largest engenvalues (occupancy) is also shown. We see that the results from AFQMC/pseudo-BCS are not exact, despite using the exact density matrix to generate the pseudo-BCS. This is not surprising, because the pseudo-BCS wave function is not the exact many-body wave function, and apparently still incurs a finite constraint error in AFQMC.

B. Hubbard model at half-filling with t
In this section we study the Hubbard model at halffilling with both near-neighbor (t) and next-nearestneighbor hopping (t ). These systems provide a good test case for us, with the presence of different magnetic orders [26][27][28][29] . We will focus on its magnetic correlations, applying pinning fields 30-32 under cylindrical boundary conditions, i.e., periodic along x-direction and open along y. The pinning fields are applied at the edges of the cylinder, adding a one-body external potential term i,σ u i,σ n i,σ to the Hamiltonian of Eq. (1), with u i,↑ = −u i,↓ = (−1) ix u 0 for i y = 1 and L y , and u i,σ = 0 for all other sites. The strength of the pinning field is fixed at u 0 = 0.25t in our calculations. The cylindrical boundary condition and pinning fields break translational symmetry along the y-direction and induces AFM correlations. Under this setting, two-body spin correlation functions in periodic systems can be probed by onebody spin densities: S z i = (n i,↑ − n i,↓ )/2 . It provides a convenient way to detect the presence and nature of long-range AFM orders including collective modes such as stripes 21,33 . We will study width-4 cylinders which can be treated very accurately by density matrix renormalization group (DMRG) calculations 34 , which we perform TABLE I. Comparison of the computed total energies from self-consistent AFQMC using trial wave functions of pseudo-BCS form (SC AFQMC/pBCS) and single Slater determinant (SC AFQMC/SD), as well as one-shot AFQMC using free-electron trial wave functions (AFQMC/FE) for various systems, compared to DMRG in 4 × 8 and 4 × 16 cylinders with pinning fields applied at the edges. All systems have U = 4t, and the details of the pinning fields in the cylindrical systems are given in the text.
Lx  Convergence of the self-consistent AFQMC calculation using a trial wave function of the pseudo-BCS form versus a single Slater determinant from natural orbitals. The system is a 4 × 8 Hubbard cylinder at half-filling, with pinning field applied on both edges, with U = 4t and t = 0.3t. The top panel shows the ground-state energy and the bottom panel shows the mean squared error δ 2 S z in the spin density with respect to DMRG. The SC procedures were repeated multiple times with different random number seeds to estimate the uncertainties in the convergence process, with the standard deviations shown by the shading on the curve.
We carry out self-consistent AFQMC calculations and benchmark the computed ground-state energy and spin densities. For the SC AFQMC/SD calculations, we keep the trial wave function in the form of a single Slater determinant, which is obtained in the self-consistent it-eration from the natural orbitals of the computed density matrix, taking the N σ leading natural orbitals with the largest occupancies 21 . For the SC AFQMC/pBCS, we compute the one-body density matrix with backpropagation 14 and then apply our density matrix decomposition procedure discussed in Sec. II B. In Fig. 3 we show the convergence of the ground-state energy and spin density as a function of self-consistency iterations. For the spin density, we measure the mean squared deviation δ 2 S z ≡ Ns i (S z i − S z i,DMRG ) 2 /N s with respect to the DMRG reference result. Both self-consistency processes were initialized using the free-electron wave function, which provides the wrong initial input in most cases as it is not magnetically ordered. Both sets of SC calculations yield improved results over the initial AFQMC/FE result, as expected. The SC AFQMC/pBCS shows both a faster convergence and better converged results over the SC AFQMC/SD.
The ground-state energies are listed in Table I for the system above and several other cylindrical systems, using the two different self-consistent approaches, as well as one-shot AFQMC calculations with the freeelectron trial wave function. The final energy for selfconsistent AFQMC with pseudo-BCS is calculated using the pseudo-BCS trial wave function with lowest variational energy after density matrix convergence. The free electron trial wave function is degenerate in 4 × 16 with t /t = 0.3. We break the degeneracy by solving the non-interacting Hamiltonian with a small twisted angle (0.01, 0.01) applied in the boundary condition. Consistent with the trend observed in Fig. 3, both SC procedures are seen to improve the energy, with the SC AFQMC/pBCS giving systematic errors of ≈ 0.2% or less compared to DMRG.
In Fig. 4, we study the spin density more systematically. As a function of the next near-neighbor hopping amplitude t , different magnetic correlations arise, which provides an excellent test ground for the self-consistency procedure via pseudo-BCS decomposition. The converged spin densities from self-consistent AFQMC with pseudo-BCS decomposition are shown and compared to the exact results from DMRG. Results from single determinant trial wave functions are also presented, including both the one-shot calculation using the free-electron trial wave function and the self-consistent constraint from truncated natural orbitals. We see that the SC AFQMC The variational energy of a pseudo-BCS wave function can be evaluated by a Monte Carlo sampling of the ket, which we have used in earlier studies 14,15 and which we briefly describe below.
To evaluate the expectation value of an observable with respect to a pseudo-BCS (or AGP) wave function, we can use the expansion in Eq. (12) to write the ket |Ψ as a linear combination of Slater determinant, |Ψ = Φ c Φ |Φ , where c Φ = d n1 ...d n Nσ and |Φ is the Slater determinant obtained by occupying the set of natural orbital pairs specified by the indices {n 1 , n 2 , · · · , n Nσ }. Equation (26) can then be written as which is similar to the "mixed estimator" in AFQMC. The sum over Φ contains a combinatorial number of terms, which makes an explicit summation impractical for larger systems. We can sample the Slater determinants |Φ according to |c 2 Φ | by, for example, a Markov chain Monte Carlo procedure in the discrete space of indices {n}, proposing to swap one of the occupied index n i with one from the unoccupied set, n i , and accepting the move based on |d n i /d ni | 2 . The mixed estimate between Ψ| and |Φ in the numerator can be computed using Green functions as shown in Eqs. (24) and (25), and applying Wick's theorem for the two-body term in V as needed 14 . (The one-body and two-body terms of a pseudo-BCS wave function could also be computed using related characteristic polynomial 36 .) The phases of {d i } only enter in the local energy in Eq. (27). They do not affect the Monte Carlo sampling or the sampled |Φ 's. The phases in |Ψ can be optimized using, for example, standard variational minimization techniques employed in variational Monte Carlo 37 . In the Hubbard model, the kinetic term K is invariant with respect to the phase; the V term easily decouples into single particle forms for the two spin sectors.

B. Generalization of the pseudo-BCS form
We briefly comment on generalization of the pseudo-BCS approach beyond systems we have discussed, which has been restricted to Hamiltonians with no spin-flip terms and singlet pairing with spin balance.
With a one-body term K that contains spin-flip terms, such as systems with spin-orbit coupling (SOC), the Slater determinant state |Φ contains spin-orbitals and is described by a 2N s × (N ↑ + N ↓ ) matrix 38 : Φ ↑ Φ ↓ . The corresponding one-body density matrix is a 2N s × 2N s matrix. The eigen decomposition of the density matrix gives 2N s natural spin-orbitals U . The pseudo-BCS wave function is given by a matrix F , in the form of: F = U DU T , where D is a skew-symmetric matrix. Similar to the spindecoupled case, the elements of D, {d i,j }, are to be determined by the eigenvalues {λ} of the density matrix, where a pair of comparable eigenvalues λ i = λ j gives d i,j = −d j,i .
The equations for overlap and the mixed estimator between a general Slater determinant Φ with spin-orbitals and a pseudo-BCS wave function Ψ from the above are given by: and where i (and j) is a general spin-orbital index (i.e., combining both the site index i and spin σ in the Hubbard Hamiltonian), and pf denotes Pfaffian 39 .

VI. SUMMARY
In summary, we introduced an approach to use the onebody density matrix to build pseudo-BCS wave functions for a many-body systems, i.e., wave functions with a similar form to number-projected BCS. They can be thought of as a generalization of BCS to systems with repulsive interactions such as in molecules and real materials. Such wave functions can provide a better mean-field ansatz than any single Slater determinants, and a natural and very powerful extension beyond UHF. The pseudo-BCS wave functions are more versatile and has more variational freedom than the usual number-projected BCS wave functions. They can be manipulated with costs similar to AGP or Slater determinants, and can be used as the reference state in a variety of many-body methods. By coupling pseudo-BCS wave functions to AFQMC self-consistently via a density matrix decomposition approach, we achieve a self-consistent AFQMC method, which yields improved results in the Hubbard model over state-of-the-art self-consistent AFQMC calculations based on single Slater determinants.

VII. ACKNOWLEDGEMENT
We thank Ettore Vitali, Mingpu Qin, Yuan-Yao He for valuable discussion. This work is supported by the Simons Foundation and the many-electron collaboration. Computing was carried out at the computational facilities at William & Mary and the Flatiron Institute. The Flatiron Institute is a division of the Simons Foundation.