# Constrained-Path Auxiliary-Field Quantum Monte Carlo for Coupled Electrons
and Phonons

Joonho Lee<EMAIL_ADDRESS>Department of Chemistry, Columbia University,
New York, New York 10027, USA  Shiwei Zhang<EMAIL_ADDRESS>Center for Computational Quantum Physics, Flatiron Institute, New York, NY
10010, USA  Department of Physics, College of William and Mary, Williamsburg,
Virginia 23187, USA  David R. Reichman<EMAIL_ADDRESS>Department of
Chemistry, Columbia University, New York, New York 10027, USA

###### Abstract

We present an extension of constrained-path auxiliary-field quantum Monte
Carlo (CP-AFQMC) for the treatment of correlated electronic systems coupled to
phonons. The algorithm follows the standard CP-AFQMC approach for description
of the electronic degrees of freedom while phonons are described in first
quantization and propagated via a diffusion Monte Carlo approach. Our method
is tested on the one- and two-dimensional Holstein and Hubbard-Holstein
models. With a simple semiclassical trial wavefunction, our approach is
remarkably accurate for $\omega/(2\text{d}t\lambda)<1$ for all parameters in
the Holstein model considered in this study. In addition, we empirically show
that the autocorrelation time scales as $1/\omega$ for $\omega/t\lesssim 1$,
which is an improvement over the $1/\omega^{2}$ scaling of the conventional
determinant quantum Monte Carlo algorithm. In the Hubbard-Holstein model, the
accuracy of our algorithm is found to be consistent with that of standard CP-
AFQMC for the Hubbard model when the Hubbard $U$ term dominates the physics of
the model, and is nearly exact when the ground state is dominated by the
electron-phonon coupling scale $\lambda$. The approach developed in this work
should be valuable for understanding the complex physics arising from the
interplay between electrons and phonons in both model lattice problems and ab-
initio systems.

##  I Introduction

The coupling of electrons to nuclear lattice distortions is responsible for
myriad important physical phenomena in bulk materials.Giustino (2017) In
particular, the thermodynamic and transport properties of solids are crucially
influenced by electron-phonon (el-ph) interactions. Perhaps the most
spectacular consequence of el-ph interactions is the emergence of
superconductivity as described by the Bardeen-Cooper-Schrieffer (BCS) theory.
Here, the el-ph interaction mediates an effective electron-electron (el-el)
attraction which results in the Cooper pairing of electrons of opposite
spin.Mahan (2013) The BCS theory provides a quantitative framework for the
description of conventional superconductivity such as that found at low
temperatures in simple metals.

A simple microscopic picture is unfortunately not available for unconventional
superconductors such as the cuprates, whose critical temperature ($T_{c}$) can
be above 90K at ambient pressure.Cava  _et al._ (1987); Wu  _et al._ (1987) It
is believed that the el-ph interaction alone cannot give rise to these high
$T_{c}$ values. Scalapino (2012) However, experimental evidence exists which
indicates that non-negligible el-ph interactions are present in these
materials.Müller (1990); Song and Annett (1995); Alexandrov (2000); Cuk  _et
al._ (2005); Khatami  _et al._ (2008); Gunnarsson and Rösch (2008); Van Heumen
_et al._ (2009); Gadermaier  _et al._ (2010); He  _et al._ (2018); Yang  _et
al._ (2019); Hu  _et al._ (2019); Rosenstein and Shapiro (2019); Grissonnanche
_et al._ (2020); Rosenstein and Shapiro (2020); Banerjee  _et al._ (2020);
Shneyder  _et al._ (2020); Mishchenko  _et al._ (2020); Sreedhar  _et al._
(2020); Peng  _et al._ (2020) It remains unclear what role el-ph interactions
play in the cuprates and related materials, and if a potentially delicate
interplay between el-el and el-ph interactions may influence their
superconducting properties.

The canonical model Hamiltonian used to capture the physics of the cuprates is
the two-dimensional (2D) repulsive Hubbard model.Hubbard (1963) The ground
state of the hole-doped 2D Hubbard model has been thought to support d-wave
superconductivity for many years.Anderson (1987) A recent joint numerical
study using two state-of-the-art approaches, density matrix renormalization
group (DMRG) and constrained path (CP) auxiliary-field quantum Monte Carlo
(AFQMC), indicates that the ground state of the standard 2D repulsive Hubbard
model with near-neighbor hopping supports modulated phases (e.g. stripes) that
are not superconducting over a range of repulsion strengths and doping
levelsQin  _et al._ (2020) expected to describe the cuprates. This suggests
that features beyond those included in the simple Hubbard model, such as the
effects of multiple bands, longer ranged Coulomb interactions, and/or the role
of el-ph interactions, may be needed to tip the balance of the ground state
towards superconductivity for realistic values of doping levels and the
magnitude of el-el repulsions.

Our work is motivated by precisely these considerations, namely the
development of a scalable and accurate numerical approach that can treat el-ph
effects on the same footing as el-el correlations. This is a challenging task,
as treating the complex electronic degrees of freedom in the pure 2D Hubbard
model is already difficult, even with state-of-the-art numerical
approaches.LeBlanc  _et al._ (2015); Zheng  _et al._ (2017); Qin  _et al._
(2020) The addition of el-ph effects, as contained in, e.g., the 2D Hubbard-
Holstein model, thus requires non-trivial extensions of these approaches in
order to treat electrons and phonons on an equal footing.

Several methods have been formulated or extended to coupled el-ph problems,
including DMRG,Jeckelmann and White (1998); Jeckelmann  _et al._ (1999);
Tezuka  _et al._ (2005, 2007); Fehske  _et al._ (2008); Ejima and Fehske
(2010) variational exact diagonalization,Bonča  _et al._ (1999) variational
Monte Carlo,Ohgoe and Imada (2014, 2017); Karakuzu  _et al._ (2017) dynamical
mean-field theory,Jeon  _et al._ (2004); Paci  _et al._ (2006); Werner and
Millis (2007); Murakami  _et al._ (2013); Li  _et al._ (2017) density matrix
embedding theory, Sandhoefer and Chan (2016); Reinhard  _et al._ (2019) and
coupled-cluster theory.Sibaev  _et al._ (2020); Dresselhaus  _et al._ (2020);
White  _et al._ (2020) There are difficulties facing each approach. For
example, large el-ph couplings and/or small phonon frequencies are challenging
to handle in most methods based on a second quantized representation of
phonons, because of the necessity of truncating the phonon Hilbert space. When
a large number of phonons per site is required, the computational cost
associated with treating them can grow prohibitively expensive. In addition to
the demand of treating the phononic Hilbert space, there is of course the
interacting many-electron problem. Clearly, the treatment of correlated el-ph
coupled systems in two and higher dimensions over a wide range of the
parameter space in an exact or near-exact manner is a forefront challenge.

The method that we propose here is an extension of the CP-AFQMC method
developed and popularized by Zhang and co-workers.Zhang  _et al._ (1995, 1997)
For purely electronic problems, the CP-AFQMC approach is similar to the
determinant quantum Monte Carlo (DQMC) methodBlankenbecler  _et al._ (1981);
Scalapino and Sugar (1981a, b); Johnston  _et al._ (2013); Mendl  _et al._
(2017); Karakuzu  _et al._ (2018); Costa  _et al._ (2020) in the sense that
the two-body propagation is aided by the Hubbard-Stratonovich transformation
Hirsch (1983) and is formulated in the space of determinants. There, however,
are several key differences. CP-AFQMC reformulates the imaginary-time
propagation by working with open-ended random walks. An exact boundary
condition is introduced in auxiliary-field space, which can be approximately
imposed using a trial wave function, to avoid the notorious fermion sign
problem. The open-ended random walk approach allows easy access to zero
temperature results, and is often much less prone to ergodicity problems in
the Monte Carlo sampling. Moreover, CP-AFQMC can be naturally extended to ab-
initio Hamiltonians while coping with the fermionic phase problem associated
with these more complex models using the phaseless approximation instead of
the constrained path approximation. Zhang and Krakauer (2003); Al-Saidi  _et
al._ (2006) Because of the constraint imposed on walker trajectories, CP-AFQMC
is no longer exact, unlike DQMC. Furthermore, due to the constraint, the
ground-state energy computed via the usual mixed estimator is not
variational.Carlson  _et al._ (1999) On the other hand, CP-AFQMC can be used
to access a wider range of interaction strengths and doping regimes in which
DQMC cannot be used due to the inherent sign problem. It should be noted that,
in addition to its flexibility, CP-AFQMC has been shown to yield excellent
accuracy for strongly correlated electrons.Zheng  _et al._ (2017); Qin  _et
al._ (2020)

In this work we devise an extension of CP-AFQMC to treat both electrons and
phonons on an equal footing, while retaining its benefits for electrons. Our
framework is similar to the extension of Green’s function Monte Carlo (GFMC)
as formulated by McKenzie and others,McKenzie  _et al._ (1996) where the
phonons are treated in a first quantized space. We present the formulation of
this new CP-AFQMC approach, provide thorough benchmark results on the 1D and
2D Holstein and Hubbard-Holstein models for various phonon frequencies and el-
ph couplings, and discuss the current scope and limitations of the proposed
approach.

The paper is organized as follows: In Section II we outline the model we study
and the important parameters that control its physics. In Section III we
outline our algorithm. Section IV is devoted to a discussion of trial wave
functions. Section V.1 and Section V.2 discuss distinct perturbative
approaches to the problem outlined in Section II. Section VI and Section VII
present results for the Holstein and Hubbard-Holstein models, respectively.
Section VIII discusses the extension of our approach to realistic ab-initio
problems. In Section IX we conclude.

##  II Model

###  II.1 The Hubbard-Holstein Hamiltonian

Although the approach we outline is general, we focus on a paradigmatic model
of a correlated system coupled to phonons, namely the Hubbard-Holstein
model.Hubbard (1963); Holstein (1959a) The Hubbard-Holstein model is defined
by the following Hamiltonian:

$\displaystyle\hat{\mathcal{H}}$
$\displaystyle=\hat{\mathcal{H}}_{\text{el}}^{(1)}+\hat{\mathcal{H}}_{\text{el}}^{\text{(2)}}+\hat{\mathcal{H}}_{\text{ph}}+\hat{\mathcal{H}}_{\text{el-
ph}},$ (1)

where

$\displaystyle\hat{\mathcal{H}}_{\text{el}}$
$\displaystyle=\hat{\mathcal{H}}_{\text{el}}^{\text{(1)}}+\hat{\mathcal{H}}_{\text{el}}^{\text{(2)}},$
(2) $\displaystyle\hat{\mathcal{H}}_{\text{el}}^{\text{(1)}}$
$\displaystyle=-t\sum_{\sigma\in\\{\uparrow,\downarrow\\}}\sum_{\langle
ij\rangle}\hat{a}_{i_{\sigma}}^{\dagger}\hat{a}_{j_{\sigma}},$ (3)
$\displaystyle\hat{\mathcal{H}}_{\text{el}}^{\text{(2)}}$
$\displaystyle=U\sum_{i}\hat{n}_{i_{\uparrow}}\hat{n}_{i_{\downarrow}},$ (4)
$\displaystyle\hat{\mathcal{H}}_{\text{ph}}$
$\displaystyle=\omega\sum_{i}\hat{b}_{i}^{\dagger}\hat{b}_{i},$ (5)

and

$\displaystyle\hat{\mathcal{H}}_{\text{el-ph}}$
$\displaystyle=-g\sum_{i}\hat{n}_{i}(\hat{b}_{i}+\hat{b}_{i}^{\dagger}).$ (7)

with

$\hat{n}_{i}=\sum_{\sigma\in\\{\uparrow,\downarrow\\}}\hat{n}_{i_{\sigma}}$
(8)

The nearest-neighbor electronic hopping is controlled by $t$ and the on-site
repulsion is characterized by the parameter $U$. The phonons are treated as
harmonic oscillators with a single frequency $\omega$. The electronic density
is coupled to the phonon degrees of freedom characterized by a coupling
constant $g$.

There are three relevant dimensionless parameters to define. The first is the
adiabaticity ratio in units of the hopping parameter

$\alpha=\frac{\omega}{t}.$ (9)

The second is the effective on-site repulsion in units of the hopping
parameter

$\frac{U}{t}.$ (10)

Lastly, we define the dimensionless el-ph coupling $\lambda$,

$\lambda=\frac{g^{2}}{2\text{d}t\omega},$ (11)

where d is the dimensionality of the system. When $U$ is the dominant
parameter, a spin density wave (SDW) phase similar to that found in the
Hubbard model is expected to arise. When $\lambda$ dominates, a charge density
wave (CDW) phase similar to that found in the Holstein model arises. A
metallic or superconducting phase can arise when the system transitions
between these two phases.Costa  _et al._ (2020)

###  II.2 Phonons in First Quantization

Since the Hamiltonian in Eq. 1 does not commute with the phonon number
operator $\hat{b}_{i}^{\dagger}\hat{b}_{i}$, the number of phonons in the
system is not conserved. Therefore, one needs to work with an infinitely large
phonon Hilbert space in order to study eigenstates of the Hubbard-Holstein
model. Methods working in a second quantized space such as DMRGJeckelmann and
White (1998) generally require a specification of the maximum number of
phonons a priori for the sake of computational tractability. Limiting the
maximum number of phonons effectively truncates the infinite Hilbert space,
which may introduce significant errors, particularly when $\alpha$ is small
and/or $\lambda$ is large.

For this reason we work within the framework of first quantization, namely
with position and momentum operators on each site $i$,

$\displaystyle\hat{X}_{i}$
$\displaystyle=\sqrt{\frac{1}{2m\omega}}(\hat{b}_{i}^{\dagger}+\hat{b}_{i}),$
(12) $\displaystyle\hat{P}_{i}$
$\displaystyle=i\sqrt{\frac{m\omega}{2}}(\hat{b}_{i}^{\dagger}-\hat{b}_{i}),$
(13)

and thus re-express

$\displaystyle\hat{\mathcal{H}}_{\text{ph}}$
$\displaystyle=\sum_{i}(\frac{m\omega^{2}}{2}\hat{X}_{i}^{2}+\frac{1}{2m}\hat{P}_{i}^{2}-\frac{\omega}{2}),$
(14) $\displaystyle\hat{\mathcal{H}}_{\text{el-ph}}$
$\displaystyle=-g\sqrt{2m\omega}\sum_{i}\hat{n}_{i}\hat{X}_{i},$ (15)

in the Hubbard-Holstein Hamiltonian in Eq. 1. Here, we introduced a fictitious
mass $m$ and throughout this work we use $m=1/\omega$. Working in a first
quantized space allows one to work directly at the complete basis set limit
for the phonons and avoids the issues posed by a truncated phonon Hilbert
space.

##  III Constrained-Path Auxiliary-Field Quantum Monte Carlo

AFQMC for mixed fermions and bosons was first formulated and studied by
Rubenstein, Zhang, and Reichman.Rubenstein  _et al._ (2012) In their
formulation, bosons are treated within a second-quantized framework.
Therefore, their approach would naturally suffer from the truncation of the
infinite bosonic Hilbert space if applied to the Hubbard-Holstein model. In
this work, we will reformulate the procedure to treat fermions in a second-
quantized space and bosons in a first-quantized space. Such a formulation is
closely related to that of Ref. 63, however our work allows the control of the
sign problem and introduces the full advantage of the CP-AFQMC approach in
treating the electronic degrees of freedom.

