# Full configuration interaction quantum Monte Carlo for coupled
electron–boson systems and infinite spaces

Robert J. Anderson  Department of Physics, King’s College London, Strand,
London WC2R 2LS, United Kingdom  Charles J. C. Scott  Department of Physics,
King’s College London, Strand, London WC2R 2LS, United Kingdom  George H.
Booth<EMAIL_ADDRESS>Department of Physics, King’s College London,
Strand, London WC2R 2LS, United Kingdom

###### Abstract

We extend the scope of full configuration interaction quantum Monte Carlo
(FCIQMC) to be applied to coupled fermion-boson hamiltonians, alleviating the
a priori truncation in boson occupation which is necessary for many other wave
function based approaches to be tractable. Detailing the required algorithmic
changes for efficient excitation generation, we apply FCIQMC in two
contrasting settings. The first is a sign-problem-free Hubbard–Holstein model
of local electron-phonon interactions, where we show that with care to control
for population bias via importance sampling and/or reweighting, the method can
achieve unbiased energies extrapolated to the thermodynamic limit, without
suffering additional computational overheads from relaxing boson occupation
constraints. Secondly, we apply the method as a “solver” within a quantum
embedding scheme which maps electronic systems to local electron-boson
auxiliary models, with the bosons representing coupling to long-range
plasmonic-like fluctuations. We are able to sample these general electron-
boson hamiltonians with ease despite a formal sign problem, including a
faithful reconstruction of converged reduced density matrices of the system.

##  I Introduction

Quantum systems of interacting electrons explicitly coupled to bosonic degrees
of freedom are ubiquitous in nature. Most notably, electrons can couple to
nuclear vibrations, characterized by collective phononic quasiparticles which
obey bosonic statistics [1, 2]. The explicit effect of these lattice
vibrations on the electronic structure is neglected in Born-Oppenheimer
approximations, but strong coupling between the electrons and these phonons
can substantially renormalize the effective interactions in the system, and
lead to not just quantitative changing of electronic expectation values, but
also qualitative new physics emerging from the system [3, 4, 5]. As a famous
example, the coupling between electrons and phonons in this way can
renormalize interactions of certain solid state systems to such a degree that
it can provide the mechanism for the formation of Cooper pairs and
superconducting order [6].

Phonons are not the only bosonic species of physical relevance with which
electrons can interact. The emerging field of tailored electronic structure
via optical cavities exposes an exciting experimental approach to designing
interactions between electrons, and inducing new decoherence pathways for
electronic systems [7, 8, 9]. In this, the photonic field is quantized via
quantum electrodynamics, with these bosonic cavity modes explicitly coupled to
the electronic degrees of freedom in the dipole limit to form mixed
“polaritonic” quasiparticles [10, 11, 12, 13, 14]. These correlated systems
can be found to engage in novel reaction pathways and even room-temperature
superfluidity [15, 16].

Beyond coupling to these “external” bosonic fields, collective composite
quasiparticles of electrons themselves can manifest with emergent
(quasi-)bosonic statistics. These encompass plasmons, magnons and other
bosonic-like collective excitations, which can mix strongly with other
electronic degrees of freedom [17, 18, 19, 11, 20, 21]. Model Hamiltonians can
be constructed whereby the physics of these long-range excitations is included
via coupling of a (perhaps local) set of correlated electronic degrees of
freedom to explicit bosons describing these physical collective modes. This
“bath” of bosons can induce novel entanglement, decoherence and screening of
the electronic interactions due to the coupling to the electrons [22, 23, 24,
25]. In this way, the bosons can represent the physical process of a large
number of external electronic degrees of freedom and their effect on a chosen
subspace.

However, we are limited in available tools for simulation of these interacting
mixed-species quantum systems, which is exceptionally challenging when the
instantaneous correlations are fully included, beyond a density functional or
Green’s function framework [26, 27, 28, 29, 30, 31, 32, 33]. From a wave
function perspective, the boson number is generally not conserved in the
hamiltonian, leading to a formally infinite Hilbert space of bosonic
occupation for any single mode. To ensure wave function approaches are still
tractable, it is common to introduce a threshold on the occupancy of any one
of these modes, truncating the high-energy parts of the bosonic Hilbert space
to ensure a return to a finite-dimensional problem. This approach has been
used to develop exact diagonalization [34, 24], density matrix renormalization
group (DMRG) [35, 36], coupled-cluster [12, 37, 38] and other wave function
approaches adapted for coupled electron-boson systems [39]. However, ensuring
convergence with respect to these truncations, especially with strongly
coupled and/or low-energy boson modes with high average occupation, can still
rapidly lead to intractable calculations compared to their fully fermionic
counterparts.

To avoid this truncation, other approaches look to either stochastically
sample over all bosonic occupations [40, 41, 42], parameterize an ansatz over
the bosonic sector (e.g. as a coherent Gaussian state) [43, 44, 45], or
implicitly renormalize over all boson occupations. This latter approach can be
achieved in a wave function context within an alternative coupled-cluster
formalism, whereby more highly occupied bosonic modes are formed from the
product of low bosonic occupation implicitly through the exponential ansatz
[46, 47, 48, 37]. Monte Carlo methods on the other hand look to stochastically
sample from this space without truncation, but while they represent the de
facto approach for pure bosonic systems, are substantially limited in their
applicability to coupled fermion-boson systems by the manifestation of the
fermion sign problem in its various guises, which can limit their efficacy
[41].

In this work, we look to adapt the modern approach to full configuration
interaction quantum Monte Carlo (FCIQMC) to tackle mixed fermion-boson models.
This stochastic approach describes an ensemble of signed “walkers” which
reside on many-body configurations, representing a coarse-grained snapshot of
the underlying wave function of interest [49, 50, 51, 52, 53, 54, 55, 56, 57,
58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75]. By
evolving the walkers according to simple rules underpinned by a master
equation derived from the imaginary-time Schrödinger equation, expectation
values can be averaged over this evolution, providing systematically
improvable estimates of correlated observables. Explicit annihilation events
between walkers of different sign, as well as “initiator” approximations can
overcome the fermion sign problem in this context and enable long-time stable
dynamics [76, 77, 78]. This efficient sampling allows for system sizes (and
Hilbert spaces) to be considered well beyond the size that can be treated with
formally exact approaches, and is emerging as a complementary method to DMRG
in fermionic systems.

The FCIQMC method has been used with success previously both for a wide array
of fermionic models [79, 80, 81, 82], as well as applied to bosons [83, 84],
however has not to date been applied to mixed-species systems. In particular,
we believe that an application of FCIQMC to coupled electron-boson systems is
particularly appealing, since it is not expected to require an explicit
truncation of the bosonic mode occupation, allowing a sampling over all
relevant bosonic degrees of freedom. An intrinsic strength of the FCIQMC
method is to efficiently seek out relevant parts of the configuration space,
and to exploit any natural emergent sparsity in the wave function
representation present in low-energy states. This will be of particular
importance in the formally infinite Hilbert space of these challenging mixed-
species systems, and we expect to provide an important tool for their
efficient simulation.

In section II we detail the specific form of the hamiltonians used in this
work, while in section III we describe the key changes required in the FCIQMC
algorithm for their efficient simulation, in particular focusing on the
excitation generation for general electron-boson hamiltonians. This has been
implemented in our FCIQMC package, M7 [85], which extends all modern aspects
of FCIQMC functionality, including excited states and semi-stochastic
adaptations to these mixed-species models. In Sec. IV we move on to an
application of electron-boson FCIQMC to the Hubbard–Holstein model, a
paradigmatic model where strongly interacting electrons and local phonon
coupling on a lattice gives rise to a number of electronically and
phononically driven phase transitions. The approach will be compared to
deterministic exact diagonalization results with an explicitly truncated
bosonic occupation. Finally, in sec. V we move towards a more general
implementation, where electron-boson FCIQMC is used to solve a model of
electrons coupled to long-range bosonic quasi-particles. This model requires
more general bosonic coupling to arbitrary density fluctuations in the
electronic states, as well as general (four-point) electronic interactions
amongst the electrons. Furthermore, we demonstrate convergence of density
matrix sampling in this general context, as well as algorithmic changes
required for excitation generation towards an ab initio setting.

##  II Electron-boson hamiltonians

The general form of the second-quantised Hamiltonian considered in this work
can be written as

$\displaystyle\hat{H}$
$\displaystyle=\hat{H}_{\mathrm{elec}}+\sum_{m}\omega_{m}\hat{a}^{\dagger}_{m}\hat{a}_{m}$
$\displaystyle+\sum_{mpq}\sum_{\sigma\in\\{\uparrow,\downarrow\\}}(V_{mp_{\sigma}q_{\sigma}}\hat{c}^{\dagger}_{p_{\sigma}}\hat{c}_{q_{\sigma}}\hat{a}_{m}+g_{m}\hat{a}_{m})+\mathrm{h.c.}$
(1)

where $\hat{c}$ symbolises the fermionic annihilation operator, and
$\hat{H}_{\mathrm{elec}}$ contains the standard one- and two-electron
(interaction) terms over the single fermion states (labelled $p,q$) and spin
labels $\sigma,\tau$:

$\displaystyle\hat{H}_{\mathrm{elec}}$
$\displaystyle=\sum_{pq}\sum_{\sigma\in\\{\uparrow,\downarrow\\}}h_{p_{\sigma}q_{\sigma}}\hat{c}^{\dagger}_{p_{\sigma}}\hat{c}_{q_{\sigma}}$
$\displaystyle+\sum_{pqrs}\sum_{\sigma\tau\in\\{\uparrow,\downarrow\\}}g_{p_{\sigma}q_{\tau}r_{\sigma}s_{\tau}}\hat{c}^{\dagger}_{p_{\sigma}}\hat{c}^{\dagger}_{q_{\tau}}\hat{c}_{s_{\tau}}\hat{c}_{r_{\sigma}}$
(2)

The bosonic annihilation operator is denoted $\hat{a}$, with indices $m$
extending over all modes. These bosons can represent phononic, photonic or
other arbitrary bosons, which in this work will interact in a linear fashion
with the quantum fluctuations in the fermionic charge density (although these
can in principle couple to e.g. spin or pairing fluctuations or other higher-
order electronic operators). The effect of this explicit coupling on the
electronic structure is a dynamic or energy-dependent change the effective
interactions between the electrons, mediated by these bosons, to screen the
electrons at different length or energy scales in the system.

A simplified paradigmatic model system which captures these effects as a
special case of Eq. II, is the Hubbard–Holstein model [86, 87, 88, 89, 90, 91,
92, 93, 94, 95, 96, 97], as

$\hat{H}=\hat{H}_{\mathrm{Hubbard}}+\sum_{m=1}^{L}g\hat{n}_{m}(\hat{a}_{m}+\hat{a}_{m}^{\dagger})+\sum_{m=1}^{L}\omega_{0}\hat{a}^{\dagger}_{m}\hat{a}_{m}$
(3)

where
$\hat{n}_{m}\equiv\hat{n}_{m_{\uparrow}}+\hat{n}_{m_{\downarrow}}\equiv\hat{c}^{\dagger}_{m_{\uparrow}}\hat{c}_{m_{\uparrow}}+\hat{c}^{\dagger}_{m_{\downarrow}}\hat{c}_{m_{\downarrow}}$
is the fermion number operator, and the electronic
$\hat{H}_{\mathrm{Hubbard}}$ is parametrised as

