Calculating ground state properties of correlated fermionic systems with BCS trial wave functions in Slater determinant path-integral approaches

TL;DR

This paper presents a stable and efficient method to incorporate BCS trial wave functions into path-integral quantum Monte Carlo calculations, enhancing the study of strongly correlated fermionic systems with pairing correlations.

Contribution

It introduces a novel technique for using BCS wave functions in AFQMC, improving efficiency and accuracy, especially in systems with pairing and sign problems.

Findings

Enhanced Monte Carlo sampling efficiency

Reduced imaginary time for projection in AFQMC

Successful benchmark results on the attractive Hubbard model

Abstract

We introduce an efficient and numerically stable technique to make use of a BCS trial wave function in the computation of correlation functions of strongly correlated quantum fermion systems. The technique is applicable to any projection approach involving paths of independent-fermion propagators, for example in mean-field or auxiliary-field quantum Monte Carlo (AFQMC) calculations. Within AFQMC, in the absence of the sign problem, the methodology allows the use of a BCS reference state which can greatly reduce the required imaginary time of projection, and improves Monte Carlo sampling efficiency and statistical accuracy for systems where pairing correlations are important. When the sign problem is present, the approach provides a powerful generalization of the constrained-path AFQMC technique which usually uses Slater determinant trial wave functions. As a demonstration of the capability of the methodology, we present benchmark results for the attractive Hubbard model, both spin-balanced (no sign problem) and with a finite spin polarization (with sign problem).

Figures

Equations

$$|\varphi ⟩={\hat{\varphi }}_{1,\uparrow }^{†}\dots {\hat{\varphi }}_{N_{\uparrow },\uparrow }^{†}{\hat{\varphi }}_{1,\downarrow }^{†}\dots {\hat{\varphi }}_{N_{\downarrow },\downarrow }^{†}|0⟩$$

$${\hat{\varphi }}_{i,\sigma }^{†}=\sum _{\mathbf{r}}{\left({\Phi }^{\sigma }\right)}_{\mathbf{r},i}{\hat{c}}_{\mathbf{r},\sigma }^{†}$$

$$|{\psi }_{BCS}⟩=\left(\prod _{i=1}^{N_{\uparrow }-N_{\downarrow }}{\hat{d}}_i^{†}\right)\frac{1}{N_{\downarrow }!}{\left(\sum _{\mathbf{r},{\mathbf{r}}^{\prime }}{F}_{\mathbf{r},{\mathbf{r}}^{\prime }}{\hat{c}}_{\mathbf{r},\uparrow }^{†}{\hat{c}}_{{\mathbf{r}}^{\prime },\downarrow }^{†}\right)}^{N_{\downarrow }}|0⟩$$

$${\hat{d}}_i^{†}=\sum _{\mathbf{r}}{\left(D\right)}_{\mathbf{r},i}{\hat{c}}_{\mathbf{r},\uparrow }^{†}$$

$$F_{\mathbf{r},{\mathbf{r}}^{\prime }}=\sum _{\alpha }{f}_{\alpha }{\mathcal{U}}_{\mathbf{r},\alpha }{\mathcal{V}}_{{\mathbf{r}}^{\prime },\alpha }$$

$$\mathcal{A}=\left(\left(\Phi^{\uparrow}\right)^{T}D^{\star}\bigg{|}\left(\Phi^{\uparrow}\right)^{T}F^{\star}\Phi^{\downarrow}\right)$$

$$S=⟨{\psi }_{BCS}|\varphi ⟩=(-1{)}^{N_{\downarrow }(N_{\downarrow }-1)/2}\det \left(\mathcal{A}\right)$$

$${\mathcal{G}}_{\mathbf{r}\sigma ,{\mathbf{r}}^{\prime }{\sigma }^{\prime }}={\delta }_{\sigma ,{\sigma }^{\prime }}\frac{⟨{\psi }_{BCS}|{\hat{c}}_{\mathbf{r},\sigma }^{†}{\hat{c}}_{{\mathbf{r}}^{\prime },\sigma }|\varphi ⟩}{⟨{\psi }_{BCS}|\varphi ⟩}$$

$${\mathcal{G}_{{\bf{r}}\uparrow,{\bf{r}}^{\prime}\uparrow}=\left(\left(D^{\star}\big{|}\left(F^{\star}\Phi^{\downarrow}\right)\right)\mathcal{A}^{-1}\left(\Phi^{\uparrow}\right)^{T}\right)_{{\bf{r}},{\bf{r}}^{\prime}}}$$

$${\mathcal{G}}_{\mathbf{r}\downarrow ,{\mathbf{r}}^{\prime }\downarrow }=\sum _{i=1}^{N_{\downarrow }}{\left(F^{†}{\Phi }^{\uparrow }{\left({\mathcal{A}}^{-1}\right)}^T\right)}_{\mathbf{r},N_u+i}{\left({\Phi }^{\downarrow }\right)}_{i,{\mathbf{r}}^{\prime }}^T,$$

$${\mathcal{F}}_{\mathbf{r}\uparrow ,{\mathbf{r}}^{\prime }\downarrow }=-\sum _{i=1}^{N_{\downarrow }}{\left({\Phi }^{\uparrow }{\left({\mathcal{A}}^{-1}\right)}^T\right)}_{\mathbf{r},N_u+i}{\left({\Phi }^{\downarrow }\right)}_{i,{\mathbf{r}}^{\prime }}^T$$

$$\begin{split}&\overline{\mathcal{F}}_{{\bf{r}}\uparrow,{\bf{r}}^{\prime}\downarrow}=F^{\star}_{{\bf{r}},{\bf{r}}^{\prime}}-\left((D^{\star}\,\big{|}\,F^{\star}\Phi_{\downarrow})\mathcal{A}^{-1}\left(\left(\Phi^{\uparrow}\right)^{T}F^{\star}\right)\right)_{{\bf{r}},{\bf{r}}^{\prime}}\,.\end{split}$$