In AFQMC, as in other projector QMC methods,Blankenbecler and Sugar (1983) we
obtain the ground state via

$|\Psi_{0}\rangle\propto\lim_{\tau\rightarrow\infty}e^{-\tau\hat{\mathcal{H}}}|\Phi_{0}\rangle,$
(16)

where $|\Psi_{0}\rangle$ is the true ground state, $\tau$ denotes imaginary
time, and $|\Phi_{0}\rangle$ is a trial wavefunction with non-zero overlap
with the true ground state. Since $\hat{\mathcal{H}}$ involves both fermions
and bosons and so do the wavefunctions $|\Psi_{0}\rangle$ and
$|\Phi_{0}\rangle$, we represent these global vibronic wavefunctions as a
function of imaginary time $\tau$ in a mixed basis

$|\Psi(\tau)\rangle=\sum_{i}\omega_{i}|\psi_{i}(\tau),\mathbf{X}_{i}(\tau)\rangle,$
(17)

where $|\psi_{i}\rangle$ is the electronic wavefunction and
$|\mathbf{X}_{i}\rangle$ is a set of coordinates that represents the phonon
degrees of freedom. In our algorithm, these basis states each take a product
form:

$|\psi_{i}(\tau),\mathbf{X}_{i}(\tau)\rangle\equiv|\psi_{i}(\tau)\rangle\otimes|\mathbf{X}_{i}(\tau)\rangle\,,$
(18)

where $|\psi_{i}\rangle$ is a single Slater determinant. We will show below
that the projection process in Eq. (16) can be turned into a random walk in
the space of product states of the form defined in Eq. (18). We note that it
is also possible to work in momentum space
($|\mathbf{P}_{i}\rangle$),Hohenadler  _et al._ (2004) however it is more
convenient to work in position space here since it makes the application of
the e-ph coupling term straightforward.

We write the propagator for a finite timestep $\Delta\tau$ as

$\displaystyle\exp(-\Delta\tau\hat{\mathcal{H}})$
$\displaystyle=e^{-\Delta\tau\hat{\mathcal{H}}_{\text{el}}^{\text{(2)}}}e^{-\Delta\tau(\hat{\mathcal{H}}_{\text{el}}^{\text{(1)}}+\hat{\mathcal{H}}_{\text{el-
ph}})}e^{-\Delta\tau\hat{\mathcal{H}}_{\text{ph}}}$
$\displaystyle+\mathcal{O}(\Delta\tau^{2})$ (19)

using the standard first-order Trotter approximation. By virtue of the
Thouless theorem,Thouless (1960) $|\psi_{i}(\tau)\rangle$ remains a single
Slater determinant after propagation via
$\hat{\mathcal{H}}_{\text{el}}^{\text{(1)}}$ and $\hat{\mathcal{H}}_{\text{el-
ph}}$ (note the latter is diagonal in $|\mathbf{X}_{i}\rangle$ space).
$\hat{\mathcal{H}}_{\text{el}}^{\text{(2)}}$ is represented as a one-body
operator coupled to Ising variables and therefore a single Slater determinant
remains in the same manifold after propagation by
$e^{-\Delta\tau\hat{\mathcal{H}}_{\text{el}}^{\text{(2)}}}$. The phonon
propagation generated by $\hat{\mathcal{H}}_{\text{ph}}$ follows a commonly
used diffusion MC (DMC) algorithm.Umrigar  _et al._ (1993); Hammond  _et al._
(1994)

Before elaborating on the propagation more concretely, let us introduce
importance sampling, using a trial vibronic wavefunction $|\Psi_{T}\rangle$.
We re-write the global vibronic wavefunction in Eq. (17) in the following
form, Zhang and Krakauer (2003)

$|\Psi(\tau)\rangle=\sum_{i}\omega_{i}\frac{|\psi_{i}(\tau),\mathbf{X}_{i}(\tau)\rangle}{\langle\Psi_{T}|\psi_{i}(\tau),\mathbf{X}_{i}(\tau)\rangle}$
(20)

to perform imaginary-time propagation, namely

$|\Psi(\tau+\Delta\tau)\rangle=e^{-\Delta\tau\hat{\mathcal{H}}_{\text{el}}^{\text{(2)}}}e^{-\Delta\tau(\hat{\mathcal{H}}_{\text{el}}^{\text{(1)}}+\hat{\mathcal{H}}_{\text{el-
ph}})}e^{-\Delta\tau\hat{\mathcal{H}}_{\text{ph}}}|\Psi(\tau)\rangle.$ (21)

In the propagation of the phonon degrees of freedom we sample from the
distribution