$\displaystyle\hat{H}_{\mathrm{Hubbard}}$ $\displaystyle=-t\sum_{\langle
i,j\rangle}\sum_{\sigma\in\\{\uparrow,\downarrow\\}}\left(\hat{c}^{\dagger}_{j_{\sigma}}\hat{c}_{i_{\sigma}}+\hat{c}^{\dagger}_{i_{\sigma}}\hat{c}_{j_{\sigma}}\right)$
$\displaystyle+U\sum_{i=1}^{L}\hat{n}_{i_{\uparrow}}\hat{n}_{i_{\downarrow}}$
(4)

by the usual hopping between nearest neighbours ($\langle i,j\rangle$) and on-
site repulsion term. In this way, the sites of the Hubbard–Holstein model are
each coupled with strength $g$ to a local boson mode of frequency
$\omega_{0}$. The effect of these local boson modes, which can drive effective
electronic pairing between the fermions, compete with both the on-site
electronic repulsion and single-particle kinetic energy terms. It is a
prototypical model for phonon-driven emergence of charge density waves in
correlated materials and Peierls instabilities, with metallic, Mott insulating
and charge density modulated phases present depending on the dominant terms,
representing an important and unsolved model system.

##  III Adaptation of FCIQMC for electron-boson systems

Due to the coupling with bosonic degrees of freedom, the many-body Hilbert
space in which the FCIQMC walker population must evolve is one formed by the
product of fermionic and bosonic configurations. While the basic ingredients
of the algorithm are largely the same, the necessities of generating boson-
coupled excitations and accommodating bosonic occupation in the many-body
basis are significant departures from the assumptions of a determinant space
implementation. While the determinant is stored as a bit string, the bosonic
part of the configuration must be stored as an integer array allowing
arbitrary occupancy [83]. Since the Hamiltonian in Eq. II does not conserve
boson number, there is in theory no limit to the number of bosons that can be
occupied in the CI space basis. Practically speaking, high boson occupation is
energetically penalised by the positive frequency $\omega_{n}$ of the modes,
with the expected distribution of bosonic occupation following a Poissonian
distribution. This implies an exponential accuracy in representation of the
wave function with increasing truncation of the maximum boson occupation. In
practice, we store the boson occupations in the FCIQMC results of this paper
as an array of unsigned characters, giving a maximum cutoff occupation of 255
bosons per mode (though we can control this in order to enable direct
comparison to approaches with lower occupation truncations). An occupancy of
255 is well beyond any reasonable occupation visited during the stochastic
dynamics on any practical timescale, and so this can be taken to be
effectively infinite within the sampling approximation.

###  III.1 Recap of FCIQMC algorithm

FCIQMC [49, 98] is a projector monte carlo method which represents the many-
body wavefunction as a population of signed walkers at each projective
iteration. The average number of walkers on each many-body basis function will
- given a total number of walkers sufficient to resolve any present sign
problem - converge to the exact FCI coefficient. Thus, expectation values with
respect to the ground state wavefunction can be extracted by averaging their
values over a number of walker distributions generated by iterative
application of a stochastically-realised integrator.

Usually, this is based on the linear order truncation of the imaginary time
evolution operator, giving

$\Delta C_{\mathbf{i}}=-\Delta\tau\sum_{\mathbf{j\neq
i}}H_{\mathbf{ij}}C_{\mathbf{j}}-\Delta\tau(H_{\mathbf{ii}}-E_{S})C_{\mathbf{i}}$
(5)

as the average update of the coefficients. Expressed in the first term are the
off-diagonal spawning contributions. These are stochastised by making on
average $|C_{\mathbf{j}}|$ attempts to draw basis functions $\mathbf{i}$ with
non-zero matrix element $H_{\mathbf{ij}}$ in a process referred to as
excitation generation. The cost of the calculation therefore scales linearly
in the 1-norm of the amplitudes, or number of walkers
$N_{W}\equiv\sum_{\mathbf{i}}|C_{\mathbf{i}}|$.

Production scale calculations with FCIQMC codes require parallelisation, with
each basis function assigned a rank index within a message passing
communicator. Spawned contributions are sent via this interface to the process
on which the sampled $C_{\mathbf{i}}$ would be accumulated. The incoming
spawns are combined with each other and any existing walker population
$C_{\mathbf{i}}$ in the annihilation procedure. Sufficient incidence of
annihilation between positive and negative walkers is understood to be the
crucial factor in suppressing the growth of the unphysical stoquastic signal,
and thus overcoming the sign problem [76, 54]. In problems with a severe sign
problem, the required number of walkers can approach the FCI Hilbert space
dimension - nullifying the advantages of the sparse representation afforded by
the FCIQMC approach. In response to this, the initiator approximation was
introduced, which imposes a dynamic truncation of the Hamiltonian prohibiting
spawns from low-weight walkers in exchange for an approximate, but
systematically-improvable sampling of a physical ground state eigenvector. The
initiator approximation is most successful when the wavefunction is localised
in orbital space, as is often the case in weak-to-moderately correlated MO-
based ab initio or reciprocal space lattice model systems, wherein systematic
errors can be reduced to within acceptable bounds well before the unmodified
FCIQMC would be able to tractably sample the physical ground state.

The second term in Eq. 5 is a diagonal scaling term which includes the shift
$E_{S}$. By modulating the shift, the total walker creation due to
spawning/annihilation can be matched by the number removed via this diagonal
death process, maintaining an approximately constant $N_{W}$ and providing in
$E_{S}$ an estimation of the energy expectation value associated with the
sampled wavefunction.

###  III.2 Generalized Excitation Generation

Excitation generation constitutes a core component of FCIQMC, whereby
individual off-diagonal terms in the hamiltonian are selected according to a
normalized probability distribution for a given configuration, providing
another “connected” configuration in the full Hilbert space [99, 100, 101,
102]. For maximum efficiency, these terms should be selected with a
probability approximately proportional to the absolute magnitude of the
hamiltonian matrix element connecting the configurations. Since this will in
general depend on the configuration in question, this can not be efficiently
achieved in practice, and so a number of approaches have been developed,
aiming at minimizing the cost of this stochastic selection algorithm, and
ensuring as optimal a probability distribution in this sampling as possible.
We initially divide these types of
$|\mathbf{i}\rangle\rightarrow|\mathbf{j}\rangle$ excitations into different
classes, assigning initial probabilities of these class types to each. These
now include both fermionic single and double excitations, as well as the
(separate) probabilities for bosonic excitation and de-excitations. The first
step involves selecting a class of excitation to generate according to these
normalized probabilities, and the probabilities of generation within each
class are continually optimized through the calculation in order to minimize
multiple walker creation events in a single attempt, and therefore allow the
largest timestep possible with stable dynamics.

The algorithm for generating excitations from the electronic part of the
hamiltonian proceeds as normal, with single excitations
$|\mathbf{i}\rangle\rightarrow\hat{c}^{\dagger}_{i}\hat{c}_{j}|\mathbf{i}\rangle$
(with $i,j$ distinct) drawn by uniform selection among occupied and unoccupied
spin orbitals of the same spin, taking into account any abelian symmetry
constraints of these orbitals. Double excitations are drawn as
$|\mathbf{i}\rangle\rightarrow\hat{c}^{\dagger}_{i}\hat{c}^{\dagger}_{j}\hat{c}_{l}\hat{c}_{k}|\mathbf{i}\rangle$
(with $i,j,k,l$ all distinct) according to the precomputed heat bath algorithm
(PCHB), which approximately draws these indices according to the magnitude of
the hamiltonian matrix elements [99, 101, 102]. We focus here on the boson
coupling term, which introduces two novel excitation generation cases, namely
fermion hopping coupled to boson excitation and de-excitation, as

$|\mathbf{i}\rangle\rightarrow\hat{c}^{\dagger}_{i}\hat{c}_{j}\hat{a}^{\dagger}_{n}|\mathbf{i}\rangle;\quad|\mathbf{i}\rangle\rightarrow\hat{c}^{\dagger}_{i}\hat{c}_{j}\hat{a}_{n}|\mathbf{i}\rangle$
(6)

(with $i,j$ distinct). Boson number-conserving terms in the Hamiltonian are
diagonal, contributing only to the death rate of walkers, and so do not need
to be considered in the excitation generation.

A stochastically efficient excitation generator for the connections of Eq. 6
is straightforward to achieve via a reapplication of the principles of the
heat bath algorithm to effectively weight the generation of these terms. To
ensure hermiticity, the coefficients $V_{ijn}$ of the boson excitations are
the same as those $V_{jin}$ for the de-excitations, and so a single PCHB
sampler suffices for both cases. When the PCHB sampler is invoked for a boson
excitation, a uniform single excitation
$|\mathbf{i}\rangle\rightarrow\hat{c}^{\dagger}_{i\sigma}\hat{c}_{j\sigma}|\mathbf{i}\rangle$
is drawn, with the $(i,j)$ compound index referencing a row in the PCHB table.
Note that this selection of the $(i,j)$ pair is not invariant to the order of
the selection of the two indices, differing from the two-electron integral
precomputed table, which uses a triangular mapping to leverage the $i<j$
permutational symmetry of the uniformly picked annihilation operator pair.

A bosonic index is then picked with probability $|V_{ijn}|/\sum_{m}|V_{ijm}|$,
whose values are precomputed in an alias table for efficient selection. This
ensures that the these bosonic excitations are drawn from a probability
distribution as close as possible to the true hamiltonian magnitude between
the states. Unlike the fermionic PCHB implementation, drawn excitations are
only discarded as null when mode $n$ is already at the imposed
$n_{\mathrm{boson}}$ cutoff (which never occurs in practical settings where
the cutoff is 255). De-excitations are drawn in exactly the same way, except
that the indices $i,j$ are exchanged to find the correct row in the sampler
table, and are discarded when mode $n$ is unoccupied.

##  IV 1D Hubbard–Holstein model

We first test the ability to faithfully sample from the boson spaces in these
models via the 1D Hubbard–Holstein model, as given in Eq. 3. Due to the form
of the model, the excitation generation is simplified from its generality
described in Sec. III.2. The electronic hopping is stochastised by simply
selecting an occupied spin orbital and then randomly choosing to move left or
right, generating a null excitation if the adjacent spin orbital is already
occupied. Furthermore, since the model only couples bosonic operators to the
electronic density (in the site basis), not the hopping, the parent
configuration bit-string is also used for the bosonic de-excitation generation
to select a boson mode. This is efficiently achieved by caching the occupied
sites coordinated to a boson mode with non-zero occupation. An analogous list
is not prepared for the boson creation excitation generator, since in general
one is interested in the limit in which no practically no cutoff is imposed on
mode occupations, and therefore all modes coupled to occupied sites can be
excited.

Even though the term in the Hubbard–Holstein model Hamiltonian which
introduces coupling between the electronic density and bosonic
(de-)excitations has a positive coefficient (therefore inducing both positive
and negative walkers into the dynamics), the overall FCIQMC simulation in the
site basis of the model is sign-problem-free. This is because these bosonic
interactions only connect adjacent boson number sectors of the Hilbert space,
and ensures that the sign on any single configuration is constrained to a
single sign, precluding annihilation events from this part of the hamiltonian.
Walkers that are positive in one boson number sector are sign-coherently
connected to negatively-signed walkers in the adjacent sectors. Furthermore,
if there is an odd (even) number of electrons in each spin sector with
periodic (anti-periodic) boundary conditions, the Hubbard hopping term is also
sign-problem-free, while the interactions are diagonal in the basis. Despite
being a sign-problem-free system, the model can still be challenging to solve
and is the subject of continued research, making it an ideal test bed for the
efficiency and convergence of the FCIQMC algorithm.