$${\mathcal{C}}_{{\mathbf{r}}_1\downarrow ,{\mathbf{r}}_2\downarrow ,{\mathbf{r}}_3\uparrow ,{\mathbf{r}}_4\uparrow }=\frac{⟨{\psi }_{BCS}|{\hat{c}}_{{\mathbf{r}}_1,\downarrow }^{†}{\hat{c}}_{{\mathbf{r}}_2,\downarrow }{\hat{c}}_{{\mathbf{r}}_3,\uparrow }^{†}{\hat{c}}_{{\mathbf{r}}_4,\uparrow }|\varphi ⟩}{⟨{\psi }_{BCS}|\varphi ⟩}$$

$$\begin{array}{rl}& {\mathcal{C}}_{{\mathbf{r}}_1\downarrow ,{\mathbf{r}}_2\downarrow ,{\mathbf{r}}_3\uparrow ,{\mathbf{r}}_4\uparrow }={\mathcal{G}}_{{\mathbf{r}}_1\downarrow ,{\mathbf{r}}_2\downarrow }{\mathcal{G}}_{{\mathbf{r}}_3\uparrow ,{\mathbf{r}}_4\uparrow }\\ & -{\overline{\mathcal{F}}}_{{\mathbf{r}}_3\uparrow ,{\mathbf{r}}_1\downarrow }{\mathcal{F}}_{{\mathbf{r}}_4\uparrow ,{\mathbf{r}}_2\downarrow },\end{array}$$

$${\mathcal{C}}_{{\mathbf{r}}_1\sigma ,{\mathbf{r}}_2\sigma ,{\mathbf{r}}_3\sigma ,{\mathbf{r}}_4\sigma }=\frac{⟨{\psi }_{BCS}|{\hat{c}}_{{\mathbf{r}}_1,\sigma }^{†}{\hat{c}}_{{\mathbf{r}}_2,\sigma }{\hat{c}}_{{\mathbf{r}}_3,\sigma }^{†}{\hat{c}}_{{\mathbf{r}}_4,\sigma }|\varphi ⟩}{⟨{\psi }_{BCS}|\varphi ⟩}$$

$$\begin{array}{rl}& {\mathcal{C}}_{{\mathbf{r}}_1\sigma ,{\mathbf{r}}_2\sigma ,{\mathbf{r}}_3\sigma ,{\mathbf{r}}_4\sigma }={\delta }_{{\mathbf{r}}_3,{\mathbf{r}}_2}{\mathcal{G}}_{{\mathbf{r}}_1\sigma ,{\mathbf{r}}_4\sigma }\\ & +{\mathcal{G}}_{{\mathbf{r}}_1\sigma ,{\mathbf{r}}_2\sigma }{\mathcal{G}}_{{\mathbf{r}}_3\sigma ,{\mathbf{r}}_4\sigma }-{\mathcal{G}}_{{\mathbf{r}}_1\sigma ,{\mathbf{r}}_4\sigma }{\mathcal{G}}_{{\mathbf{r}}_3\sigma ,{\mathbf{r}}_2\sigma }\end{array}$$

$$\hat{H}=\hat{K}+\hat{V}$$

$$\hat{K}=-t\sum _{⟨\mathbf{r},{\mathbf{r}}^{\prime }⟩\sigma }{\hat{c}}_{\mathbf{r},\sigma }^{†}{\hat{c}}_{{\mathbf{r}}^{\prime },\sigma }+h.c.,\hat{V}=U\sum _{\mathbf{r}}{\hat{n}}_{\mathbf{r},\uparrow }{\hat{n}}_{\mathbf{r},\downarrow }$$

$$\underset{\tau \to +\infty }{\lim }{e}^{-\tau \left(\hat{H}-E_0\right)}|{\varphi }_0⟩\propto |{\Psi }_0⟩.$$

$$e^{-\tau \left(\hat{H}-E_0\right)}\simeq {\left(e^{-\delta \tau \left(\hat{H}-E_0\right)}\right)}^n\simeq {\left(\int d\mathbf{x}p(\mathbf{x})\hat{B}(\mathbf{x})\right)}^n$$

$$e^{-\delta \tau U{\hat{n}}_{\mathbf{r},\uparrow }{\hat{n}}_{\mathbf{r},\downarrow }}=e^{-\delta \tau U({\hat{n}}_{\mathbf{r}}-1)/2}\sum _{x=\pm 1}\frac{1}{2}{e}^{\gamma x({\hat{n}}_{\mathbf{r}}-1)}$$

$$\hat{B}(\mathbf{x})=e^{\tau E_0}{e}^{-\delta \tau \hat{K}/2}\prod _{\mathbf{r}}{\hat{b}}_{\mathbf{r}}(x(\mathbf{r}))e^{-\delta \tau \hat{K}/2}$$

$${\hat{b}}_{\mathbf{r}}(x)=e^{-\delta \tau U({\hat{n}}_{\mathbf{r}}-1)/2}{e}^{\gamma x({\hat{n}}_{\mathbf{r}}-1)}$$

$$|\Psi (n)⟩=\sum _w{\mathcal{W}}_w(n)\frac{|{\varphi }^w(n)⟩}{⟨{\psi }_T|{\varphi }^w(n)⟩}$$

$$e^{-\delta \tau \hat{K}/2}{\mathcal{W}}_w(n)\frac{|{\varphi }^w(n)⟩}{⟨{\psi }_T|{\varphi }^w(n)⟩},$$

$${\mathcal{W}}_w\frac{⟨{\psi }_T|e^{-\delta \tau \hat{K}/2}{\varphi }^w⟩}{⟨{\psi }_T|{\varphi }^w⟩}\frac{e^{-\delta \tau \hat{K}/2}|{\varphi }^w⟩}{⟨{\psi }_T|e^{-\delta \tau \hat{K}/2}{\varphi }^w⟩},$$

$${\mathcal{W}}_w\to {\mathcal{W}}_w\frac{⟨{\psi }_T|e^{-\delta \tau \hat{K}/2}{\varphi }^w⟩}{⟨{\psi }_T|{\varphi }^w⟩}$$

$$|{\varphi }^w⟩\to e^{-\delta \tau \hat{K}/2}|{\varphi }^w⟩.$$

$$\left(\int d\mathbf{x}p(\mathbf{x})\prod _{\mathbf{r}}{\hat{b}}_{\mathbf{r}}(x(\mathbf{r}))\right){\mathcal{W}}_w\frac{|{\varphi }^w⟩}{⟨{\psi }_T|{\varphi }^w⟩},$$

