Many-body computations by stochastic sampling in Hartree-Fock-Bogoliubov space

We describe the computational ingredients for an approach to treat interacting fermion systems with pairing fields, based on path integrals in the space of Hartree-Fock-Bogoliubov (HFB) wave functions. The path integrals can be evaluated by Monte Carlo, via random walks of HFB wave functions whose orbitals evolve stochastically. The approach combines the advantage of HFB theory in paired fermion systems and many-body quantum Monte Carlo techniques. The properties of HFB states, written in the form of either product states or Thouless states, are discussed. The states preserve forms when propagated by generalized one-body operators. They can be stabilized for numerical iteration. Overlaps and one-body Green's functions between two such states can be computed. A constrained-path or phaseless approximation can be applied to the random walks of the HFB states if a sign problem or phase problem is present. The method is illustrated with an exact numerical projection in the Kitaev model, and in the Hubbard model with attractive interaction under an external pairing field.

Figure 1

Energy vs imaginary time during projection in the Kitaev model. The lattice size $L_1$ is 100, and the model parameters are $t=1.0,\mathrm{\Delta }=2.0$, and $\mu =-3.2$. A time step $\mathrm{\Delta }\tau =0.01$ was used. Results from propagating product states are numerically the same as those from propagating Thouless states. The mixed estimator and the growth estimator are consistent with each other, and they converge to the exact answer for sufficiently large $\tau$. The inset shows results for $\tau$ from 2 to 10, with log-scale of the energy.

Figure 2

Pairing order $\langle c_r^{†}{c}_{r+1}^{†}\rangle$ vs lattice position $r$ in the Kitaev model computed from $\left|\psi \right(\tau )\rangle$ at different projection times $\tau$, with the same parameters as in Fig. 1. The order parameter converges to the exact solution at the large imaginary-time limit. For clarity, data in the middle of the lattice are shown at every third value of $r$ for $\tau =10$.

Figure 3

QMC calculations by projecting BCS random walkers. Average particle number (for $\uparrow$ electrons) is shown vs chemical potential. The lattice size is $4\times 4$, and model parameters are $t=1.0,U=-12.0$. An imaginary-time step of $\mathrm{\Delta }\tau =0.01$ was chosen, with projection time $\beta =64t$. Our BCS initial wave function has $\langle M_{\uparrow }\rangle =2.0$. The algorithm converges to different densities as $\mu$ is varied and gives accurate results. The plateaus indicate integer particle numbers.

Figure 4

Pairing correlation functions vs distance $r$ (between site 0 and i) computed from QMC and ED. The chemical potential is tuned in QMC to match particle numbers in the ED calculations. The same run parameters are used as in Fig. 3 and Table 1. QMC statistical error bars are smaller than symbol size.

Figure 5

Pairing order vs distance from site (0,0). The lattice sizes are $8\times 8$ and $16\times 16$, with the total number of particles tuned to 10 and 40, respectively. The model parameters are $t=1.0$ and $U=-8.0$. We choose a time step $\mathrm{\Delta }\tau =0.01$ and a projection time $\beta =64$. The pinning field is put on two neighboring sites (0,0) and (1,0), with $h_i=1$.