Figure 1: Correlation energies for the Hubbard–Holstein model ($L=8$,
$n_{\mathrm{elec}}=8$, $U=2.0t$, $\omega_{0}=5.0t$, $g=\sqrt{10.0}t$) with
maximum boson occupation per mode of $n_{\mathrm{boson}}=1$, compared to the
exact diagonalisation value. These energy estimates are unbiased for
population control bias by the reweighting of shift estimates of the energy.
This correction is performed across the time series of shift measurements
($E_{S}$) and instantaneous $N_{W}$ for different correction orders ($n$).
Consecutive points in the time series from which the corrected energies are
computed were 10 Monte Carlo iterations apart, and the shift updates occur
with the same frequency and with damping parameter $\gamma=5$. The FCIQMC was
performed for $16\times 10^{6}$ MC cycles with a mean walker number of
$\bar{N}_{W}=188$.

Due to this lack of sign problem, the exact ground state eigenvalue can in
principle be estimated from the shift estimator ($E_{S}$) with an arbitrarily
small $N_{W}$, given a long enough evolution of the walker population and a
correction for the population control bias. It has been shown in recent work
[103] that this correction can be made without alterations to the FCIQMC
algorithm, and that given only the $N_{W}$ and shift time series the
systematic bias due to the statistical covariance between the shift and the
stochastic wave function can be substantially removed [104, 83, 105, 106]. The
essence of this approach is to retrospectively “undo” the dynamical scaling of
the wavefunction due to the shift-dependent death step, so that the covariance
between the wavefunction and the shift can be reduced. The crucial factor to
be determined is the optimal span of time series points $n$ on which these
“reweightings” should be performed so as to effectively remove the bias
without introducing excess random error due to the reduction in number of
effectively sampled points. We call $n$ the “correction order” in Fig. 1.

We first test this approach for a small lattice where we can perform exact
diagonalization calculations (8 sites). Obtaining the exact diagonalization
data also necessitates a severe truncation in the maximum boson occupation,
which we limit to a single boson occupancy per mode. The shift-bias corrected
FCIQMC data with the same equivalent artificial boson truncation is shown in
Fig. 1. Note that the figure shows the expected unbiased energy as a function
of the length of history of the walker population used in order to control for
population bias, not iteration number. This shows similar characteristics to
the equivalent plot in Ref. 103 for a single calculation. Initially there is
convergence towards the exact eigenvalue, then we observe the expected
diminishing returns of attempting to further reduce systematic error by
elongating the history over which the correction is computed, as the
statistical uncertainty in the estimator grows. For this model, a correction
order around 2000 appears to be a good compromise. This however serves as a
proof-of-principle of the exactness of the implementation, and ability to
compute unbiased energies in these (albeit small) systems compared to exact
results.

###  IV.1 Importance Sampling

The bias in the uncorrected shift estimator for sign-problem-free systems has
also been shown to reduce significantly with the introduction of importance
sampling, the basic premise of which is to sample not the ground state CI
coefficients $C_{\mathbf{i}}$, but the modified wave function
$\tilde{C}_{\mathbf{i}}=\langle\phi_{\mathrm{guide}}|\mathbf{i}\rangle
C_{\mathbf{i}}$ where $|\phi_{\mathrm{guide}}\rangle$ is a “guiding” wave
function. Equivalently viewed as a similarity transform of the Hamiltonian,
this approach effectively modifies the spawning probabilities such that moves
to configurations with larger diagonal matrix element are probabilistically
suppressed. Expectation values can then be unbiased for this modification to
the dynamic. [107, 108] Even with a very simple form for the guiding wave
function, such as that of the Gutzwiller-like form

$|\phi_{\mathrm{guide}}\rangle=\sum_{\mathbf{i}}e^{-\alpha
H_{\mathbf{\mathbf{ii}}}}|\mathbf{i}\rangle,$ (7)

significant improvements have been found in fermionic model systems. We expect
this guiding wave function to be even more powerful in the case of boson
number non-conserving Hamiltonians where inhibition of spawning to more
energetic configurations with highly occupied bosonic modes is an obvious
route to achieving a more compressed stochastic representation of the ground
state.

(a) 8-site chain with walker number of $N_{W}=2\times 10^{4}$ and
$n_{\mathrm{boson}}=2$, an exactly diagonalizable system.

(b) Reduction in population control bias for increasing Hubbard–Holstein chain
lengths. Each series was obtained with $n_{\mathrm{boson}}=255$ and a target
$N_{W}$ of $10^{6}$.

Figure 2: Population control bias amelioration by importance sampling of the
Hubbard–Holstein ground state ($U=2.0t$, $\omega_{0}=5.0t$, $g=\sqrt{10.0}t$)
using a guiding wave function of the type defined in Eq. 7, with variation in
the single adjustable parameter, $\alpha$. All energies are shown relative to
the average shift with no importance sampling ($\alpha=0$).

Figure 2(a) illustrates an example application of this adaptation where an
$\alpha$ value of $\sim 0.2$ is seen to remove the bias almost completely,
even without the explicit reweighting algorithm of Fig. 1, though a
significant bosonic occupation cutoff is applied ($n_{\textrm{boson}}=2$) in
order to enable comparison to exact diagonalization. The plot also shows the
expected trend of increasing random error in the shift estimator due to the
importance sampling [103]. Intuitively, this is due to the reciprocal
relationship between the importance sampling factors for proposed moves
$\mathbf{i}\rightarrow\mathbf{j}$ and $\mathbf{j}\rightarrow\mathbf{i}$. As a
move $\mathbf{i}\rightarrow\mathbf{j}$ becomes more suppressed, the variance
of the walker population (and true CI coefficient) on $\mathbf{j}$ increases,
and the back-spawning magnitudes also increase. This boosted reproduction rate
of ill-defined walker occupations propagates more error. Fortunately, it seems
that the increases in noise are modest in comparison to the dramatic reduction
in the systematic sampling error due to population control bias.

###  IV.2 Relaxing Bosonic Occupation Constraints

Having established the ability to practically obtain unbiased ground state
energetics via electron-boson FCIQMC in small systems with boson occupation
constraints, we now turn to the practicality of relaxing these boson
occupation constraints of any mode. This requires the FCIQMC to effectively
sample from an infinite Hilbert space of configurations comprising all
possible combinations of boson occupations of each mode. We will consider the
convergence with respect to relaxing the bosonic threshold on both the
(population bias corrected and uncorrected) ground state energies, as well as
the convergence of these FCIQMC estimates with respect to the number of
walkers.

We start by considering the Mott regime of the 8-site Hubbard–Holstein model,
with $U=4t$, $\omega_{0}=0.5t$ and $g=\sqrt{0.15}t$. In this hamiltonian, the
phononic renormalization is not sufficient to overcome the strong Coulomb
repulsion, and Mott order remains in the system. While strongly correlated,
the relatively weak boson coupling should allow for consideration of the
convergence of the energies with number of walkers compared to exact
diagonalization as the occupation is increased. Furthermore, we can transform
the hamiltonian, by performing a shifting of the boson operators to remove the
coupling of the phonons to the static charge density of the system (which is
known exactly by symmetry). This zero phonon mode removal (ZPMR)
transformation can leave the exact eigenstates of the full Hamiltonian
unchanged (as long as it is accounted for in the expectation values), but can
dramatically increase the speed of convergence with respect to average bosonic
occupation. This is because the bosons now only couple to charge fluctuations,
rather than the mean-field static component. This allows exact diagonalization
to converge faster with boson occupation, and we can compare this to the
unconstrained FCIQMC as walker numbers increase. We note that while this
transformation of the hamiltonian can in principle also be performed for
FCIQMC, it is not in this case, since it would induce a sign problem in the
resulting dynamics. The details of this transformation are given in Appendix
A.

Figure 3: Approximately linear convergence of systematic energy errors with
respect to $1/N_{W}$ in the Mott insulator phase of the Hubbard–Holstein model
at half filling with $L=8$, $U=4t$, $\omega_{0}=0.5t$, and $g=\sqrt{0.15}t$,
as well as energies explicitly corrected for this bias showing agreement with
exact diagonalization results. We show uncorrected ($E_{S}$) and explicitly
reweighted (corr. $E_{S}$) average energy (shift) estimators from the FCIQMC,
for the model with $n_{\textrm{boson}}$ truncated to 1, 2 and untruncated
($n_{\textrm{boson}}=255$ for practical purposes). These first two results are
compared to exact diagonalization in the same Hilbert space, while we compare
the untruncated FCIQMC results to exact diagonalization with the shifted
bosonic operators and $n_{\textrm{boson}}=2$, which should be well converged
in this regime (see Appendix A), showing good agreement of the corrected
results at all walker numbers.

Fig. 3 shows the convergence of the ground state FCIQMC energy estimates with
respect to reciprocal walker number, for both the uncorrected and a posteriori
corrected shift energy estimates. As has been noted elsewhere (although still
somewhat debated), we find a linear convergence of the population control bias
with respect to reciprocal walker number. While this error is reduced with the
importance sampling of the dynamic via the trial wave function, this doesn’t
sufficiently remove the bias, and the reweighting approach is still required.

Once this reweighting is applied, the energies at all walker numbers are
unbiased, with the exact diagonalization results for a boson occupation
truncation of one and two within the statistical errors of the FCIQMC.
Relaxing the boson occupation entirely ($n_{\textrm{boson}}=255$ for technical
reasons, but this threshold is never reached) the energies are relaxed by a
further $\sim 0.2t$. These agree well with the exact diagonalization result in
the presence of the ZPMR transformation of the boson operators, despite a
boson occupation truncated to $n_{\textrm{boson}}=2$, which we believe to also
be well converged in this regime with weak boson coupling. We find our best
FCIQMC result for the model to be $-7.1787(17)t$ at $\sim 7\times 10^{3}$
walkers, with the ZPMR exact diagonalization result with
$n_{\textrm{boson}}=2$ to be $-7.18408t$.

(a) Convergence of the uncorrected and corrected shift energies for the ground
state with respect to $n_{\mathrm{boson}}$ truncation. The shaded area
represents the random error in the corrected shift energy with
$n_{\mathrm{boson}}=255$ (in practice, an untruncated bosonic occupation). All
means and standard errors were obtained by blocking [109] of the shift time
series with $\bar{N}_{W}\sim 800$.

(b) Histograms of FCIQMC walker weight according to either the largest single
mode boson occupation of the configuration (blue), or the sum of the boson
occupation over all bosonic modes (orange). Note that appreciable amplitude of
the walkers in this system reside on configurations with up to ten bosons in
any one mode, reflecting the energetic impact of the bosonic truncation seen
in (a). The calculation was performed with $n_{\mathrm{boson}}=255$ and a
large enough $\bar{N}_{W}$ ($\sim 2\times 10^{5}$) to obtain the correct shift
without the use of reweighting.

Figure 4: Demonstration of the effect of many-body wavefunction truncation due
to the imposition of boson mode occupation cutoffs for a small
Hubbard–Holstein chain without importance sampling ($L=8$,
$n_{\mathrm{elec}}=8$, $U=2.0t$, $\omega_{0}=5.0t$, $g=\sqrt{10.0}t$).