$$\int d\mathbf{x}\pi (\mathbf{x}){\mathcal{W}}_w\frac{{\prod }_{\mathbf{r}}{\hat{b}}_{\mathbf{r}}(x(\mathbf{r}))|{\varphi }^w⟩}{⟨{\psi }_T|{\prod }_{\mathbf{r}}{\hat{b}}_{\mathbf{r}}(x(\mathbf{r})){\varphi }^w⟩},$$

Full text

I Introduction

Strongly correlated many-body systems are a central challenge of modern physics. These systems exhibit a variety of exotic phenomena in a wide array of physical contexts, from high-temperature superconductors to the cores of neutron stars. Despite the fundamental importance of strongly correlated systems, there are at present relatively few quantitative theoretical treatments of these systems at the many-body level. One cutting-edge method that has demonstrated strong capabilities in the study of ground states and excitations of strongly correlated many-body systems is auxiliary-field quantum Monte Carlo (AFQMC) Zhang et al. (1997); Zhang (1999); Zhang and Krakauer (2003). This method formulates the calculation of ground-state properties as an imaginary-time propagation of an entangled ensemble of independent-particle solutions in auxiliary-fields. The framework united key ingredients of determinantal Monte Carlo Blankenbecler et al. (1981); Sugiyama and Koonin (1986); Assaad and real-space diffusion Monte Carlo (DMC) Ceperley and Alder (1986); Ceperley (1995), casting the solution of the many-body ground state problem as an iterative process involving mean-field states, as in standard density-functional theory (DFT) calculations, living in fluctuating external fields. It is efficiently realized computationally as open-ended, importance-sampled random walks in a manifold of Slater determinants.

In certain situations the AFQMC technique is sign-problem-free (as is the standard determinantal Monte Carlo approach), and numerically exact results can be obtained. There are many examples of important applications along these lines Qin et al. (2016); Vitali et al. (2016); Karakuzu et al. (2018); Otsuka et al. (2018); Hohenadler et al. (2012, 2014); Li et al. (2015, 2018); Lee et al. (2009), with the frontier of interesting sign-problem-free situations still being actively expanded. A key class of applications with direct experimental and theoretical importance are in systems involving ultracold Fermi atoms Carlson et al. (2011); Shi et al. (2015); Anderson and Drut (2015); Bulgac et al. (2008); Paiva et al. (2011); Gilbreth and Alhassid (2013); Alhassid, Y. et al. (2014); Paiva et al. (2004) and more recently, in cold atom systems with spin-orbit coupling Shi et al. (2016); Rosenberg et al. (2017). It was also shown that the use of a BCS trial wave function in the unitary Fermi gas can both reduce projection time (to reach the ground state) and dramatically improve the statistical accuracy of the computed ground-state energy Carlson et al. (2011). However, with BCS trial wave functions, the computation of observables that do not commute with the Hamiltonian, and correlation functions, requires additional methodological steps, which we address in the present paper.

When the sign problem is present, a constraint can be applied to remove the exponential increase of the statistical noise. The constraint is implemented by introducing a trial wave function, whose overlap with each random walker is used to define a gauge condition on the sign or phase of that walker. In such a way the random walks of each walker are constrained to observe the same sign Zhang et al. (1997); Zhang (1999) or gauge Zhang and Krakauer (2003). This is an approximation, which becomes exact if the trial wave function becomes exact. In a large variety of problems, the constraint has been shown to be a rather loose one which yields accurate results and whose accuracy is quite insensitive to the details of the trial wave function Motta and Zhang (2018b); Shi and Zhang (2013). Typically, optimized single Slater determinants (SD) from Hartree-Fock theory or DFT are used as trial wave functions. In systems where strong pairing is present, it can be expected that the use of a BSC trial wave function would improve the results. This was demonstrated in the Fermi gas Carlson et al. (2011) for the calculation of the ground-state energy. (Indeed even in molecular systems where the electron-electron interaction is repulsive, the BCS form can give an improved ansatz Casula and Sorella (2003), certainly more general than the SD.) Thus, in constrained-path AFQMC, the implementation of a BCS trial wave function is potentially more valuable, but the same obstacles in terms of computation of observables must be removed as in the sign-problem-free cases.

In this paper we introduce a general approach which removes these obstacles. The approach allows one to explicitly use the BCS wave function as a trial wave function, whether the system is free of or has a sign/phase problem. The method we discuss can in principle be extended to even more general wave functions. Here we focus on lattice models, although, from the methodological point of view, the extension to quantum chemical systems or realistic solids is straightforward. The BCS wave function, as mentioned, provides a very powerful generalization with respect to single Slater determinants. The applicability and accuracy of this approach make it an essential step towards the quantitative description of many strongly-correlated many-body phenomena, including such exotic behaviors as finite-momentum (FFLO) pairing states and high-temperature superconductivity.

The remainder of the paper is organized as follows. In Section I we introduce the formal definitions of Slater determinants and BCS wave functions, and present some algebraic results that are crucial to the implementation of the methodology. In Section II we briefly review the AFQMC algorithm and discuss the force bias technique, which allows us to perform a very efficient sampling of auxiliary-field configurations. In Section III, we introduce our new technique for computing physical properties. In Section IV, we discuss the choice of BCS wave functions and several algorithmic issues. Finally, in Section V, we present benchmark results including comparisons with exact diagonalization, before drawing our conclusions in Section VI.

II Slater determinants and projected BCS wave functions

Slater determinants and BCS pairing wave functions are essential ingredients in the study of strongly-correlated systems. They serve as the building blocks of several theoretical and computational approaches. Since our purpose is to use BCS wave functions as trial wave functions in AFQMC simulations, we first establish our notation and provide the formal definitions of the wave functions that will be used in the remainder of the paper.

We will limit our attention to the sector of the Hilbert space with fixed number of spin up fermions, $N_{\uparrow}$, and spin down fermions, $N_{\downarrow}$. We assume $N_{\uparrow}\geq N_{\downarrow}$ for clarity. The most general Slater determinant under these assumptions has the form:

[TABLE]

where:

[TABLE]