$\displaystyle f[(\mathbf{X}(\tau+\Delta\tau)]$
$\displaystyle\equiv\frac{\langle\Psi_{T}|\psi(\tau),\mathbf{X}(\tau+\Delta\tau)\rangle}{\langle\Psi_{T}|\psi(\tau),\mathbf{X}(\tau)\rangle}$
$\displaystyle\times\langle
X(\tau+\Delta\tau)|e^{-\Delta\tau\hat{\mathcal{H}}_{\text{ph}}}|X(\tau)\rangle\,,$
(22)

where we have omitted the walker index $i$ in the subscript. One can derive
the following MC move for the updating of the variable
$\mathbf{X}(\tau)\equiv\mathbf{X}$:

$\mathbf{X}(\tau+\Delta\tau)=\mathbf{X}+\mathcal{N}(\mu=0,\sigma=\sqrt{m\Delta\tau})+\frac{\Delta\tau}{m}\frac{\nabla_{\mathbf{X}}\Psi_{T}(\mathbf{X})}{\Psi_{T}(\mathbf{X})}.$
(23)

Here $\mathcal{N}(\mu,\sigma)$ is a normally distributed random number with
mean $\mu$ and variance $\sigma^{2}$, and the last term is the so-called drift
term. Updates for the walker weights are carried out as

$w(\tau+\Delta\tau)=w(\tau)e^{{-\frac{\Delta\tau}{2}(E_{\text{ph}}(\tau+\Delta\tau)+E_{\text{ph}}(\tau)-2E_{\text{shift}})}},$
(24)

where $E_{\text{shift}}$ is a constant shift that can be adjusted to control
walker weight fluctuations, and we define

$E_{\text{ph}}(\tau)=\frac{\langle\Psi_{T}|\hat{H}_{\text{ph}}|\mathbf{X}(\tau)\rangle}{\Psi_{T}(\mathbf{X}(\tau))}.$
(25)

This algorithm is the same as the standard diffusion Monte Carlo
algorithm.Moskowitz  _et al._ (1982); Reynolds  _et al._ (1982); Umrigar  _et
al._ (1993); Hammond  _et al._ (1994)

Propagation arising from $\hat{\mathcal{H}}_{\text{el-ph}}$ is straightforward
to implement, since

$\sum_{i}\hat{n}_{i}\hat{X}_{i}|\psi(\tau),\mathbf{X}(\tau)\rangle=\sum_{i}\hat{n}_{i}X_{i}(\tau)|\psi(\tau),\mathbf{X}(\tau)\rangle.$
(26)

$\sum_{i}\hat{n}_{i}X_{i}(\tau)$ is thus a diagonal matrix in the single-
particle space with its $i$-th entry being $X_{i}(\tau)$. It is then
straightforward to exponentiate this matrix and apply it along with
$\hat{\mathcal{H}}_{\text{el}}^{(1)}$ to the Slater determinant.

Lastly, propagation generated by $\hat{\mathcal{H}}_{\text{el}}^{(2)}$ is the
same as that for the standard AFQMC algorithm for the Hubbard model. We employ
the discrete Hirsch spin decomposition for the two-body propagator Hirsch
(1983):

$e^{-\Delta\tau
U\hat{n}_{i_{\uparrow}}\hat{n}_{i_{\downarrow}}}=\frac{1}{2}e^{-\Delta\tau
U(\hat{n}_{i_{\uparrow}}+\hat{n}_{i_{\downarrow}})/2}\sum_{x_{i}=\pm
1}e^{\gamma x_{i}(\hat{n}_{i_{\uparrow}}-\hat{n}_{i_{\downarrow})}},$ (27)

where the constant $\gamma$ is determined by

$\cosh(\gamma)=e^{-\Delta\tau U/2}.$ (28)

For a given $\mathbf{x}$, the action of Eq. 27 on a single Slater determinant
keeps the Slater determinant in the single determinant manifold. In AFQMC we
keep track of the overlap between the walker wavefunction and a chosen trial
wavefunction. More specifically, we measure the overlap ratio of the $i$-th
walker,

$r_{i}=\frac{\langle\Psi_{T}|\psi_{i}(\tau+\Delta\tau),\mathbf{X}_{i}(\tau+\Delta\tau)\rangle}{\langle\Psi_{T}|\psi_{i}(\tau),\mathbf{X}_{i}(\tau+\Delta\tau)\rangle}.$
(29)

If $r_{i}$ is negative, the constraint condition is invoked and we set the
weight $w_{i}$ to zero, which then causes the walker to be removed from the
simulation. Furthermore, we apply heat-bath samplingZhang  _et al._ (1997)
using this ratio to importance sample the Ising variables. This completes the
description of our algorithm for the Hubbard-Holstein Hamiltonian.

The local energy evaluation at $\tau$ with the Hubbard-Holstein model is
straightforward via the one-body walker Green’s function

$G_{i_{\sigma}j_{\sigma}}(\tau)\equiv\frac{\langle\Psi_{T}|\hat{a}^{\dagger}_{i_{\sigma}}\hat{a}_{j_{\sigma}}|\psi(\tau),\mathbf{X}(\tau)\rangle}{\langle\Psi_{T}|\psi(\tau),\mathbf{X}(\tau)\rangle},$
(30)

and the two-body walker Green’s function,

$\Gamma_{i_{\uparrow}i_{\downarrow}}=\frac{\langle\Psi_{T}|\hat{a}^{\dagger}_{i_{\uparrow}}\hat{a}_{i_{\uparrow}}\hat{a}^{\dagger}_{i_{\downarrow}}\hat{a}_{i_{\downarrow}}|\psi(\tau),\mathbf{X}(\tau)\rangle}{\langle\Psi_{T}|\psi(\tau),\mathbf{X}(\tau)\rangle}.$
(31)

We will also need the mixed estimator for the phonon displacement,

$\langle\hat{X}_{i}\rangle(\tau)\equiv\frac{\langle\Psi_{T}|\hat{X}_{i}|\psi(\tau),\mathbf{X}(\tau)\rangle}{\langle\Psi_{T}|\psi(\tau),\mathbf{X}(\tau)\rangle}=X_{i}(\tau),$
(32)

and for the squared phonon momentum,

$\langle\hat{P}_{i}^{2}\rangle(\tau)=-\frac{\langle\Psi_{T}|\nabla_{X_{i}}^{2}|\psi(\tau),\mathbf{X}(\tau)\rangle}{\langle\Psi_{T}|\psi(\tau),\mathbf{X}(\tau)\rangle}\,,$
(33)

where $\nabla_{X_{i}}^{2}$ can be applied to the left on $|\Psi_{T}\rangle$.
Using these mixed estimators, the local energy can be evaluated as

$\displaystyle E_{L}$ $\displaystyle=-t\sum_{\sigma}\sum_{\langle
ij\rangle}G_{i_{\sigma}j_{\sigma}}+U\sum_{i}{\Gamma}_{i_{\uparrow}i_{\downarrow}}$
$\displaystyle+\sum_{i}(\frac{m\omega^{2}}{2}X_{i}^{2}+\frac{1}{2m}\langle\hat{P}_{i}^{2}\rangle-\frac{\omega}{2})$
$\displaystyle-g\sqrt{2m\omega}\sum_{i}(G_{i_{\uparrow}i_{\uparrow}}+G_{i_{\downarrow}i_{\downarrow}})X_{i}.$
(34)

##  IV Trial Wavefunctions

The choice of the trial wavefunction can affect the quality of the CP
approximation in treating the electronic degree of freedom. It can also affect
the computational efficiency in treating the electronic and especially the
phononic degrees of freedom; in particular, a poor choice of the importance
function can magnify or even introduce additional ergodicity issues,
especially in an el-ph system when multiple phonon modes are pronounced. It is
highly advantageous if an accurate trial wavefunction allows the overlap ratio
in Eq. 29, and the local energy in Eq. 34, to be efficiently evaluated.

###  IV.1 Semiclassical state

The simplest variational trial wavefunction that we employ in this work takes
a simple product form between electronic and bosonic degrees of freedom,

$|\Psi_{T}\rangle=|\psi_{T}\rangle\otimes|\phi_{T}\rangle,$ (35)

where $|\psi_{T}\rangle$ is a single determinant and $|\phi_{T}\rangle$ is a
coherent state (or a shifted harmonic oscillator state). This wave function
has been referred to as a “semiclassical state” in literature.Kalosakas  _et
al._ (1998); Romero  _et al._ (1998) Due to its simple product form, there is
no explicit entanglement between electrons and phonons. The electronic trial
wavefunction, $|\psi_{T}\rangle$, is parametrized by orbital rotation
$\boldsymbol{\theta}$,

$|\psi_{T}(\boldsymbol{\theta})\rangle=e^{\hat{\kappa}}|\psi_{0}\rangle,$ (36)

where

$\hat{\kappa}=\sum_{ij}(\theta_{ij}-\theta_{ji})\hat{a}_{i}^{\dagger}\hat{a}_{j},$
(37)

and $|\psi_{0}\rangle$ is some initial determinant (normally obtained by
diagonalizing the one-body electronic Hamiltonian). Single determinant trial
wavefunctions have been widely used in previous AFQMC studies of the Hubbard
model.Shi and Zhang (2013); LeBlanc  _et al._ (2015); Qin  _et al._ (2016a);
Zheng  _et al._ (2017)

The phonon trial wavefunction, $|\phi_{T}\rangle$, is parametrized by coherent
state displacements $\boldsymbol{\beta}$,

$|\phi_{T}(\boldsymbol{\beta})\rangle=e^{\sum_{i}\beta_{i}(\hat{b}_{i}^{\dagger}-\hat{b}_{i})}|0\rangle\equiv\hat{\mathcal{D}}(\vec{\beta})|0\rangle,$
(38)

where $\hat{\mathcal{D}}(\vec{\beta})$ is the displacement operator. We
optimize the energy of $|\Psi_{T}\rangle$ in Eq. 35 variationally over
$\boldsymbol{\theta}$ and $\boldsymbol{\beta}$ and use this as the final trial
wavefunction. $|\phi_{T}(\boldsymbol{\beta})\rangle$ technically contains
infinitely many bosons, but it has a convenient property which allows for an
efficient AFQMC algorithm

$\hat{b}_{i}|\phi_{T}(\boldsymbol{\beta})\rangle=\beta_{i}|\phi_{T}(\boldsymbol{\beta})\rangle.$
(39)

Using this fact, one can show that the projection of
$\langle\phi_{T}(\boldsymbol{\beta})|$ on to $|\mathbf{X}\rangle$ is

$\langle\phi_{T}(\boldsymbol{\beta})|\mathbf{X}\rangle=\prod_{i}\left(\frac{m\omega}{\pi}\right)^{\frac{1}{4}}e^{-\frac{m\omega}{2}(X_{i}-\sqrt{\frac{2}{m\omega}}\beta_{i})^{2}}.$
(40)

Similarly, the numerator of Eq. 33 is straightforward to evaluate as well
using

$-\langle\phi_{T}(\boldsymbol{\beta})|\nabla_{X_{i}}^{2}|\mathbf{X}\rangle=-\nabla_{X_{i}}^{2}\langle\phi_{T}(\boldsymbol{\beta})|\mathbf{X}\rangle.$
(41)

This semiclassical trial wavefunction therefore can be efficiently combined
with the AFQMC algorithm.

The variational energy of the semiclassical state can be obtained within the
Born-Oppenheimer (BO) approximation. After some algebra, it can be shown that
the lowest energy of a semiclassical state can be obtained by minimizing

$\langle\hat{H}_{\text{el}}\rangle-\frac{g^{2}}{\omega}\langle\sum_{i}\hat{n}_{i}^{2}\rangle$
(42)

over the variational parameters in $|\psi_{T}\rangle$. For a fixed $\lambda$,
variations in $\alpha$ do not change the energy of the semiclassical state.

The semiclassical state is exact in (1) the limit $g\rightarrow\infty$, (2)
the adiabatic limit $\omega\rightarrow 0$ (for a fixed $\lambda$) with
$U\rightarrow 0$, and (3) the atomic limit $U\rightarrow\infty$. When
$\lambda\rightarrow\infty$ (or $g\rightarrow\infty$ for a fixed $\omega$), the
use of a single semiclassical state is not problematic even though the BO
potential strongly bifurcates. This situation is physically similar to that of
the atomic limit of the Hubbard model ($U\rightarrow\infty$) where spin flips
do not cause an energy penalty and a degeneracy occurs amongst all possible
$2^{N}$ spin flips. Similarly, in the Holstein model, charge swapping does not
cause an energy penalty and the same macroscopic degeneracy occurs. In other
words, any one of the degenerate semiclassical states is equally well-suited
as an approximate wave function.

Aside from these limits, the semiclassical state itself can be inaccurate, but
we find that the subsequent AFQMC calculation with the semiclassical trial
wave function is often numerically exact. The most difficult parameter regime
for our AFQMC framework is when the Holstein coupling strength $g$ takes an
intermediate value. That is, $g$ is large enough that the el-ph correlation is
strong but is small enough that the macroscopic degeneracy does not occur. A
straightforward way to probe these situations is to increase $\omega$ for a
fixed $\lambda$ value. In this case, $g$ can be much larger than $t$ but is
always smaller than $\omega$ as long as $2\text{d}t\lambda<\omega$. In these
situations, the semiclassical state can be a poor choice of a trial wave
function in AFQMC, as we shall see. This is because no correlation between
electronic and bosonic degrees of freedom is built into this trial wave
function.

From a different point of view, the difficulty of semiclassical states was
noted in the work of Proville and Aubry, who defined the “quantumness” of the
phonons as Proville and Aubry (2000)

$\zeta=\frac{\omega^{2}}{g^{2}}=\frac{\omega}{2\text{d}t\lambda}.$ (43)

As $\zeta$ increases, the semiclassical state qualitatively fails.Proville and
Aubry (2000) This is consistent with the picture described above in that for a
given $\lambda$, both $g$ and $\zeta$ increase as $\omega$ increases. We
attribute the difficulties associated with semiclassical states to the
increase in correlation between electrons and phonons instead of the quantal
effect associated with the phonons alone.

###  IV.2 Multivibronic state

A linear combination of multiple semiclassical states can be used to correlate
electrons and phonons

$|\Psi_{\text{MS}}\rangle=\sum_{i}c_{i}|\psi_{T}(\boldsymbol{\theta}_{i})\rangle\otimes|\phi_{T}(\boldsymbol{\beta}_{i})\rangle,$
(44)

where one may determine $\mathbf{c}$, $\boldsymbol{\theta}$, and
$\boldsymbol{\beta}$ variationally. We refer to this wavefunction as a
multivibronic (MV) wavefunction. Similar to the multi-Slater determinant trial
states employed to study purely electronic systems.Shi and Zhang (2013) The MV
wavefunctions of this form would need exponentially many states for large
systems. Nevertheless, due to their simplicity, multivibronic states can be
valuable for the study of small systems.

Figure 1:  Born-Oppenheimer potential energy surfaces in units of $t$ for the
2-electron and 2-site Holstein model: (a) $\omega=t$, $\lambda=0.1t$, and
$g=0.447t$ and (b) $\omega=t$, $\lambda=t$, and $g=1.414t$. The minimum in (a)
is $E=-2.40t$ at ($X_{1}$=0.63, $X_{2}$=0.63) while the two minima in (b) are
$E=-8.25t$ at ($X_{1}$=3.90, $X_{2}$=0.10) and ($X_{1}$=0.10, $X_{2}$=3.90),
respectively.

A particular flavor of MV wavefunction that we focus on in this work is
closely tied to the underlying order encoded in the semiclassical states
themselves. Let us consider the two-electron, two-site Holstein model. It is
well known that for the Holstein model at a large coupling $\lambda$, the BO
surface bifurcates and develops into a double well potential.Holstein (1959a,
b) For the two-site problem, the BO potential energy surfaces (PESs) are
characterized by

$\displaystyle\hat{H}_{\text{BO}}(X_{1},X_{2})$
$\displaystyle=\hat{H}_{\text{el}}-g\sqrt{2m\omega}(\hat{n}_{1}X_{1}+\hat{n}_{2}X_{2})$
$\displaystyle+\frac{m\omega^{2}}{2}(X_{1}^{2}+X_{2}^{2}),$ (45)

where $X_{1}$ and $X_{2}$ are constant scalars denoting the coordinates of the
classical phonons. We can find the ground-state electronic wavefunction of the
Hamiltonian in Eq. 45 by exactly diagonalizing it and forming a potential
energy surface for each combination of ($X_{1}$,$X_{2}$). In Fig. 1, a
representative example of the BO PESs is given. Fig. 1(a) illustrates an
example of the weak coupling case, where the coherent states have the same
centers for all sites and no charge modulation occurs. In such cases, the
minimum BO state (i.e. the semiclassical state) is an excellent variational
wavefunction. In Fig. 1(b) there are two distinct minima with equal BO
energies. Here, a wave function of a single semiclassical state with a
Gaussian function centered at one of the two BO minima in position space would
not provide a good description of the system. When used as an importance
function, it can introduce or exacerbate ergodicity problems in the Monte
Carlo sampling and induce a large or even infinite variance in the energy
fluctuations.

We propose the following improved variational wavefunction in this situation
for the Holstein model on a bipartite lattice. At half-filling there are two
exactly degenerate semiclassical states. In particular, one state is
characterized by

$\beta_{i}^{(1)}=\begin{cases}\beta_{e},&\text{if}\ i\ \text{on A
sublattice}\\\ \beta_{o},&\text{if}\ i\ \text{on B sublattice},\end{cases}$
(46)

where $i$ is a site index. The other degenerate solution is given by switching
the A and B sublattices. One can smoothly interpolate between the two states
by defining a convex combination

$\vec{\beta}(\alpha)=\alpha\vec{\beta}^{(1)}+(1-\alpha)\vec{\beta}^{(2)}$ (47)

for $\alpha\in[0,1]$. For each $\vec{\beta}(\alpha)$, we find a single
determinant that minimizes the energy of a single semiclassical state. One can
take a linear combination of all of these semiclassical states along the line
that interpolates two solutions to form a MV wavefunction. We refer this to as
the Thouless path (TP) wavefunction, $|\Psi_{\text{TP}}\rangle$,

$|\Psi_{\text{TP}}\rangle=\sum_{\alpha}c_{\alpha}|\Psi_{T}(\beta(\alpha))\rangle,$
(48)

where $c_{\alpha}$ is determined by variationally minimizing the energy. The
construction of the TP wavefunction can be generalized to arbitrary filling
fractions and number of sites because different filling fractions simply give
rise to ordered states with different wavelengths. The cost for its
construction is negligible compared to the optimization of a semiclassical
state. Its use in AFQMC as a trial wave function simply introduces a prefactor
depending on the number of states included in Eq. 48. We will refer to a TP
wave function with $n$ interpolation values of $\alpha$ as TP($n$). While TP
wave functions provide a simple and accurate importance function for
bifurcated potential energy surfaces, they also become inaccurate when the
correlation between electrons and phonons becomes strong.

###  IV.3 Variational Lang-Firsov trial wavefunctions

A simple, widely used way to incorporate correlation effects between electrons
and phonons is to use the polaron transformation or the Lang-Firsov (LF)
transformation,Lang and Firsov (1963)

$|\Psi_{\text{LF}}\rangle=\hat{U}_{\text{LF}}(\boldsymbol{\xi})|\psi_{T}(\boldsymbol{\theta})\rangle\otimes|\phi_{T}(\boldsymbol{\beta})\rangle,$
(49)

where

$\hat{U}_{\text{LF}}(\boldsymbol{\xi})=e^{\frac{1}{\sqrt{2}}\sum_{i}\xi_{i}\hat{n}_{i}(\hat{b}^{\dagger}_{i}-\hat{b}_{i})},$
(50)

and the set $\boldsymbol{\xi}$ are referred to as the LF amplitudes which are
variational parameters along with $\boldsymbol{\theta}$ and
$\boldsymbol{\beta}$. From the wavefunction viewpoint, Eq. 49 provides a way
to explicitly build a wavefunction with non-perturbative el-ph correlation on
top of semiclassical states via a unitary transformation. Typical LF
implementations involve the phonon vacuum state as opposed to the coherent
state in Eq. 49. We find that having the coherent state provides additional
variational flexibility and thereby yields lower energies compared to those
that use the vacuum state. Since it does not complicate the underlying
optimization problem, we use the coherent state as written in Eq. 49.

While the details of the LF transformation and its variational optimization
have been well-documented,Lang and Firsov (1963); Silbey and Harris (1984);
Harris and Silbey (1985); Pouthier (2013) we briefly summarize them to provide
a self-contained description. Our goal is to simultaneously determine
$\boldsymbol{\xi}$, $\boldsymbol{\theta}$, and $\boldsymbol{\beta}$
variationally. To carry this out, we find that it is simpler to work with the
unitary-transformed Hamiltonian, $\hat{H}^{\text{LF}}$, based on
$\hat{U}_{\text{LF}}$, and optimize the variational energy of
$\hat{H}^{\text{LF}}$ evaluated with the semiclassical wavefunction. We start
from

$\hat{U}_{\text{LF}}(\boldsymbol{\xi})^{\dagger}a_{i_{\sigma}}^{\dagger}a_{j_{\sigma}}\hat{U}_{\text{LF}}(\boldsymbol{\xi})=a_{i_{\sigma}}^{\dagger}a_{j_{\sigma}}e^{(-\xi_{i}(b^{\dagger}_{i}-b_{i})+\xi_{j}(b^{\dagger}_{j}-b_{j}))}\,$
(51)

and

$\displaystyle\hat{U}_{\text{LF}}(\boldsymbol{\xi})^{\dagger}b_{i}^{\dagger}\hat{U}(\boldsymbol{\xi})$
$\displaystyle=b_{i}^{\dagger}-\xi_{i}\hat{n}_{i},$ (52)
$\displaystyle\hat{U}_{\text{LF}}(\boldsymbol{\xi})^{\dagger}b_{i}\hat{U}(\boldsymbol{\xi})$
$\displaystyle=b_{i}-\xi_{i}\hat{n}_{i}.$ (53)

The LF transformed Hamiltonian reads

$\displaystyle\hat{{\mathcal{H}}}_{\text{el}}^{\text{LF}}$
$\displaystyle=-t\sum_{\sigma}\sum_{\langle
ij\rangle}a_{i_{\sigma}}^{\dagger}a_{j_{\sigma}}e^{(-\frac{\xi_{i}}{\sqrt{2}}(b^{\dagger}_{i}-b_{i})+\frac{\xi_{j}}{{\sqrt{2}}}(b^{\dagger}_{j}-b_{j}))}$
$\displaystyle+U\sum_{i}\hat{n}_{i_{\uparrow}}\hat{n}_{i_{\downarrow}},$ (54)
$\displaystyle\hat{{\mathcal{H}}}_{\text{ph}}^{\text{LF}}$
$\displaystyle=\omega\sum_{i}(\hat{b}^{\dagger}_{i}+\frac{\xi_{i}}{\sqrt{2}}\hat{n}_{i})(\hat{b}_{i}+\frac{\xi_{i}}{\sqrt{2}}\hat{n}_{i}),$
(55)

and

$\displaystyle\hat{{\mathcal{H}}}_{\text{el-ph}}^{\text{LF}}$
$\displaystyle=-g\sum_{i}\hat{n}_{i}(\hat{b}_{i}+\hat{b}_{i}^{\dagger}+\sqrt{2}\xi_{i}\hat{n}_{i}).$
(57)

All of the energy terms are straightforward to evaluate with semiclassical
trial wavefunctions. The electronic kinetic energy is more complex than its
bare Hamiltonian counterpart due to the presence of exponential bosonic
operators, so we provide more details here. To utilize Eq. 39, we write the
exponential term in the kinetic energy operator as

$\displaystyle
e^{\frac{1}{\sqrt{2}}(t_{j}(\hat{b}_{j}^{\dagger}-\hat{b}_{j})-t_{i}(\hat{b}_{i}^{\dagger}-\hat{b}_{i}))}$
$\displaystyle=e^{\frac{1}{\sqrt{2}}(t_{j}\hat{b}_{j}^{\dagger}-t_{i}\hat{b}_{i}^{\dagger})}e^{\frac{1}{\sqrt{2}}(t_{i}\hat{b}_{i}-t_{j}\hat{b}_{j})}e^{-\frac{1}{4}(t_{i}^{2}+t_{j}^{2})},$
(58)

where we have used

$e^{\frac{1}{\sqrt{2}}t_{j}(\hat{b}_{j}^{\dagger}-\hat{b}_{j})}=e^{\frac{1}{\sqrt{2}}t_{j}\hat{b}_{j}^{\dagger}}e^{-\frac{1}{\sqrt{2}}t_{j}\hat{b}_{j}}e^{-\frac{1}{4}t_{j}^{2}}.$
(59)

The expectation value of Eq. 58 is simple to evaluate with the semiclassical
state of Eq. 35

$\displaystyle\langle\Psi_{T}|e^{\frac{1}{\sqrt{2}}(t_{j}\hat{b}_{j}^{\dagger}-t_{i}\hat{b}_{i}^{\dagger})}e^{\frac{1}{\sqrt{2}}(t_{i}\hat{b}_{i}-t_{j}\hat{b}_{j})}|\Psi_{T}\rangle
e^{-\frac{1}{4}(t_{i}^{2}+t_{j}^{2})}$
$\displaystyle=e^{\frac{1}{\sqrt{2}}(t_{j}\beta_{j}-t_{i}\beta_{i})}e^{\frac{1}{\sqrt{2}}(t_{i}\beta_{i}-t_{j}\beta_{j})}e^{-\frac{1}{4}(t_{i}^{2}+t_{j}^{2})}.$
(60)

The variational LF wavefunction is expected to be more accurate than the
semiclassical state due to the explicit correlation between electrons and
phonons. Furthermore, the limit of $\omega\rightarrow\infty$ which is
difficult for simple semiclassical wavefunctions to treat, can be exactly
treated by the LF wavefunction, because the el-ph coupling term in
$\hat{H}^{\text{LF}}$ can be removed by setting $\xi_{i}=\sqrt{2}g/\omega$.
Due to the fact that phonon displacements are significantly penalized in this
limit, the variational optimization over $\boldsymbol{\beta}$ naturally yields
$\boldsymbol{\beta}=\boldsymbol{0}$. Therefore, the bosonic operators in the
hopping amplitude in Eq. 54 all vanish. Provided that one can handle the
remaining electronic Hamiltonian terms exactly, the variational LF
wavefunction should be exact in this limit. We note that for many-electron
systems in the $\omega\rightarrow\infty$ limit, the LF Hamiltonian takes the
same form as the attractive Hubbard model, which is another sign-free lattice
model that can be efficiently simulated in AFQMC.Shi and Zhang (2016, 2017)

Despite these desirable properties, there seems to be no simple and general
way to use this wavefunction in AFQMC without invoking a major increase in
scaling. As an exception to this, we mention here the work of Hohenadler and
co-workersHohenadler  _et al._ (2004) where a QMC algorithm with the LF
Hamiltonian was formulated for single electron problems. It was demonstrated,
however, that the transformed electronic Hamiltonian in Eq. 54 creates a
complex phase problem.

Therefore, we briefly investigate a simpler linearized LF (LLF) wavefunction
of the form,

$|\Psi_{\text{LLF}}\rangle=(1+\frac{1}{\sqrt{2}}\sum_{i}\xi_{i}\hat{n}_{i}\hat{b}_{i}^{\dagger})|\psi_{T}(\boldsymbol{\theta})\rangle\otimes|\phi_{T}(\boldsymbol{\beta})\rangle,$
(61)

where we have omitted a term that is proportional to $\hat{n}_{i}\hat{b}_{i}$
since the action of $\hat{b}_{i}$ on $|\phi_{T}(\boldsymbol{\beta})\rangle$ is
trivial due to Eq. 39. We variationally optimize $\boldsymbol{\xi}$ in Eq. 61
to maximize the accuracy of the LLF trial wavefunction. The AFQMC algorithm
presented in Section III can be efficiently implemented for Eq. 61.

It is possible to formulate a simple extension of the LLF wavefunction in the
spirit of the TP wavefunction:

$|\Psi_{\text{TP-
LLF}}\rangle=\sum_{\alpha=1}^{n}c_{\alpha}|\Psi_{\text{LLF}}^{(\alpha)}\rangle,$
(62)

where each of the $|\Psi_{\text{LLF}}^{(\alpha)}\rangle$ terms has its own
variational parameters. Following the discussion of the TP trial wavefunction,
it may be possible to determine these variational parameters via a convex
interpolation of $\boldsymbol{\beta}$ and $\boldsymbol{\xi}$ as in Eq. 47. We
refer to this wavefunction as the TP-LLF($n$) wavefunction which goes beyond
both the TP$(n)$ and the LLF wavefunctions in sophistication.

In contrast with the LF form, a trial wave function with an el-ph Jastrow
factor can be used more straightforwardly in AFQMC, since $\hat{X}_{i}$
operators are involved in the exponent instead of $\hat{P}_{i}$ as in LF. The
unitary transformation in the LF wave function can be thought of as a simple
Jastrow factor that encodes correlation between the electronic density and the
phonon momentum on a site. However, $\hat{U}_{\text{LF}}(\boldsymbol{\xi})$ is
unitary and we thus expect this transformation to behave differently from
Jastrow factors in el-ph problems.Ohgoe and Imada (2014, 2017); Karakuzu  _et
al._ (2017) (It is also different from the coupled-cluster operators
considered in recent studies of el-ph problems .Sibaev  _et al._ (2020);
Dresselhaus  _et al._ (2020); White  _et al._ (2020)) Given the performance
improvement with the LLF trial wave function (as discussed below), we expect
an el-ph Jastrow trial wave function will greatly reduce the difficulties in
parameter regimes with strong el-ph coupling, and result in a major
improvement in our AFQMC approach. We leave the implementation and systematic
studies using an el-ph Jastrow trial wavefunction in AFQMC for future work.

###  IV.4 Additional details

The semiclassical state in Eq. 35 can describe two competing phases, SDW and
CDW. To obtain the variational wavefunction for these two distinct states, we
employ the following protocol:

  1. 1.

For a CDW state, we perform a variational optimization of a semiclassical
state with spin-restriction. Due to the spin-restriction, any states that
arise from minimization are not capable of describing SDW order.

  2. 2.

For an SDW state, we perform a variational optimization of a spin-unrestricted
Hartree-Fock (UHF) wavefunction to minimize the electronic energy. Once a UHF
state is obtained, we determine the shift vector $\boldsymbol{\beta}$
variationally while fixing the electronic degrees of freedom. As long as the
UHF state exhibits SDW order, such a coupled el-ph semiclassical state with
exhibit the same SDW order. We have used an ad hoc effective repulsion
strength ($U_{\text{eff}}/t$) of 0.5 Shi and Zhang (2013) in our UHF
calculations to obtain SDW trial wavefunctions for the Hubbard-Holstein model
in this work. The CP-AFQMC results are not sensitive to this particular
choice. (We note that it is possible to determine this effective repulsion
strength via a self-consistent procedure with CP-AFQMC Qin  _et al._ (2016b).)

##  V Perturbation Theory

###  V.1 Coherent State Møller-Plesset Perturbation Theory

It is instructive to consider low-order perturbation theory for comparison to
numerically exact approaches. We first note that

$\displaystyle\hat{\mathcal{D}}(\vec{\beta})^{\dagger}b_{i}^{\dagger}\hat{\mathcal{D}}(\vec{\beta})$
$\displaystyle=b_{i}^{\dagger}+\beta_{i},$ (63)
$\displaystyle\hat{\mathcal{D}}(\vec{\beta})^{\dagger}b_{i}\hat{\mathcal{D}}(\vec{\beta})$
$\displaystyle=b_{i}+\beta_{i}.$ (64)

Using this property, we write

$\displaystyle\hat{\bar{\mathcal{H}}}_{\text{ph}}\equiv\hat{\mathcal{D}}(\vec{\beta})^{\dagger}\hat{\mathcal{H}}_{\text{ph}}\hat{\mathcal{D}}(\vec{\beta})$
$\displaystyle=\omega\sum_{i}(b_{i}^{\dagger}+\beta_{i})(b_{i}+\beta_{i}),$
(65)

and

$\displaystyle\hat{\bar{\mathcal{H}}}_{\text{el-
ph}}\equiv\hat{\mathcal{D}}(\vec{\beta})^{\dagger}\hat{\mathcal{H}}_{\text{el-
ph}}\hat{\mathcal{D}}(\vec{\beta})$
$\displaystyle=-g\sum_{i}\hat{n}_{i}(\hat{b}_{i}+\hat{b}_{i}^{\dagger}+2\beta_{i}).$
(67)

Thus we have

$\displaystyle\hat{\bar{\mathcal{H}}}=\hat{\mathcal{H}}_{\text{el}}+\hat{\bar{\mathcal{H}}}_{\text{ph}}+\hat{\bar{\mathcal{H}}}_{\text{el-
ph}}.$ (68)

We note that the following zeroth-order Hamiltonian naturally has the
semiclassical state of Eq. 35 as its ground state

$\hat{\mathcal{H}}_{0}=\hat{\mathcal{F}}+\omega\sum_{i}(\beta_{i}^{2}+\hat{b}_{i}^{\dagger}\hat{b}_{i}),$
(69)

where $\hat{\mathcal{F}}$ is the Fock operator defined as (for spin
$\sigma=\uparrow$ or $\downarrow$)

$\displaystyle\hat{\mathcal{F}}_{\sigma}$
$\displaystyle=\hat{\mathcal{F}}_{\sigma}^{\text{el}}-2g\sum_{i}\hat{n}_{i_{\sigma}}\beta_{i},$
(70)

with the electronic Fock operators:

$\displaystyle\hat{\mathcal{F}}_{\uparrow}^{\text{el}}$
$\displaystyle=-t\sum_{\langle
ij\rangle}\hat{a}_{i_{\uparrow}}^{\dagger}\hat{a}_{j_{\uparrow}}+U\sum_{i}\hat{n}_{i_{\uparrow}}\langle\hat{n}_{i_{\downarrow}}\rangle_{\psi_{T}},$
(71) $\displaystyle\hat{\mathcal{F}}_{\downarrow}^{\text{el}}$
$\displaystyle=-t\sum_{\langle
ij\rangle}\hat{a}_{i_{\downarrow}}^{\dagger}\hat{a}_{j_{\downarrow}}+U\sum_{i}\hat{n}_{i_{\downarrow}}\langle\hat{n}_{i_{\uparrow}}\rangle_{\psi_{T}}.$
(72)

where

$\langle\hat{n}_{i_{\sigma}}\rangle_{\psi_{T}}=\frac{\langle\psi_{T}|\hat{n}_{i_{\sigma}}|\psi_{T}\rangle}{\langle\psi_{T}|\psi_{T}\rangle}$
(73)

It is straightforward to show that $|\Psi_{T}\rangle$ is an eigenstate of
$\hat{\mathcal{H}}_{0}$. From this starting point, one can develop an order-
by-order perturbation theory to capture all of the correlation effects among
electrons and between electrons and phonons built through

$\displaystyle\hat{V}$
$\displaystyle=\hat{\bar{\mathcal{H}}}-\hat{\mathcal{H}}_{0}$
$\displaystyle=(\hat{\mathcal{H}}_{\text{el}}-\hat{\mathcal{F}}_{\text{el}})+\sum_{i}(\omega\beta_{i}-g\hat{n}_{i})(\hat{b}_{i}+\hat{b}_{i}^{\dagger}).$
(74)

We note that such a partitioning of the Hamiltonian resembles the widely used
Møller-Plesset (MP) perturbation theory in quantum chemistry.Shavitt and
Bartlett (2009) We refer this perturbation theory to as “coherent state
Møller-Plesset perturbation theory” (CSMP) since a coherent state (or a
semiclassical state) is an eigenstate of the zeroth-order Hamiltonian. This
was also recently discussed in the work of White and co-workers in the context
of coupled-cluster theory.White  _et al._ (2020)

Similar to MP, CSMP recovers the energy of the semiclassical state with the
first-order perturbation correction,

$\displaystyle
E^{(0)}+E^{(1)}=\langle\hat{\mathcal{H}}_{\text{el}}\rangle_{\psi_{T}}+\omega\sum_{i}\beta_{i}^{2}-2g\sum_{i}\langle\hat{n}_{i}\rangle_{\psi_{T}}\beta_{i},$
(75)

where we have defined

$\langle\hat{n}_{i}\rangle_{\psi_{T}}\equiv\langle\psi_{T}|\hat{n}_{i_{\uparrow}}+\hat{n}_{i_{\downarrow}}|\psi_{T}\rangle.$
(76)

In this work, we are interested in comparing the second-order perturbation
theory (CSMP2) with AFQMC. The evaluation of the CSMP2 energy is most natural
in the molecular orbital (MO) basis rather than in the site basis. The MO
basis is defined by a set of orbitals, $\\{\psi_{i_{\sigma}}\\}$, that satisfy

$\hat{\mathcal{F}}_{\sigma}\psi_{p_{\sigma}}=\epsilon_{p_{\sigma}}\psi_{p_{\sigma}},$
(77)

where $\epsilon_{p_{\sigma}}$ is the $p$-th MO energy and the $p$-th MO,
$\psi_{p_{\sigma}}$, is expanded via a set of site orbitals,
$\\{\phi_{\mu_{\sigma}}\\}$,

$\psi_{p_{\sigma}}=\sum_{\mu}C_{\mu p_{\sigma}}\phi_{\mu_{\sigma}}.$ (78)

We then transform Eq. 74 from the site basis to the MO basis using the
coefficient matrix $\mathbf{C}$ for each spin

$\displaystyle\hat{\bar{V}}$
$\displaystyle=\sum_{\sigma\in\\{\uparrow,\downarrow\\}}\sum_{pq}\sum_{\mu}\left((\omega\beta_{\mu}-g_{pq})\hat{a}_{p_{\sigma}}^{\dagger}\hat{a}_{q_{\sigma}}\hat{b}_{\mu}^{\dagger}+\text{h.c.}\right)$
$\displaystyle+\sum_{pqrs}U_{p_{\uparrow}q_{\downarrow}r_{\downarrow}s_{\uparrow}}\hat{a}_{p_{\uparrow}}^{\dagger}\hat{a}_{q_{\downarrow}}^{\dagger}\hat{a}_{r_{\downarrow}}\hat{a}_{s_{\uparrow}},$
(79)

where

$\displaystyle g_{pq}$ $\displaystyle=g\sum_{\mu}(C_{\mu p})^{*}C_{\mu q},$
(80) $\displaystyle U_{p_{\uparrow}q_{\downarrow}r_{\downarrow}s_{\uparrow}}$
$\displaystyle=U\sum_{\mu}(C_{\mu p_{\uparrow}})^{*}(C_{\mu
q_{\downarrow}})^{*}C_{\mu r_{\downarrow}}C_{\mu s_{\uparrow}}.$ (81)

The CSMP2 energy expression follows in a spin-orbital MO basis,

$\displaystyle E^{(2)}$
$\displaystyle=-\sum_{\mu=1}^{M}\frac{(\omega\beta_{\mu}-g\sum_{i}(C_{\mu
i})^{*}C_{\mu i})^{2}}{\omega}$
$\displaystyle-\sum_{ia}\frac{-|g_{ai}|^{2}}{\omega+\epsilon_{a}-\epsilon_{i}}$
$\displaystyle-\sum_{i_{\uparrow}a_{\uparrow}}\sum_{j_{\downarrow}b_{\downarrow}}\frac{|U_{i_{\uparrow}j_{\downarrow}b_{\downarrow}a_{\uparrow}}|^{2}}{\epsilon_{a_{\uparrow}}+\epsilon_{b_{\downarrow}}-\epsilon_{i_{\uparrow}}-\epsilon_{j_{\downarrow}}},$
(82)

where the orbital energies $\\{\epsilon_{p}\\}$ are eigenvalues of the Fock
operator $\hat{\mathcal{F}}$. We note that the first term in Eq. 82 is zero if
the semiclassical reference state is fully optimized.

###  V.2 Lang-Firsov Perturbation Theory

It may be useful to develop a second-order perturbation theory from a
reference state given by the Lang-Firsov transformationLang and Firsov (1963)
in Eq. 49. In the spirit of the original LF transformation,Lang and Firsov
(1963) we set

$\xi_{i}=\frac{\sqrt{2}g}{\omega},$ (83)

which removes the Holstein coupling in the transformed framework. With this
choice of the LF amplitudes instead of the variational LF amplitudes, the
transformed Hamiltonian is simplified to

$\displaystyle\hat{\mathcal{H}}^{\text{LF}}$
$\displaystyle=\omega\sum_{i}\hat{b}_{i}^{\dagger}\hat{b}_{i}-\frac{g^{2}}{\omega}\sum_{i}\hat{n}_{i}+(U-\frac{2g^{2}}{\omega})\sum_{i}\hat{n}_{i_{\uparrow}}\hat{n}_{i_{\downarrow}}$
$\displaystyle-t\sum_{\sigma}\sum_{\langle
ij\rangle}a_{i_{\sigma}}^{\dagger}a_{j_{\sigma}}e^{(-\frac{g}{\omega}(b^{\dagger}_{i}-b_{i})+\frac{g}{\omega}(b^{\dagger}_{j}-b_{j}))}.$
(84)

For convenience, we rewrite

$e^{(-\frac{g}{\omega}(b^{\dagger}_{i}-b_{i})+\frac{g}{\omega}(b^{\dagger}_{j}-b_{j}))}=e^{-\frac{g^{2}}{\omega^{2}}}e^{-\frac{g}{\omega}(b^{\dagger}_{i}-b^{\dagger}_{j})}e^{\frac{g}{\omega}(b_{i}-b_{j})}.$
(85)

Based on the approach of Bonca, Trugman, and co-workers,Bonča  _et al._
(2000); Bonča and Trugman (2000) we choose the zeroth order Hamiltonian as

$\hat{\mathcal{H}}_{0}^{\text{LF}}=\omega\sum_{i}\hat{b}_{i}^{\dagger}\hat{b}_{i}-\frac{g^{2}}{\omega}\sum_{i}\hat{n}_{i}+(U-\frac{2g^{2}}{\omega})\sum_{i}\hat{n}_{i_{\uparrow}}\hat{n}_{i_{\downarrow}},$
(86)

with the perturbing Hamiltonian

$\hat{V}^{\text{LF}}=-te^{-\frac{g^{2}}{\omega^{2}}}\sum_{\sigma}\sum_{\langle
ij\rangle}a_{i_{\sigma}}^{\dagger}a_{j_{\sigma}}e^{-\frac{g}{\omega}(b^{\dagger}_{i}-b^{\dagger}_{j})}e^{\frac{g}{\omega}(b_{i}-b_{j})}.$
(87)

For concreteness, we consider the specific case of second-order perturbation
theory for a two-electron, two-site problem. The ground state of
$\hat{\mathcal{H}}_{0}^{\text{LF}}$ is either
$|\uparrow_{1}\downarrow_{1}\rangle\otimes|0,0\rangle$ or
$|\uparrow_{2}\downarrow_{2}\rangle\otimes|0,0\rangle$ with an energy
$U-\frac{4g^{2}}{\omega}$. We start from an unperturbed state,

$|\Psi^{(0)}_{0}\rangle=\frac{1}{\sqrt{2}}(|\uparrow_{1}\downarrow_{1}\rangle\otimes|0,0\rangle,+|\uparrow_{2}\downarrow_{2}\rangle\otimes|0,0\rangle),$
(88)

noting that an excited state which can be connected to the unperturbed ground
state via $\hat{V}^{\text{LF}}$ takes the form

$|\Psi_{mn}^{(0)}\rangle=\frac{1}{\sqrt{2}}(|\uparrow_{1}\downarrow_{2}\rangle-|\downarrow_{1}\uparrow_{2}\rangle)\otimes|m,n\rangle.$
(89)

It is then easy to show that the first-order energy correction to the
unperturbed state is zero. The second-order energy reads

$\displaystyle E^{(2)}_{\text{LFPT}}=-2t^{2}e^{-\frac{2g^{2}}{\omega^{2}}}$
$\displaystyle\times\sum_{m=0}^{\infty}\sum_{n=0}^{\infty}\frac{(g/\omega)^{2(m+n)}}{m!n!}\frac{1+(-1)^{m+n}}{(m+n)\omega-U+2\frac{g^{2}}{\omega}}.$
(90)

We numerically evaluate this expression in a brute-force manner, observing
that a maximum $n$ of 200 is sufficient to converge the energy to machine
precision. We note that this expression differs from that of Bonca, Trugman
and co-workers Bonča  _et al._ (2000); Bonča and Trugman (2000) since in their
work a single bipolaron was considered in the continuum limit (infinite
lattice) whereas in our work we focus on a 2-site problem.

##  VI The Holstein Model

To study the behavior of the proposed AFQMC algorithm with simple trial
wavefunctions such as the semiclassical, LLF, and TP wavefunctions, we shall
investigate the 1D and 2D Holstein models first, namely we set $U=0$ in Eq. 1.

For the Holstein model, the sign problem is absent, as is well-known in the
determinant quantum Monte Carlo (DQMC) approach Batrouni and Scalettar (2019).
Similar to DQMC, the overlap function
$\langle\Psi_{T}|\psi_{i}(\tau),\mathbf{X}_{i}(\tau)\rangle$ in Eq. 20 remains
non-negative throughout the imaginary-time propagation, since the phonon
component, $\phi_{T}(\mathbf{X}(\tau))$ (omitting the walker index again), is
positive everywhere, and the electronic component,
$\langle\psi_{T}|\psi\rangle=|\langle\psi_{T,\uparrow}|\psi_{\uparrow}\rangle|^{2}$
with a spin-restricted form, is also non-negative. Thus, in the Holstein model
the difference between our approach and DQMC is primarily in the way the Monte
Carlo sampling is conducted. AFQMC uses a branching random walk with a
population of walkers to construct the imaginary-time path iteratively, as we
have described, while DQMC treats the entire path as a path integral or
worldline, and updates it by sweeping different imaginary-time locations via a
Metropolis-like algorithm. A second difference is the introduction of an
importance functions in our approach via the similarity transformation, as
indicated in Eqs. (22) and (29). These factors can affect the behavior of the
Monte Carlo sampling, and yield different performances in different regimes of
the parameter space, including efficiency, auto-correlation time, and possibly
different levels of difficulty with ergodicity. The examples in the Holstein
model below serve as a first test of the AFQMC method in this context.

###  VI.1 2-Electron 2-site model

We start with this small problem where we easily can compare results against
exact diagonalization (ED). Since there are only two sites in our model, we
compute energies with open boundary conditions (OBCs).

Figure 2:  Error in the total energy in units of $t$ for the 2-site 2-electron
Holstein model as a function of $\omega$ for various $\lambda$ values: (a)
CSMP2, (b) LFPT2 and (c) AFQMC results. For $\lambda=0.1$, LFPT2 energy errors
lie outside the plotted range. In (c), for $\lambda\geq 0.5$ AFQMC/TP(11)
results are shown while for $\lambda=0.1$ we present AFQMC/S results. Note the
different vertical scales in panel (c).

In Fig. 2, we present the error in the total energy of CSMP2, LFPT2 and AFQMC
compared to ED. Understanding the behaviors of the two flavors of perturbation
theory helps gauge non-perturbative effects in our system. In Fig. 2(a), it is
clear that the CSMP2 energy becomes more inaccurate as we increase $\omega$.
This is because the zeroth order wavefunction, a semiclassical trial
wavefunction, starts to degrade when increasing el-ph correlation. Perhaps the
most striking behavior to note concerning CSMP2 is that this approach performs
worst for intermediate $\lambda$ values, (e.g. $\lambda=0.5$) and is in fact
more accurate for larger $\lambda$ values such as $\lambda=2.0$. This can also
be understood in terms of the increase in el-ph correlation as explained in
Section IV.1. On the other hand, LFPT2 in Fig. 2(b) is comparatively more
accurate than CSMP2 for $\lambda$ values larger than $\lambda=0.5$. The LF
reference state (namely two electrons occupying one site and with a phonon
vacuum state) is qualitatively incorrect when the el-ph coupling is small. In
such cases, we cannot treat the the hopping term perturbatively. This is
clearly reflected in Fig. 2(b) as LFPT2 exhibits large errors for small
$\lambda$ values. As LFPT2 is well known to produce accurate results for
strong coupling, it is remarkable that a weak coupling perturbation theory,
CSMP2, performs equally well even at $\lambda=2.0$.

We also show the performance of AFQMC for the same 2-site Holstein dimer. The
error of AFQMC is shown on a much smaller (100 times) scale. We have tested
both the single semiclassical state and TP($11$), i.e., the TP state with a
superposition of 11 semiclassical states, as trial wave function. These are
referred to as AFQMC/S and AFQMC/TP($11$), respectively. Results are shown in
the figure, with AFQMC/S for smaller $\lambda$ and AFQMC/TP($11$) for
$\lambda\geq 0.5$. Near-exact energies are obtained for all parameters
examined here. We observed that results can become severely biased with
AFQMC/S for large $\lambda$, as a consequence of a poor importance function
causing large, or even diverging, variances. Even with an improved importance
function TP($11$), small residual effects can be present (via underestimation
of the statistical error, or bias from population size). We also note a large
statistical error at $\omega=2t$ which is maximized at an intermediate value
of $\lambda=0.5$ (or $g=\sqrt{2}t$). Nevertheless, with TP(11), the bias (if
any) is smaller than 0.001 $t$ for the Holstein dimer, which highlights the
accuracy and sampling efficiency of AFQMC/TP(11). We discuss the issue of bias
in AFQMC in sign problem-free models further in Section VI.2 and Section VI.5.

###  VI.2 1D 4-site model at half-filling

To further investigate the effect of the importance function on the sampling
result and any potential bias, we consider a 1D 4-site Holstein model
employing periodic boundary conditions (PBCs), at half-filling with
$\lambda=0.5$, $\omega=4t$, and $g=2t$. We compute the ground state energies
with AFQMC using the following trial wave functions: a single semiclassical
state (S); the TP wavefunction with 13 semiclassical states (TP(13)); the LLF
wavefunction; and the TP-LLF wavefunction with two LLF states (TP-LLF(2)).
With AFQMC/TP-LLF(2), a ground-state energy of $-10.293(2)$ is obtained,
compared to the exact result of $-10.292$ (obtained from DMRG using
iTensorFishman  _et al._ (2020), although ED can also be done here). In
contrast, a biased result is seen with each of the other forms of the trial
wave function. The bias is about $0.4$% relative to the exact result, using
the computational parameters specified in Sec. X, and is essentially
independent of whether bifurcation is accounted for or not in the trial wave
function. These results suggest that, to remove the sampling bias in this
parameter regime, it is critical to capture in the importance function both a
means to overcome the adiabatic potential bifurcation issue and treat el-ph
correlation.

###  VI.3 1D 20-site model at half-filling

Next, we discuss a 20-site 1D Holstein model at half-filling employing PBCs.
ED is no longer feasible for systems of this size so we used the
iTensorFishman  _et al._ (2020) package to perform DMRG calculationsJeckelmann
and White (1998). The DMRG calculations were carried out by placing
alternating fermionic and bosnic sites on a 1D lattice so that overall there
are twice the number of sites compared to the physical lattice problem. While
it is possible to use an optimized phonon basisZhang  _et al._ (1998) to
handle larger el-ph coupling cases, here we employed the most primitive
version of DMRG for simple comparisons. The bond dimension we used was fixed
at 1000 and the maximum number of bosons for each site was taken to be 60.

Figure 3:  Total energy per site in units of $t$ for the 20-site 20-electron
1D Holstein model as a function of $\omega$ for various values of $\lambda$
values: (a) $\lambda=0.1$ results, (b) $\lambda=0.3$ results, (c)
$\lambda=0.8$ results, and (d) $\lambda=2.0$ results. Note that all error bars
for AFQMC are too small to be seen on the plotted scales, DMRG reference
values are unavailable for $\omega<0.8$ when $\lambda=0.8$ and $\lambda=2.0$,
and AFQMC/TP(11) is not presented for $\lambda=0.1$ because the results for
this trial wavefunction are nearly identical to those of AFQMC/S.

We compare the total energy per site within DMRG, CSMP2, variational LF,
AFQMC/S, and AFQMC/TP($11$) in Fig. 3 for various $\lambda$ and $\omega$.
Given the discussion of Section VI.2, it is desirable to employ the TP-
LLF($n$) wavefunctions in general, but we leave a more detailed study with
this trial wavefunction for a future study. Here, we focus on AFQMC with
simpler and less accurate trial wavefunctions (AFQMC/S and AFQMC/TP(11)).

Similarly to the two-site problem, we observe that CSMP2 follows the (near-
exact) answers given by DMRG and AFQMC closely at small $\lambda$ (e.g.
$\lambda=0.1$) as the frequency of the phonon mode is varied. However, we see
a clear quantitative deviation of CSMP2 from the other curves as $\omega$
increases. The deviation is again maximized at an intermediate coupling
$\lambda=0.8$ and is smaller at weak and strong couplings. The variational LF
wavefunction works better than CSMP2 for all $\lambda>0.1$. Its strength over
CSMP2 is highlighted as $\omega$ increases. This clearly suggests that the
variational LF wavefunction includes el-ph correlation beyond the second-order
contribution provided in CSMP2.

The performance of AFQMC/S is very good at all coupling strengths considered
here for $\omega\leq 2.0$. Similar issues with biased final estimates from
poor importance functions are seen at $\lambda=0.3$ and $\lambda=0.8$.
AFQMC/TP($11$) shows improvement over the simplest semiclassical importance
function in the case of $\lambda=0.3$. However its improvement for
$\lambda=0.8$ as $\omega$ becomes larger is very small. For example, the
residual bias is still visible at $\omega=2.0$. This points to the need to
improve the importance function over the forms we have used. We expect that
incorporating el-ph correlation directly, as in the LF-type wavefunctions,
will ameliorate this sampling bias greatly, as observed in Section VI.2. For
the rest of the paper, we focus on AFQMC with the simplest trial wavefunction,
namely a single semiclassical trial wavefunction, because the observed
sampling bias is small enough that it does not affect the conclusions of this
work.

###  VI.4 2D 4x4 model at half-filling

We have established the expected behavior of AFQMC with semicalssical
wavefunctions as importance functions from studying one-dimensional problems
such as the Holstein dimer and a 1D chain. Here, we explore higher dimensions
by investigating a 2D square lattice problem with a 4x4 geometry. We employed
PBCs along x-axis and OBCs along y-axis. The main reason for choosing this
boundary condition is to ease the convergence of the DMRG calculations. We
were able to converge DMRG calculations only for $0.8\geq\omega$ and
$\lambda\leq 0.5$ where we used a bond dimension of 2500 and a maximum number
of bosons of 25.

In Fig. 4, the energy per site as a function of $\omega$ for various $\lambda$
values is presented for this 2D model. We observe conclusions similar to our
previous one-dimensional examples. CSMP2 quantitatively fails as $\omega$
increases. Furthermore, for a fixed $\omega$, CSMP2 performs worst for
intermediate $\lambda$ values and is more accurate for small and large
$\lambda$ values. Similarly to the 1D 20-site case, the variational LF energy
is more accurate than CSMP2 for $\lambda>0.1$ and its improvement over CSMP2
becomes larger as $\omega$ increases. AFQMC/S is well behaved in the range
$\omega\in[0.1,2.0]$. Its maximum error occurs at $\lambda=0.3$ and
$\omega=2.0$, where clear indications of sampling bias arises. Nevertheless,
the range of parameters where AFQMC/S can be reliably performed with the
simplest possible semiclassical trial wavefunction is quite broad even in 2D,
highlighting the utility and potential of this approach.

Figure 4:  Total energy per site in units of $t$ for the 4x4 2D Holstein model
at half-filling as a function of $\omega$ for various $\lambda$ values: (a)
$\lambda=0.1$ results, (b) $\lambda=0.3$ results, (c) $\lambda=0.5$ results,
and (d) $\lambda=2.0$ results. Note that all error bars of CP-AFQMC/S are too
small to be seen on the plotted scales, DMRG reference values are only
available for $\omega\geq 0.8$ when $\lambda<2.0$.

###  VI.5 Autocorrelation time and variance control

As we mentioned at the beginning of this section, in the Holstein model the
difference between our AFQMC approach and DQMC is mainly in the details of the
Monte Carlo sampling algorithm. The two methods can thus have different
behavior in terms of efficiency in different regimes of the parameter space.
Here we look into this to help understand the domain of applicability. We note
that this is not the focus of our study, since in the most general case where
electron interactions are present, the branching random walk approach must be
adopted in order to control the sign or phase problem.

The standard DQMC algorithm based on local updates exhibits a long
autocorrelation time in the Holstein model for $\omega<0.5$ and for low
temperatures. It has been found that this is a consequence of an ergodicity
problem. A careful mathematical analysis of the causes of this problem can be
found in the work of Hohenadler and co-workers.Hohenadler  _et al._ (2004,
2005); Hohenadler and von der Linden (2007) These authors have shown that the
condition number of the bosonic action sampled in DQMC for small values of
$\Delta\tau$ scales as $1/(\omega\Delta\tau)^{2}$. This poorly conditioned
action leads to a long autocorrelation time that scales quadratically with
increasing $\omega^{-1}$. There have been attempts to ameliorate this problem
based on global moves such as the Langevin dynamics approachKarakuzu  _et al._
(2018); Batrouni and Scalettar (2019); Hébert  _et al._ (2019) and the self-
learning Monte Carlo approach.Chen  _et al._ (2018) We also mention that the
work of Hohenadler and co-workers removed the autocorrelation problem using
the Lang-Firsov transformation along with a principal component
analysis.Hohenadler  _et al._ (2004)

Figure 5:  Log-log plot of autocorrelation time $\tau_{\text{ac}}/\Delta\tau$
and phonon frequency $\omega/t$ for $\lambda=0.1$ and $\lambda=1.0$. The black
dotted lines are linear fits for each curve. For $\lambda=0.1$, the slope is
-0.9876 with $R^{2}=0.9993$ and for $\lambda=1.0$, the slope is -0.8617 with
$R^{2}=0.9942$.

Since AFQMC applies a projector to the entire set of electronic and phonon
degrees of freedom using a population of random walkers, it is less prone to
ergodicity problems. The Monte Carlo time coincides with the imaginary-time
direction, with open-ended random walks evolving along the worldines, which
makes them less likely to become trapped in particular configurations of the
phonon paths. There is a deep connection between this fact and the necessity
to resort to this sampling approach in order to impose a CP or phaseless gauge
condition Zhang  _et al._ (1997); Zhang and Krakauer (2003). To quantify this,
we directly compute the energy autocorrelation function,

$c_{E}(\tau)=\frac{1}{N-\tau}\sum_{n=1}^{N-\tau}(E_{n}-\langle
E\rangle)(E_{n+\tau}-\langle E\rangle),$ (91)

which gives an estimate of the integrated autocorrelation time via

$\tau_{\text{ac}}=\sum_{\tau=-\infty}^{\infty}\frac{c_{E}(\tau)}{c_{E}(0)}.$
(92)

The summation in Eq. 92 needs to be performed within some window instead of
over the entire set of samples since, for $\tau\gg\tau_{\text{ac}}$, the
summed noise becomes comparable to the actual signal. We follow Sokal’s
prescription of the automated windowing procedure to handle this issue.mac ;
Sokal (1997) In Fig. 5, we present estimates of autocorrelation times for
$\lambda=0.1$ and $\lambda=1.0$ for the 4x4 2D Holstein model. For both values
of $\lambda$, we observe a near linear behavior in the log-log scale
correlation between $\tau_{\text{ac}}$ and $\omega$ as in Fig. 5. Empirically,
we find that the autocorrelation time scales as $1/\omega^{0.9876}$ for
$\lambda=0.1$ and $1/\omega^{0.8617}$ for $\lambda=1.0$. This scaling is a
significant improvement over that of the standard DQMC algorithms where
$\tau_{\text{ac}}$ scales as $1/\omega^{2}$.Hohenadler  _et al._ (2004)

On the other hand, in AFQMC we use an importance function to guide the random
walks. If the quality of the importance function is very poor, the variance
can grow and even become infinite.Shi and Zhang (2016) In cases where the
trial wave function suppresses certain regions of the Hilbert space being
sampled with a qualitatively incorrect functional form, the autocorrelation
time and thus the variance can diverge, as mentioned earlier in this section.
This situation was seen in the examples with the semiclassical wave function
where there is a strong bifurcation of the adiabatic potential in the Holstein
dimer. Another example occurs with the use of semiclasscial trial wave
functions where the lack of explicit el-ph correlation leads to large sampling
biases. In extreme cases, calculations will be seemingly well-behaved in
“normal-sized” runs, as the Monte Carlo sampling is strongly biased by the
wrong importance function and the auto-correlation time grows exponentially.
These situations require careful analysis of the variance and study of the
dependence on the details of the importance function to reveal the problem.Shi
and Zhang (2016) Separate but related to the quality of the importance
function is the issue of efficiently sampling of multi-modal landscapes in the
el-ph models, as we have only incorporated local moves in our random walks. In
the worst cases, AFQMC can, even with reasonable choices of importance
functions, experience difficulties with long autocorrelations as occurs in
DQMC. In AFQMC the use of a population of open-ended random walkers with
branching can help avoid the sampling being stuck.

##  VII The Hubbard-Holstein Model

The focus of our CP-AFQMC method is on doped systems and more realistic
Hamiltonians, where the sign problem or phase problem will be present. In the
previous section, we studied the 1D and 2D Holstein models with simple trial
wavefunctions. In this section, we present benchmark data using the same
approach for the Hubbard-Holstein model with $U/t=4$. Because of the
competition between $U$ and $g$, we carefully study trial wavefunctions with
both CDW and SDW order. Note that the CP-AFQMC algorithm is no longer exact
because the el-el repulsion will lead to a sign problem. Karakuzu, Seki, and
Sorella have presented an efficient QMC algorithm which is free of the sign-
problem at half-filling as long as $U>2g^{2}/\omega$.Karakuzu  _et al._ (2018)
Within CP-AFQMC, we can simulate any parameter regime without the sign problem
at the expense of the constraint bias. Since CP-AFQMC has been extensively
benchmarked for electronic systems in the past,Qin  _et al._ (2016a) we focus
on any additional biases that may arise from the interplay between electrons
and phonons in this section. We use an ad hoc $U_{\text{eff}}=0.5\,t$ in the
electronic mean-field part to generate all the SDW trial wavefunctions for CP-
AFQMC in this section.

Below, we first examine the behavior of our algorithm in 1D, and then in 2D
both at half-filling and $1/8$-doping. We will focus on benchmarking the
accuracy of the computed ground-state energies. We note that for the purely
electronic cases with $\lambda=0$, all of our models in 1D and at half-filling
in 2D are sign-problem-free. However, CP-AFQMC can incur a systematic error in
the energy in these cases, because of an “artificial node” in auxiliary-field
space; Zhang  _et al._ (1997) this error can be removed straightforwardly. Shi
and Zhang (2013); Qin  _et al._ (2016a) Instead of invoking the scheme to
remove this artificial bias, we will perform the CP-AFQMC calculation in the
generic way as described above, since our main focus in this work is the most
general situation of a doped Hubbard-Holstein model where the sign problem is
present.

###  VII.1 1D 20-site model at half-filling

Figure 6:  The total energy per site in units of $t$ for the 20-site 1D
Hubbard-Holstein model at half-filling with $U=4t$ as a function of $\lambda$
for various $\omega$ values: (a) $\omega=0.1$, (b) $\omega=0.4$, (c)
$\omega=1.0$, and (d) $\omega=2.0$. Note that all error bars of CP-AFQMC are
too small to be seen on the plotted scales. Note that the DMRG results are
unavailable for $\omega=0.1$ and $\omega=0.4$ as well as for $\lambda>2$.

We benchmark CP-AFQMC against DMRG for the 20-site 1D Hubbard-Holstein model
at half-filling with PBCs. Unlike for the case of the pure Holstein model,
CSMP2 and the variational LF approach are quantitatively and qualitatively
inaccurate for all parameters examined here. This is not surprising because
the on-site repulsion term for $U/t=4$ is not small, so the failure of mean-
field theories and a low-order perturbation theory on the el-el interaction is
expected. For this reason we do not discuss CSMP2 and variational LF results
here.

To study the CDW and SDW phases and the possibility of a phase transition
between them, we carry out CP-AFQMC calculations using two different mean-
field trial wave functions with the corresponding broken symmetry. Comparison
of the computed energies indicates which one is the ground state at each
Hamiltonian parameter choice as well as the existence and location of a
transition, although the fact that our CP-AFQMC energies computed from the
mixed estimate are not variational Carlson  _et al._ (1999) adds a subtlety to
this procedure. Here the calculation leading to the higher energy can be
thought of as the constraint acting to “hold” the projection to an excited
state compatible with the broken symmetry of the trial wave function. In
actual applications, we could use a self-consistent CP-AFQMC procedure Qin
_et al._ (2016b) to tune the trial wave function and reduce its effect on the
result, but for the purpose of benchmarks we will only perform one-shot
calculations here using UHF trial wave functions generated with a fixed
$U_{\rm eff}$, and rely on comparison with DMRG results to gauge the accuracy.

In Fig. 6, we compare two sets of CP-AFQMC results at different $\omega$
values, one set with SDW trial wavefunctions (denoted by CP-AFQMC/SDW) and
another with CDW trial wavefunctions (denoted CP-AFQMC/CDW). The trial wave
functions themselves show a SDW to CDW transition at $\lambda\sim 0.5$, given
our ad hoc choice of $U_{\text{eff}}=0.5$. We see that CP-AFQMC/CDW leads to
higher energies than CP-AFQMC/SDW for $\lambda\leq 1.0$. The energy
differences are large enough to make it straightforward to identify the
correct phase.

Figure 7:  Error in total energy per site compared to DMRG in units of
$0.001t$ for the 20-site 1D Hubbard-Holstein model at half-filling with
$U/t=4$ as a function of $\omega$ for various $\lambda$ values: (a)
$\lambda=0.1$ results and (b) $\lambda=1.0$ results. The black dotted line in
(a) indicates the electronic CP-AFQMC energy error (i.e., $\lambda=0.0$). In
(b), DMRG energies are only available for $\omega=1.0,1.5,2.0$. Note that for
(b) the reference DMRG values are potentially unconverged with a maximum
phonon number of 40 and the remaining time step error of CP-AFQMC is not
negligible on the plotted energy scale.

We next make more quantitative comparison of the total energy per site
obtained from CP-AFQMC and DMRG in Fig. 6. DMRG calculations are performed
with a bond dimension of 1000 and with the maximum number of bosons of 40. In
Fig. 6 (c) and (d), we observe that CP-AFQMC/SDW closely follows the DMRG
energies from $\lambda=0.1$ to $\lambda=1.0$. At $\lambda=2.0$, the energy
obtained from CP-AFQMC/SDW is significantly higher than that from CP-
AFQMC/CDW, with the latter in good agreement with DMRG. The procedure
described above of combining the lowest energy curves between CP-AFQMC/SDW and
CP-AFQMC/CDW produces quantitatively accurate results across the full range of
parameters that we study. Variations in the value of the phonon frequency do
not change the qualitative conclusions.

In Fig. 7 we show a magnified view of the absolute discrepancies between CP-
AFQMC and DMRG energies. As a comparison, for the purely electronic Hubbard
model ($\lambda=0$), CP-AFQMC exhibits an error per site of 0.00245(8) $t$
with respect to the DMRG reference values. Similarly in the Hubbard-Holstein
model, CP-AFQMC/SDW energies exhibit an error per site of approximately
0.002–0.003 $t$ for $\lambda=0.1,0.3,0.8$. When the system reaches values as
large as $\lambda=1.0$, we observe a small increase in the error as $\omega$
increases. The largest error found for $\lambda=1.0$ is 0.0064(2) $t$ at
$\omega=2.0$, which is slightly larger than the constraint bias found in the
purely electronic problem. The point at which the largest error was observed
coincides with the expected phase transition point between SDW and CDW (see
below). At $\lambda=2.0$, the discrepancy between CP-AFQMC and DMRG is an
order of magnitude smaller, with a maximum deviation of CP-AFQMC per site of
-0.00040(2) $t$ at $\omega=2.0$. Possible reasons for this small discrepancy
are that the DMRG energies with a maximum boson number of 40 are not fully
converged and the residual time step error in CP-AFQMC is not negligible on
the plotted energy scale.

It has been shown by several methodsOhgoe and Imada (2014); Karakuzu  _et al._
(2018); Reinhard  _et al._ (2019); Costa  _et al._ (2020) that in the
thermodynamic limit the Hubbard-Holstein model undergoes a transition between
SDW and CDW at

$U\approx\frac{2g^{2}}{\omega}=4\text{d}t\lambda.$ (93)

This value of $U$ is where the effective on-site interaction changes sign as
shown in Eq. 84. For $U/t$ = 4 and $\text{d}=1$, we expect the phase
transition to occur at approximately $\lambda=1.0$. We find that in the
20-site model, despite the expected finite size effects, the onset of the
phase transition is captured quite well. In particular, the crossover between
CP-AFQMC/SDW and CP-AFQMC/CDW occurs roughly at $\lambda=1.0$ in Fig. 6. While
this is encouraging, detailed phase diagram studies with CP-AFQMC should be
carried out in the future. We note that CP-AFQMC often restores the symmetry
breaking of the underlying mean-field trial wave function, Purwanto  _et al._
(2008); Lee  _et al._ (2020) as would be expected of an exact many-body
computation. Therefore, a proper phase diagram study with CP-AFQMC should
involve a direct measurement of correlation functions Chang and Zhang (2010)
or order parameters with explicit symmetry-breaking induced.Qin  _et al._
(2020) Furthermore, there may be intermediate phases such as metallic or
superconducting phases near the onset of the phase transition between the SDW
and CDW phases. Studying these putative intermediate phases is of great
interest.Costa  _et al._ (2020)

###  VII.2 2D 4x4 model at half-filling and 1/8 hole-doping

Instead of comparing CP-AFQMC with other methods for 2D Hubbard-Holstein
systems, we simply report the computed total energy per site using PBCs along
both the $x$ and $y$ directions, as shown in Fig. 8. Based on the benchmark
studies in the previous sections and on experience from the purely electronic
model, we expect that our results will be of similar accuracy (or better
because of effective reduction of the el-el interaction from the el-ph
coupling) to that in the Hubbard model for most parameters considered in this
work. In Fig. 8, we see that CP-AFQMC/CDW has a lower energy for $\lambda>0.5$
at both frequencies ($\omega=0.1$ and $\omega=2.0$) at both half-filling and
1/8 hole-doping. Based on Eq. 93, it is expected that the onset of the
crossover occurs around $\lambda=0.5$ in 2D, consistent with our numerical
results. (Note that our definition of $\lambda$ includes dimensionality, hence
the change in the crossover value from 1D to 2D).

Consistent with our previous results, for $\lambda=0.1$ and $\lambda=0.3$ at
all frequencies up to $\omega=2.0$, the CP-AFQMC error per site is
approximately 0.009 $t$ or slightly larger. At the onset of the crossover
between CP-AFQMC/CDW and CP-AFQMC/SDW ($\lambda\sim 0.5$), we expect the error
to be maximized and larger than that of the CP-AFQMC bias for the electronic
problem, similar to the 1D 20-site model at half-filling. For $\lambda=0.8$
and $\lambda=2.0$, we expect that our results will be nearly exact for the
$\omega$ values studied. For the purely electronic Hubbard model
($\lambda=0$), CP-AFQMC exhibits an error per site of about 0.00901(9) $t$ at
half-filling and 0.00469(4) $t$ at 1/8 hole-doping, using a UHF trial wave
function. As mentioned, the error at half-filling is “artificial” and can be
removed, Shi and Zhang (2013); Qin  _et al._ (2016a) but this is not done
here. For $\lambda\leq 1$ at half-filling we expect an error of comparable
size. Comparing two different fillings, we do not see qualitative differences
in physical behavior in our finite-sized lattice, and the value of the phonon
frequency does not appear to make qualitative differences as well. We note
that the energy difference between CP-AFQMC/SDW and CP-AFQMC/CDW noticeably
shrinks as the phonon frequency $\omega$ increases.

Figure 8:  Total energy per site in units of $t$ for the 4x4 2D Hubbard-
Holstein model with $U/t=4$ as a function of $\omega$: (a) $\omega=0.1$
results at half-filling, (b) $\omega=2.0$ results at half-filling, (c)
$\omega=0.1$ results at 1/8 hole-doping, and (d) $\omega=2.0$ results at 1/8
hole-doping. Note that all error bars of CP-AFQMC are too small to be seen on
the plotted scales.

##  VIII Towards ab-initio Hamiltonians

We briefly discuss the extension of the presented algorithm for general ab-
initio Hamiltonians. The ab-initio Hamiltonian that describes el-ph problems
typically involves linear el-ph coupling. Therefore, the most widely used ab-
initio Hamiltonian has the same form as Eq. 1 with more general Hamiltonian
matrix elements,

$\displaystyle\hat{\mathcal{H}}_{\text{el}}^{\text{(1)}}$
$\displaystyle=\sum_{\sigma\in\\{\uparrow,\downarrow\\}}\sum_{pq}h_{p_{\sigma}q_{\sigma}}\hat{a}_{p_{\sigma}}^{\dagger}\hat{a}_{q_{\sigma}},$
(94) $\displaystyle\hat{\mathcal{H}}_{\text{el}}^{\text{(2)}}$
$\displaystyle=\frac{1}{2}\sum_{\sigma,\sigma^{\prime}\in\\{\uparrow,\downarrow\\}}\sum_{pqrs}\left(p_{\sigma}r_{\sigma}|q_{\sigma^{\prime}}s_{\sigma^{\prime}}\right)\hat{a}_{p_{\sigma}}^{\dagger}\hat{a}_{q_{\sigma^{\prime}}}^{\dagger}\hat{a}_{s_{\sigma^{\prime}}}\hat{a}_{r_{\sigma}},$
(95) $\displaystyle\hat{\mathcal{H}}_{\text{ph}}$
$\displaystyle=\sum_{I}\omega_{I}\hat{b}_{I}^{\dagger}\hat{b}_{I},$ (96)

and

$\displaystyle\hat{\mathcal{H}}_{\text{el-ph}}$
$\displaystyle=\sum_{\sigma\in\\{\uparrow,\downarrow\\}}\sum_{pqI}g_{p_{\sigma}q_{\sigma}I}\hat{a}_{p_{\sigma}}^{\dagger}\hat{a}_{q_{\sigma}}(\hat{b}_{I}^{\dagger}+\hat{b}_{I}),$
(98)

where we have suppressed other quantum numbers such as $\mathbf{k}$-point
dependencies and have expressed everything in terms of the electronic
($\\{p,q,r,s,\cdot\cdot\cdot\\}$) and phononic
bands($\\{I,J,K,L,\cdot\cdot\cdot\\}$). The computation of these matrix
elements at the level of density functional theory has been well-documented
Giustino (2017); Poncé  _et al._ (2016); Zhou  _et al._ (2020) so here we
focus on briefly describing the phaseless AFQMC (ph-AFQMC) algorithmZhang and
Krakauer (2003) for these realistic el-ph problems.

The walkers take the same form as in Eq. 20. Therefore, the essence of the
propagation algorithm remains unchanged. The only complication arises from the
generalized form of $\hat{\mathcal{H}}_{\text{el}}^{\text{(2)}}$ which
necessitates the use of a continuous Hubbard-Stratonovich
transformation.Stratonovich (1957); Hubbard (1959) The continuous
transformation leads to the fermionic phase problem which can be removed via
the phaseless constraint. Zhang and Krakauer (2003) The propagation is carried
out the same way with appropriate modifications to the constraint to account
for the phase problem. The ab-initio generalization of semiclassical states
used in this work is also straightforward. The trial wavefunction still takes
the form of Eq. 35, Eq. 36, and Eq. 38. A variational minimization of the
total energy then leads to a trial wavefunction that can be used in ph-AFQMC.
The projection of the trial wavefunction onto phonon displacements
$|\mathbf{X}\rangle$ is identical to Eq. 40 except that the phonon mass and
frequency now depend on band indices $\\{I\\}$. The ab-initio generalization
of the LF wavefunction may also be carried out straightforwardly by extending
the LF generator in Eq. 50. We expect that the ab-initio ph-AFQMC approach
will become a valuable tool for understanding polaronic physics in realistic
correlated materials in the future.

##  IX Conclusions and outlook

In this work, we have introduced an extension of CP-AFQMC to describe
correlated systems with el-ph coupling. Our approach utilizes a mixed
first/second-quantized representation where the phonons are propagated in
first quantization following the commonly used diffusion MC algorithm, and the
electronic degrees of freedom are handled in second quantization via AFQMC.
The resulting algorithm is compared with numerically-exact DMRG and low-order
perturbation theories for the Holstein model as a first test of the basic
algorithm. We have demonstrated that the autocorrelation time problems that
arises in the commonly used DQMC methods is greatly ameliorated in AFQMC, with
autocorrelation time that scales roughly as $1/\omega$.

While the Holstein model is sign-problem free, AFQMC with the simplest trial
wavefunctions, namely semiclassical states, is found to introduce a small bias
when the underlying adiabatic surface bifurcates and/or $g$ is larger than $t$
but smaller than $\omega$ (e.g., $\omega\rightarrow\infty$ for a fixed
$\lambda$). Based on a 4-site model, we have shown that this bias can be
removed by using an improved trial wavefunction where both bifurcations and
increased el-ph correlation are accounted for. We have demonstrated the
remarkable accuracy of AFQMC for both 1D and 2D Holstein models over a
reasonably broad set of coupling and phonon frequency parameters via direct
comparison with DMRG.

We have tested CP-AFQMC on the finite sized versions of the 1D and 2D Hubbard-
Holstein models with $U/t=4$, using the simplest form of trial wave functions
consisting of a semiclassical state with a single Slater determinant. For the
1D Hubbard-Holstein model, we have compared CP-AFQMC against numerically exact
DMRG results. When $\lambda$ is small and the ground state is dominated by the
Hubbard $U$ term, we find that the error of our algorithm is roughly the same
as that expected from standard CP-AFQMC applied to the purely electronic
Hubbard model. Furthermore, when the ground state is dominated by the el-ph
coupling term and exhibits charge density wave order, we find that the overall
error becomes remarkably small (smaller than that expected in purely
electronic systems). These facts have motivated the production of what we
believe are benchmark results for the finite sized 2D Hubbard-Holstein model
for various values of $\lambda$ and $\omega$ at half-filling and 1/8 hole-
doping.

For $U<2g^{2}/\omega$ at half-filling and for all parameter regimes at any
hole-doping, standard QMC approaches suffer from the sign problem.Karakuzu
_et al._ (2018) Therefore, our AFQMC approach should become an essential tool
for producing accurate results scalable to large system sizes for this model.
We have investigated the competing spin- and charge-density wave orders in the
Hubbard-Holstein model. At the onset of the phase transition between these
phases, we observe a crossover in the energies between two AFQMC calculations
targeting the two phases. Lastly, we have briefly discussed the extension of
the presented algorithm to ab-initio Hamiltonians that can be easily
formulated based on the phaseless AFQMC method for general electronic
Hamiltonians. Zhang and Krakauer (2003)

Some immediate future directions include using this framework to provide a
detailed study of the phase diagram of the Hubbard-Holstein model and other
lattice models, and extending this framework to finite-temperature problems
based on the constrained path approximation.Zhang (1999) As mentioned, a trial
wave function with an el-ph Jastrow factor can be implemented
straightforwardly in AFQMC, which is expected to further reduce the bias and
improve sampling efficiency in anti-adiabatic regimes with
$2\text{d}t\lambda<\omega$. It will also be valuable to further investigate
the implementation of the full LF trial wave function. Furthermore,
application of the proposed AFQMC approach to ab-initio systems will be of
great interest as well. While there are several algorithmic aspects that can
be further improved, including improved forms of importance functions and
better sampling in large phonon frequency regimes, we believe that the
algorithms and insights presented in this work will serve as stepping stones
towards simulating model as well as ab-initio systems with a non-trivial
interplay between electronic correlation and el-ph couplings, which continue
to be of great importance in modern condensed matter physics.

##  X Computational Details

Our algorithm was implemented in a public open-source auxiliary-field quantum
Monte Carlo package called PAUXY.pau  The blocking analysis was performed with
pyblock.pyb  The pair branching algorithm was used for population
control.Wagner  _et al._ (2009) Variational calculations were aided by
automatic differentiation using JAX.Bradbury  _et al._ (2018) A total of 640
walkers and a time step of 0.005 $t^{-1}$ were used in all calculations except
for the Holstein dimer and the 4-site 1D Holstein model. For the Holstein
dimer, we used a time step of 0.0005 $t^{-1}$ for $\omega\leq 1.6$ and 0.00025
$t^{-1}$ for $\omega>1.6$ with 6400 walkers. For the 4-site 1D Holstein model,
we used 6400 walkers and a time step of 0.0005 $t^{-1}$. The population
control bias and time step error were found to be smaller than 0.001 $t$ in
the absolute total energy per site.

##  XI Acknowledgement

We are grateful to Zihang Li and Hao Shi for their contributions during early
stages of this work. We thank Hao Shi, Michael Lindsey, Fionn Malone for
fruitful discussions and Miles Stoudenmire for help with iTensor calculations.
DRR acknowledges support from grant NSF-CHE 1954791. This work was conducted
using computational resources and services at the Flatiron Institute. The
Flatiron Institute is a division of the Simons Foundation.

## References

  * Giustino (2017) F. Giustino, Rev. Mod. Phys. 89, 015003 (2017).
  * Mahan (2013) G. D. Mahan,  _Many-particle physics_ (Springer Science & Business Media, 2013).
  * Cava  _et al._ (1987) R. Cava, R. Van Dover, B. Batlogg, and E. Rietman, Phys. Rev. Lett. 58, 408 (1987).
  * Wu  _et al._ (1987) M.-K. Wu, J. R. Ashburn, C. Torng, P. H. Hor, R. L. Meng, L. Gao, Z. J. Huang, Y. Wang, and a. Chu, Phys. Rev. Lett. 58, 908 (1987).
  * Scalapino (2012) D. J. Scalapino, Rev. Mod. Phys. 84, 1383 (2012).
  * Müller (1990) K. A. Müller, Z. Phys. B. Con. Mat. 80, 193 (1990).
  * Song and Annett (1995) J. Song and J. F. Annett, Phys. Rev. B 51, 3840 (1995).
  * Alexandrov (2000) A. S. Alexandrov, J. Supercond. Nov. Magn. 13, 985 (2000).
  * Cuk  _et al._ (2005) T. Cuk, D. H. Lu, X. J. Zhou, Z. X. Shen, T. P. Devereaux, and N. Nagaosa, Phys. Status Solidi B Basic Res. 242, 11 (2005).
  * Khatami  _et al._ (2008) E. Khatami, A. MacRidin, and M. Jarrell, Phys. Rev. B 78, 1 (2008).
  * Gunnarsson and Rösch (2008) O. Gunnarsson and O. Rösch, J. Condens. Matter Phys. 20 (2008).
  * Van Heumen  _et al._ (2009) E. Van Heumen, E. Muhlethaler, A. B. Kuzmenko, H. Eisaki, W. Meevasana, M. Greven, and D. Van Der Marel, Phys. Rev. B 79, 1 (2009).
  * Gadermaier  _et al._ (2010) C. Gadermaier, A. S. Alexandrov, V. V. Kabanov, P. Kusar, T. Mertelj, X. Yao, C. Manzoni, D. Brida, G. Cerullo, and D. Mihailovic, Phys. Rev. Lett. 105, 257001 (2010).
  * He  _et al._ (2018) Y. He, M. Hashimoto, D. Song, S. D. Chen, J. He, I. M. Vishik, B. Moritz, D. H. Lee, N. Nagaosa, J. Zaanen, T. P. Devereaux, Y. Yoshida, H. Eisaki, D. H. Lu, and Z. X. Shen, Science 362, 62 (2018).
  * Yang  _et al._ (2019) S. L. Yang, J. A. Sobota, Y. He, D. Leuenberger, H. Soifer, H. Eisaki, P. S. Kirchmann, and Z. X. Shen, Phys. Rev. Lett. 122, 1 (2019).
  * Hu  _et al._ (2019) Y. Hu, X. Chen, S. T. Peng, C. Lane, M. Matzelle, Z. L. Sun, M. Hashimoto, D. H. Lu, E. F. Schwier, M. Arita, T. Wu, R. S. Markiewicz, K. Shimada, X. H. Chen, Z. X. Shen, A. Bansil, S. D. Wilson, and J. F. He, Phys. Rev. Lett. 123, 1 (2019).
  * Rosenstein and Shapiro (2019) B. Rosenstein and B. Y. Shapiro, Phys. Rev. B 100 (2019).
  * Grissonnanche  _et al._ (2020) G. Grissonnanche, S. Thériault, A. Gourgout, M. E. Boulanger, E. Lefrançois, A. Ataei, F. Laliberté, M. Dion, J. S. Zhou, S. Pyon, T. Takayama, H. Takagi, N. Doiron-Leyraud, and L. Taillefer, Nat. Phys. (2020).
  * Rosenstein and Shapiro (2020) B. Rosenstein and B. Y. Shapiro, (2020), arXiv:2003.03521 .
  * Banerjee  _et al._ (2020) S. Banerjee, W. A. Atkinson, and A. P. Kampf, , 1 (2020), arXiv:2008.01401 .
  * Shneyder  _et al._ (2020) E. I. Shneyder, E. I. Shneyder, S. V. Nikolaev, S. V. Nikolaev, M. V. Zotova, M. V. Zotova, R. A. Kaldin, S. G. Ovchinnikov, and S. G. Ovchinnikov, Phys. Rev. B 101, 235114 (2020).
  * Mishchenko  _et al._ (2020) A. S. Mishchenko, N. Nagaosa, and N. Prokof’ev, , 1 (2020), arXiv:2007.09888 .
  * Sreedhar  _et al._ (2020) S. A. Sreedhar, A. Rossi, J. Nayak, Z. Anderson, Y. Tang, B. Gregory, M. Hashimoto, D. H. Lu, E. Rotenberg, R. J. Birgeneau, M. Greven, M. Yi, and I. M. Vishik, (2020), arXiv:2005.07233 .
  * Peng  _et al._ (2020) Y. Peng, A. Husain, M. Mitrano, S.-L. Sun, T. Johnson, A. Zakrzewski, G. MacDougall, A. Barbour, I. Jarrige, V. Bisogni, and P. Abbamonte, Phys. Rev. Lett. 097002, 1 (2020).
  * Hubbard (1963) J. Hubbard, Proc. R. Soc. A 276, 238 (1963).
  * Anderson (1987) P. W. Anderson, Science 235, 1196 (1987).
  * Qin  _et al._ (2020) M. Qin, C.-M. Chung, H. Shi, E. Vitali, C. Hubig, U. Schollwöck, S. R. White, S. Zhang, _et al._ , Phys. Rev. X 10, 031016 (2020).
  * LeBlanc  _et al._ (2015) P. F. LeBlanc, A. E. Antipov, F. Becca, I. W. Bulik, G. K. L. Chan, C. M. Chung, Y. Deng, M. Ferrero, T. M. Henderson, C. A. Jiménez-Hoyos, E. Kozik, X. W. Liu, A. J. Millis, N. V. Prokof’ev, M. Qin, G. E. Scuseria, H. Shi, B. V. Svistunov, L. F. Tocchio, I. S. Tupitsyn, S. R. White, S. Zhang, B. X. Zheng, Z. Zhu, and E. Gull, Phys. Rev. X 5, 1 (2015).
  * Zheng  _et al._ (2017) B.-X. Zheng, C.-M. Chung, P. Corboz, G. Ehlers, M.-P. Qin, R. M. Noack, H. Shi, S. R. White, S. Zhang, and G. K.-L. Chan, Science 358, 1155 (2017).
  * Jeckelmann and White (1998) E. Jeckelmann and S. R. White, Phys. Rev. B 57, 6376 (1998).
  * Jeckelmann  _et al._ (1999) E. Jeckelmann, C. Zhang, and S. R. White, Phys. Rev. B 60, 7950 (1999).
  * Tezuka  _et al._ (2005) M. Tezuka, R. Arita, and H. Aoki, Physica B: Condensed Matter 359, 708 (2005).
  * Tezuka  _et al._ (2007) M. Tezuka, R. Arita, and H. Aoki, Phys. Rev. B 76, 155114 (2007).
  * Fehske  _et al._ (2008) H. Fehske, G. Hager, and E. Jeckelmann, EPL 84, 57001 (2008).
  * Ejima and Fehske (2010) S. Ejima and H. Fehske, J. Phys.: Conf. Ser 200, 012031 (2010).
  * Bonča  _et al._ (1999) J. Bonča, S. Trugman, and I. Batistić, Phys. Rev. B 60, 1633 (1999).
  * Ohgoe and Imada (2014) T. Ohgoe and M. Imada, Phys. Rev. B 89, 195139 (2014).
  * Ohgoe and Imada (2017) T. Ohgoe and M. Imada, Phys. Rev. Lett. 119, 197001 (2017).
  * Karakuzu  _et al._ (2017) S. Karakuzu, L. F. Tocchio, S. Sorella, and F. Becca, Phys. Rev. B 96, 205145 (2017).
  * Jeon  _et al._ (2004) G. S. Jeon, T. H. Park, J. H. Han, H. C. Lee, and H. Y. Choi, Phys. Rev. B 70, 1 (2004).
  * Paci  _et al._ (2006) P. Paci, M. Capone, E. Cappelluti, S. Ciuchi, and C. Grimaldi, Phys. Rev. B 74, 1 (2006).
  * Werner and Millis (2007) P. Werner and A. J. Millis, Phys. Rev. Lett. 99, 1 (2007).
  * Murakami  _et al._ (2013) Y. Murakami, P. Werner, N. Tsuji, and H. Aoki, Phys. Rev. B 88, 1 (2013), arXiv:1305.5771 .
  * Li  _et al._ (2017) S. Li, E. Khatami, and S. Johnston, Phys. Rev. B 95, 1 (2017).
  * Sandhoefer and Chan (2016) B. Sandhoefer and G. K. L. Chan, Phys. Rev. B 94 (2016).
  * Reinhard  _et al._ (2019) T. E. Reinhard, U. Mordovina, C. Hubig, J. S. Kretchmer, U. Schollwöck, H. Appel, M. A. Sentef, and A. Rubio, J. Chem. Theory Comput. 15, 2221 (2019).
  * Sibaev  _et al._ (2020) M. Sibaev, I. Polyak, F. R. Manby, and P. J. Knowles, J. Chem. Phys. 153, 124102 (2020).
  * Dresselhaus  _et al._ (2020) T. Dresselhaus, C. Bungey, P. J. Knowles, and F. R. Manby, arXiv preprint arXiv:2010.04654 (2020).
  * White  _et al._ (2020) A. F. White, Y. Gao, A. J. Minnich, and G. K. Chan, arXiv preprint arXiv:2009.13568 (2020).
  * Zhang  _et al._ (1995) S. Zhang, J. Carlson, and J. E. Gubernatis, Phys. Rev. Lett. 74, 3652 (1995).
  * Zhang  _et al._ (1997) S. Zhang, J. Carlson, and J. E. Gubernatis, Phys. Rev. B 55, 7464 (1997).
  * Blankenbecler  _et al._ (1981) R. Blankenbecler, D. J. Scalapino, and R. L. Sugar, Phys. Rev. D 24, 2278 (1981).
  * Scalapino and Sugar (1981a) D. J. Scalapino and R. L. Sugar, Phys. Rev. Lett. 46, 519 (1981a).
  * Scalapino and Sugar (1981b) D. J. Scalapino and R. L. Sugar, Phys. Rev. B 24, 4295 (1981b).
  * Johnston  _et al._ (2013) S. Johnston, E. A. Nowadnick, Y. F. Kung, B. Moritz, R. T. Scalettar, and T. P. Devereaux, Phys. Rev. B 87, 235133 (2013).
  * Mendl  _et al._ (2017) C. B. Mendl, E. A. Nowadnick, E. W. Huang, S. Johnston, B. Moritz, and T. P. Devereaux, Phys. Rev. B 96, 1 (2017).
  * Karakuzu  _et al._ (2018) S. Karakuzu, K. Seki, and S. Sorella, Phys. Rev. B 98, 1 (2018).
  * Costa  _et al._ (2020) N. C. Costa, K. Seki, S. Yunoki, and S. Sorella, Commun. Phys. 3, 1 (2020).
  * Hirsch (1983) J. E. Hirsch, Phys. Rev. B 28, 4059 (1983).
  * Zhang and Krakauer (2003) S. Zhang and H. Krakauer, Phys. Rev. Lett. 90, 136401 (2003).
  * Al-Saidi  _et al._ (2006) W. Al-Saidi, S. Zhang, and H. Krakauer, J. Chem. Phys. 124, 224101 (2006).
  * Carlson  _et al._ (1999) J. Carlson, J. E. Gubernatis, G. Ortiz, and S. Zhang, Phys. Rev. B 59, 12788 (1999).
  * McKenzie  _et al._ (1996) R. H. McKenzie, C. Hamer, and D. Murray, Phys. Rev. B 53, 9676 (1996).
  * Holstein (1959a) T. Holstein, Ann. Phys. 8, 325 (1959a).
  * Rubenstein  _et al._ (2012) B. M. Rubenstein, S. Zhang, and D. R. Reichman, Phys. Rev. A 86, 053606 (2012).
  * Blankenbecler and Sugar (1983) R. Blankenbecler and R. Sugar, Phys. Rev. D 27, 1304 (1983).
  * Hohenadler  _et al._ (2004) M. Hohenadler, H. G. Evertz, and W. Von der Linden, Phys. Rev. B 69 (2004).
  * Thouless (1960) D. J. Thouless, Nucl. Phys. 21, 225 (1960).
  * Umrigar  _et al._ (1993) C. Umrigar, M. Nightingale, and K. Runge, J. Chem. Phys. 99, 2865 (1993).
  * Hammond  _et al._ (1994) B. L. Hammond, W. A. Lester, and P. J. Reynolds,  _Monte Carlo methods in ab initio quantum chemistry_ , Vol. 1 (World Scientific, 1994).
  * Moskowitz  _et al._ (1982) J. W. Moskowitz, K. Schmidt, M. A. Lee, and M. H. Kalos, J. Chem. Phys. 77, 349 (1982).
  * Reynolds  _et al._ (1982) P. J. Reynolds, D. M. Ceperley, B. J. Alder, and W. A. Lester Jr, J. Chem. Phys. 77, 5593 (1982).
  * Kalosakas  _et al._ (1998) G. Kalosakas, S. Aubry, and G. Tsironis, Phys. Rev. B 58, 3094 (1998).
  * Romero  _et al._ (1998) A. H. Romero, D. W. Brown, and K. Lindenberg, J. Chem. Phys. 109, 6540 (1998).
  * Shi and Zhang (2013) H. Shi and S. Zhang, Phys. Rev. B 88, 125132 (2013).
  * Qin  _et al._ (2016a) M. Qin, H. Shi, and S. Zhang, Phys. Rev. B 94, 085103 (2016a).
  * Proville and Aubry (2000) L. Proville and S. Aubry, The European Physical Journal B-Condensed Matter and Complex Systems 15, 405 (2000).
  * Holstein (1959b) T. Holstein, Ann. Phys. 8, 343 (1959b).
  * Lang and Firsov (1963) I. Lang and Y. A. Firsov, Sov. Phys. JETP 16, 1301 (1963).
  * Silbey and Harris (1984) R. Silbey and R. A. Harris, J. Chem. Phys. 80, 2615 (1984).
  * Harris and Silbey (1985) R. A. Harris and R. Silbey, J. Chem. Phys. 83, 1069 (1985).
  * Pouthier (2013) V. Pouthier, J. Chem. Phys. 138, 044108 (2013).
  * Shi and Zhang (2016) H. Shi and S. Zhang, Physical Review E 93, 033303 (2016).
  * Shi and Zhang (2017) H. Shi and S. Zhang, Phys. Rev. B 95, 045144 (2017).
  * Qin  _et al._ (2016b) M. Qin, H. Shi, and S. Zhang, Phys. Rev. B 94, 235119 (2016b).
  * Shavitt and Bartlett (2009) I. Shavitt and R. J. Bartlett,  _Many-body methods in chemistry and physics: MBPT and coupled-cluster theory_ (Cambridge university press, 2009).
  * Bonča  _et al._ (2000) J. Bonča, T. Katrasnik, and S. Trugman, Phys. Rev. Lett. 84, 3153 (2000).
  * Bonča and Trugman (2000) J. Bonča and S. Trugman, J. Supercond. 13, 999 (2000).
  * Batrouni and Scalettar (2019) G. Batrouni and R. T. Scalettar, Phys. Rev. B 99, 035114 (2019).
  * Fishman  _et al._ (2020) M. Fishman, S. R. White, and E. M. Stoudenmire, arXiv preprint arXiv:2007.14822 (2020).
  * Zhang  _et al._ (1998) C. Zhang, E. Jeckelmann, and S. R. White, Phys. Rev. Lett. 80, 2661 (1998).
  * Hohenadler  _et al._ (2005) M. Hohenadler, H. G. Evertz, and W. Von Der Linden, Phys. Status Solidi B Basic Res. 242, 1406 (2005).
  * Hohenadler and von der Linden (2007) M. Hohenadler and W. von der Linden, Springer Series in Materials Science 103, 463 (2007).
  * Hébert  _et al._ (2019) F. Hébert, B. Xiao, V. G. Rousseau, R. T. Scalettar, and G. G. Batrouni, Phys. Rev. B 99, 075108 (2019).
  * Chen  _et al._ (2018) C. Chen, X. Y. Xu, J. Liu, G. Batrouni, R. Scalettar, and Z. Y. Meng, Phys. Rev. B 98, 1 (2018).
  * (96) “Autocorrelation time estimation,” https://dfm.io/posts/autocorr/, accessed: 2020-09-30.
  * Sokal (1997) A. Sokal, “Monte carlo methods in statistical mechanics: Foundations and new algorithms,” in  _Functional Integration: Basics and Applications_, edited by C. DeWitt-Morette, P. Cartier, and A. Folacci (Springer US, Boston, MA, 1997) pp. 131–192.
  * Purwanto  _et al._ (2008) W. Purwanto, W. A. Al-Saidi, H. Krakauer, and S. Zhang, J. Chem. Phys. 128, 114309 (2008).
  * Lee  _et al._ (2020) J. Lee, F. D. Malone, and M. A. Morales, J. Chem. Theory Comput. 16, 3019 (2020).
  * Chang and Zhang (2010) C.-C. Chang and S. Zhang, Physical review letters 104, 116402 (2010).
  * Poncé  _et al._ (2016) S. Poncé, E. Margine, C. Verdi, and F. Giustino, Comput. Phys. Commun. 209, 116 (2016).
  * Zhou  _et al._ (2020) J.-J. Zhou, J. Park, I.-T. Lu, I. Maliyov, X. Tong, and M. Bernardi, “Perturbo: a software package for ab initio electron-phonon interactions, charge transport and ultrafast dynamics,” (2020), arXiv:2002.02045 .
  * Stratonovich (1957) R. L. Stratonovich, Sov. Phys. Dokl. 2, 416 (1957).
  * Hubbard (1959) J. Hubbard, Phys. Rev. Lett. 3, 77 (1959).
  * Zhang (1999) S. Zhang, Phys. Rev. Lett. 83, 2777 (1999).
  * (106) See https://github.com/pauxy-qmc/pauxy for details on how to obtain the source code.
  * (107) See https://github.com/jsspencer/pyblock for details on how to obtain the source code.
  * Wagner  _et al._ (2009) L. K. Wagner, M. Bajdich, and L. Mitas, J. Comput. Phys. 228, 3390 (2009).
  * Bradbury  _et al._ (2018) J. Bradbury, R. Frostig, P. Hawkins, M. J. Johnson, C. Leary, D. Maclaurin, and S. Wanderman-Milne, “JAX: composable transformations of Python+NumPy programs,”  (2018).