We now turn to a parameter regime of the 1D Hubbard–Holstein model where the
bosonic renormalization effects are strong enough to induce a Peierls charge
ordering, breaking the small Mott order present. This anti-adiabatic phase
will be characterized by $U=2t$, $\omega_{0}=5.0t$, and $g=\sqrt{10.0}t$. In
this regime, we expect a much slower convergence of the energy with respect to
increasing boson occupation, and we can explicitly verify that this
convergence correctly reproduces the infinite boson occupation limit, without
any increase in bias or statistical error in the FCIQMC simulations. These
results are shown in Fig. 4(a), where explicit boson occupations up to eight
are required for convergence to within the small FCIQMC errorbars, with the
remaining population control bias accounting for $\sim 0.15t$. These FCIQMC
calculations were performed without importance sampling with $\sim 800$
walkers at all values of $n_{\textrm{boson}}$, including the infinite boson
limit, which agrees with the $n_{\textrm{boson}}=10$ result both in value and
magnitude of stochastic error. This indicates that no additional difficulties
arise in the FCIQMC by fully relaxing this boson occupation constraint, which
would not be possible to obtain with exact diagonalisation. In Fig. 4(b), the
fraction of wavefunction amplitude histogrammed according to either the sum of
boson occupation, or the largest single mode occupation can be seen,
illustrating non-negligable walker weight for this system on configurations
with up to ten bosons in any single mode. We find a best FCIQMC estimate of
this point in the parameter space of the 8-site 1D Hubbard–Holstein model of
$-28.612(6)t$.

###  IV.3 Approach to the Thermodynamic Limit

With the ability of FCIQMC to faithfully and efficiently sample from
unconstrained bosonic occupations, we now turn to increasing the physical
system size, and demonstrate convergence to the thermodynamic limit of the
model. We remain in the anti-adiabatic Peierls regime where strong bosonic
coupling distorts the charge density and large bosonic occupations are
required to be sampled as shown previously in Fig. 4(a).

Having recognised the benefits of the use of importance sampling, we return to
the matter of optimising a guiding wavefunction. Figure 2 shows the shift
energy improvement with respect to $\alpha$ for a selection of HH chains. In
each of these systems the shift converges from above, and so it is valid to
take the $\alpha$ responsible for the most negative energy difference to be
the optimal value.

In Figure 2(b), the absolute energy density correction due to importance
sampling for a selection of chain lengths for a constant computational cost
$N_{W}$ is evident, with an $\alpha$ of $\sim 0.05$ conferring optimal bias
reduction regardless of chain length. These data demonstrate that optimisation
of the $\alpha$ parameter can be done on cheaper, low $N_{W}$ FCIQMC
evolutions before the resulting optimal value is used with larger $N_{W}$
values, and in conjunction with the a posteriori reweighted shift correction.
Neither of these measures are capable of eliminating the bias completely, but
they are effective means of accelerating the systematic improvement with
respect to $N_{W}$.

Fig. 5 shows the convergence of the FCIQMC energy per site as the chain length
increases to 48 sites, with $\sim 10^{7}$ walkers in each simulation. A power
law extrapolation of the energies achieves a highly accurate unbiased FCIQMC
prediction for the thermodynamic limit energy density of $-3.55790(13)t$ per
site for this point in the phase diagram. The fitting was repeated with the
omission of the longer chains which are associated with the largest random
error; the resulting energy densities of $-3.55782(18)t$ for $L\leq 42$ and
$-3.55812(4)t$ for $L\leq 24$ are found to agree with the $L\leq 48$ estimate
within standard errors. The latter result indicates a $\sim 2\times 10^{-4}t$
under-correction of the bias for $L>24$, but this is within standard error and
expected given that $N_{W}$ was not increased with $L$. Also shown is the
increase in the random error (as estimated by blocking analysis) of the
corrected shift with respect to $L$. Since each system was tackled with an
equal $N_{W}$ and in the same amount of CPU time, this error can be identified
as directly corresponding to the scaling in the effective cost of the
calculations for the unbiased energy to a given statistical uncertainty as a
function of increasing chain length.

Figure 5: Convergence of the FCIQMC energy per site of the Peierls regime 1D
Hubbard–Holstein model ($U=2.0t$, $\omega_{0}=5.0t$, $g=\sqrt{10.0}t$).
Approximately $10^{7}$ walkers were used in each calculation, with importance
sampling and the reweighting correction both applied to obtain unbiased
results. The curve follows a power law fit to these data giving a
thermodynamic limit energy density of $-3.55790(13)t$. The inset shows the
random error in the corrected shift as a function of the chain length.

##  V FCIQMC for Electron-Boson Embedded Models

While a more extensive investigation of the Hubbard–Holstein model and its
phase diagram with FCIQMC is planned for future work, we want to also stress
the generality of the electron-boson FCIQMC approach and implementation,
including to more complex interactions and couplings, and in the presence of
sign problems. Rather than an ab initio electron-phonon hamiltonian, we will
consider the challenge of solving an auxiliary effective electron-boson
hamiltonian that arises within the context of a recently-developed quantum
embedding theory. Here, the bosons in the embedded hamiltonian are fictitious,
but their presence mimics the effect of interactions between an electronic
fragment of a system with its (potentially large) electronic environment. In
this way, the bosons describe the coupling of local electronic degrees of
freedom to long-range collective electronic excitations in its environment.
While these collective excitations represent composite fermionic
quasiparticles, they obey (quasi-)bosonic statistics. Bosons therefore
naturally arise in the description of wholly electronic systems. The
hamiltonian of the embedded cluster in this context has a “dense” coupling of
the bosons to arbitrary fermionic hopping in the fragment electronic
subsystem, as well as a dense four-point electronic interaction and presence
of a sign problem. Furthermore, reproducing the properties of the full system
in this approach requires a sampling of reduced density matrices of this model
hamiltonian, and we will investigate the fidelity by which these important
quantities (including the coupled electron-boson reduced density matrix) can
be sampled in FCIQMC.

This recent quantum embedding theory, which maps an electronic system to an
auxiliary electron-boson coupled problem describing a small fragment of the
original system, is dubbed the “extended density matrix embedding theory”
(EDMET)[25], building on a “parent” mean-field quantum embedding formalism.
The resulting local cluster hamiltonian has the form of Eq. II. The fermionic
degrees of freedom are comprised of a small fragment of the sites of the full
system, coupled to a fermionic bath of the same dimensionality. The
hamiltonian in this fermionic space is algebraically constructed to ensure
that the mean field properties of the fragment match between the cluster model
and full system. In a similar spirit, the bosonic modes are constructed
algebraically from the random phase approximation (RPA) solution over the full
system, in order to require that the RPA description of the two-body fragment
properties match those of the RPA over the full system. In this way, a simple
construction rigorously maps from a large electronic system to a coupled
electron-boson model for a fragment of the solution, with the properties of
the whole system accessible from the reduced density matrices of this embedded
problem. Full details of the method and the specifics of the cluster
hamiltonian construction can be found in Ref. 25.

We now consider the ability of FCIQMC to provide the required accuracy in the
solution of this embedded cluster hamiltonian. Two differences from the form
of Eq. II exist in the resulting cluster model; $g_{n}=0$ (i.e. there are no
single boson terms not coupled to the charge density), and the electron-boson
coupling term is zero when the fermionic indicies $i$ and $j$ co-incide, if
the fermionic cluster basis is canonicalized to a mean-field representation
(i.e. the bosons represent couplings to the fragment density fluctuations,
rather than the mean-field density, and so ZPMR is not effective). To date,
only exact diagonalization has been used as a solver of these embedded
hamiltonians, which puts a severe constraint on both the size of the fragment
spaces which can be chosen, as well as the boson occupation truncation to
ensure computational tractability. Physically, this truncation corresponds to
limiting the energy-scales of long-range plasmonic excitations that the
fragment can couple to, and it would be advantageous to not truncate these
effects, and so alternative solvers for these auxiliary model hamiltonians is
urgently sought. Using electron-boson FCIQMC will however require an ability
to sample the reduced density matrices (RDM) of the coupled electron-boson
model in FCIQMC, including electron-boson terms (which will also be useful for
other property estimates in a more general context).

###  V.1 Reduced Density Matrix Sampling

The standard way to stochastically estimate (off-diagonal) elements of RDMs in
FCIQMC is to repurpose the spawning step of the algorithm. This step can be
considered the act of randomly drawing a configuration pair, and a product of
two instantaneous walker occupations taken from statistically independent
replica simulations, which can then be used to estimate their contribution to
unbiased RDMs [110, 111]. This has proven a successful means of sampling one-
and two-body RDMs in fermionic systems for a diversity of applications, e.g.
molecular properties [112, 113], MCSCF orbital optimisations [59, 114, 115,
116] and contributions from internally contracted perturbers in multireference
perturbation theories [117, 118]. However, this approach makes an
approximation that the set of configurational pairs, $(\mathbf{i},\mathbf{j})$
with non-zero coefficient product which contribute to the RDMs of interest
($\Gamma$-connections) is equal to the set of $(\mathbf{i},\mathbf{j})$ pairs
connected by a non-zero Hamiltonian matrix element ($H$-connections). This is
required as it is impractical to essentially take the outer product of a
sparse, computationally distributed representation of the FCI wavefunction so
that all $\Gamma$-connections are explicitly included, and so the feasibility
of accurate, fully stochastic RDM estimates within an FCIQMC computational
framework is predicated on the success of the assumed equivalence of $H$\- and
$\Gamma$-connections.

A stark example of the breakdown of this assumption is in the case of
Hartree–Fock canonical orbital representations, in which single excitations of
the Hartree–Fock determinant can have substantial coefficients, but do not
interact with the reference through the hamiltonian, $H$, so relying on
spawning alone to sample the 1RDM would neglect some large contributions. This
class of $\Gamma$-connections is therefore explicitly included as a separate
contribution. An approach to including these contributions more widely, is to
ensure that important $\Gamma-$connected configuration pairs (which may not be
significantly weighted $H-$connected configurations) are included in a
“deterministic” space, which is already widely used in FCIQMC propagation[119,
120]. Here, coefficients on configurations designated to comprise the
deterministic subspace are collected from across all processes and multiplied
exactly by a sparse representation of the Hamiltonian projected into that
subspace. In a semi-stochastic RDM calculation, another sparse matrix object
can be initialised which only stores those $\Gamma$-connections which are not
also $H-$connections so that all RDM contributions between configurations in
the deterministic subspace are made without approximation. As the size of the
deterministic subspace increases, any error of the kind described above is
rigorously removed.

However, ensuring all $\Gamma-$connected configurations are included in a
deterministic subspace may not always be efficient or practical. Sampling of
the remaining $\Gamma$-connections (those which do not involve pairs of
configurations in the deterministic subspace) can be achieved by the
implementation of special “ghost” excitation generators, which are propagated
separately to the spawning events of FCIQMC and only carry the instantaneous
population of the parent configuration (ignoring generated connections that
would already be accounted for by the deterministic subspace), specifically
for RDM accumulation. Such a situation has already been covered in detail in
the accumulation of the three-body RDM and four-body contracted auxiliary
matrices for MRPT2[117], where these higher-rank contributions are explicitly
outside the rank of the standard electronic hamiltonian. The necessity to
confront the discrepancy between $H-$ and $\Gamma-$connected configuration
pairs manifests again in the present context, since the fermion-boson RDM
$\langle\hat{c}^{\dagger}_{i}\hat{c}_{j}\hat{a}^{\dagger}_{n}\rangle$ would
not receive contributions to any elements of the form
$\langle\hat{c}^{\dagger}_{i}\hat{c}_{i}\hat{a}^{\dagger}_{n}\rangle$ since
there are no direct interactions ($H-$connections) between the mean-field
electronic densities and boson (de-)excitations. Thus, a whole class of
configuration pairs is omitted by assuming $H-$ and $\Gamma$-connection
equivalence when propagating the EDMET cluster Hamiltonian. The same would be
true not just for the EDMET hamiltonian, but also general hamiltonians where
the ZPMR transformation of the bosons has been applied (generalized to remove
coupling to a mean-field density).