with ${\bf{r}}$ denoting lattice sites or basis index, whose dimension will be denoted by $L$. The manifold of Slater determinants can thus be parametrized by $2$ complex $L\times N_{\sigma}$ matrices $\Phi^{\sigma}$, for $\sigma=\uparrow,\downarrow$, which are made of the orbitals occupied by the independent fermions.

The most general particle-number-projected singlet BCS wave function has the form:

[TABLE]

It is made of a set of unpaired orbitals:

[TABLE]

and a pairing part which describes pairs of fermions in a singlet-state, with the function $F_{{\bf{r}},{\bf{r}}^{\prime}}$ the two-body wave function of the pair. The wave function in (3) is parametrized by an $L\times(N_{\uparrow}-N_{\downarrow})$ complex matrix $D$ and a $L\times L$ complex matrix $F$. Performing a singular value decomposition,

[TABLE]

we see that (1) is actually a special case of (3) and corresponds to the situation when only $N_{\downarrow}$ of the singular values in $f_{\alpha}$ are non-zero. More generally, (3) is a linear combination of Slater determinants. However, it is much more convenient and computationally efficient compared to a generic linear combination, i.e. a multi-determinant, within the AFQMC approach, since the complexity remains comparable to the single determinant case, as we will discuss below. Physically, the pairing wave function is clearly the most natural choice where fermion pairing correlation is expected.

We next outline some algebraic results for suitable matrix elements of operators between a BCS wave function as in (3) and a Slater determinant as in (1). Some of the results have been derived before Carlson et al. (2011); Shi and Zhang (2017) but we include them here to facilitate ensuing discussions. The central object is the overlap matrix:

[TABLE]

where the vertical bar means that $\mathcal{A}$ is obtained by horizontally stacking the $N_{\uparrow}\times(N_{\uparrow}-N_{\downarrow})$ matrix $\left(\Phi^{\uparrow}\right)^{T}D^{\star}$ and the $N_{\uparrow}\times N_{\downarrow}$ matrix $\left(\Phi^{\uparrow}\right)^{T}F^{\star}\Phi^{\downarrow}$, which results in an $N_{\uparrow}\times N_{\uparrow}$ matrix. We observe that the actual computation of $\mathcal{A}$, if we store the full matrix $F$ and do not assume any additional property, for example translational symmetry, requires three matrix multiplications, with complexity $O(L^{2}N)$. We will need to calculate the overlap between the BCS wave function and a Slater determinant:

[TABLE]

where the sign factor is due to the convention we use in (1), writing the Slater determinants with all the spin up first and the spin down later. The Green function matrix elements:

[TABLE]

also have simple expressions:

[TABLE]

and

[TABLE]

where $N_{u}=N_{\uparrow}-N_{\downarrow}$ is the number of unpaired orbitals. We note that the Green functions are $L\times L$ complex matrices and their computation requires inversion of the matrix $\mathcal{A}$ ($\mathcal{O}(N_{\uparrow}^{3})$ operations) and matrix multiplications ($\mathcal{O}(L^{2}N)$ operations).

Another important component of these simulations is the calculation of two-body correlation functions, i.e., four point correlators. These can be constructed from the Green functions above and the anomalous correlators:

[TABLE]

and

[TABLE]

These computations share the same complexity as the Green functions, and have similar definitions. For example, $\mathcal{F}_{{\bf{r}},{\bf{r}}^{\prime}}$ gives the matrix element of the operator $\hat{c}_{{\bf{r}},\uparrow}\hat{c}_{{\bf{r}}^{\prime},\downarrow}$ between the Slater determinant $|\phi\rangle$ and a BCS wave function which does not conserve the number of particles, but which is defined by the same pairing matrix $F$ and the same set of unpaired orbitals (i.e., the parent wave function from which the $\psi_{BCS}$ is derived via number-projection). With some care we can apply Wick’s theorem to obtain the expressions for the two-body correlations. For example,

[TABLE]

is given by

[TABLE]

and

[TABLE]

is given by

[TABLE]

III The AFQMC algorithm with force bias

With the formalism introduced in the previous section, we now describe the technique. We will focus on lattice models and present the algorithm in the constrained-path formulation, commonly referred to as constrained-path auxiliary-field quantum Monte Carlo (CP-AFQMC). However, as mentioned, the algorithm we introduce can be generalized to the phaseless AFQMC approach for molecules and solids. It is also straightforward to apply it to standard determinantal Monte Carlo.

For the purpose of illustration and concreteness, we start with a single-band Hubbard model:

[TABLE]

with:

[TABLE]

where, as usual, $\hat{n}_{{\bf{r}},\sigma}=\hat{c}^{\dagger}_{{\bf{r}},\sigma}\hat{c}_{{\bf{r}},\sigma}$ and the brackets denote nearest-neighbor sites. The approach relies on the imaginary time evolution operator $\exp({-\tau\hat{H}})$ which, for large $\tau$, projects onto the ground state wave function $\left|\,\Psi_{0}\,\right\rangle$ of the model. For any chosen initial wave function $\left|\,\phi_{0}\,\right\rangle$, not orthogonal to the ground state, we have:

[TABLE]

For simplicity, we assume $\left|\,\phi_{0}\,\right\rangle$ to be a single Slater determinant, for example, the non-interacting Fermi sea. (It is straightforward to use a linear combination of Slater determinants. Moreover, it may become important to use a BCS initial wave function, especially in the path-integral formalism typically used for sign-problem-free situations. This will be discussed in the next section.)

The relation (19) defines a path integral in the Hilbert space of the system which connects the initial wave function to the ground state wave function. The essence of CP-AFQMC is to map (19) onto a random walk in the manifold of Slater determinants, which can be efficiently sampled through Monte Carlo methods. The key tool to achieve this is a combination of the Trotter decomposition and the Hubbard-Stratonovich (HS) transformation:

[TABLE]

In (20) the integration is carried over the configurations ${\bf{x}}$ of an auxiliary-field defined on the lattice, each configuration being weighted by a probability density $p({\bf{x}})$. The operator $\hat{B}({\bf{x}})$ is a one-body propagator, describing the evolution of a non-interacting Fermi gas embedded in an external random field. The crucial advantage of using the HS transformation is that the operator $\hat{B}({\bf{x}})$ transforms Slater determinants into Slater determinants: whenever $|\phi\rangle$ is a Slater determinant, $\hat{B}({\bf{x}})\,|\phi\rangle$ is also a Slater determinant. More explicitly, for the attractive Hubbard model, the auxiliary-field is an Ising field on the lattice: ${\bf{x}}({\bf{r}})=\pm 1$, with uniform probability density $p({\bf{x}})=1/2^{L}$, and the propagator $\hat{B}({\bf{x}})$ can be constructed using the operator identity:

[TABLE]

where $\hat{n}_{{\bf{r}}}=\hat{n}_{{\bf{r}},\uparrow}+\hat{n}_{{\bf{r}},\downarrow}$ is the local particle density operator. Note that there are variants of this form which can affect the Monte Carlo sampling efficiency and the systematic accuracy of the constraints (see, e.g., Ref. Shi and Zhang (2013)), but we will not distinguish them here and will focus on the general formalism instead. The one body propagator takes the form:

[TABLE]

with:

[TABLE]

Using these ingredients, Eq. (19) is mapped onto an ensemble of paths in the manifold of Slater determinants and the ensemble average recovers the fully correlated problem. This average corresponds to the multidimensional integral over the space-time dependent auxiliary field configurations in (20). In order to compute this integral, an importance sampling scheme is implemented through the introduction of a trial wave function $\left|\,\psi_{T}\,\right\rangle$, which is assumed to be a good approximation to the ground state wave function. At imaginary time $\tau=n\delta\tau$, a stochastic linear combination of the form:

[TABLE]

is generated, where $|\phi^{w}(n)\rangle$ are Slater determinants, i.e., the walkers, which are labeled by $w$, and $\mathcal{W}_{w}(n)$ are their weights.

At the initial time $n=0$, we let $|\phi^{w}(n=0)\rangle=|\phi_{0}\rangle$, making all the walkers start from the initial wave function, and the weights are $\mathcal{W}_{w}(n=0)=1$, such that $\left|\,\Psi(n)\,\right\rangle=|\phi_{0}\rangle$, apart from an irrelevant normalization factor.

For $n>0$, to build $\left|\,\Psi(n+1)\,\right\rangle$ starting from $\left|\,\Psi(n)\,\right\rangle$, we proceed as follows. For each walker $w$, to apply the kinetic energy (or one-body part of the Hamiltonian) term:

[TABLE]

we have

[TABLE]

where we have omitted the dependence on $n$ to keep the notations simple. This implies that the kinetic contribution to the propagation in imaginary time leads to the simple updates:

[TABLE]

and

[TABLE]

Let us now turn to the interaction part:

[TABLE]

which can be rewritten as:

[TABLE]

where

[TABLE]

The idea of importance sampling is to sample ${\bf{x}}$ not from the bare $p({\bf{x}})$, but according to a suitable approximation for the function $\pi({\bf{x}})$, which favors walkers with larger overlap with the trial wave function. Such importance sampling can improve the efficiency dramatically even with a modest trial wave function, because of the large dimensionality involved (i.e., of ${\bf x}$ and the fact that this is repeatedly applied over many iterations $n$).

In the present paper we implement importance sampling in the above by using a force bias Zhang and Krakauer (2003). We introduce an additional (not normalized) probability density $\tilde{p}({\bf{x}})$, with $\tilde{\mathcal{N}}=\int d{\bf{x}}\,\tilde{p}({\bf{x}})$ and rewrite Eq. (30) as

[TABLE]

and choose:

[TABLE]

The force bias above is a discrete version of the typical form in a shifted Gaussian probability Zhang and Krakauer (2003), and allows us to continue to use the Ising auxiliary-fields from Eq. (21). It is closely related to the form used in Ref. Shi et al. (2015), obtained from $\pi({\bf{x}})$ by expanding the exponential up to order $\sqrt{\delta\tau}$, as an approximation for $\pi({\bf{x}})$ in the small time step limit.

In summary, to implement the interaction part of the propagation, we first compute the mixed estimator of the density:

[TABLE]

sample a configuration ${\bf{x}}$ drawn from $\tilde{p}({\bf{x}})$, which consists of sampling independent random variables on each lattice site, and update:

[TABLE]

and

[TABLE]

Iterating this procedure and extrapolating to infinite total imaginary time, infinite number of walkers, and zero time step, the stochastic linear combination $\left|\,\Psi(n)\,\right\rangle$ converges to the ground state wave function. Note that at the limit of $\delta\tau\rightarrow 0$, the constraint is automatically imposed by the importance sampling in the case of a sign problem. With a twist boundary condition or when a magnetic field is imposed, a “weak” phase problem arises for which a straightforward generalization can be applied Chang and Zhang (2008). When a more intrinsic phase problem (with realistic electron-electron interactions, for example) is present, a phaseless approximation needs to be imposed which requires an additional projection beyond the importance sampling Zhang and Krakauer (2003).

IV Computation of observables

Suppose we choose a BCS wave function, $|\psi_{BCS}\rangle$, as the trial wave function for a CP-AFQMC calculation, or use $|\psi_{BCS}\rangle$ as the initial wave wavefunction $|\phi_{0}\rangle$ in a path-integral formulation as is adopted in more standard sign-problem-free calculations. To compute a general expectation value of a physical observable in either case, the propagation of $|\psi_{BCS}\rangle$ by the one-body propagator in Eq. (22) is necessary. This causes a computational difficulty as we discuss below, followed by a proposed solution.

The simplest kind of estimator that we can build is a mixed estimator, defined as:

[TABLE]

Using the stochastic linear combination (24) as an approximation to the ground state $|\,\Psi_{0}\rangle$ we get immediately:

[TABLE]

This estimator can be readily computed using the sampling scheme described above, combined with the algebraic relations that we introduced in Sec. II. This is the approach used to compute the Bertsch parameter in Ref. Carlson et al. (2011) and the equation of state in the two-dimensional Fermi gas Shi et al. (2015).

For a general observable, or correlation functions, the mixed estimator is insufficient and we need to compute the pure estimator:

[TABLE]

In the special case when the operator $\hat{O}$ commutes with the Hamiltonian operator, and thus with the imaginary time propagator $\exp(-\tau\hat{H})$, the mixed estimator in Eq. (37) becomes equivalent to the pure estimator in Eq. (39). In general, however, the two are different and the mixed estimator is biased.