A possible solution is to implement a uniform excitation generator to draw
“ghost” connections of the form
$|\mathbf{i}\rangle\rightarrow\hat{a}^{\dagger}_{n}|\mathbf{i}\rangle$,
$|\mathbf{i}\rangle\rightarrow\hat{a}_{n}|\mathbf{i}\rangle$, purely for the
purpose of sampling these neglected contributions to the fermion-boson coupled
RDM, as was done for triple excitations in the generation of 3-body RDMs.
However in the present work, this class of $\Gamma$-connection was taken into
account by the deterministic space alone, and so any RDM sampling must be
checked for convergence with respect to the number of configurations in this
subspace, denoted $N_{\text{SS}}$. It should be stressed that a more general
hamiltonian which also coupled to the full quadratic fermion terms would not
encounter this problem.

###  V.2 Extended Hubbard Model

We consider FCIQMC as the auxiliary hamiltonian “solver” within an EDMET
framework, applied to the 1D Hubbard model with extended interactions. This
purely fermionic model can be written as

${\hat{H}}_{\mathrm{Ex-Hubbard}}={\hat{H}}_{\mathrm{Hubbard}}+V\sum_{\langle
i,j\rangle}\hat{n}_{i}\hat{n}_{j},$ (8)

and we aim to converge to the thermodynamic limit of the model. The additional
nearest neighbour repulsion fundamentally changes the phase diagram from the
simple Mott insulating phase of the 1D ${\hat{H}}_{\textrm{Hubbard}}$, to one
where the competing $V$ term can drive the system into a charge ordered phase
to minimize nearest neighbour density-density repulsion. The phase diagram
therefore roughly consists of Mott insulating phases for $V\lesssim U/2$ and a
charge ordered phase for $V\gtrsim U/2$ [121, 122, 123, 124]. Rather than
solve this directly in the large system limit, we aim to solve a local model
of a cluster of $n_{\textrm{frag}}$ local sites, coupled to fermionic and
bosonic degrees of freedom approximating the coupling and entanglement with
the wider system (a bath space). The long range charge fluctuations induced by
the nearest neighbour interactions will then be modelled by the random phase
approximation, and used to map to a local model of bosons coupled to the
interacting local fermionic space. Charge ordering in the model can then be
induced by these bosons in the auxiliary model, rather than collective long-
range fermionic fluctuations. The energy density of the model over the
fragment space can then be computed from the one- and two-body fermionic
density matrices of the model, as well as the boson-fermion coupled density
matrix, as $\langle\hat{c}_{i}^{\dagger}{\hat{c}}_{j}\hat{a}_{n}\rangle$ and
$\langle\hat{c}_{i}^{\dagger}{\hat{c}}_{j}\hat{a}^{\dagger}_{n}\rangle$.

We consider a specific point in the phase diagram with $U=2t$ and $V=1t$,
which is very close to the phase transition point between Mott insulator and
charge ordered phases [124]. We expect therefore that the model will feature
significant long range two-point fluctuations, which are hard to model in a
local model and will have to be mimicked via strong coupling to the bosonic
bath. We first consider a four-site fragment space, to ensure convergence of
the system properties with FCIQMC. Mapping a 42 site chain (close to the
thermodynamic limit for the energy density) to the auxiliary hamiltonian of
Eq. II, with eight fermionic degrees of freedom and 10 boson modes, we
consider the convergence with number of walkers in Fig. 6. Additionally, we
check that changes to the deterministic space does not change the fidelity of
the RDMs, from which the energies are computed.

Figure 6: The FCIQMC-EDMET energy per site for an extended Hubbard model
($L=42,U=2.0,V=1.0$) as determined from multiple independent FCIQMC
estimations of the required RDMs of a 4-site fragment size. This leads to a
auxiliary hamiltonian consisting of 8 fermions and 10 bosons, with no cutoff
on bosonic occupation. The calculations were carried out for two different
sizes of semi-stochastic subspace to ensure convergence with respect to RDM
contributions not generated through spawning

We find a very rapid convergence of the FCIQMC energies with respect to
increasing walker number, with all values between $N_{W}=10^{5}$ and
$N_{W}=10^{7}$ statistically indistinguishable within their random error bars,
estimated from multiple calculations with different random number seeds.
Convergence with walker number is expected due to the use of the “initiator”
approximation in the FCIQMC dynamics to control the sign problem. In these
simulations, the initiator approximation is employed with the threshold walker
number of $n_{a}=3$ in order to define an initiator configuration. This is
without any restriction of bosonic occupation, and shows that even in the
presence of a sign problem in this model, these systems can be faithfully
converged, relying on a combination of annihilation events and the initiator
approximation to converge. Furthermore, changing the size of the deterministic
space ($N_{\mathrm{SS}}$) does not further change the quality of the results,
but does reduce the statistical fluctuations. The scale of these fluctuations
is however already very small, with the EDMET energy density for this fragment
size reliably found within $10^{-4}t$ by $\sim 10^{6}$ walkers, demonstrating
the applicability of fermion-boson FCIQMC as a cluster solver for this
embedding method.

Figure 7: Comparison of DMRG and EDMET using the FCIQMC cluster solver for a
half-filled extended Hubbard model ($L=42,U=2.0,V=1.0$) as the system size
increases. The inset shows the systematic improvement of the EDMET per-site
energy with respect to the size of the number of sites in the fragment.
Convergence with fragment size is imperceptible on the scale of the main plot.

We then consider the convergence of the FCIQMC-EDMET energies as the fragment
size increases. This model can be effectively solved within the density matrix
renormalization group (DMRG) given the one-dimensional nature of the model,
for which we use the block code, with 1000 renormalized states retained in the
matrix product state of each site. In Fig. 7 we show the convergence of the
energy density to the thermodynamic limit for this parameter regime as the
chain length increases. Taking the largest chain length of 42 sites, we also
map to the local coupled electron-boson cluster model with EDMET for different
numbers of fragment sites, and solve with FCIQMC. The convergence with
$n_{\mathrm{frag}}$ is shown in the inset, and demonstrates the monotonic
convergence, to an accuracy in the ground state energy of $\sim 10^{-3}t$,
with both the EDMET mapping and the FCIQMC solver giving systematic errors
smaller than the residual finite size errors. The $n_{\textrm{frag}}=6$
cluster problem consists of a 12 orbital, 12 electron fermionic space, coupled
to 21 auxiliary bosonic modes, which is easily treated, converging within
$\sim 10^{-4}t$ between $10^{6}$ and $3\times 10^{6}$ walkers after 3000 RDM
accumulation cycles.

Further improvements in the FCIQMC-EDMET approach can be found via self-
consistency, which we neglect in this demonstration. In this, the RDMs are
used to update both the mean-field and interaction kernel of the RPA,
subsequently modifying the resulting cluster hamiltonian in the bath space. We
will explore the potential of FCIQMC-EDMET in future work, with the main
conclusion here being the applicability of FCIQMC to electron-boson systems
for more general hamiltonians and as a solver in this workflow, with the
initiator approximation and annihilation events of the sparse representation
allowing for faithful convergence of both energies and RDMs even in the
presence of a sign problem in these infinite Hilbert spaces.

##  VI Conclusions and Outlook

We have developed an extension of the FCIQMC method to enable application to
coupled electron-boson systems, avoiding the truncation in bosonic occupancy
which can limit the efficiency of other wave function based approaches. We
describe the changes for efficient excitation generation in these systems, as
well as considerations in sampling reduced density matrices. We demonstrate
the convergence and applicability of the method in two very different physical
contexts. Firstly it was applied to a sign-problem-free system (the 1D
Hubbard–Holstein model) where sparse sampling and care to control for
population bias was shown to enable unbiased results in the thermodynamic
limit. Furthermore, the method was used as a solver of an embedded hamiltonian
within the extended density matrix embedding method, which maps large
electronic systems to local interacting electron-boson auxiliary models to
solve. These generalized electron-boson models induce sign problems in the
FCIQMC dynamics, which we show can be easily controlled via standard FCIQMC
techniques, without having to resort to truncation of the boson occupancy.
Furthermore, the extraction of the full system energy density requires a
faithful sampling of the reduced density matrices (including the electron-
boson RDM), which we demonstrate are faithfully sampled and compare to DMRG
reference values in the full model.

Although the results presented here are for the ground state, FCIQMC’s
dynamical adaptations are in principle completely compatible with mixed
fermion-boson Hamiltonians. In this regard however, it should be noted that
immunity to the sign problem — which is key to the method’s efficacy in the
Hubbard–Holstein model — is a state-specific attribute, and that any attempt
to sample excited states within the walker dynamic, either by explicit
projection [63] or real-time propagation [67] will necessitate annihilations,
thereby reintroducing a sign problem. Krylov projection [56, 65] on the other
hand, does not interfere with the sign-problem-free propagation towards the
ground state; but attempting to sample matrix elements, particularly those of
the overlap matrix between Krylov space-spanning replica walker populations
would be marred by the highly delocalised nature of the wavefunction in the
site basis. None of these caveats however are especially applicable to the
more localised and sign-problematic low-lying eigenfunctions encountered in
the systems described in Section V.

Having established the applicability of FCIQMC in a new domain, a wide array
of physical phenomena are now available for study, and future work must apply
the approach to open problems to assess the merits and limits of the method in
the context of complementary methods. We expect particular efficiency of
FCIQMC in cases of general interacting hamiltonian forms, with many bosons
that are relatively weakly coupled. This case allows for a sparse
representation of the wave function to be efficiently sampled within FCIQMC,
and enables the annihilation and initiator approximations to work at their
best, without requiring a priori truncations on bosonic occupation. The use of
FCIQMC as a solver within EDMET is an example of this type of setting, and
will be further explored, including self-consistency in future work, where the
ability to obtain unbiased RDMs is a particular strength of the method
compared to more traditional QMC approaches to electron-boson problems.

###### Acknowledgements.

G.H.B. gratefully acknowledges support from the Royal Society via a University
Research Fellowship, as well as funding from the European Union’s Horizon 2020
research and innovation programme under grant agreement No. 759063. We are
grateful to the UK Materials and Molecular Modelling Hub for computational
resources, which is partially funded by EPSRC (EP/P020194/1).