We introduce here a new technique to compute the pure expectation value of Eq. (39) within the auxiliary-field framework when a BCS wave function is involved. Below we describe the technique in a CP-AFQMC calculation. It is straightforward to generalize it to the path-integral formulation for sign-problem-free calculations, or indeed in several other contexts (for example mean-field calculations by projection), which we will discuss in Sec. VI.2.

We make use of the following:

[TABLE]

for large $\overline{\tau}=m\delta\tau$, which will be accomplished by carrying out many discrete time-steps. We first consider, for the purpose of illustration, a single time-step:

[TABLE]

where we have reinstated the time label $n$ which is useful here and we have used the same manipulations described in the previous section.

In the usual approach to (41), referred to as back propagation Zhang et al. (1997); Motta and Zhang (2018a), a walker $|\phi^{w}(n)\rangle$ is first propagated in the forward direction as described in the previous section, and the path ${\bf{x}}$ that is sampled is then used to explicitly transform:

[TABLE]

building a new bra which is used in combination with $|\phi^{w}(n)\rangle$ to compute the estimator. The update (42) can be implemented for a BCS wave function, resulting in another BCS wave function with a different set of unpaired orbitals and a different pairing matrix. More specifically in the expression in Eq. (5), the unpaired orbitals and $\mathcal{U}_{{\bf{r}},\alpha}$ are multiplied by the spin-$\uparrow$ part of the propagator, while $\mathcal{V}_{{\bf{r}}^{\prime},\alpha}$ is multiplied by the spin-$\downarrow$ part of the propagator, similar to how the $\uparrow$- and $\downarrow$-spin parts of a single Slater determinant are propagated. However, when iterated many times (to reach large $m$), the orbitals will become numerically contaminated. With Slater determinants, this can be controlled by Gram-Schmidt re-orthonormalization, in which all but the resulting orthonormal orbitals can be discarded without affecting the algorithm. With $\langle\psi_{BCS}|$, this is not the case, and we did not find a satisfactory way to keep an arbitrarily high level of numerical stability.

Guided by the desire to have more general and numerically stable calculations in the BCS case, we devise a different approach, which we present in the remainder of this section. Let us consider the case:

[TABLE]

with $|\mu\rangle=\sum_{{\bf{r}}}\langle{\bf{r}}|\mu\rangle|{\bf{r}}\rangle$ and similarly for $|\nu\rangle$. We use the operator identity:

[TABLE]

which rests on the fact that $\hat{B}({\bf{x}})$ is a single particle propagator. Intuitively, if we want $\hat{B}({\bf{x}})$ to “jump” over a $\hat{c}^{\dagger}_{\mu}$ or $\hat{c}_{\nu}$ operator, we need to let the orbitals $|\mu\rangle$ and $|\nu\rangle$ evolve under the action of $\hat{B}({\bf{x}})$, in the forward and backward time direction, respectively:

[TABLE]

We can thus write:

[TABLE]

The right-hand side of (46) closely resembles the usual mixed estimator and has the further advantage that we can apply Gram-Schmidt decomposition to the Slater determinant as it is propagated:

[TABLE]

without changing the result. The difference is that we need to follow the evolution of the orbitals (45), which is straightforward to implement.

The generalization of the above single time-step procedure to more time steps $m>1$ follows immediately by iteration. In particular, if one wishes to compute the full one-body density matrix, it is sufficient to compute the evolution matrix:

[TABLE]

as well as the corresponding one for the backward evolution involving the product of the inverse adjoint of the propagator matrices.

We observe that, although we are able to stabilize the Slater determinant as it is propagated forward, we still have to deal with matrices of the form in Eq. (48), which unavoidably leads to numerical instabilities. To control these instabilities, we assume that:

[TABLE]

is composed of orthonormal orbitals, which can be achieved easily via the forward stabilization already in place in AFQMC. We use the following properties of creation and destruction operators:

[TABLE]

where:

[TABLE]

In (50), the first equality follows from canonical anticommutation relations, while the other relations are closely related to the Pauli exclusion principle. Intuitively, when we destroy an orbital in $|\phi^{\prime}\rangle$, only the component in the linear span of the orthonormal orbitals defining $|\phi^{\prime}\rangle$ contributes; for the creation operator, the “opposite” is true, where we have to project out all linear dependencies Vitali et al. (2016).

We thus obtain for the Green function:

[TABLE]

where, in the last line, we have used the fact that $\langle\tilde{\nu}^{\prime}|\tilde{\mu}^{\prime}\rangle=0$ by construction. In this way, instead of (45), we have the modified evolution:

[TABLE]

and its obvious iteration for $m>1$. The key advantage is that, thanks to (51), the projection on the linear span of the orbitals defining $|\phi^{\prime}\rangle$ and on its orthogonal complement, we can keep the instability under control. Empirically, we always find that we can achieve a robust extrapolation to $m\to+\infty$, while keeping the statistical noise level small enough. This completes the full CP-AFQMC technique which uses a BCS wave function as a trial wave function and enables the computation of the ground-state expectation of any physical property.

V Results

In this section we present two sets of benchmark results obtained using the method we have described. We compare these results with exact diagonalization (ED) calculations in small lattices and standard CP-AFQMC with Slater determinants in larger lattices.

V.1 Spin-balanced system

As a demonstration of the accuracy of the methodology, we have performed several benchmark calculations, which we compare with ED for a system of $N_{\uparrow}=N_{\downarrow}=5$ fermions moving on a two-dimensional lattice with $3\times 4$ sites in periodic boundary conditions and an interaction strength $U/t=-8$. In Fig. 1, we show the comparison for three different correlations functions, which are important for the study of correlated systems. The spin correlation function $\left\langle{\bf{\hat{S}}}({\bf{r}})\cdot{\bf{\hat{S}}}({\bf{r}^{\prime}})\right\rangle$, the density correlation function $\left\langle{\hat{n}}({\bf{r}}){\hat{n}}({\bf{r}^{\prime}})\right\rangle$, and the on-site s-wave pairing correlation function $\left\langle{\hat{\Delta}}^{\dagger}({\bf{r}}){\hat{\Delta}}({\bf{r}^{\prime}})\right\rangle$. The spin operator is defined as