## References

  * Bernardi [2016] M. Bernardi, First-principles dynamics of electrons and phonons*, The European Physical Journal B 89, 239 (2016).
  * Giustino [2017] F. Giustino, Electron-phonon interactions from first principles, Rev. Mod. Phys. 89, 015003 (2017).
  * Hasan and Kane [2010] M. Z. Hasan and C. L. Kane, Colloquium: Topological insulators, Rev. Mod. Phys. 82, 3045 (2010).
  * Fröhlich [1954] H. Fröhlich, Electrons in lattice fields, Advances in Physics 3, 325 (1954), https://doi.org/10.1080/00018735400101213 .
  * Holstein [1959] T. Holstein, Studies of polaron motion: Part i. the molecular-crystal model, Annals of Physics 8, 325 (1959).
  * Bardeen [1973] J. Bardeen, Electron-phonon interactions and superconductivity, Science 181, 1209 (1973), https://www.science.org/doi/pdf/10.1126/science.181.4106.1209 .
  * von der Linden  _et al._ [1995a] W. von der Linden, E. Berger, and P. Valášek, The hubbard-holstein model, Journal of Low Temperature Physics 99, 517 (1995a).
  * Ribeiro  _et al._ [2018] R. F. Ribeiro, L. A. Martínez-Martínez, M. Du, J. Campos-Gonzalez-Angulo, and J. Yuen-Zhou, Polariton chemistry: controlling molecular dynamics with optical cavities, Chem. Sci. 9, 6325 (2018).
  * Kéna-Cohen and Forrest [2010] S. Kéna-Cohen and S. R. Forrest, Room-temperature polariton lasing in an organic single-crystal microcavity, Nature Photonics 4, 371 (2010).
  * Ruggenthaler  _et al._ [2018] M. Ruggenthaler, N. Tancogne-Dejean, J. Flick, H. Appel, and A. Rubio, From a quantum-electrodynamical light–matter description to novel spectroscopies, Nature Reviews Chemistry 2, 0118 (2018).
  * Pellegrini  _et al._ [2015] C. Pellegrini, J. Flick, I. V. Tokatly, H. Appel, and A. Rubio, Optimized effective potential for quantum electrodynamical time-dependent density functional theory, Phys. Rev. Lett. 115, 093001 (2015).
  * Mordovina  _et al._ [2020] U. Mordovina, C. Bungey, H. Appel, P. J. Knowles, A. Rubio, and F. R. Manby, Polaritonic coupled-cluster theory, Phys. Rev. Research 2, 023262 (2020).
  * Rokaj  _et al._ [2018] V. Rokaj, D. M. Welakuh, M. Ruggenthaler, and A. Rubio, Light–matter interaction in the long-wavelength limit: no ground-state without dipole self-energy, Journal of Physics B: Atomic, Molecular and Optical Physics 51, 034005 (2018).
  * Rivera  _et al._ [2019] N. Rivera, J. Flick, and P. Narang, Variational theory of nonrelativistic quantum electrodynamics, Phys. Rev. Lett. 122, 193603 (2019).
  * Lerario  _et al._ [2017] G. Lerario, A. Fieramosca, F. Barachati, D. Ballarini, K. S. Daskalakis, L. Dominici, M. De Giorgi, S. A. Maier, G. Gigli, S. Kéna-Cohen, and D. Sanvitto, Room-temperature superfluidity in a polariton condensate, Nature Physics 13, 837 (2017).
  * Stranius  _et al._ [2018] K. Stranius, M. Hertzog, and K. Börjesson, Selective manipulation of electronically excited states through strong light–matter interactions, Nature Communications 9, 2273 (2018).
  * Trügler  _et al._ [2017] A. Trügler, U. Hohenester, and F. J. García de Abajo, Plasmonics simulations including nonlocal effects using a boundary element method approach, International Journal of Modern Physics B 31, 1740007 (2017), https://doi.org/10.1142/S0217979217400070 .
  * Ruggenthaler  _et al._ [2014] M. Ruggenthaler, J. Flick, C. Pellegrini, H. Appel, I. V. Tokatly, and A. Rubio, Quantum-electrodynamical density-functional theory: Bridging quantum optics and electronic-structure theory, Phys. Rev. A 90, 012508 (2014).
  * Flick  _et al._ [2019] J. Flick, D. M. Welakuh, M. Ruggenthaler, H. Appel, and A. Rubio, Light–matter response in nonrelativistic quantum electrodynamics, ACS Photonics 6, 2757–2778 (2019), pMID: 31788500, https://doi.org/10.1021/acsphotonics.9b00768 .
  * Sjöstrand and Aryasetiawan [2022] T. J. Sjöstrand and F. Aryasetiawan, Magnon-phonon interaction and the underlying role of the pauli exclusion principle, Phys. Rev. B 105, 014433 (2022).
  * Wölfle [2018] P. Wölfle, Quasiparticles in condensed matter systems, Reports on Progress in Physics 81, 032501 (2018).
  * Rohringer  _et al._ [2018] G. Rohringer, H. Hafermann, A. Toschi, A. A. Katanin, A. E. Antipov, M. I. Katsnelson, A. I. Lichtenstein, A. N. Rubtsov, and K. Held, Diagrammatic routes to nonlocal correlations beyond dynamical mean field theory, Rev. Mod. Phys. 90, 025003 (2018).
  * Bickers  _et al._ [1989] N. E. Bickers, D. J. Scalapino, and S. R. White, Conserving approximations for strongly correlated electron systems: Bethe-salpeter equation and dynamics for the two-dimensional hubbard model, Phys. Rev. Lett. 62, 961 (1989).
  * Medvedeva  _et al._ [2017] D. Medvedeva, S. Iskakov, F. Krien, V. V. Mazurenko, and A. I. Lichtenstein, Exact diagonalization solver for extended dynamical mean-field theory, Phys. Rev. B 96, 235149 (2017).
  * Scott and Booth [2021] C. J. C. Scott and G. H. Booth, Extending density matrix embedding: A static two-particle theory, Phys. Rev. B 104, 245114 (2021).
  * Lazzeri  _et al._ [2008] M. Lazzeri, C. Attaccalite, L. Wirtz, and F. Mauri, Impact of the electron-electron correlation on phonon dispersion: Failure of lda and gga dft functionals in graphene and graphite, Phys. Rev. B 78, 081406 (2008).
  * Flick and Narang [2018] J. Flick and P. Narang, Cavity-correlated electron-nuclear dynamics from first principles, Phys. Rev. Lett. 121, 113002 (2018).
  * de Melo and Marini [2016] P. M. M. C. de Melo and A. Marini, Unified theory of quantized electrons, phonons, and photons out of equilibrium: A simplified ab initio approach based on the generalized baym-kadanoff ansatz, Phys. Rev. B 93, 155102 (2016).
  * Faber  _et al._ [2011] C. Faber, J. L. Janssen, M. Côté, E. Runge, and X. Blase, Electron-phonon coupling in the c60 fullerene within the many-body $gw$ approach, Phys. Rev. B 84, 155104 (2011).
  * Yin  _et al._ [2013] Z. P. Yin, A. Kutepov, and G. Kotliar, Correlation-enhanced electron-phonon coupling: Applications of $gw$ and screened hybrid functional to bismuthates, chloronitrides, and other high-${T}_{c}$ superconductors, Phys. Rev. X 3, 021011 (2013).
  * Antonius  _et al._ [2014] G. Antonius, S. Poncé, P. Boulanger, M. Côté, and X. Gonze, Many-body effects on the zero-point renormalization of the band structure, Phys. Rev. Lett. 112, 215501 (2014).
  * Faber  _et al._ [2015] C. Faber, P. Boulanger, C. Attaccalite, E. Cannuccia, I. Duchemin, T. Deutsch, and X. Blase, Exploring approximations to the $gw$ self-energy ionic gradients, Phys. Rev. B 91, 155109 (2015).
  * Monserrat [2016] B. Monserrat, Correlation effects on electron-phonon coupling in semiconductors: Many-body theory along thermal lines, Phys. Rev. B 93, 100301 (2016).
  * Cuoco  _et al._ [2003] M. Cuoco, C. Noce, J. Ranninger, and A. Romano, Boson-fermion model: An exact diagonalization study, Phys. Rev. B 67, 224504 (2003).
  * Jeckelmann and White [1998] E. Jeckelmann and S. R. White, Density-matrix renormalization-group study of the polaron problem in the holstein model, Phys. Rev. B 57, 6376 (1998).
  * Baiardi  _et al._ [2017] A. Baiardi, C. J. Stein, V. Barone, and M. Reiher, Vibrational density matrix renormalization group, Journal of Chemical Theory and Computation 13, 3764–3777 (2017), pMID: 28679054, https://doi.org/10.1021/acs.jctc.7b00329 .
  * White  _et al._ [2020] A. F. White, Y. Gao, A. J. Minnich, and G. K.-L. Chan, A coupled cluster framework for electrons and phonons, The Journal of Chemical Physics 153, 224112 (2020), https://doi.org/10.1063/5.0033132 .
  * Haugland  _et al._ [2020] T. S. Haugland, E. Ronca, E. F. Kjønstad, A. Rubio, and H. Koch, Coupled cluster theory for molecular polaritons: Changing ground and excited states, Phys. Rev. X 10, 041043 (2020).
  * Wang  _et al._ [2020] Y. Wang, I. Esterlis, T. Shi, J. I. Cirac, and E. Demler, Zero-temperature phases of the two-dimensional hubbard-holstein model: A non-gaussian exact diagonalization study, Phys. Rev. Research 2, 043258 (2020).
  * Costa  _et al._ [2020] N. C. Costa, K. Seki, S. Yunoki, and S. Sorella, Phase diagram of the two-dimensional hubbard-holstein model, Communications Physics 3, 80 (2020).
  * [41] Quantum monte carlo methods for fermion-boson problems, https://www.cond-mat.de/events/correl18/manuscripts/assaad.pdf, accessed: 2022-07-20.
  * Weber [2022] M. Weber, Quantum monte carlo simulation of spin-boson models using wormhole updates, Phys. Rev. B 105, 165129 (2022).
  * Bonča  _et al._ [1999] J. Bonča, S. A. Trugman, and I. Batistić, Holstein polaron, Phys. Rev. B 60, 1633 (1999).
  * Karakuzu  _et al._ [2017] S. Karakuzu, L. F. Tocchio, S. Sorella, and F. Becca, Superconductivity, charge-density waves, antiferromagnetism, and phase separation in the hubbard-holstein model, Phys. Rev. B 96, 205145 (2017).
  * Shi  _et al._ [2018] T. Shi, E. Demler, and J. Ignacio Cirac, Variational study of fermionic and bosonic systems with non-gaussian states: Theory and applications, Annals of Physics 390, 245–302 (2018).
  * Christiansen [2004a] O. Christiansen, A second quantization formulation of multimode dynamics, The Journal of Chemical Physics 120, 2140 (2004a), https://doi.org/10.1063/1.1637578 .
  * Faucheaux and Hirata [2015] J. A. Faucheaux and S. Hirata, Higher-order diagrammatic vibrational coupled-cluster theory, The Journal of Chemical Physics 143, 134105 (2015), https://aip.scitation.org/doi/pdf/10.1063/1.4931472 .
  * Christiansen [2004b] O. Christiansen, Vibrational coupled cluster theory, The Journal of Chemical Physics 120, 2149 (2004b), https://doi.org/10.1063/1.1637579 .
  * Booth  _et al._ [2009] G. H. Booth, A. J. W. Thom, and A. Alavi, Fermion monte carlo without fixed nodes: A game of life, death, and annihilation in slater determinant space, The Journal of Chemical Physics 131, 054106 (2009).
  * Cleland  _et al._ [2010a] D. Cleland, G. H. Booth, and A. Alavi, Communications: Survival of the fittest: Accelerating convergence in full configuration-interaction quantum monte carlo, The Journal of Chemical Physics 132, 041103 (2010a).
  * Shepherd  _et al._ [2012a] J. J. Shepherd, G. Booth, A. Grüneis, and A. Alavi, Full configuration interaction perspective on the homogeneous electron gas, Phys. Rev. B 85, 081103 (2012a).
  * Booth  _et al._ [2013] G. H. Booth, A. Grüneis, G. Kresse, and A. Alavi, Towards an exact description of electronic wavefunctions in real solids, Nature 493, 365 (2013).
  * Overy  _et al._ [2014a] C. Overy, G. H. Booth, N. S. Blunt, J. J. Shepherd, D. Cleland, and A. Alavi, Unbiased reduced density matrices and electronic properties from full configuration interaction quantum monte carlo, The Journal of Chemical Physics 141, 244117 (2014a).
  * Shepherd  _et al._ [2014] J. J. Shepherd, G. E. Scuseria, and J. S. Spencer, Sign problem in full configuration interaction quantum monte carlo: Linear and sublinear representation regimes for the exact wave function, Phys. Rev. B 90, 155130 (2014).
  * Blunt  _et al._ [2015a] N. S. Blunt, S. D. Smart, G. H. Booth, and A. Alavi, An excited-state approach within full configuration interaction quantum monte carlo, The Journal of Chemical Physics 143, 134117 (2015a).
  * Blunt  _et al._ [2015b] N. S. Blunt, A. Alavi, and G. H. Booth, Krylov-projected quantum monte carlo method, Phys. Rev. Lett. 115, 050603 (2015b).
  * Booth  _et al._ [2012] G. H. Booth, D. Cleland, A. Alavi, and D. P. Tew, An explicitly correlated approach to basis set incompleteness in full configuration interaction quantum monte carlo, The Journal of Chemical Physics 137, 164112 (2012).
  * Thomas  _et al._ [2015a] R. E. Thomas, D. Opalka, C. Overy, P. J. Knowles, A. Alavi, and G. H. Booth, Analytic nuclear forces and molecular properties from full configuration interaction quantum monte carlo, The Journal of Chemical Physics 143, 054108 (2015a).
  * Thomas  _et al._ [2015b] R. E. Thomas, Q. Sun, A. Alavi, and G. H. Booth, Stochastic multiconfigurational self-consistent field theory, Journal of Chemical Theory and Computation 11, 5316 (2015b).
  * Li Manni  _et al._ [2016] G. Li Manni, S. D. Smart, and A. Alavi, Combining the complete active space self-consistent field method and the full configuration interaction quantum monte carlo within a super-ci framework, with application to challenging metal-porphyrins, Journal of Chemical Theory and Computation 12, 1245 (2016).
  * Kersten  _et al._ [2016] J. A. F. Kersten, G. H. Booth, and A. Alavi, Assessment of multireference approaches to explicitly correlated full configuration interaction quantum monte carlo, The Journal of Chemical Physics 145, 054117 (2016), https://doi.org/10.1063/1.4959245 .
  * Booth and Chan [2012] G. H. Booth and G. K.-L. Chan, Communication: Excited states, dynamic correlation functions and spectral properties from full configuration interaction quantum monte carlo, The Journal of Chemical Physics 137, 191102 (2012).
  * Blunt  _et al._ [2017a] N. S. Blunt, G. H. Booth, and A. Alavi, Density matrices in full configuration interaction quantum monte carlo: Excited states, transition dipole moments, and parallel distribution, The Journal of Chemical Physics 146, 244105 (2017a).
  * Bogdanov  _et al._ [2018] N. A. Bogdanov, G. Li Manni, S. Sharma, O. Gunnarsson, and A. Alavi, New superexchange paths due to breathing-enhanced hopping in corner-sharing cuprates, Arxiv arXiv:1803.07026 (2018).
  * Blunt  _et al._ [2018] N. S. Blunt, A. Alavi, and G. H. Booth, Nonlinear biases, stochastically sampled effective hamiltonians, and spectral functions in quantum monte carlo methods, Phys. Rev. B 98, 085118 (2018).
  * Blunt [2018] N. S. Blunt, Communication: An efficient and accurate perturbative correction to initiator full configuration interaction quantum monte carlo, The Journal of Chemical Physics 148, 221101 (2018).
  * Guther  _et al._ [2018] K. Guther, W. Dobrautz, O. Gunnarsson, and A. Alavi, Time propagation and spectroscopy of fermionic systems using a stochastic technique, Phys. Rev. Lett. 121, 056401 (2018).
  * Samanta  _et al._ [2018a] P. K. Samanta, N. S. Blunt, and G. H. Booth, Response formalism within full configuration interaction quantum monte carlo: Static properties and electrical response, Journal of Chemical Theory and Computation 14, 3532 (2018a).
  * Li Manni and Alavi [2018] G. Li Manni and A. Alavi, Understanding the mechanism stabilizing intermediate spin states in fe(ii)-porphyrin, Journal of Physical Chemistry A 122, 4935 (2018).
  * Luo and Alavi [2018] H. Luo and A. Alavi, Combining the transcorrelated method with full configuration interaction quantum monte carlo: Application to the homogeneous electron gas, Journal of Chemical Theory and Computation 14, 1403 (2018), pMID: 29431996, https://doi.org/10.1021/acs.jctc.7b01257 .
  * Li Manni  _et al._ [2020] G. Li Manni, W. Dobrautz, and A. Alavi, Compression of spin-adapted multiconfigurational wave functions in exchange-coupled polynuclear spin systems, Journal of Chemical Theory and Computation 16, 2202 (2020), pMID: 32053374.
  * Dobrautz  _et al._ [2019a] W. Dobrautz, S. D. Smart, and A. Alavi, Efficient formulation of full configuration interaction quantum monte carlo in a spin eigenbasis via the graphical unitary group approach, The Journal of Chemical Physics 151, 094104 (2019a), https://doi.org/10.1063/1.5108908 .
  * Li Manni  _et al._ [2019] G. Li Manni, D. Kats, D. P. Tew, and A. Alavi, Role of valence and semicore electron correlation on spin gaps in fe(ii)-porphyrins, Journal of Chemical Theory and Computation 15, 1492 (2019).
  * Anderson  _et al._ [2020a] R. J. Anderson, T. Shiozaki, and G. H. Booth, Efficient and stochastic multireference perturbation theory for large active spaces within a full configuration interaction quantum monte carlo framework, The Journal of Chemical Physics 152, 054101 (2020a), https://doi.org/10.1063/1.5140086 .
  * Ghanem  _et al._ [2019] K. Ghanem, A. Y. Lozovoi, and A. Alavi, Unbiasing the initiator approximation in full configuration interaction quantum monte carlo, The Journal of Chemical Physics 151, 224108 (2019), https://doi.org/10.1063/1.5134006 .
  * Spencer  _et al._ [2012] J. S. Spencer, N. S. Blunt, and W. M. Foulkes, The sign problem and population dynamics in the full configuration interaction quantum monte carlo method, The Journal of Chemical Physics 136, 054110 (2012), https://doi.org/10.1063/1.3681396 .
  * Cleland  _et al._ [2010b] D. Cleland, G. H. Booth, and A. Alavi, Communications: Survival of the fittest: Accelerating convergence in full configuration-interaction quantum monte carlo, The Journal of Chemical Physics 132, 041103 (2010b), https://doi.org/10.1063/1.3302277 .
  * Cleland  _et al._ [2011] D. M. Cleland, G. H. Booth, and A. Alavi, A study of electron affinities using the initiator approach to full configuration interaction quantum monte carlo, The Journal of Chemical Physics 134, 024112 (2011), https://doi.org/10.1063/1.3525712 .
  * Dobrautz  _et al._ [2022] W. Dobrautz, V. M. Katukuri, N. A. Bogdanov, D. Kats, G. Li Manni, and A. Alavi, Combined unitary and symmetric group approach applied to low-dimensional heisenberg spin systems, Phys. Rev. B 105, 195123 (2022).
  * Dobrautz  _et al._ [2019b] W. Dobrautz, H. Luo, and A. Alavi, Compact numerical solutions to the two-dimensional repulsive hubbard model obtained via nonunitary similarity transformations, Phys. Rev. B 99, 075119 (2019b).
  * Schwarz  _et al._ [2015] L. R. Schwarz, G. H. Booth, and A. Alavi, Insights into the structure of many-electron wave functions of mott-insulating antiferromagnets: The three-band hubbard model in full configuration interaction quantum monte carlo, Phys. Rev. B 91, 045139 (2015).
  * Shepherd  _et al._ [2012b] J. J. Shepherd, G. Booth, A. Grüneis, and A. Alavi, Full configuration interaction perspective on the homogeneous electron gas, Phys. Rev. B 85, 081103 (2012b).
  * Brand  _et al._ [2022] J. Brand, M. Yang, and E. Pahl, Stochastic differential equation approach to understanding the population control bias in full configuration interaction quantum monte carlo, Phys. Rev. B 105, 235144 (2022).
  * Yang  _et al._ [2022] M. Yang, M. Čufar, E. Pahl, and J. Brand, Polaron-depleton transition in the yrast excitations of a one-dimensional bose gas with a mobile impurity, Condensed Matter 7, 10.3390/condmat7010015 (2022).
  * Anderson [2020] R. J. Anderson, M7, https://github.com/robert-anderson/M7 (2020).
  * von der Linden  _et al._ [1995b] W. von der Linden, E. Berger, and P. Valášek, The hubbard-holstein model, Journal of Low Temperature Physics 99, 517–525 (1995b).
  * Berger  _et al._ [1995] E. Berger, P. Valášek, and W. von der Linden, Two-dimensional hubbard-holstein model, Phys. Rev. B 52, 4806 (1995).
  * Hirsch [1985] J. E. Hirsch, Phase diagram of the one-dimensional molecular-crystal model with coulomb interactions: Half-filled-band sector, Phys. Rev. B 31, 6022 (1985).
  * Zheng  _et al._ [1990] H. Zheng, D. Feinberg, and M. Avignon, Quantum lattice fluctuations in the one-dimensional peierls-hubbard model, Phys. Rev. B 41, 11557 (1990).
  * Yonemitsu and Imada [1996] K. Yonemitsu and M. Imada, Spin-gap phase in nearly half-filled one-dimensional conductors coupled with phonons, Phys. Rev. B 54, 2410 (1996).
  * Pao and Schüttler [1998] C.-H. Pao and H.-B. Schüttler, Superconducting instability in the holstein-hubbard model: A numerical renormalization-group study, Phys. Rev. B 57, 5051 (1998).
  * Clay and Hardikar [2005] R. T. Clay and R. P. Hardikar, Intermediate phase of the one dimensional half-filled hubbard-holstein model, Phys. Rev. Lett. 95, 096401 (2005).
  * Koller  _et al._ [2005] W. Koller, D. Meyer, A. Hewson, and Y. Ōno, Phase diagram and dynamic response functions of the holstein–hubbard model, Physica B: Condensed Matter 359-361, 795–797 (2005), proceedings of the International Conference on Strongly Correlated Electron Systems.
  * Tezuka  _et al._ [2007] M. Tezuka, R. Arita, and H. Aoki, Phase diagram for the one-dimensional hubbard-holstein model: A density-matrix renormalization group study, Phys. Rev. B 76, 155114 (2007).
  * Murakami  _et al._ [2013] Y. Murakami, P. Werner, N. Tsuji, and H. Aoki, Ordered phases in the holstein-hubbard model: Interplay of strong coulomb interaction and electron-phonon coupling, Phys. Rev. B 88, 125126 (2013).
  * 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, Density-matrix embedding theory study of the one-dimensional hubbard–holstein model, Journal of Chemical Theory and Computation 15, 2221 (2019), pMID: 30807149, https://doi.org/10.1021/acs.jctc.8b01116 .
  * Sandhoefer and Chan [2016] B. Sandhoefer and G. K.-L. Chan, Density matrix embedding theory for interacting electron-phonon systems, Phys. Rev. B 94, 085115 (2016).
  * Booth  _et al._ [2014a] G. H. Booth, S. D. Smart, and A. Alavi, Linear-scaling and parallelisable algorithms for stochastic quantum chemistry, Mol. Phys. 112, 1855 (2014a), https://doi.org/10.1080/00268976.2013.877165 .
  * Neufeld and Thom [2019] V. A. Neufeld and A. J. W. Thom, Exciting determinants in quantum monte carlo: Loading the dice with fast, low-memory weights, Journal of Chemical Theory and Computation 15, 127 (2019), https://doi.org/10.1021/acs.jctc.8b00844 .
  * Booth  _et al._ [2014b] G. H. Booth, S. D. Smart, and A. Alavi, Linear-scaling and parallelisable algorithms for stochastic quantum chemistry, Molecular Physics 112, 1855 (2014b), https://doi.org/10.1080/00268976.2013.877165 .
  * Holmes  _et al._ [2016] A. A. Holmes, H. J. Changlani, and C. J. Umrigar, Efficient heat-bath sampling in fock space, Journal of Chemical Theory and Computation 12, 1561–1571 (2016), pMID: 26959242, https://doi.org/10.1021/acs.jctc.5b01170 .
  * Guther  _et al._ [2020] K. Guther, R. J. Anderson, N. S. Blunt, N. A. Bogdanov, D. Cleland, N. Dattani, W. Dobrautz, K. Ghanem, P. Jeszenszki, N. Liebermann, G. L. Manni, A. Y. Lozovoi, H. Luo, D. Ma, F. Merz, C. Overy, M. Rampp, P. K. Samanta, L. R. Schwarz, J. J. Shepherd, S. D. Smart, E. Vitale, O. Weser, G. H. Booth, and A. Alavi, Neci: N-electron configuration interaction with an emphasis on state-of-the-art stochastic methods, The Journal of Chemical Physics 153, 034107 (2020), https://doi.org/10.1063/5.0005754 .
  * Ghanem  _et al._ [2021] K. Ghanem, N. Liebermann, and A. Alavi, Population control bias and importance sampling in full configuration interaction quantum monte carlo, Phys. Rev. B 103, 155135 (2021).
  * Nightingale and Blöte [1988] M. P. Nightingale and H. W. J. Blöte, Monte carlo calculation of free energy, critical point, and surface critical behavior of three-dimensional heisenberg ferromagnets, Phys. Rev. Lett. 60, 1562 (1988).
  * Cerf and Martin [1995] N. Cerf and O. C. Martin, Finite population-size effects in projection monte carlo methods, Phys. Rev. E 51, 3679 (1995).
  * Umrigar  _et al._ [1993] C. J. Umrigar, M. P. Nightingale, and K. J. Runge, A diffusion monte carlo algorithm with very small time‐step errors, The Journal of Chemical Physics 99, 2865 (1993), https://doi.org/10.1063/1.465195 .
  * Ceperley and Alder [1986] D. Ceperley and B. Alder, Quantum monte carlo, Science 231, 555 (1986), https://www.science.org/doi/pdf/10.1126/science.231.4738.555 .
  * Kolodrubetz and Clark [2012] M. Kolodrubetz and B. K. Clark, Partial node configuration-interaction monte carlo as applied to the fermi polaron, Phys. Rev. B 86, 075109 (2012).
  * Flyvbjerg and Petersen [1989] H. Flyvbjerg and H. G. Petersen, Error estimates on averages of correlated data, The Journal of Chemical Physics 91, 461 (1989), https://doi.org/10.1063/1.457480 .
  * Blunt  _et al._ [2017b] N. S. Blunt, G. H. Booth, and A. Alavi, Density matrices in full configuration interaction quantum monte carlo: Excited states, transition dipole moments, and parallel distribution, The Journal of Chemical Physics 146, 244105 (2017b), https://doi.org/10.1063/1.4986963 .
  * Overy  _et al._ [2014b] C. Overy, G. H. Booth, N. S. Blunt, J. J. Shepherd, D. Cleland, and A. Alavi, Unbiased reduced density matrices and electronic properties from full configuration interaction quantum monte carlo, The Journal of Chemical Physics 141, 244117 (2014b), https://doi.org/10.1063/1.4904313 .
  * Thomas  _et al._ [2015c] R. E. Thomas, D. Opalka, C. Overy, P. J. Knowles, A. Alavi, and G. H. Booth, Analytic nuclear forces and molecular properties from full configuration interaction quantum monte carlo, The Journal of Chemical Physics 143, 054108 (2015c), https://doi.org/10.1063/1.4927594 .
  * Samanta  _et al._ [2018b] P. K. Samanta, N. S. Blunt, and G. H. Booth, Response formalism within full configuration interaction quantum monte carlo: Static properties and electrical response, Journal of Chemical Theory and Computation 14, 3532 (2018b), pMID: 29897746, https://doi.org/10.1021/acs.jctc.8b00454 .
  * Thomas  _et al._ [2015d] R. E. Thomas, Q. Sun, A. Alavi, and G. H. Booth, Stochastic multiconfigurational self-consistent field theory, Journal of Chemical Theory and Computation 11, 5316–5325 (2015d), https://doi.org/10.1021/acs.jctc.5b00917 .
  * Li Manni  _et al._ [2016] G. Li Manni, S. D. Smart, and A. Alavi, Combining the complete active space self-consistent field method and the full configuration interaction quantum monte carlo within a super-ci framework, with application to challenging metal-porphyrins, Journal of Chemical Theory and Computation 12, 1245–1258 (2016), pMID: 26808894, https://doi.org/10.1021/acs.jctc.5b01190 .
  * Weser  _et al._ [2022] O. Weser, K. Guther, K. Ghanem, and G. Li Manni, Stochastic generalized active space self-consistent field: Theory and application, Journal of Chemical Theory and Computation 18, 251–272 (2022), pMID: 34898215, https://doi.org/10.1021/acs.jctc.1c00936 .
  * Anderson  _et al._ [2020b] R. J. Anderson, T. Shiozaki, and G. H. Booth, Efficient and stochastic multireference perturbation theory for large active spaces within a full configuration interaction quantum monte carlo framework, The Journal of Chemical Physics 152, 054101 (2020b), https://doi.org/10.1063/1.5140086 .
  * Halson  _et al._ [2020] J. J. Halson, R. J. Anderson, and G. H. Booth, Improved stochastic multireference perturbation theory for correlated systems with large active spaces, Molecular Physics 118, e1802072 (2020), https://doi.org/10.1080/00268976.2020.1802072 .
  * Blunt  _et al._ [2015c] N. S. Blunt, S. D. Smart, J. A. F. Kersten, J. S. Spencer, G. H. Booth, and A. Alavi, Semi-stochastic full configuration interaction quantum monte carlo: Developments and application, The Journal of Chemical Physics 142, 184107 (2015c), https://doi.org/10.1063/1.4920975 .
  * Petruzielo  _et al._ [2012] F. R. Petruzielo, A. A. Holmes, H. J. Changlani, M. P. Nightingale, and C. J. Umrigar, Semistochastic projector monte carlo method, Phys. Rev. Lett. 109, 230201 (2012).
  * Sengupta  _et al._ [2002] P. Sengupta, A. W. Sandvik, and D. K. Campbell, Bond-order-wave phase and quantum phase transitions in the one-dimensional extended hubbard model, Phys. Rev. B 65, 155113 (2002).
  * Tsuchiizu and Furusaki [2002] M. Tsuchiizu and A. Furusaki, Phase diagram of the one-dimensional extended hubbard model at half filling, Phys. Rev. Lett. 88, 056402 (2002).
  * Jeckelmann [2002] E. Jeckelmann, Ground-state phase diagram of a half-filled one-dimensional extended hubbard model, Phys. Rev. Lett. 89, 236401 (2002).
  * Ejima and Nishimoto [2007] S. Ejima and S. Nishimoto, Phase diagram of the one-dimensional half-filled extended hubbard model, Phys. Rev. Lett. 99, 216403 (2007).

##  Appendix A Zero Phonon Mode Removal

When studying the Hubbard–Holstein model with periodic boundary conditions, it
is routine to bring about a reduction in the coupling between the electron and
boson sector by applying the zero phonon mode removal (ZPMR) transformation to
the bosonic second quantised operators. This can be written as

$\hat{a}_{m}\rightarrow\hat{a}_{m}-\frac{g}{\omega_{0}}\langle n\rangle$ (9)

where $\langle n\rangle=N_{\mathrm{elec}}/L$ stands for the average electronic
density per site, Under this transformation, the electron-boson coupling
becomes

$g\hat{n}_{m}\left(\hat{a}_{m}^{\dagger}+\hat{a}_{m}-\frac{2g}{\omega_{0}}\langle
n\rangle\right).$ (10)

Since the fermion number is fixed, the average fermion occupation per site is
known by symmetry, and the resultant on-site one-electron contribution amounts
to an energy shift of $-2g^{2}N_{\mathrm{elec}}{}^{2}/(\omega_{0}L)$.
Likewise, the boson energy contributions transform to

$\omega_{0}\left(\hat{a}_{m}^{\dagger}\hat{a}_{m}-\frac{g}{\omega_{0}}\langle
n\rangle(\hat{a}_{m}^{\dagger}+\hat{a}_{m})+\frac{g^{2}}{\omega_{0}^{2}}\langle
n\rangle^{2}\right)$ (11)