[TABLE]

with $\vec{\sigma}_{\sigma,\sigma^{\prime}}$ denoting the elements of the Pauli matrices. The density operator is ${\hat{n}}({\bf{r}})=\sum_{\sigma}\hat{c}^{\dagger}_{{\bf{r}},\sigma}\hat{c}_{{\bf{r}},\sigma}$, while the pairing operator is ${\hat{\Delta}}({\bf{r}})=\hat{c}_{{\bf{r}},\downarrow}\hat{c}_{{\bf{r}},\uparrow}$.

As seen from the figure, the method is numerically exact, even at large interaction strengths, in sign-problem-free systems.

V.2 Spin-imbalanced system

We have demonstrated that a BCS trial wave function can yield numerically exact results in spin-balanced, sign-problem-free systems. In these systems the use of a BCS trial wave function can improve the statistical efficiency, and reduce the projection length in imaginary-time for reaching the ground state. We now consider the more computationally challenging case of non-zero spin polarization, which leads to the emergence of a sign-problem.

When the sign-problem is present, the systematic accuracy, as well as the efficiency of the simulation, can be affected by the choice of trial wave function. Here we show that, for spin-imbalanced systems with a sign-problem, BCS trial wave functions obtained via the procedure outlined in Sec. VI, which makes use of a Hartree-Fock-Bogoliubov (HFB) transformation, offer a significant improvement over the standard choice of Slater determinants.

We provide an example of this improvement in Figs. 2,  3, and  4, which compare the accuracy of results for different observables calculated using AFQMC with a Slater determinant or BCS trial wave function. Results are shown for the case of $N_{\uparrow}=7$, $N_{\downarrow}=5$ on a two-dimensional lattice with $4\times 4$ sites and periodic boundary conditions, at an interaction strength of $U/t=-4$, in which we can perform exact diagonalization for comparison. For each observable, although AFQMC yields good results with qualitatively correct predictions of the correlation functions using either a single determinant (SD) trial wave function or a projected BCS trial wave function derived from HFB, quantitative differences are evident. The results obtained using the BCS trial wave function show better agreement with the exact result. The insets in each figure, which plot the magnitude of the error relative to the exact result for either choice of trial wave function, show generally smaller errors for all three observables with the choice of BCS trial wave function.

V.3 Large systems

As a final example, we apply the new method to a larger supercell to illustrate that it can scale straightforwardly and that it retains complete numerical stability for large system size and long imaginary-time projection.

Figure 5 plots the on-site pairing correlation function for a $12\times 12$ lattice, hosting $122$ particles with $U/t=-4$. Though this system is spin-balanced and therefore sign-problem-free, it can be used to test the constrained-path approximation, by allowing the force bias to naturally take effect and impose an artificial constraint. As we mentioned at the end of Sec. III, when there is a sign problem (as opposed to a phase problem), the use of the force bias imposes a barrier at $\langle\psi_{T}|\phi^{w}\rangle=0$ such that, in the limit of $\delta\tau\rightarrow 0$, no walker would be able to cross. In an actual application, once the absence of the sign problem is established (either with analytical arguments or empirical testing), we can simply add a small shift in the importance function, similar to constraint release Shi and Zhang (2013), to remove the barrier. But if we do not introduce this step, the CP calculation would proceed, and a constraint bias can be present even if there is no sign problem. Here we take advantage of this feature as a way to both benchmark and illustrate the new method with BCS trial wave function.

We use the path-integral formalism, as is commonly employed in sign-problem-free cases, which samples fixed length paths in auxiliary-field space using the Metropolis algorithm. The Metropolis sampling procedure can be accelerated using the idea of force bias, as introduced in Ref. Shi et al. (2015), and the initial wave functions are chosen as single determinants. Exact results are obtained with this approach to use as a benchmark here. We then carry out two sets of CP calculations using the method discussed above, with either a single Slater determinant or a projected BCS trial wave function.

As Fig. 5 reveals, the CP-AFQMC calculation with Slater determinant trial wave function shows a bias, with the pairing correlation function strength considerably under-estimated. The BCS trial wave function, on the other hand, enables the CP calculation to recover the exact result. The tiny statistical error bars seen in the inset of Fig 5 are a reflection of the numerical stability and good efficiency of the new method as discussed earlier.

These results illustrate the potential of the new method for treating a variety of interesting systems where pairing correlation is important, for example spin-imbalanced atomic gas systems and optical lattices, which are experimentally accessible and which may host exotic phases like FFLO states, as well as various model systems including those relevant to unconventional superconductivity. Almost all of these systems would incur a sign problem. The method we have introduced will allow us to use CP-AFQMC with BCS trial wave functions to study these systems.

VI Discussion

In the preceding sections we outlined a method for incorporating BCS trial wave functions into AFQMC simulations for computing observables and correlation functions, and demonstrated its implementation with a set of benchmark calculations. Here we address some additional technical aspects and comment on possible extensions and applications.

VI.1 The choice of the BCS wave function

We begin by considering the construction of the BCS state to be used as trial wave function. The method we have introduced is independent of this discussion. However, we expect these considerations to be useful in applications for constructing the optimal trial wave functions.

In the simple case of translationally invariant systems where singlet pairing with zero total momentum of the pair is expected, the formalism of standard BCS theory can be applied. The ”textbook” definition of the fully-paired BCS wave function is:

[TABLE]

where $u_{{\bf{k}}}$ and $v_{{\bf{k}}}$ are the coefficients of the celebrated Bogoliubov transformation. When projected onto a sector of the Hilbert space with a fixed number of particles, the wave function (55) can easily be recast in the form (3), with:

[TABLE]

where some care has to be taken if some of the coefficients $u_{{\bf{k}}}$ vanish. If there are unpaired fermions, their orbitals will be simple plane-waves.

More generally, the Hartree-Fock-Bogoliubov (HFB) methodology Ring and Schuck (1980) can be used in order to obtain a generalized Bogoliubov transformation of the form:

[TABLE]

where the matrix of the transformation has $4L\times 4L$ elements. Here we show how to connect the transformation (57) to a wave function of the form (3). We build the unpaired orbitals, i.e. the matrix $D$, and the pairing matrix $F$ in such a way that the wave function (3) is the vacuum of the algebra of operators $\left\{\hat{\gamma}^{\dagger}\,,\,\hat{\gamma}\right\}$, that is:

[TABLE]

We first perform a singular value decomposition of the matrix $V$ in (57):

[TABLE]

where $\mathcal{A}$ and $\mathcal{B}$ are unitary complex matrices, while $\mathcal{Z}$ is real and diagonal, and introduce the $2L\times 2L$ complex matrix:

[TABLE]

We next consider the density matrix and the pairing tensor, defined on the vacuum of the operators $\left\{\hat{\gamma}^{\dagger}\,,\,\hat{\gamma}\right\}$:

[TABLE]

It is easy to see,

[TABLE]

meaning that the matrix $\mathcal{A}$ is made of natural orbitals, which we express as $\{|\alpha\rangle,|\beta\rangle\dots\}$. We have:

[TABLE]

Similarly,

[TABLE]

and we have:

[TABLE]

providing the interpretation of the elements of the matrix $\mathcal{W}$ introduced in Eq. (60).

Now if we introduce the new operators

[TABLE]

and

[TABLE]

we see that the transformation in Eq. (57) is equivalent to the following:

[TABLE]

The unitarity of the matrix on the right-hand side above implies that

[TABLE]

considering that $\mathcal{Z}$ is real and diagonal. In other words, the rows and columns of the matrix $\mathcal{W}$ are orthogonal vectors, their squared norm is $\left(1-\zeta_{\alpha}^{2}\right)$, and they are identically zero whenever $\zeta_{\alpha}^{2}=1$.

We write explicitly:

[TABLE]

We have two possible situations that lead us to classify the states as occupied (”o”) or paired (”p”):

[TABLE]

We define $\tilde{\mathcal{W}}$ to be the matrix obtained from $\mathcal{W}$ keeping only the rows and columns corresponding to the $p$ orbitals, for which $\zeta^{2}_{p}\in[0,1)$. We similarly define $\tilde{\mathcal{Z}}$. With the preceding definitions we construct the antisymmetric matrix:

[TABLE]

and the wave function:

[TABLE]

where $o$ runs only over the occupied states, while $p$ and $p^{\prime}$ run over the paired states. The pre-factor $\mathcal{N}$ is included for normalization. This wave function is the BCS wave function we set out to construct.

It is a simple exercise to verify that

[TABLE]

In order to go back to the original basis of lattice sites, we simply need to use the transformation

[TABLE]

In the case of singlet pairing, $\left|\psi_{BCS}\right\rangle$ can be immediately recast in the form (3), with:

[TABLE]

and:

[TABLE]

The relations (76) and (77) outline as algorithm that provides an interface between the HFB calculation and the QMC simulation. This interface allows us to use a wave function obtained from an HFB transformation as a trial wave function for a CP-AFQMC calculation. We also note that the same HFB transformation can be used to construct an initial wave function that has good overlap with the trial wave function.

VI.2 Possible extensions

While the methodology we have described focuses on the use of BCS states in the constrained-path formalism, the same technique can be applied to wave functions of different forms. This includes multi-determinant trial wave functions, especially when they are composed of non-orthogonal Slater determinants. Multi-determinant wave functions are of considerable importance in quantum chemistry calculations Motta and Zhang (2018b); Landinez Borda et al. (2019), and have been shown to provide high-accuracy results.

This technique can be implemented in the path-integral AFQMC framework, in which fixed-length paths in auxiliary-field space are sampled using the Metropolis algorithm. Path-integral AFQMC has a long record of success in the treatment of various model Hamiltonians, including those with exotic pairing behaviors, such as the attractive Fermi gas Carlson et al. (2011); Bulgac et al. (2008); Shi et al. (2015); Anderson and Drut (2015); Shi et al. (2016); Rosenberg et al. (2017). The use of BCS wave functions would help extend the reach of the method and enable the simulation of even larger systems, with important pairing properties, in both two and three dimensions. One can use a BCS trial wave function at one end of the path, while using either a single-/multi-determinant trial wave function or a BCS trial wave function on the other. (If we use BCS trial wave function on both ends of the path, we could alternatively view the formalism as propagating in HFB space, which has been discussed in Ref. Shi and Zhang (2017).)

Similarly, the same technique we have discussed can be used in a mean-field context, which can be considered a specialized case of the AFQMC, with only a single path instead of the path integral. For example, in a mean-field calculation formulated as an imaginary-time projection Purwanto and Zhang (2004), exactly the same formalism can be adopted. The technique may also be useful for time-dependent mean-field calculation for dynamical properties.

VII Conclusions

We have presented a method for computations of observables and correlation functions using a BCS state as a trial wave function in many-body computations involving path integrals in Slater determinant space. We illustrated the method with a set of benchmark calculations comparing with existing technology and also in situations where a sign problem is present, by comparison with exact diagonalization. We demonstrated that the method removes any numerical instabilities in propagating BCS states. The method controls the sign problem and has computational cost scaling as a low power with system size.

For attractive interactions, the methodology provides a clear improvement over the cutting-edge CP-AFQMC technique relying on Slater determinant trial wave functions. This development will enable high-accuracy computations of exotic superfluid phases, like FFLO states. In the realm of repulsive models and molecular/solid systems, this approach will allow the direct use of a pairing trial wave function in AFQMC calculations, which can be advantageous in, for example, the study of models for interacting electrons where pairing is expected.

The work we have presented here can serve as a general framework for incorporating BCS and other beyond-Slater determinant wave functions into quantum Monte Carlo calculations working in second quantization. The formalism introduced, namely to replace back-propagation in the computation of observables and correlation functions by forward propagating the corresponding Green function matrix, can be useful in other contexts. We hope these developments will enable many applications in a variety of problems, and also stimulate further methodological improvements in the study of strongly correlated models.

We thank J. Carlson, A. Gezerlis, Lianyi He, and Hao Shi for helpful discussions. This work was supported by NSF (Grant No. DMR-1409510) and by the Simons Foundation. Computing was carried out at the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation grant number ACI-1053575, and the computational facilities at William and Mary. The Flatiron Institute is a division of the Simons Foundation.