from which the second term shows that there are now single boson
(de-)excitations that are uncoupled to the electronic modes, and an additional
constant. After summation over sites, this gives another constant energy shift
of $g^{2}N_{\mathrm{elec}}{}^{2}/(\omega_{0}L)$, resulting in an overall
energy shift of $-g^{2}N_{\mathrm{elec}}{}^{2}/(\omega_{0}L)$ to the ZPMR
Hamiltonian spectrum relative to the untransformed definition.

This transformation redefines the vacuum state in the boson modes in such a
way that the boson occupation-truncated wave function converges faster to the
untruncated wave function with respect to $n_{\mathrm{boson}}$, making larger
systems tractable by exact eigensolvers provided that the coupling to the
bosons is weak. Physically, this amounts to decoupling the bosons from the
static charge density of the system, ensuring only coupling to the
instantaneous quantum density fluctuations (higher moments of the density
distribution) remain. It is tempting to bring the ZPMR transformation into
effect in the hope of improving the performance of FCIQMC by condensing the
walker distribution around configurations with lower-occupied modes. Indeed,
this is likely to be advantageous for general hamiltonians. However, for the
Hubbard–Holstein model, it must be noted that the density-uncoupled
(de-)excitations between boson number sectors would always carry the opposite
sign to the density-coupled (de-)excitations. This introduces a sign problem,
which would negate the advantages in wave function compactification achieved
by the transformation. Therefore, use of the ZPMR transformation in the
present work is restricted to obtaining faster converging FCI energies with
respect to boson occupation cutoff for small systems, to allow comparison to
FCIQMC results.