# Diffusion Monte Carlo using domains in configuration space

Roland Assaraf  Laboratoire de Chimie Théorique, Sorbonne-Université, Paris,
France  Emmanuel Giner  Laboratoire de Chimie Théorique, Sorbonne-Université,
Paris, France  Vijay Gopal Chilkuri  Laboratoire de Chimie et Physique
Quantiques (UMR 5626), Université de Toulouse, CNRS, UPS, France  Pierre-
François Loos  Laboratoire de Chimie et Physique Quantiques (UMR 5626),
Université de Toulouse, CNRS, UPS, France  Anthony Scemama  Laboratoire de
Chimie et Physique Quantiques (UMR 5626), Université de Toulouse, CNRS, UPS,
France  Michel Caffarel  Corresponding author<EMAIL_ADDRESS>Laboratoire de Chimie et Physique Quantiques (UMR 5626), Université de
Toulouse, CNRS, UPS, France

###### Abstract

The sampling of the configuration space in diffusion Monte Carlo (DMC) is done
using walkers moving randomly. In a previous work on the Hubbard model
[Assaraf et al. Phys. Rev. B 60, 2299 (1999)], it was shown that the
probability for a walker to stay a certain amount of time in the same state
obeys a Poisson law and that the on-state dynamics can be integrated out
exactly, leading to an effective dynamics connecting only different states.
Here, we extend this idea to the general case of a walker trapped within
domains of arbitrary shape and size. The equations of the resulting effective
stochastic dynamics are derived. The larger the average (trapping) time spent
by the walker within the domains, the greater the reduction in statistical
fluctuations. A numerical application to the Hubbard model is presented.
Although this work presents the method for finite linear spaces, it can be
generalized without fundamental difficulties to continuous configuration
spaces.

##  I Introduction

Diffusion Monte Carlo (DMC) is a class of stochastic methods for evaluating
the ground-state properties of quantum systems. They have been extensively
used in virtually all domains of physics and chemistry where the many-body
quantum problem plays a central role (condensed-matter physics,Foulkes  _et
al._ (2001); Kolorenc and Mitas (2011) quantum liquids,Holzmann  _et al._
(2006) nuclear physics,Carlson  _et al._ (2015); Carlson (2007) theoretical
chemistry,Austin  _et al._ (2012) etc). DMC can be used either for systems
defined in a continuous configuration space (typically, a set of particles
moving in space) for which the Hamiltonian operator is defined in a Hilbert
space of infinite dimension or systems defined in a discrete configuration
space where the Hamiltonian reduces to a matrix. Here, we shall consider only
the discrete case, that is, the general problem of calculating the lowest
eigenvalue and/or eigenstate of a (very large) matrix. The generalization to
continuous configuration spaces presents no fundamental difficulty.

In essence, DMC is based on stochastic power methods, a family of well
established numerical approaches able to extract the largest or smallest
eigenvalues of a matrix (see, e.g., Ref. Golub and Van Loan, 2012). These
approaches are particularly simple as they merely consist in applying a given
matrix (or some simple function of it) as many times as required on some
arbitrary vector belonging to the linear space. Thus, the basic step of the
corresponding algorithm essentially reduces to successive matrix-vector
multiplications. In practice, power methods are employed under more
sophisticated implementations, such as, e.g., the Lanczòs algorithm (based on
Krylov subspaces) Golub and Van Loan (2012) or Davidson’s method where a
diagonal preconditioning is performed. Davidson (1975) When the size of the
matrix is too large, matrix-vector multiplications become unfeasible and
probabilistic techniques to sample only the most important contributions of
the matrix-vector product are required. This is the basic idea of DMC. There
exist several variants of DMC known under various names: pure DMC, Caffarel
and Claverie (1988) DMC with branching, Reynolds  _et al._ (1982) reptation
Monte Carlo, Baroni and Moroni (1999) stochastic reconfiguration Monte Carlo,
Sorella (1998); Assaraf  _et al._ (2000) etc. Here, we shall place ourselves
within the framework of pure DMC whose mathematical simplicity is particularly
appealing when developing new ideas. However, all the ideas presented in this
work can be adapted without too much difficulty to the other variants, so the
denomination DMC must ultimately be understood here as a generic name for this
broad class of methods.

Without entering into the mathematical details (which are presented below),
the main ingredient of DMC in order to perform the matrix-vector
multiplications probabilistically is to introduce a stochastic matrix (or
transition probability matrix) that generates stepwise a series of states over
which statistical averages are evaluated. The critical aspect of any Monte
Carlo scheme is the amount of computational effort required to reach a given
statistical error. Two important avenues to decrease the error are the use of
variance reduction techniques (for example, by introducing improved estimators
Assaraf and Caffarel (1999)) or to improve the quality of the sampling
(minimization of the correlation time between states). Another possibility, at
the heart of the present work, is to integrate out exactly some parts of the
dynamics, thus reducing the number of degrees of freedom and, hence, the
amount of statistical fluctuations.

In previous works,Assaraf  _et al._ (1999); Caffarel and Assaraf (2000) it has
been shown that the probability for a walker to stay a certain amount of time
in the same state obeys a Poisson law and that the on-state dynamics can be
integrated out to generate an effective dynamics connecting only different
states with some renormalized estimators for the properties. Numerical
applications have shown that the statistical errors can be very significantly
decreased. Here, we extend this idea to the general case where a walker
remains a certain amount of time in a finite domain no longer restricted to a
single state. It is shown how to define an effective stochastic dynamics
describing walkers moving from one domain to another. The equations of the
effective dynamics are derived and a numerical application for a model (one-
dimensional) problem is presented. In particular, it shows that the
statistical convergence of the energy can be greatly enhanced when domains
associated with large average trapping times are considered.

It should be noted that the use of domains in quantum Monte Carlo is not new.
Domains have been introduced within the context of the Green’s function Monte
Carlo (GFMC) method pioneered by KalosKalos (1962, 1970) and later developped
and applied by Kalos and others.Kalos  _et al._ (1974); Ceperley and Kalos
(1979); Ceperley (1983); Moskowitz and Schmidt (1986) In GFMC, an approximate
Green’s function that can be sampled is needed to propagate stochastically the
wavefunction. In the so-called Domain’s GFMC version of GFMC introduced in
Ref.Kalos, 1970 and Kalos  _et al._ , 1974 the sampling is realized by using
the restriction of the Green’s function to a small domain consisting of the
Cartesian product of small spheres around each particle, the potential being
considered as constant within the domain. Fundamentally, the method presented
in this work is closely related to the Domain’s GFMC, although the way we
present the formalism in terms of walkers trapped within domains and derive
the equations may appear as different. However, a key difference here is that
we show how to use domains of arbitrary size, a point which can greatly
enhance the efficiency of the simulations when domains are suitably chosen, as
illustrated in our numerical application.

Finally, from a general perspective, it is interesting to emphasize that the
present method illustrates how suitable combinations of stochastic and
deterministic techniques lead to a more efficient and valuable method. In
recent years, a number of works have exploited this idea and proposed hybrid
stochastic/deterministic schemes. Let us mention, for example, the semi-
stochastic approach of Petruzielo et al., Petruzielo  _et al._ (2012) two
different hybrid algorithms for evaluating the second-order perturbation
energy in selected configuration interaction methods, Garniron  _et al._
(2017); Sharma  _et al._ (2017) the approach of Willow et al. for computing
stochastically second-order many-body perturbation energies, Willow  _et al._
(2012) or the zero variance Monte Carlo scheme for evaluating two-electron
integrals in quantum chemistry. Caffarel (2019)

The paper is organized as follows. Section II presents the basic equations and
notations of DMC. First, the path integral representation of the Green’s
function is given in Subsec. II.1. Second, the probabilistic framework
allowing the Monte Carlo calculation of the Green’s function is presented in
Subsec. II.2. Section III is devoted to the use of domains in DMC. We recall
in Subsec. III.1 the case of a domain consisting of a single state. The
general case is then treated in Subsec. III.2. In Subsec. III.3, both the
time- and energy-dependent Green’s function using domains are derived. Section
IV presents the application of the approach to the one-dimensional Hubbard
model. Finally, in Sec.V, some conclusions and perspectives are given. Atomic
units are used throughout.

##  II Diffusion Monte Carlo

###  II.1 Path-integral representation

As previously mentioned, DMC is a stochastic implementation of the power
method defined by the following operator:

$T=\mathds{1}-\tau(H-E\mathds{1}),$ (1)

where $\mathds{1}$ is the identity operator, $\tau$ a small positive parameter
playing the role of a time step, $E$ some arbitrary reference energy, and $H$
the Hamiltonian operator. For any initial vector $\ket{\Psi_{0}}$ provided
that $\innerproduct{\Phi_{0}}{\Psi_{0}}\neq 0$ and for $\tau$ sufficiently
small, we have

$\lim_{N\to\infty}T^{N}\ket{\Psi_{0}}=\ket{\Phi_{0}},$ (2)

where $\ket{\Phi_{0}}$ is the ground-state wave function, i.e.,
$H\ket{\Phi_{0}}=E_{0}\ket{\Phi_{0}}$. The equality in Eq. (2) holds up to a
global phase factor playing no role in physical quantum averages. At large but
finite $N$, the vector $T^{N}\ket{\Psi_{0}}$ differs from $\ket{\Phi_{0}}$
only by an exponentially small correction, making it straightforward to
extrapolate the finite-$N$ results to $N\to\infty$.

Likewise, ground-state properties may be obtained at large $N$. For example,
in the important case of the energy, one can project out the vector
$T^{N}\ket{\Psi_{0}}$ on some approximate vector, $\ket{\Psi_{\text{T}}}$, as
follows

$E_{0}=\lim_{N\to\infty}\frac{\matrixelement{\Psi_{0}}{T^{N}}{H\Psi_{T}}}{\matrixelement{\Psi_{0}}{T^{N}}{\Psi_{T}}}.$
(3)

$\ket{\Psi_{\text{T}}}$ is known as the trial wave vector (function), and is
chosen as an approximation of the true ground-state vector.

To proceed further we introduce the time-dependent Green’s matrix $G^{(N)}$
defined as

$G^{(N)}_{ij}=\matrixelement{i}{T^{N}}{j}.$ (4)

where $\ket{i}$ and $\ket{j}$ are basis vectors. The denomination “time-
dependent Green’s matrix” is used here since $G$ may be viewed as a short-time
approximation of the (time-imaginary) evolution operator $e^{-N\tau H}$, which
is usually referred to as the imaginary-time dependent Green’s function.

Introducing the set of $N-1$ intermediate states, $\\{\ket{i_{k}}\\}_{1\leq
k\leq N-1}$, between each $T$ in $T^{N}$, $G^{(N)}$ can be written in the
following expanded form

$G^{(N)}_{i_{0}i_{N}}=\sum_{i_{1}}\sum_{i_{2}}\cdots\sum_{i_{N-1}}\prod_{k=0}^{N-1}T_{i_{k}i_{k+1}},$
(5)

where $T_{ij}=\matrixelement{i}{T}{j}$. Here, each index $i_{k}$ runs over all
basis vectors.

In quantum physics, Eq. (5) is referred to as the path-integral representation
of the Green’s matrix (or function). The series of states
$\ket{i_{0}},\ldots,\ket{i_{N}}$ is interpreted as a “path” in the Hilbert
space starting at vector $\ket{i_{0}}$ and ending at vector $\ket{i_{N}}$,
where $k$ plays the role of a time index. Each path is associated with a
weight $\prod_{k=0}^{N-1}T_{i_{k}i_{k+1}}$ and the path integral expression of
$G$ can be recast in the more suggestive form as follows:

$G^{(N)}_{i_{0}i_{N}}=\sum_{\text{all paths
$\ket{i_{1}},\ldots,\ket{i_{N-1}}$}}\prod_{k=0}^{N-1}T_{i_{k}i_{k+1}}.$ (6)

This expression allows a simple and vivid interpretation of the solution. In
the limit $N\to\infty$, the $i_{N}$th component of the ground-state wave
function (obtained as $\lim_{N\to\infty}G^{(N)}_{i_{0}i_{N}})$ is the weighted
sum over all possible paths arriving at vector $\ket{i_{N}}$. This result is
independent of the initial vector $\ket{i_{0}}$, apart from some irrelevant
global phase factor. We illustrate this fundamental property pictorially in
Fig. 1. When the size of the linear space is small, the explicit calculation
of the full sums involving $M^{N}$ terms (where $M$ is the size of the Hilbert
space) can be performed. In such a case, we are in the realm of what one would
call “deterministic” power methods, such as the Lanczòs or Davidson
approaches. If not, probabilistic techniques for generating only the paths
contributing significantly to the sums are to be used. This is the central
theme of the present work.

Figure 1:  Path integral representation of the exact coefficient
$c_{i}=\innerproduct{i}{\Phi_{0}}$ of the ground-state wave function
$\ket{\Phi_{0}}$ obtained as an infinite sum of paths starting from
$\ket{i_{0}}$ and ending at $\ket{i}$ [see Eq.(6)]. Each path carries a weight
$\prod_{k}T_{i_{k}i_{k+1}}$ computed along it. The result is independent of
the choice of the initial state $\ket{i_{0}}$, provided that
$\innerproduct{i_{0}}{\Phi_{0}}\neq 0$. Here, only four paths of infinite
length have been represented.

###  II.2 Probabilistic framework

To derive a probabilistic expression for the Green’s matrix, we introduce a
guiding wave function $\ket{\Psi^{+}}$ having strictly positive components,
i.e., $\Psi^{+}_{i}>0$, in order to perform a similarity transformation of the
operators $G^{(N)}$ and $T$,

$\displaystyle\bar{T}_{ij}$
$\displaystyle=\frac{\Psi^{+}_{j}}{\Psi^{+}_{i}}T_{ij},$
$\displaystyle\bar{G}^{(N)}_{ij}$
$\displaystyle=\frac{\Psi^{+}_{j}}{\Psi^{+}_{i}}G^{(N)}_{ij}.$ (7)

Note that, thanks to the properties of similarity transformations, the path
integral expression relating $G^{(N)}$ and $T$ [see Eq. (6)] remains unchanged
for $\bar{G}^{(N)}$ and $\bar{T}$.

Now, the key idea to take advantage of probabilistic techniques is to rewrite
the matrix elements of $\bar{T}$ as those of a stochastic matrix multiplied by
some residual weights (here, not necessarily positive), namely

$\bar{T}_{ij}=p_{i\to j}w_{ij}.$ (8)

Here, we recall that a stochastic matrix is defined as a matrix with positive
entries that obeys

$\sum_{j}p_{i\to j}=1.$ (9)

Using this representation for $\bar{T}_{ij}$ the similarity-transformed
Green’s matrix components can be rewritten as

$\bar{G}^{(N)}_{i_{0}i_{N}}=\sum_{i_{1},\ldots,i_{N-1}}\quantity(\prod_{k=0}^{N-1}p_{i_{k}\to
i_{k+1}})\prod_{k=0}^{N-1}w_{i_{k}i_{k+1}},$ (10)

which is amenable to Monte Carlo calculations by generating paths using the
transition probability matrix $p_{i\to j}$.

Let us illustrate this in the case of the energy as given by Eq. (3). Taking
$\ket{\Psi_{0}}=\ket{i_{0}}$ as initial state, we have

$E_{0}=\lim_{N\to\infty}\frac{\sum_{i_{N}}G^{(N)}_{i_{0}i_{N}}{(H\Psi_{\text{T}})}_{i_{N}}}{\sum_{i_{N}}{G}^{(N)}_{i_{0}i_{N}}{\Psi_{\text{T}}}_{i_{N}}}.$
(11)

which can be rewritten probabilistically as

$E_{0}=\lim_{N\to\infty}\frac{\expectationvalue{\prod_{k=0}^{N-1}w_{i_{k}i_{k+1}}\frac{{(H\Psi_{\text{T}})}_{i_{N}}}{\Psi^{+}_{i_{N}}}}}{\expectationvalue{\prod_{k=0}^{N-1}w_{i_{k}i_{k+1}}\frac{{\Psi_{\text{T}}}_{i_{N}}}{\Psi^{+}_{i_{N}}}}},$
(12)

where $\expectationvalue{...}$ is the probabilistic average defined over the
set of paths $\ket{i_{1}},\ldots,\ket{i_{N}}$ occurring with probability

$\text{Prob}_{i_{0}}(i_{1},\ldots,i_{N})=\prod_{k=0}^{N-1}p_{i_{k}\to
i_{k+1}}.$ (13)

Using Eq. (9) and the fact that $p_{i\to j}\geq 0$, one can easily verify that
$\text{Prob}_{i_{0}}$ is positive and obeys

$\sum_{i_{1},\ldots,i_{N}}\text{Prob}_{i_{0}}(i_{1},\ldots,i_{N})=1,$ (14)

as it should.

The rewriting of $\bar{T}_{ij}$ as a product of a stochastic matrix times some
general real weight does not introduce any constraint on the choice of the
stochastic matrix, so that, in theory, any stochastic matrix could be used.
However, in practice, it is highly desirable that the magnitude of the
fluctuations of the weight during the Monte Carlo simulation be as small as
possible. A natural solution is to choose a stochastic matrix equal to a good
approximation of $\bar{T}_{ij}$. This is done as follows.

Let us introduce the following operator

$T^{+}=\mathds{1}-\tau\quantity(H^{+}-E_{\text{L}}^{+}\mathds{1}),$ (15)

where $H^{+}$ is the matrix obtained from $H$ by imposing the off-diagonal
elements to be negative

$H^{+}_{ij}=\begin{cases}\phantom{-}H_{ij},&\text{if $i=j$},\\\
-\absolutevalue{H_{ij}},&\text{if $i\neq j$}.\end{cases}$ (16)

Here, $E_{\text{L}}^{+}\mathds{1}$ is the diagonal matrix whose diagonal
elements are defined as

$(E_{\text{L}}^{+})_{i}=\frac{\sum_{j}H^{+}_{ij}\Psi^{+}_{j}}{\Psi^{+}_{i}}.$
(17)

The vector $E_{\text{L}}^{+}$ is known as the local energy vector associated
with $\Psi^{+}$. By construction, the operator
$H^{+}-E_{\text{L}}^{+}\mathds{1}$ in the definition of $T^{+}$ [see Eq. (15)]
has been chosen to admit $\ket{\Psi^{+}}$ as a ground-state wave function with
zero eigenvalue, i.e., $\quantity(H^{+}-E_{L}^{+}\mathds{1})\ket{\Psi^{+}}=0$,
leading to the relation

$T^{+}\ket{\Psi^{+}}=\ket{\Psi^{+}}.$ (18)

We are now in the position to define the stochastic matrix as

$p_{i\to
j}=\frac{\Psi^{+}_{j}}{\Psi^{+}_{i}}T^{+}_{ij}=\begin{cases}1-\tau\quantity[H^{+}_{ii}-(E_{\text{L}}^{+})_{i}],&\mbox{if\quad}i=j,\\\
\tau\frac{\Psi^{+}_{j}}{\Psi^{+}_{i}}\absolutevalue{H_{ij}}\geq
0,&\mbox{if\quad}i\neq j.\end{cases}$ (19)

As readily seen in Eq. (19), the off-diagonal terms of the stochastic matrix
are positive, while the diagonal terms can be made positive if $\tau$ is
chosen sufficiently small via the condition

$\tau\leq\frac{1}{\max_{i}\absolutevalue{H^{+}_{ii}-(E_{\text{L}}^{+})_{i}}}.$
(20)

The sum-over-states condition [see Eq. (9)],

$\sum_{j}p_{i\to
j}=\frac{\matrixelement{i}{T^{+}}{\Psi^{+}}}{\Psi^{+}_{i}}=1,$ (21)

follows from the fact that $\ket{\Psi^{+}}$ is eigenvector of $T^{+}$, as
evidenced by Eq. (18). This ensures that $p_{i\to j}$ is indeed a stochastic
matrix.

At first sight, the condition defining the maximum value of $\tau$ [see Eq.
(20)] may appear rather tight since, for very large matrices, it may impose an
extremely small value of the time step. However, in practice, during the
simulation only a (tiny) fraction of the linear space is sampled, and the
maximum absolute value of $H^{+}_{ii}-(E_{\text{L}}^{+})_{i}$ for the sampled
states turns out to be not too large. Hence, reasonable values of $\tau$ can
be selected without violating the positivity of the transition probability
matrix. Note that one can eschew this condition via a simple generalization of
the transition probability matrix:

$p_{i\to
j}=\frac{\frac{\Psi^{+}_{j}}{\Psi^{+}_{i}}\absolutevalue{\matrixelement{i}{T^{+}}{j}}}{\sum_{j}\frac{\Psi^{+}_{j}}{\Psi^{+}_{i}}\absolutevalue{\matrixelement{i}{T^{+}}{j}}}=\frac{\Psi^{+}_{j}\absolutevalue*{T^{+}_{ij}}}{\sum_{j}\Psi^{+}_{j}\absolutevalue*{T^{+}_{ij}}}.$
(22)

This new transition probability matrix with positive entries reduces to Eq.
(19) when $T^{+}_{ij}$ is positive as $\sum_{j}\Psi^{+}_{j}T^{+}_{ij}=1$.

Now, we need to make the connection between the transition probability matrix,
$p_{i\to j}$, defined from the approximate Hamiltonian $H^{+}$ via $T^{+}$ and
the operator $T$ associated with the exact Hamiltonian $H$. This is done
thanks to Eq. (8) that connects $p_{i\to j}$ and $T_{ij}$ through the weight

$w_{ij}=\frac{T_{ij}}{T^{+}_{ij}},$ (23)

derived from Eqs. (7) and (19).

To calculate the probabilistic averages, an artificial (mathematical)
“particle” called walker (or psi-particle) is introduced. During the Monte
Carlo simulation, the walker moves in configuration space by drawing new
states with probability $p_{i_{k}\to i_{k+1}}$, thus realizing the path of
probability $\text{Prob}_{i_{0}}$. Note that, instead of using a single
walker, it is common to introduce a population of independent walkers and to
calculate the averages over this population. In addition, thanks to the
ergodicity property of the stochastic matrix (see, for example, Ref. Caffarel
and Claverie, 1988), a unique path of infinite length from which sub-paths of
length $N$ can be extracted may also be used. We shall not insist here on
these practical details that are discussed, for example, in Refs. Foulkes  _et
al._ , 2001; Kolorenc and Mitas, 2011.

##  III Diffusion Monte Carlo with domains

###  III.1 Single-state domains

During the simulation, walkers move from state to state with the possibility
of being trapped a certain number of times on the same state before exiting to
a different state. This fact can be exploited in order to integrate out some
part of the dynamics, thus leading to a reduction of the statistical
fluctuations. This idea was proposed some time ago and applied to the SU($N$)
one-dimensional Hubbard model.Assaraf and Caffarel (1999); Assaraf  _et al._
(1999); Caffarel and Assaraf (2000)

Considering a given state $\ket{i}$, the probability that a walker remains
exactly $n$ times in $\ket{i}$ (with $n\geq 1$) and then exits to a different
state $j$ (with $j\neq i$) is

$\mathcal{P}_{i\to j}(n)=\quantity(p_{i\to i})^{n-1}p_{i\to j}.$ (24)

Using the relation

$\sum_{n=1}^{\infty}(p_{i\to i})^{n-1}=\frac{1}{1-p_{i\to i}}$ (25)

and the normalization of the $p_{i\to j}$’s [see Eq. (9)], one can check that
the probability is properly normalized, i.e.,

$\sum_{j\neq i}\sum_{n=1}^{\infty}\mathcal{P}_{i\to j}(n)=1.$ (26)

Naturally, the probability of being trapped during $n$ steps is obtained by
summing over all possible exit states

$P_{i}(n)=\sum_{j\neq i}\mathcal{P}_{i\to j}(n)=\quantity(p_{i\to
i})^{n-1}\quantity(1-p_{i\to i}),$ (27)

and this defines a Poisson law with an average number of trapping events

$\bar{n}_{i}=\sum_{n=1}^{\infty}nP_{i}(n)=\frac{1}{1-p_{i\to i}}.$ (28)

Introducing the continuous time $t_{i}=n_{i}\tau$, the average trapping time
is thus given by

$\bar{t_{i}}=\quantity[H^{+}_{ii}-(E_{\text{L}}^{+})_{i}]^{-1},$ (29)

and, in the limit $\tau\to 0$, the Poisson probability takes the usual form

$P_{i}(t)=\frac{1}{\bar{t}_{i}}\exp(-\frac{t}{\bar{t}_{i}}).$ (30)

The time-averaged contribution of the on-state weight can then be easily
calculated to be

$\bar{w}_{i}=\sum_{n=1}^{\infty}w^{n}_{ii}P_{i}(n)=\frac{T_{ii}}{T^{+}_{ii}}\frac{1-T^{+}_{ii}}{1-T_{ii}}$
(31)

Details of the implementation of this effective dynamics can be in found in
Refs. Assaraf  _et al._ , 1999 and Caffarel and Assaraf, 2000.

###  III.2 Multi-state domains

Let us now extend the results of Sec. III.1 to a general domain. To do so, we
associate to each state $\ket{i}$ a set of states, called the domain of
$\ket{i}$ denoted $\mathcal{D}_{i}$, consisting of the state $\ket{i}$ plus a
certain number of states. No particular constraints on the type of domains are
imposed. For example, domains associated with different states may be
identical, and they may or may not have common states. The only important
condition is that the set of all domains ensures the ergodicity property of
the effective stochastic dynamics, that is, starting from any state, there is
a non-zero probability to reach any other state in a finite number of steps.
In practice, it is not difficult to impose such a condition.

Let us write an arbitrary path of length $N$ as

$\ket{i_{0}}\to\ket{i_{1}}\to\cdots\to\ket{i_{N}},$ (32)

where the successive states are drawn using the transition probability matrix
$p_{i\to j}$. This path belongs to the set of paths that can be represented as
follows

$(\ket*{I_{0}},n_{0})\to(\ket*{I_{1}},n_{1})\to\cdots\to(\ket*{I_{p}},n_{p}),$
(33)

where $\ket{I_{0}}=\ket{i_{0}}$ is the initial state, $n_{0}$ is the number of
times the walker remains in the domain of $\ket{i_{0}}$ (with $1\leq n_{0}\leq
N+1$), $\ket{I_{1}}$ is the first exit state that does not belong to
$\mathcal{D}_{i_{0}}$, $n_{1}$ is the number of times the walker remains in
$\mathcal{D}_{i_{1}}$ (with $1\leq n_{1}\leq N+1-n_{0}$), $\ket{I_{2}}$ is the
second exit state, and so on. Here, the integer $p$ (with $0\leq p\leq N$)
indicates the number of exit events occurring along the path. The two extreme
values, $p=0$ and $p=N$, correspond to the cases where the walker remains in
the initial domain during the entire path, and where the walker exits a domain
at each step, respectively. In what follows, we shall systematically label
exit states with upper-case letters $\ket{I_{k}}$, while lower-case letters
denote elementary states $\ket{i_{k}}$. Making this distinction is important
since the effective stochastic dynamics used in practical Monte Carlo
calculations only involve exit states $\ket{I_{k}}$, the contribution from the
elementary states $\ket{i_{k}}$ being exactly integrated out.

Figure 2 exemplifies how a path can be decomposed as proposed in Eq. (33). To
make things as clear as possible, let us describe how the path drawn in Fig. 2
evolves in time. The walker starts at $\ket{i_{0}}$ in the domain
$\mathcal{D}_{i_{0}}$. Then, it performs two steps from $\ket{i_{0}}$ to
$\ket{i_{2}}$ and has left the domain at $\ket{i_{3}}$, which is thus the
first exit state, i.e., $\ket{I_{1}}=\ket{i_{3}}$. The trapping time in
$\mathcal{D}_{i_{0}}$ is $n_{0}=3$ since three states in $\mathcal{D}_{i_{0}}$
have been visited (namely, $\ket{i_{0}}$, $\ket{i_{1}}$, and $\ket{i_{2}}$).
During the next steps, the domains $\mathcal{D}_{I_{1}}$,
$\mathcal{D}_{I_{2}}$, and $\mathcal{D}_{I_{3}}$ are successively visited with
$n_{1}=2$, $n_{2}=3$, and $n_{3}=1$, respectively. The corresponding exit
states are $\ket{I_{2}}=\ket{i_{5}}$, $\ket{I_{3}}=\ket{i_{8}}$, and
$\ket{I_{4}}=\ket{i_{9}}$, respectively.

Figure 2: Representation of a path in terms of exit states $\ket{I_{k}}$ and
trapping times $\ket{n_{k}}$. The states $\ket{i_{k}}$ along the path are
represented by small black circles and the exit states, $\ket{I_{k}}$, by
larger black squares. By convention, the initial state is denoted using a
capital letter, i.e., $\ket{i_{0}}=\ket{I_{0}}$, since it is the first state
of the effective dynamics involving only exit states. See text for additional
comments on the time evolution of the path.

Generalizing the single-state case treated previously, let us define the
probability of remaining $n$ times in the domain of $\ket{I_{0}}$ and to exit
at $\ket{I}\notin\mathcal{D}_{I_{0}}$

$\mathcal{P}_{I_{0}\to
I}(n)=\sum_{\ket{i_{1}}\in\mathcal{D}_{I_{0}}}\cdots\sum_{\ket{i_{n-1}}\in\mathcal{D}_{I_{0}}}p_{I_{0}\to
i_{1}}\ldots p_{i_{n-2}\to i_{n-1}}p_{i_{n-1}\to I}.$ (34)

Since the sums are restricted to states belonging to the domain, it is
convenient to introduce a projector over each domain

$P_{I}=\sum_{\ket{i}\in\mathcal{D}_{I}}\outerproduct{i}{i},$ (35)

as well as the projection of $T^{+}$ over $\mathcal{D}_{I}$, i.e.,
$T^{+}_{I}=P_{I}T^{+}P_{I}$, which governs the dynamics of the walkers trapped
in this domain. Using Eqs. (19) and (34), the probability can be rewritten as

$\mathcal{P}_{I_{0}\to
I}(n)=\frac{1}{\Psi^{+}_{I_{0}}}\matrixelement{I_{0}}{\quantity(T^{+}_{I_{0}})^{n-1}F^{+}_{I_{0}}}{I}\Psi^{+}_{I},$
(36)

where the operator $F^{+}_{I}=P_{I}T^{+}(1-P_{I})$, corresponding to the last
move connecting the inside and outside regions of the domain, has the
following matrix elements:

$(F^{+}_{I})_{ij}=\begin{cases}T^{+}_{ij},&\mbox{if\quad}\ket{i}\in\mathcal{D}_{I}\;\;{\rm
and}\;\;\ket{j}\notin\mathcal{D}_{I},\\\ 0,&\text{otherwise}.\end{cases}$ (37)

Physically, $F$ may be seen as a flux operator through the boundary of
$\mathcal{D}_{I}$.

Knowing the probability of remaining $n$ times in the domain and, then, to
exit to some state, it is possible to obtain the probability of being trapped
$n$ times in $\mathcal{D}_{I}$, just by summing over all possible exit states:

$P_{I}(n)=\frac{1}{\Psi^{+}_{I}}\matrixelement{I}{\quantity(T^{+}_{I})^{n-1}F^{+}_{I}}{\Psi^{+}}.$
(38)

The normalization of this probability can be verified using the fact that

$\quantity(T^{+}_{I})^{n-1}F^{+}_{I}=\quantity(T^{+}_{I})^{n-1}T^{+}-\quantity(T^{+}_{I})^{n},$
(39)

leading to111 The property results from the fact that the series is a
telescoping series and that the general term
$\matrixelement{I}{\quantity(T^{+}_{I})^{n}}{\Psi^{+}}$ goes to zero as $n$
goes to infinity.

$\sum_{n=0}^{\infty}P_{I}(n)=\frac{1}{\Psi^{+}_{I}}\sum_{n=1}^{\infty}\quantity[\matrixelement{I}{\quantity(T^{+}_{I})^{n-1}}{\Psi^{+}}-\matrixelement{I}{\quantity(T^{+}_{I})^{n}}{\Psi^{+}}]=1.$
(40)

The average trapping time defined as ${\bar{t}}_{I}={\bar{n}}_{I}\tau$ where
${\bar{n}}_{I}=\sum_{n}nP_{I}(n)$ is calculated to be

${\bar{t}}_{I}=\frac{1}{\Psi^{+}_{I}}\matrixelement{I}{{\quantity[P_{I}\quantity(H^{+}-E_{\text{L}}^{+}\mathds{1})P_{I}]}^{-1}}{\Psi^{+}}.$
(41)

In practice, the various quantities restricted to the domain will be computed
by inverting the matrix $(H^{+}-E_{\text{L}}^{+}\mathds{1})$ in
$\mathcal{D}_{I}$. Note that it is possible only if the dimension of the
domains is not too large (say, less than a few thousands).

###  III.3 Time-dependent Green’s matrix using domains

In this section we generalize the path-integral expression of the Green’s
matrix, Eq. (6), to the case where domains are used. To do so, we introduce
the Green’s matrix associated with each domain as follows:

$G^{(N),\mathcal{D}}_{ij}=\matrixelement{i}{T_{i}^{N}}{j}.$ (42)

Starting from Eq. (5) and using the representation of the paths in terms of
exit states and trapping times [see Eq. (33)], we get

$G^{(N)}_{i_{0}i_{N}}=\sum_{p=0}^{N}\sum_{\mathcal{C}_{p}}\sum_{(i_{1},...,i_{N-1})\in\mathcal{C}_{p}}\prod_{k=0}^{N-1}\matrixelement{i_{k}}{T}{i_{k+1}},$
(43)

where $\mathcal{C}_{p}$ is the set of paths with $p$ exit states,
$\ket{I_{k}}$, and trapping times $n_{k}$ with the constraints that
$\ket{I_{k}}\notin\mathcal{D}_{k-1}$ (with $1\leq n_{k}\leq N+1$ and
$\sum_{k=0}^{p}n_{k}=N+1$). It follows that

$G^{(N)}_{i_{0}i_{N}}=G^{(N),\mathcal{D}}_{i_{0}i_{N}}+\sum_{p=1}^{N}\sum_{\ket{I_{1}}\notin\mathcal{D}_{I_{0}},\ldots,\ket{I_{p}}\notin\mathcal{D}_{I_{p-1}}}\sum_{n_{0}\geq
1}\cdots\sum_{n_{p}\geq 1}\delta_{\sum_{k=0}^{p}n_{k},N+1}\\\
\times\quantity[\prod_{k=0}^{p-1}\matrixelement{I_{k}}{\quantity(T_{I_{k}})^{n_{k}-1}F_{I_{k}}}{I_{k+1}}]G^{(n_{p}-1),\mathcal{D}}_{I_{p}i_{N}},$
(44)

where $\delta_{i,j}$ is a Kronecker delta. This expression is the path-
integral representation of the Green’s matrix using only the variables
$(\ket{I_{k}},n_{k})$ of the effective dynamics defined over the set of
domains. The standard formula for $G^{(N)}_{i_{0}i_{N}}$ derived above [see
Eq. (6)] may be considered as the particular case where the walker exits of
the current state $\ket{i_{k}}$ at each step (no domains are introduced),
leading to a number of exit events $p$ equal to $N$. In this case, we have
$\ket{I_{k}}=\ket{i_{k}}$, $n_{k}=1$ (for $0\leq k\leq N$), and we are left
only with the $p$th component of the sum, that is,
$G^{(N)}_{i_{0}i_{N}}=\prod_{k=0}^{N-1}\matrixelement{I_{k}}{F_{I_{k}}}{I_{k+1}}$,
with $F=T$, thus recovering Eq. (6). Note that the first contribution
$G^{(N),\mathcal{D}}_{i_{0}i_{N}}$ corresponds to the case $p=0$ and collects
all contributions to the Green’s matrix coming from paths remaining
indefinitely in the domain of $\ket{i_{0}}$ (no exit event). This contribution
is isolated from the sum in Eq. (44) since, as a domain Green’s matrix, it is
calculated exactly and is not subject to a stochastic treatment. Note also
that the last state $\ket{i_{N}}$ is never an exit state because of the very
definition of our path representation.

In order to compute $G^{(N)}_{i_{0}i_{N}}$ by resorting to Monte Carlo
techniques, let us reformulate Eq. (44) using the transition probability
$\mathcal{P}_{I\to J}(n)$ introduced in Eq. (36). We first rewrite Eq. (44)
under the form

${G}^{(N)}_{i_{0}i_{N}}={G}^{(N),\mathcal{D}}_{i_{0}i_{N}}+{\Psi^{+}_{i_{0}}}\sum_{p=1}^{N}\sum_{\ket{I_{1}}\notin\mathcal{D}_{I_{0}},\ldots,\ket{I_{p}}\notin\mathcal{D}_{I_{p-1}}}\sum_{n_{0}\geq
1}\cdots\sum_{n_{p}\geq 1}\\\
\times\delta_{\sum_{k=0}^{p}n_{k},N+1}\quantity{\prod_{k=0}^{p-1}\quantity[\frac{\Psi^{+}_{I_{k+1}}}{\Psi^{+}_{I_{k}}}\matrixelement{I_{k}}{\quantity(T_{I_{k}})^{n_{k}-1}F_{I_{k}}}{I_{k+1}}]}\frac{{G}^{(n_{p}-1),\mathcal{D}}_{I_{p}i_{N}}}{\Psi^{+}_{I_{p}}}.$
(45)

Introducing the weights

$W_{I_{k}I_{k+1}}=\frac{\matrixelement{I_{k}}{\quantity(T_{I_{k}})^{n_{k}-1}F_{I_{k}}}{I_{k+1}}}{\matrixelement{I_{k}}{\quantity(T^{+}_{I_{k}})^{n_{k}-1}F^{+}_{I_{k}}}{I_{k+1}}},$
(46)

and using the effective transition probability, Eq. (36), we get

${G}^{(N)}_{i_{0}i_{N}}={G}^{(N),\mathcal{D}}_{i_{0}i_{N}}\\\
+{\Psi^{+}_{i_{0}}}\sum_{p=1}^{N}\sum_{{(I,n)}_{p,N}}{\quantity(\prod_{k=0}^{p-1}W_{I_{k}I_{k+1}})\quantity(\prod_{k=0}^{p-1}\mathcal{P}_{I_{k}\to
I_{k+1}}(n_{k}))\frac{1}{\Psi^{+}_{I_{p}}}{G}^{(n_{p}-1),\mathcal{D}}_{I_{p}i_{N}}},$
(47)

where, for clarity, $\sum_{{(I,n)}_{p,N}}$ is used as a short-hand notation
for
$\sum_{\ket{I_{1}}\notin\mathcal{D}_{I_{0}},\ldots,\ket{I_{p}}\notin\mathcal{D}_{I_{p-1}}}\sum_{n_{0}\geq
1}\cdots\sum_{n_{p}\geq 1}$ with the constraint $\sum_{k=0}^{p}n_{k}=N+1$.

Under this form, ${G}^{(N)}_{i_{0}i_{N}}$ is now amenable to Monte Carlo
calculations by generating paths using the transition probability matrix
$\mathcal{P}_{I\to J}(n)$. For example, in the case of the energy, we start
from

$E_{0}=\lim_{N\to\infty}\frac{\sum_{i_{N}}{G}^{(N)}_{i_{0}i_{N}}{(H\Psi_{\text{T}})}_{i_{N}}}{\sum_{i_{N}}{G}^{(N)}_{i_{0}i_{N}}{\Psi_{\text{T}}}_{i_{N}}},$
(48)

which can be rewritten probabilistically as

$E_{0}=\lim_{N\to\infty}\frac{{G}^{(N),\mathcal{D}}_{i_{0}i_{N}}+{\Psi^{+}_{i_{0}}}\sum_{p=1}^{N}\expectationvalue{\quantity(\prod_{k=0}^{p-1}W_{I_{k}I_{k+1}}){\cal
H}_{n_{p},I_{p}}}_{p}}{{G}^{(N),\mathcal{D}}_{i_{0}i_{N}}+{\Psi^{+}_{i_{0}}}\sum_{p=1}^{N}\expectationvalue{\quantity(\prod_{k=0}^{p-1}W_{I_{k}I_{k+1}}){\cal
S}_{n_{p},I_{p}}}_{p}},$ (49)

where $\expectationvalue{\cdots}_{p}$ is the probabilistic average defined
over the set of paths with $p$ exit events of probability
$\prod_{k=0}^{p-1}\mathcal{P}_{I_{k}\to I_{k+1}}(n_{k})$, and

$\displaystyle\mathcal{H}_{n_{p},I_{p}}$
$\displaystyle=\frac{1}{\Psi^{+}_{I_{p}}}\sum_{i_{N}}{G}^{(n_{p}-1),\mathcal{D}}_{I_{p}i_{N}}(H\Psi_{T})_{i_{N}},$
(50) $\displaystyle\mathcal{S}_{n_{p},I_{p}}$
$\displaystyle=\frac{1}{\Psi^{+}_{I_{p}}}\sum_{i_{N}}{G}^{(n_{p}-1),\mathcal{D}}_{I_{p}i_{N}}(\Psi_{T})_{i_{N}},$
(51)

two quantities taking into account the contribution of the trial wave function
at the end of the path.

In practice, the present Monte Carlo algorithm is a simple generalization of
the standard diffusion Monte Carlo algorithm. Stochastic paths starting at
$\ket{I_{0}}$ are generated using the probability $\mathcal{P}_{I_{k}\to
I_{k+1}}(n_{k})$ and are stopped when $\sum_{k}n_{k}$ is greater than some
target value $N$. Averages of the weight products
$\prod_{k=0}^{p-1}W_{I_{k}I_{k+1}}$ times the end-point contributions
${(\mathcal{H}/\mathcal{S})}_{n_{p},I_{p}}$ are computed for each $p$. The
generation of the paths can be performed using a two-step process. First, the
integer $n_{k}$ is drawn using the probability $P_{I_{k}}(n)$ [see Eq. (38)]
and, then, the exit state $\ket{I_{k+1}}$ is drawn using the conditional
probability $\mathcal{P}_{I_{k}\to I_{k+1}}(n_{k})/P_{I_{k}}(n_{k})$.

###  III.4 Domain Green’s Function Monte Carlo

The aim of this section is to show that it is possible to go further by
integrating out the trapping times, $n_{k}$, of the preceding expressions,
thus defining a new effective stochastic dynamics involving now only the exit
states. Physically, it means that we are going to compute exactly the time
evolution of all stochastic paths trapped in each domain. We shall present two
different ways to derive the new dynamics and renormalized probabilistic
averages. The first one, called the pedestrian way, consists in starting from
the preceding time expression for $G$ and make the explicit integration over
the $n_{k}$’s. The second, more direct and elegant, is based on the Dyson
equation.

####  III.4.1 The pedestrian way

Let us define the energy-dependent Green’s matrix

$G^{E}_{ij}=\tau\sum_{N=0}^{\infty}\matrixelement{i}{T^{N}}{j}=\matrixelement{i}{\quantity(H-E\mathds{1})^{-1}}{j}.$
(52)

The denomination “energy-dependent” is chosen here since this quantity is the
discrete version of the Laplace transform of the time-dependent Green’s
function in a continuous space, usually known under this name.222As
$\tau\rightarrow 0$ and $N\rightarrow\infty$ with $N\tau=t$, the operator
$T^{N}$ converges to $e^{-t(H-E\mathds{1})}$. We then have
$G^{E}_{ij}\rightarrow\int_{0}^{\infty}dt\matrixelement{i}{e^{-t(H-E\mathds{1})}}{j}$,
which is the Laplace transform of the time-dependent Green’s function
$\matrixelement{i}{e^{-t(H-E\mathds{1})}}{j}$. The remarkable property is
that, thanks to the summation over $N$ up to infinity, the constrained
multiple sums appearing in Eq. (44) can be factorized in terms of a product of
unconstrained sums, as follows

$\sum_{N=1}^{\infty}\sum_{p=1}^{N}\sum_{n_{0}\geq 1}\cdots\sum_{n_{p}\geq
1}\delta_{\sum_{k=0}^{p}n_{k},N+1}F(n_{0},\ldots,n_{N})\\\
=\sum_{p=1}^{\infty}\sum_{n_{0}=1}^{\infty}\cdots\sum_{n_{p}=1}^{\infty}F(n_{0},\ldots,n_{N}),$
(53)

where $F$ is some arbitrary function of the trapping times. Using the fact
that $G^{E}_{ij}=\tau\sum_{N=0}^{\infty}G^{(N)}_{ij}$, where $G^{(N)}_{ij}$ is
given by Eq. (44), and summing over the variables $n_{k}$, we get

${G}^{E}_{i_{0}i_{N}}={G}^{E,\mathcal{D}}_{i_{0}i_{N}}+\sum_{p=1}^{\infty}\sum_{I_{1}\notin\mathcal{D}_{0},\ldots,I_{p}\notin\mathcal{D}_{p-1}}\\\
\quantity[\prod_{k=0}^{p-1}\matrixelement{I_{k}}{{\quantity[P_{k}\quantity(H-E\mathds{1})P_{k}]}^{-1}(-H)(\mathds{1}-P_{k})}{I_{k+1}}]{G}^{E,\mathcal{D}}_{I_{p}i_{N}},$
(54)

where, ${G}^{E,\mathcal{D}}$ is the energy-dependent domain’s Green matrix
defined as
${G}^{E,\mathcal{D}}_{ij}=\tau\sum_{N=0}^{\infty}\matrixelement{i}{T^{N}_{i}}{j}$.

As a didactical example, Appendix A reports the exact derivation of this
formula in the case of a two-state system.

####  III.4.2 Dyson equation

In fact, there is a more direct way to derive the same equation by resorting
to the Dyson equation. Starting from the well-known equality

$\quantity(H-E\mathds{1})^{-1}=\quantity(H_{0}-E\mathds{1})^{-1}+\quantity(H_{0}-E\mathds{1})^{-1}(H_{0}-H)\quantity(H-E\mathds{1})^{-1},$
(55)

where $H_{0}$ is some arbitrary reference Hamiltonian, we have the Dyson
equation

$G^{E}_{ij}=G^{E}_{0,ij}+\sum_{kl}G^{E}_{0,ik}(H_{0}-H)_{kl}G^{E}_{lj},$ (56)

with $G^{E}_{0,ij}=\matrixelement{i}{\quantity(H_{0}-E\mathds{1})^{-1}}{j}$.
Let us choose $H_{0}$ such that
$\matrixelement{i}{H_{0}}{j}=\matrixelement{i}{P_{i}HP_{i}}{j}$ for all $i$
and $j$. Then, the Dyson equation (56) becomes

${G}^{E}_{ij}={G}^{E,\mathcal{D}}_{ij}+\sum_{k}\matrixelement{i}{{\quantity[P_{i}\quantity(H-E\mathds{1})P_{i}]}^{-1}(H_{0}-H)}{k}{G}^{E}_{kj}.$
(57)

Using the following identity

$\begin{split}{\quantity[P_{i}\quantity(H-E\mathds{1})P_{i}]}^{-1}(H_{0}-H)&={\quantity[P_{i}\quantity(H-E\mathds{1})P_{i}]}^{-1}(P_{i}HP_{i}-H)\\\
&={\quantity[P_{i}\quantity(H-E\mathds{1})P_{i}]}^{-1}(-H)(\mathds{1}-P_{i}),\end{split}$
(58)

the Dyson equation may be written under the form

$G^{E}_{ij}={G}^{E,\mathcal{D}}_{ij}+\sum_{k\notin\mathcal{D}_{i}}\matrixelement{i}{{\quantity[P_{i}\quantity(H-E\mathds{1})P_{i}]}^{-1}(-H)(\mathds{1}-P_{i})}{k}{G}^{E}_{kj},$
(59)

which is identical to Eq. (54) when $G^{E}_{ij}$ is expanded iteratively. Let
us use as effective transition probability density

$P_{I\to
J}=\frac{1}{\Psi^{+}(I)}\matrixelement{I}{{\quantity[P_{I}\quantity(H^{+}-E_{\text{L}}^{+}\mathds{1})P_{I}]}^{-1}(-H^{+})(\mathds{1}-P_{I})}{J}\Psi^{+}(J),$
(60)

and the weight

$W^{E}_{IJ}=\frac{\matrixelement{I}{{\quantity[P_{I}\quantity(H-E\mathds{1})P_{I}]}^{-1}(-H)(\mathds{1}-P_{I})}{J}}{\matrixelement{I}{{\quantity[P_{I}\quantity(H^{+}-E_{\text{L}}^{+}\mathds{1})P_{I}]}^{-1}(-H^{+})(\mathds{1}-P_{I})}{J}}.$
(61)

Using Eqs. (34), (36) and (39), one can easily verify that $P_{I\to J}\geq 0$
and $\sum_{J}P_{I\to J}=1$. Finally, the probabilistic expression reads

$G^{E}_{i_{0}i_{N}}={G}^{E,\mathcal{D}}_{i_{0}i_{N}}+\sum_{p=1}^{\infty}\Psi^{+}_{i_{0}}\expectationvalue{\quantity(\prod_{k=0}^{p-1}W^{E}_{I_{k}I_{k+1}})\frac{{G}^{E,\mathcal{D}}_{I_{p}i_{N}}}{\Psi^{+}_{I_{p}}}}.$
(62)

####  III.4.3 Energy estimator

To calculate the energy, we introduce the following estimator

$\mathcal{E}(E)=\frac{\matrixelement{I_{0}}{\quantity(H-E\mathds{1})^{-1}}{H\Psi_{\text{T}}}}{\matrixelement{I_{0}}{\quantity(H-E\mathds{1})^{-1}}{\Psi_{\text{T}}}},$
(63)

which is the counterpart of the quantity ${\cal
E}_{N}=\frac{\matrixelement{I_{0}}{T^{N}}{H\Psi_{\text{T}}}}{\matrixelement{I_{0}}{T^{N}}{\Psi_{\text{T}}}}$
used in the time-dependent formalism. In this case, the energy was easily
obtained by taking the large $N$-limit of ${\cal E}_{N}$, see Eq.(3). Here,
the situation is not as simple and we must find a way to extract the energy
from $\mathcal{E}(E)$.

Using the spectral decomposition of $H$, we have

$\mathcal{E}(E)=\frac{\sum_{i}\frac{E_{i}c_{i}}{E_{i}-E}}{\sum_{i}\frac{c_{i}}{E_{i}-E}},$
(64)

where
$c_{i}=\innerproduct{I_{0}}{\Phi_{i}}\innerproduct{\Phi_{i}}{\Psi_{\text{T}}}$
and $\Phi_{i}$ are the eigenstates of $H$, i.e.,
$H\ket{\Phi_{i}}=E_{i}\ket{\Phi_{i}}$. The important observation is that for
all eigenstates we have $\mathcal{E}(E_{i})=E_{i}$. Thus, to get the ground-
state energy we propose to search for the solution of the non-linear equation
$\mathcal{E}(E)=E$ in the vicinity of $E_{0}$.

In practical Monte Carlo calculations, the DMC energy is obtained by computing
a finite number of components $H_{p}$ and $S_{p}$ defined as follows

$\mathcal{E}^{\text{DMC}}(E,p_{max})=\frac{H_{0}+\sum_{p=1}^{p_{\text{max}}}H^{\text{DMC}}_{p}}{S_{0}+\sum_{p=1}^{p_{\text{max}}}S^{\text{DMC}}_{p}}.$
(65)

For $p\geq 1$, Eq. (62) gives

$\displaystyle H^{\text{DMC}}_{p}$
$\displaystyle=\Psi^{+}_{i_{0}}\expectationvalue{\quantity(\prod_{k=0}^{p-1}W^{E}_{I_{k}I_{k+1}}){\cal
H}_{I_{p}}},$ (66) $\displaystyle S^{\text{DMC}}_{p}$
$\displaystyle=\Psi^{+}_{i_{0}}\expectationvalue{\quantity(\prod_{k=0}^{p-1}W^{E}_{I_{k}I_{k+1}}){\cal
S}_{I_{p}}},$ (67)

where

$\displaystyle\mathcal{H}_{I_{p}}$
$\displaystyle=\frac{1}{\Psi^{+}_{I_{p}}}\sum_{i_{N}}{G}^{E,\mathcal{D}}_{I_{p}i_{N}}(H\Psi_{T})_{i_{N}},$
(68) $\displaystyle\mathcal{S}_{I_{p}}$
$\displaystyle=\frac{1}{\Psi^{+}_{I_{p}}}\sum_{i_{N}}{G}^{E,\mathcal{D}}_{I_{p}i_{N}}(\Psi_{T})_{i_{N}}.$
(69)

For $p=0$, the two components are exactly evaluated as

$\displaystyle H_{0}$
$\displaystyle=\matrixelement{I_{0}}{{\quantity[P_{I_{0}}\quantity(H-E\mathds{1})P_{I_{0}}]}^{-1}}{H\Psi_{\text{T}}},$
(70) $\displaystyle S_{0}$
$\displaystyle=\matrixelement{I_{0}}{{\quantity[P_{I_{0}}\quantity(H-E\mathds{1})P_{I_{0}}]}^{-1}}{\Psi_{\text{T}}}.$
(71)

Note that the evaluation of $(H_{0},S_{0})$ is possible as long as the size of
the domain is small enough to allow the calculation of the inverse matrix.
Avoiding the stochastic calculation of quantities, such as $H_{0}$ or $S_{0}$,
that can be evaluated analytically is computationally very appealing as the
statistical error associated with these quantities is completely removed. We
thus suggest to extend the exact calculation of $H_{p}$ and $S_{p}$ to higher
$p$ values, up to the point where the exponential increase of the number of
intermediate states involved in the explicit sums makes the calculation
unfeasible.

Finally, $\mathcal{E}^{\text{DMC}}(E,p_{\text{ex}},p_{\text{max}})$ is
evaluated in practice as follows

$\mathcal{E}^{\text{DMC}}(E,p_{\text{ex}},p_{\text{max}})=\frac{\sum_{p=0}^{p_{\text{ex}}-1}H_{p}+\sum_{p=p_{\text{ex}}}^{p_{\text{max}}}H^{\text{DMC}}_{p}}{\sum_{p=0}^{p_{\text{ex}}-1}S_{p}+\sum_{p=p_{\text{ex}}}^{p_{\text{max}}}S^{\text{DMC}}_{p}},$
(72)

where $p_{\text{ex}}$ is the number of components computed exactly. Let us
emphasize that, since the magnitude of $H_{p}$ and $S_{p}$ decreases as a
function of $p$, most of the statistical error is removed by computing the
dominant terms analytically. This will be illustrated in the numerical
application presented below.

It is easy to check that, in the vicinity of the exact energy $E_{0}$,
$\mathcal{E}(E)$ is a linear function of $E-E_{0}$. Therefore, in practice, we
compute the value of $\mathcal{E}(E)$ for several values of $E$, and fit these
data using a linear, quadratic or a more complicated function of $E$ in order
to obtain, via extrapolation, an estimate of $E_{0}$. In order to have a
precise extrapolation of the energy, it is best to compute $\mathcal{E}(E)$
for values of $E$ as close as possible to the exact energy. However, as $E\to
E_{0}$, both the numerators and denominators of Eq. (64) diverge. This is
reflected by the fact that one needs to compute more and more $p$-components
with an important increase of statistical fluctuations. Thus, from a practical
point of view, a trade-off has to be found between the quality of the
extrapolation and the amount of statistical fluctuations.

##  IV Numerical application to the Hubbard model

###  IV.1 Hamiltonian and trial wave function

Let us consider the one-dimensional Hubbard Hamiltonian for a chain of $N$
sites

$H=-t\sum_{\expectationvalue{ij}\sigma}\hat{a}^{+}_{i\sigma}\hat{a}_{j\sigma}+U\sum_{i}\hat{n}_{i\uparrow}\hat{n}_{i\downarrow},$
(73)

where $\expectationvalue{ij}$ denotes the summation over two neighboring
sites, $\hat{a}^{+}_{i\sigma}$ ($\hat{a}_{i\sigma}$) is the fermionic creation
(annihilation) operator of an spin-$\sigma$ electron (with $\sigma$ =
$\uparrow$ or $\downarrow$) on site $i$,
$\hat{n}_{i\sigma}=\hat{a}^{+}_{i\sigma}\hat{a}_{i\sigma}$ the number
operator, $t$ the hopping amplitude, and $U$ the on-site Coulomb repulsion. We
consider a chain with an even number of sites and open boundary conditions at
half-filling, that is, $N_{\uparrow}=N_{\downarrow}=N/2$. In the site
representation, a general vector of the Hilbert space can be written as

$\ket{n}=\ket{n_{1\uparrow},\ldots,n_{N\uparrow},n_{1\downarrow},\ldots,n_{N\downarrow}},$
(74)

where $n_{i\sigma}=0$ or $1$ is the number of electrons of spin $\sigma$ on
site $i$.

For the one-dimensional Hubbard model with open boundary conditions, the
components of the ground-state vector have the same sign (say, $c_{i}>0$). It
is then possible to equate the guiding and trial vectors, i.e.,
$\ket{c^{+}}=\ket{c_{\text{T}}}$. As trial wave function, we shall employ a
generalization of the Gutzwiller wave function Gutzwiller (1963)

$\innerproduct{n}{c_{\text{T}}}=e^{-\alpha n_{D}(n)-\beta n_{A}(n)},$ (75)

where $n_{D}(n)$ is the number of doubly occupied sites for the configuration
$\ket{n}$ and $n_{A}(n)$ the number of nearest-neighbor antiparallel pairs
defined as

$n_{A}(n)=\sum_{\expectationvalue{ij}}n_{i\uparrow}n_{j\downarrow}\delta_{n_{i\uparrow},1}\delta_{n_{j\downarrow},1}.$
(76)

The parameters $\alpha$ and $\beta$ are optimized by minimizing the
variational energy of the trial wave function, i.e.,

$E_{\text{T}}(\alpha,\beta)=\frac{\matrixelement{c_{\text{T}}}{H}{c_{\text{T}}}}{\innerproduct{c_{\text{T}}}{c_{\text{T}}}}.$
(77)

###  IV.2 Domains

As discussed above, the efficiency of the method depends on the choice of
states forming each domain. As a general guiding principle, it is advantageous
to build domains associated with a large average trapping time in order to
integrate out the most important part of the Green’s matrix.

Here, as a first illustration of the method, we shall consider the large-$U$
regime of the Hubbard model where the construction of such domains is rather
natural. Indeed, at large $U$ and half-filling, the Hubbard model approaches
the Heisenberg limit where only the $2^{N}$ states with no double occupancy,
$n_{D}(n)=0$, have a significant weight in the wave function. The contribution
of the other states vanishes as $U$ increases with a rate increasing sharply
with $n_{D}(n)$. In addition, for a given number of double occupations,
configurations with large values of $n_{A}(n)$ are favored due to their high
electronic mobility. Therefore, we build domains associated with small $n_{D}$
and large $n_{A}$ in a hierarchical way as described below.

For simplicity and reducing the number of matrix inversions to perform, we
shall consider only one non-trivial domain called here the main domain and
denoted as $\mathcal{D}$. This domain will be chosen common to all states
belonging to it, that is,

$\mathcal{D}_{i}=\mathcal{D}\mbox{\quad for\quad}\ket{i}\in\mathcal{D}.$ (78)

For the other states, we choose a single-state domain as

$\mathcal{D}_{i}=\quantity{\ket{i}}\mbox{\quad
for\quad}\ket{i}\notin\mathcal{D}.$ (79)

To define $\mathcal{D}$, let us introduce the following set of states

${\mathcal{S}_{ij}=\quantity{\ket{n}:n_{D}(n)=i\land n_{A}(n)=j}}.$ (80)

$\mathcal{D}$ is defined as containing the set of states having up to
$n_{D}^{\text{max}}$ double occupations and, for each state with a number of
double occupations equal to $m$, a number of nearest-neighbor antiparallel
pairs between $n_{A}^{\text{min}}(m)$ and $n_{A}^{\text{max}}(m)$. Here,
$n_{A}^{\text{max}}(m)$ is fixed and set at its maximum value for a given $m$,
i.e., $n_{A}^{\text{max}}(m)=\max(N-1-2m,0)$. Using these definitions, the
main domain is taken as the union of some elementary domains

$\mathcal{D}=\bigcup_{n_{D}=0}^{n_{D}^{\text{max}}}\mathcal{D}(n_{D},n_{A}^{\text{min}}(n_{D})),$
(81)

where the elementary domains are defined as

$\mathcal{D}(n_{D},n_{A}^{\text{min}}(n_{D}))=\bigcup_{j=n_{A}^{\text{min}}(n_{D})}^{n_{A}^{\text{max}}(n_{D})}\mathcal{S}_{n_{D}j}.$
(82)

The two quantities defining the main domain are thus $n_{D}^{\text{max}}$ and
$n_{A}^{\text{min}}(m)$. To give an illustrative example, let us consider the
4-site case. There are 6 possible elementary domains:

$\displaystyle\mathcal{D}(0,3)$ $\displaystyle=\mathcal{S}_{03},$
$\displaystyle\mathcal{D}(0,2)$
$\displaystyle=\mathcal{S}_{03}\cup\mathcal{S}_{02},$
$\displaystyle\mathcal{D}(0,1)$
$\displaystyle=\mathcal{S}_{03}\cup\mathcal{S}_{02}\cup\mathcal{S}_{01},$
$\displaystyle\mathcal{D}(1,1)$ $\displaystyle=\mathcal{S}_{11},$
$\displaystyle\mathcal{D}(1,0)$
$\displaystyle=\mathcal{S}_{11}\cup\mathcal{S}_{10},$
$\displaystyle\mathcal{D}(2,0)$ $\displaystyle=\mathcal{S}_{20},$

where

$\displaystyle\mathcal{S}_{03}$
$\displaystyle=\quantity{\ket{\uparrow,\downarrow,\uparrow,\downarrow},\ket{\downarrow,\uparrow,\downarrow,\uparrow}},\mbox{\quad(the
two N\'{e}el states)\quad}$ $\displaystyle\mathcal{S}_{02}$
$\displaystyle=\quantity{\ket{\uparrow,\downarrow,\downarrow,\uparrow},\ket{\downarrow,\uparrow,\uparrow,\downarrow}},$
$\displaystyle\mathcal{S}_{01}$
$\displaystyle=\quantity{\ket{\uparrow,\uparrow,\downarrow,\downarrow},\ket{\downarrow,\downarrow,\uparrow,\uparrow}},$
$\displaystyle\mathcal{S}_{11}$
$\displaystyle=\quantity{\ket{\uparrow\downarrow,\uparrow,\downarrow,0},\ket{\uparrow\downarrow,0,\uparrow,\uparrow}+\ldots},$
$\displaystyle\mathcal{S}_{10}$
$\displaystyle=\quantity{\ket{\uparrow\downarrow,\uparrow,0,\downarrow},\ket{\uparrow\downarrow,0,\uparrow,\downarrow}+\ldots},$
$\displaystyle\mathcal{S}_{20}$
$\displaystyle=\quantity{\ket{\uparrow\downarrow,\uparrow\downarrow,0,0}+\ldots}.$

For the three last cases, the dots indicate that one must also consider the
remaining states obtained by permuting the position of the pairs.

###  IV.3 DMC simulations

Let us now present our DMC calculations for the Hubbard model. In what
follows, we shall restrict ourselves to the case of the Green’s function Monte
Carlo approach where trapping times are integrated out exactly.

Let us begin with a small chain of 4 sites with $U=12$. From now on, we shall
take $t=1$. The size of the linear space is ${\binom{4}{2}}^{2}=36$ and the
ground-state energy obtained by exact diagonalization is
$E_{0}=-0.768068\ldots$. The two variational parameters of the trial vector
have been optimized and fixed at the values of $\alpha=1.292$ and
$\beta=0.552$ with a variational energy $E_{\text{T}}=-0.495361\ldots$. In
what follows, $\ket{I_{0}}$ is systematically chosen as one of the two Néel
states, e.g., $\ket{I_{0}}=\ket{\uparrow,\downarrow,\uparrow,\ldots}$.

Figure 3 shows the convergence of $H_{p}$ as a function of $p$ for different
values of the reference energy $E$. We consider the simplest case where a
single-state domain is associated to each state. Five different values of $E$
have been chosen, namely $E=-1.6$, $-1.2$, $-1.0$, $-0.9$, and $-0.8$. Only
$H_{0}$ is computed analytically ($p_{\text{ex}}=0$). At the scale of the
figure, error bars are too small to be seen.

When $E$ is far from the exact value, the convergence is very rapid and only a
few terms of the $p$-expansion are necessary. In contrast, when $E$ approaches
the exact energy, a slower convergence is observed, as expected from the
divergence of the matrix elements of the Green’s matrix at $E=E_{0}$ where the
expansion does not converge at all. The oscillations of the curves as a
function of $p$ are due to a parity effect specific to this system. In
practice, it is not too much of a problem since a smoothly convergent behavior
is nevertheless observed for the even- and odd-parity curves. The ratio,
$\mathcal{E}^{\text{DMC}}(E,p_{\text{ex}}=1,p_{\text{max}})$ as a function of
$E$ is presented in Fig. 4. Here, $p_{\text{max}}$ is taken sufficiently large
so that the convergence at large $p$ is reached. The values of $E$ are
$-0.780$, $-0.790$, $-0.785$, $-0.780$, and $-0.775$. For small $E$, the curve
is extrapolated using the so-called two-component expression:

$\mathcal{E}(E)=\frac{\frac{\epsilon_{0}c_{0}}{\epsilon_{0}-E}+\frac{\epsilon_{1}c_{1}}{\epsilon_{1}-E}}{\frac{c_{0}}{\epsilon_{0}-E}+\frac{c_{1}}{\epsilon_{1}-E}},$
(83)

which considers only the first two terms of the exact expression of
$\mathcal{E}(E)$ [see Eq. (64)]. Here, the fitting parameters that needs to be
determined are $c_{0}$, $\epsilon_{0}$, $c_{1}$, and $\epsilon_{1}$. The
estimate of the energy obtained from $\mathcal{E}(E)=E$ is $-0.76807(5)$ in
full agreement with the exact value of $-0.768068\ldots$.

Figure 3: One-dimensional Hubbard model with $N=4$ and $U=12$. $H_{p}$ as a
function of $p$ for $E=-1.6$, $-1.2$, $-1.0$, $-0.9$, and $-0.8$. $H_{0}$ is
computed analytically and $H_{p}$ ($p>0$) is computed stochastically. Error
bars are smaller than the symbol size. Figure 4: One-dimensional Hubbard model
with $N=4$ and $U=12$. $\mathcal{E}(E)$ as a function of $E$. The horizontal
and vertical lines are at $\mathcal{E}(E_{0})=E_{0}$ and $E=E_{0}$,
respectively. $E_{0}=-0.768068\ldots$ is the exact energy. The dotted line is
the two-component extrapolation defined in Eq. (83). Error bars are smaller
than the symbol size.

Table 1 illustrates the dependence of the Monte Carlo results upon the choice
of the domain. The reference energy is $E=-1$. The first column indicates the
various domains consisting of the union of some elementary domains as
explained above. The first line of the table gives the results when one uses a
single-state domain for all states, and the last one for the maximal domain
containing the full linear space. The size of the various domains is given in
the second column, the average trapping time for the state $\ket{I_{0}}$ in
the third column, and an estimate of the speed of convergence of the
$p$-expansion for the energy in the fourth column. To quantify the rate of
convergence, we report the quantity, $p_{\text{conv}}$, defined as the
smallest value of $p$ for which the energy is converged to six decimal places.
The smaller $p_{\text{conv}}$, the better the convergence is. Although this is
a rough estimate, it is sufficient here for our purpose. As clearly seen, the
speed of convergence is directly related to the magnitude of
$\bar{t}_{I_{0}}$. The longer the stochastic trajectories remain trapped
within the domain, the better the convergence. Of course, when the domain is
chosen to be the full space, the average trapping time becomes infinite. Let
us emphasize that the rate of convergence has no reason to be related to the
size of the domain. For example, the domain
$\mathcal{D}(0,3)\cup\mathcal{D}(1,0)$ has a trapping time for the Néel state
of $6.2$, while the domain $\mathcal{D}(0,3)\cup\mathcal{D}(1,1)$, having
almost the same number of states (28 states), has an average trapping time
about 6 times longer. Finally, the last column gives the energy obtained for
$E=-1$. The energy is expected to be independent of the domain and to converge
to a common value, which is indeed the case here. The exact value,
$\mathcal{E}(E=-1)=-0.75272390\ldots$, can be found in the last row of Table 1
for the case of a domain corresponding to the full space. In sharp contrast,
the statistical error depends strongly on the type of domains used. As
expected, the largest error of $3\times 10^{-5}$ is obtained in the case of a
single-state domain for all states. The smallest statistical error is obtained
for the “best” domain having the largest average trapping time. Using this
domain leads to a reduction in the statistical error as large as about three
orders of magnitude, which nicely illustrates the critical importance of the
domains employed.

Table 1: One-dimensional Hubbard model with $N=4$, $U=12$, $E=-1$, $\alpha=1.292$, $\beta=0.552$, and $p_{\text{ex}}=4$. The simulation is performed with 20 independent blocks and $10^{5}$ stochastic paths starting from the Néel state. $\bar{t}_{I_{0}}$ is the average trapping time for the Néel state. $p_{\text{conv}}$ is a measure of the convergence of $\mathcal{E}^{\text{DMC}}(p)$ as a function of $p$. See text for more details. Domain | Size | $\bar{t}_{I_{0}}$ | $p_{\text{conv}}$ | $\mathcal{E}^{\text{DMC}}$
---|---|---|---|---
Single | 1 | 0.026 | 88 | $-0.75276(3)$
$\mathcal{D}(0,3)$ | 2 | 2.1 | 110 | $-0.75276(3)$
$\mathcal{D}(0,2)$ | 4 | 2.1 | 106 | $-0.75275(2)$
$\mathcal{D}(0,1)$ | 6 | 2.1 | 82 | $-0.75274(3)$
$\mathcal{D}(0,3)\cup\mathcal{D}(1,1)$ | 14 | 4.0 | 60 | $-0.75270(2)$
$\mathcal{D}(0,3)\cup\mathcal{D}(1,0)$ | 26 | 6.2 | 45 | $-0.752730(7)$
$\mathcal{D}(0,2)\cup\mathcal{D}(1,1)$ | 16 | 10.1 | 36 | $-0.75269(1)$
$\mathcal{D}(0,2)\cup\mathcal{D}(1,0)$ | 28 | 34.7 | 14 | $-0.7527240(6)$
$\mathcal{D}(0,1)\cup\mathcal{D}(1,1)$ | 18 | 10.1 | 28 | $-0.75269(1)$
$\mathcal{D}(0,1)\cup\mathcal{D}(1,0)$ | 30 | 108.7 | 11 | $-0.75272400(5)$
$\mathcal{D}(0,3)\cup\mathcal{D}(1,1)\cup\mathcal{D}$(2,0) | 20 | 4.1 | 47 | $-0.75271(2)$
$\mathcal{D}(0,3)\cup\mathcal{D}(1,0)\cup\mathcal{D}$(2,0) | 32 | 6.5 | 39 | $-0.752725(3)$
$\mathcal{D}(0,2)\cup\mathcal{D}(1,1)\cup\mathcal{D}$(2,0) | 22 | 10.8 | 30 | $-0.75270(1)$
$\mathcal{D}(0,2)\cup\mathcal{D}(1,0)\cup\mathcal{D}$(2,0) | 34 | 52.5 | 13 | $-0.7527236(2)$
$\mathcal{D}(0,1)\cup\mathcal{D}(1,1)\cup\mathcal{D}$(2,0) | 24 | 10.8 | 26 | $-0.75270(1)$
$\mathcal{D}(0,1)\cup\mathcal{D}(1,0)\cup\mathcal{D}$(2,0) | 36 | $\infty$ | 1 | $-0.75272390$

As explained above, it is very advantageous to calculate exactly as many
$(H_{p},S_{p})$ as possible in order to avoid the statistical error on the
largest components. Table 2 shows the results both for the case of a single-
state domain and for the domain having the largest average trapping time,
namely $\mathcal{D}(0,1)\cup\mathcal{D}(1,1)$ (see Table 1). Table 2 reports
the statistical fluctuations of the energy for the simulation of Table 1.
Results show that it is indeed worth computing exactly as many components as
possible. For the single-state domain, the statistical error is reduced by a
factor two when passing from no analytical computation, $p_{\text{ex}}=0$, to
the case where eight components for $H_{p}$ and $S_{p}$ are computed exactly.
For the best domain, the impact is much more important with a huge reduction
of about three orders of magnitude in the statistical error.

Table 3 reports the energies convergence as a function of $p$ alongside their
statistical error on the last digit for $E=-0.8$, $-0.795$, $-0.79$, $-0.785$,
and $-0.78$. The values are displayed in Fig. 4. As seen, the behavior of
$\mathcal{E}$ as a function of $E$ is almost perfectly linear. The
extrapolated values obtained from the five values of the energy with three
different fitting functions are reported. Using the linear fitting function
leads to an energy of $-0.7680282(5)$ (to be compared with the exact value of
$-0.768068\ldots$). A small bias of about $4\times 10^{-5}$ is observed. This
bias vanishes within the statistical error when one relies on the quadratic
and two-component fitting functions. Our final value is in full agreement with
the exact value with about six decimal places.

Table 4 shows the evolution of the average trapping times and extrapolated
energies as a function of $U$ when using
$\mathcal{D}(0,1)\cup\mathcal{D}(1,0)$ as the main domain. We also report the
variational and exact energies together with the values of the optimized
parameters of the trial wave function. As $U$ increases the configurations
with zero or one double occupation become more and more predominant and the
average trapping time increases. For very large values of $U$ (say, $U\geq
12$) the increase of $\bar{t}_{I_{0}}$ becomes particularly steep.

Finally, in Table 5, we report the results obtained for larger systems at
$U=12$ for a chain size ranging from $N=4$ (36 states) to $N=12$ ($\sim
10^{6}$ states). No careful construction of domains maximizing the average
trapping time has been performed, we have merely chosen domains of reasonable
size (no more than 2682) by taking not too large number of double occupations
(only, $n_{D}=0,1$) and not too small number of nearest-neighbor antiparallel
pairs. As seen, as the number of sites increases, the average trapping time
for the chosen domains decreases. This point is, of course, undesirable since
the efficiency of the approach may gradually deteriorate when considering
large systems. A more elaborate way of constructing domains is clearly
desirable in this case. The exact DMC energies extrapolated using the two-
component function are also reported. Similarly to what has been done for
$N=4$, the extrapolation is performed using about five values of the reference
energy. The extrapolated DMC energies are in full agreement with the exact
value within error bars. However, an increase of the statistical error is
observed when the system size increases. To get lower error bars, more
accurate trial wave functions may be considered, better domains, and also
larger simulation times. Of course, it will also be particularly interesting
to take advantage of the fully parallelizable character of the algorithm to
get much lower error bars. All these aspects will be considered in a
forthcoming work.

Table 2: One-dimensional Hubbard model with $N=4$, $U=12$, and $E=-1$. Dependence of the statistical error on the energy with the number of $p$-components calculated analytically. The simulation is performed the same way as in Table 1. Results are presented when a single-state domain is used for all states and when $\mathcal{D}(0,1)\cup\mathcal{D}(1,0)$ is chosen as the main domain. $p_{ex}$ | single-state | $\mathcal{D}(0,1)\cup\mathcal{D}(1,0)$
---|---|---
$0$ | $4.3\times 10^{-5}$ | $347\times 10^{-8}$
$1$ | $4.0\times 10^{-5}$ | $377\times 10^{-8}$
$2$ | $3.7\times 10^{-5}$ | $43\times 10^{-8}$
$3$ | $3.3\times 10^{-5}$ | $46\times 10^{-8}$
$4$ | $2.6\times 10^{-5}$ | $5.6\times 10^{-8}$
$5$ | $2.5\times 10^{-5}$ | $6.0\times 10^{-8}$
$6$ | $2.3\times 10^{-5}$ | $0.7\times 10^{-8}$
$7$ | $2.2\times 10^{-5}$ | $0.6\times 10^{-8}$
$8$ | $2.2\times 10^{-5}$ | $0.05\times 10^{-8}$
Table 3: One-dimensional Hubbard model with $N=4$, $U=12$, $\alpha=1.292$, $\beta=0.552$, and $p_{\text{ex}}=4$. The main domain is $\mathcal{D}(0,1)\cup\mathcal{D}(1,0)$. The simulation is performed with 20 independent blocks and $10^{6}$ stochastic paths starting from the Néel state. The fits are performed with the five values of $E$. $E$ | $E^{\text{DMC}}$
---|---
$-0.8$ | $-0.7654686(2)$
$-0.795$ | $-0.7658622(2)$
$-0.79$ | $-0.7662607(3)$
$-0.785$ | $-0.7666642(4)$
$-0.78$ | $-0.7670729(5)$
$E_{0}$ linear fit | $-0.7680282(5)$
$E_{0}$ quadratic fit | $-0.7680684(5)$
$E_{0}$ two-component fit | $-0.7680676(5)$
$E_{0}$ exact | $-0.768068\ldots$
Table 4: One-dimensional Hubbard model with $N=4$. Average trapping time as a function of $U$. The main domain is $\mathcal{D}(0,1)\cup\mathcal{D}(1,0)$. $U$ | $(\alpha,\beta)$ | $E_{\text{T}}$ | $E_{\text{ex}}$ | $\bar{t}_{I_{0}}$
---|---|---|---|---
8 | (0.908, 0.520) | $-0.770342\ldots$ | $-1.117172\ldots$ | 33.5
10 | (1.116, 0.539) | $-0.604162\ldots$ | $-0.911497\ldots$ | 63.3
12 | (1.292, 0.552) | $-0.495361\ldots$ | $-0.768068\ldots$ | 108.7
14 | (1.438, 0.563) | $-0.419163\ldots$ | $-0.662871\ldots$ | 171.7
20 | (1.786, 0.582) | $-0.286044\ldots$ | $-0.468619\ldots$ | 504.5
50 | (2.690, 0.609) | $-0.110013\ldots$ | $-0.188984\ldots$ | 8040.2
200 | (4.070, 0.624) | $-0.026940\ldots$ | $-0.047315\ldots$ | 523836.0
Table 5: One-dimensional Hubbard model with $U=12$. The extrapolated DMC energies are obtained using the two-component fitting function. $N$ | Hilbert space size | Domain | Domain size | $(\alpha,\beta)$ | $\bar{t}_{I_{0}}$ | $E_{\text{T}}$ | $E_{\text{ex}}$ | $E^{\text{DMC}}$
---|---|---|---|---|---|---|---|---
4 | 36 | $\mathcal{D}(0,1)\cup\mathcal{D}(1,0)$ | 30 | (1.292, 0.552) | 108.7 | $-0.495361$ | $-0.768068$ | $-0.7680676(5)$
6 | 400 | $\mathcal{D}(0,1)\cup\mathcal{D}(1,0)$ | 200 | (1.124, 0.689) | 57.8 | $-0.633297$ | $-1.215395$ | $-1.215389(9)$
8 | 4 900 | $\mathcal{D}(0,1)\cup\mathcal{D}(1,0)$ | 1 190 | (0.984, 0.788) | 42.8 | $-0.750995$ | $-1.66395$ | $-1.6637(2)$
10 | 63 504 | $\mathcal{D}(0,5)\cup\mathcal{D}(1,4)$ | 2 682 | (0.856, 0.869) | 31.0 | $-0.855958$ | $-2.113089$ | $-2.1120(7)$
12 | 853 776 | $\mathcal{D}(0,8)\cup\mathcal{D}(1,7)$ | 1 674 | (0.739, 0.938) | 16.7 | $-0.952127$ | $-2.562529$ | $-2.560(6)$

##  V Summary and perspectives

In this work, it has been shown how to integrate out exactly — within a DMC
framework — the contribution of all stochastic trajectories trapped in some
given domains of the Hilbert space. The corresponding equations have been
derived in a general context. In such a way, a new effective stochastic
dynamics connecting only domains (and not the individual states) is defined. A
key property of this effective dynamics is that it does not depend on the
“shape” of the domains used for each state. Therefore, rather general domains
(with or without overlap) can be considered.

To obtain the effective transition probability (which provides the probability
of going from one domain to another) and the corresponding renormalized
estimators, the Green’s functions restricted to each domain must be computed
analytically, that is, in practice, matrices of the size of the number of
states for the sampled domains have to be inverted. This is the main
computationally intensive step of the present approach. The efficiency of the
method is directly related to the importance of the average time spent by the
stochastic trajectories in each domain.

Being able to define domains with large average trapping times is the key
aspect of the method since it may lead to some important reduction of the
statistical error, as illustrated in our numerical applications. Therefore, a
trade-off has to be found between maximizing the average trapping time and
minimizing the cost of computing the domain Green’s functions. In practice,
there is no general rule to construct such domains. For each system at hand,
one needs to determine, on physical grounds, which regions of the
configuration space are preferentially sampled by the stochastic trajectories
and to build domains of minimal size enclosing such regions.

In the first application presented here on the one-dimensional Hubbard model,
we exploit the physics of the large-$U$ regime that is known to approach the
Heisenberg limit where double occupations have small weights. This simple
example has been chosen to illustrate the various aspects of the approach.

Our goal is, of course, to tackle much larger systems, like those treated by
state-of-the-art methods, such as selected CI, Huron  _et al._ (1973);
Harrison (1991); Giner  _et al._ (2013); Holmes  _et al._ (2016); Schriber and
Evangelista (2016); Tubman  _et al._ (2020), FCIQMC, Booth  _et al._ (2009);
Cleland  _et al._ (2010), AFQMC, Zhang and Krakauer (2003) or DMRG. White and
Martin (1999); Chan and Sharma (2011) Here, we have mainly focused on the
theoretical aspects of the approach. In order to consider larger systems, an
elaborate implementation of the present method is necessary in order to keep
under control the cost of the simulation. This is outside the scope of the
present study and will be presented in a forthcoming work.

###### Acknowledgements.

P.F.L., A.S., and M.C. acknowledge funding from the European Research Council
(ERC) under the European Union’s Horizon 2020 research and innovation
programme (Grant agreement No. 863481). This work was supported by the
European Centre of Excellence in Exascale Computing TREX — Targeting Real
Chemical Accuracy at the Exascale. This project has received funding from the
European Union’s Horizon 2020 — Research and Innovation program — under grant
agreement no. 952165.

##  Appendix A Particular case of a $2\times 2$ matrix

For the simplest case of a system containing only two states, $\ket{1}$ and
$\ket{2}$, the fundamental equation given in Eq. (54) simplifies to

$\mathcal{I}=\matrixelement{I_{0}}{\quantity(H-E\mathds{1})^{-1}}{\Psi}=\matrixelement{I_{0}}{{\quantity[P_{0}\quantity(H-E\mathds{1})P_{0}]}^{-1}}{\Psi}+\sum_{p=1}^{\infty}\mathcal{I}_{p},$
(84)

with

$\mathcal{I}_{p}=\sum_{I_{1}\notin\mathcal{D}_{0},\ldots,I_{p}\notin\mathcal{D}_{p-1}}\quantity[\prod_{k=0}^{p-1}\matrixelement{I_{k}}{{\quantity[P_{k}\quantity(H-E\mathds{1})P_{k}]}^{-1}(-H)(1-P_{k})}{I_{k+1}}]\\\
\times\matrixelement{I_{p}}{{\quantity[P_{p}\quantity(H-E\mathds{1})P_{p}]}^{-1}}{\Psi}.$
(85)

To treat simultaneously the two possibilities for the final state, i.e.,
$\ket{i_{N}}=\ket{1}$ or $\ket{2}$, Eq. (54) has been slightly generalized to
the case of a general vector for the final state written as

$\ket{\Psi}=\Psi_{1}\ket{1}+\Psi_{2}\ket{2}.$ (86)

Let us choose a single-state domain for both states, namely
$\mathcal{D}_{1}=\quantity{\ket{1}}$ and $\mathcal{D}_{2}=\quantity{\ket{2}}$.
Note that, due to the simplicity of the present two-state model, there are
only two possible deterministic “alternating” paths, namely,
$\ket{1}\to\ket{2}\to\ket{1},\ldots$ and $\ket{2}\to\ket{1}\to\ket{2},\ldots$.

For the sake of convenience, we introduce the following quantities

$\displaystyle A_{1}$
$\displaystyle=\matrixelement{1}{{\quantity[P_{1}\quantity(H-E\mathds{1})P_{1}]}^{-1}(-H)(1-P_{1})}{2},$
(87) $\displaystyle A_{2}$
$\displaystyle=\matrixelement{2}{{\quantity[P_{2}\quantity(H-E\mathds{1})P_{2}]}^{-1}(-H)(1-P_{2})}{1},$
(88)

and

$\displaystyle C_{1}$
$\displaystyle=\matrixelement{1}{{\quantity[P_{1}\quantity(H-E\mathds{1})P_{1}]}^{-1}}{\Psi},$
(89) $\displaystyle C_{2}$
$\displaystyle=\matrixelement{2}{{\quantity[P_{2}\quantity(H-E\mathds{1})P_{2}]}^{-1}}{\Psi}.$
(90)

Without loss of generality, let us choose, for example, $\ket{I_{0}}=\ket{1}$.
Then, it is easy to show that

$\displaystyle\mathcal{I}_{2k+1}$
$\displaystyle=\frac{C_{2}}{A_{2}}(A_{1}A_{2})^{2k+1},$
$\displaystyle\mathcal{I}_{2k}$ $\displaystyle=C_{1}(A_{1}A_{2})^{2k},$ (91)

which yields

$\sum_{p=1}^{\infty}\mathcal{I}_{p}=\frac{C_{2}}{A_{2}}\sum_{p=1}^{\infty}(A_{1}A_{2})^{p}+C_{1}\sum_{p=1}^{\infty}(A_{1}A_{2})^{p}=A_{1}\frac{C_{2}+C_{1}A_{2}}{1-A_{1}A_{2}},$
(92)

and

$\matrixelement{1}{\quantity(H-E\mathds{1})^{-1}}{\Psi}=C_{1}+A_{1}\frac{C_{2}+C_{1}A_{2}}{1-A_{1}A_{2}}.$
(93)

For a $2\times 2$ matrix of the form

$H=\begin{pmatrix}H_{11}&H_{12}\\\ H_{12}&H_{22}\\\ \end{pmatrix},$ (94)

it is easy to evaluate the $A_{i}$’s. Using Eqs. (87) and (88), one gets, for
$i=1$ or $2$,

$\displaystyle A_{i}$ $\displaystyle=-\frac{H_{12}}{H_{ii}-E},$ $\displaystyle
C_{i}$ $\displaystyle=\frac{1}{H_{ii}-E}\Psi_{i},$ (95)

which finally yields

$\matrixelement{1}{\quantity(H-E\mathds{1})^{-1}}{\Psi}=\frac{H_{22}-E}{\Delta}\Psi_{1}-\frac{H_{12}}{\Delta}\Psi_{2},$
(96)

where $\Delta$ is the determinant of $H$.

Alternatively, the quantity
$\matrixelement{1}{\quantity(H-E\mathds{1})^{-1}}{\Psi}$ in the left-hand-side
of Eq. (96) can be directly calculated using the inverse of the matrix defined
in Eq. (94), yielding

$\begin{split}\quantity(H-E\mathds{1})^{-1}\ket{\Psi}&=\frac{1}{H-E\mathds{1}}\begin{pmatrix}\Psi_{1}\\\
\Psi_{2}\\\ \end{pmatrix}\\\
&=\frac{1}{\Delta}\begin{pmatrix}H_{22}-E&-H_{21}\\\ -H_{21}&H_{11}-E\\\
\end{pmatrix}\begin{pmatrix}\Psi_{1}\\\ \Psi_{2}\\\ \end{pmatrix}.\end{split}$
(97)

As readily seen, the first component of the vector
$\quantity(H-E\mathds{1})^{-1}\ket{\Psi}$ is identical to the one given in Eq.
(96), thus confirming independently the validity of this equation.

## References

  * Foulkes  _et al._ (2001) W. M. C. Foulkes, L. Mitas, R. J. Needs, and G. Rajagopal, Rev. Mod. Phys. 73, 33 (2001).
  * Kolorenc and Mitas (2011) J. Kolorenc and L. Mitas, Rep. Prog. Phys. 74, 026502 (2011).
  * Holzmann  _et al._ (2006) M. Holzmann, B. Bernu, and D. M. Ceperley, Phys. Rev. B 74, 104510 (2006).
  * Carlson  _et al._ (2015) J. Carlson, S. Gandolfi, F. Pederiva, S. C. Pieper, R. Schiavilla, K. E. Schmidt, and R. B. Wiringa, Rev. Mod. Phys. 87, 1067 (2015).
  * Carlson (2007) J. Carlson, Nucl. Physics. A 787, 516 (2007).
  * Austin  _et al._ (2012) B. M. Austin, D. Y. Zubarev, and W. A. Lester, Chem. Rev. 112, 263 (2012).
  * Golub and Van Loan (2012) G. H. Golub and C. F. Van Loan,  _Matrix Computations_ , 4th ed. (The Johns Hopkins University Press, 2012).
  * Davidson (1975) E. R. Davidson, J. Comput. Phys. 17, 87 (1975).
  * Caffarel and Claverie (1988) M. Caffarel and P. Claverie, J. Chem. Phys. 88, 1088 (1988).
  * Reynolds  _et al._ (1982) P. J. Reynolds, D. M. Ceperley, B. J. Alder, and W. A. Lester, J. Chem. Phys. 77, 5593 (1982).
  * Baroni and Moroni (1999) S. Baroni and S. Moroni, Phys. Rev. Lett. 82, 4745 (1999).
  * Sorella (1998) S. Sorella, Phys. Rev. Lett. 80, 4558 (1998).
  * Assaraf  _et al._ (2000) R. Assaraf, M. Caffarel, and A. Khelif, Phsys. Rev. E 61, 4566 (2000).
  * Assaraf and Caffarel (1999) R. Assaraf and M. Caffarel, Phys. Rev. Lett. 83, 4682 (1999).
  * Assaraf  _et al._ (1999) R. Assaraf, P. Azaria, M. Caffarel, and P. Lecheminant, Phys. Rev. B 60, 2299 (1999).
  * Caffarel and Assaraf (2000) M. Caffarel and R. Assaraf, in  _Lecture Notes in Chemistry_ , edited by M. Defranceschi and C. L. Bris (Springer, 2000) p. 45.
  * Kalos (1962) M. H. Kalos, Phys. Rev. 128, 1791 (1962).
  * Kalos (1970) M. H. Kalos, Phys. Rev. A 2, 250 (1970).
  * Kalos  _et al._ (1974) M. H. Kalos, D. Levesque, and L. Verlet, Phys. Rev. A 9, 2178 (1974).
  * Ceperley and Kalos (1979) D. Ceperley and M. Kalos,  _Monte Carlo Methods in Statistical Physics_ , edited by K.Binder (Springer, Berlin, 1979) Chap. 4.
  * Ceperley (1983) D. Ceperley, Journal of Computational Physics 51, 404 (1983).
  * Moskowitz and Schmidt (1986) J. W. Moskowitz and K. E. Schmidt, The Journal of Chemical Physics 85, 2868 (1986), https://doi.org/10.1063/1.451046 .
  * Petruzielo  _et al._ (2012) F. R. Petruzielo, A. A. Holmes, H. J. Changlani, M. P. Nightingale, and C. J. Umrigar, Phys. Rev. Lett. 109, 230201 (2012).
  * Garniron  _et al._ (2017) Y. Garniron, A. Scemama, P.-F. Loos, and M. Caffarel, J. Chem. Phys. 147, 034101 (2017).
  * Sharma  _et al._ (2017) S. Sharma, A. A. Holmes, G. Jeanmairet, A. Alavi, and C. J. Umrigar, J. Chem. Theory Comput. 13, 1595 (2017).
  * Willow  _et al._ (2012) S. Y. Willow, K. S. Kim, and S. Hirata, J. Chem. Phys. 137, 204122 (2012).
  * Caffarel (2019) M. Caffarel, J. Chem. Phys. 151, 064101 (2019).
  * Note (1) The property results from the fact that the series is a telescoping series and that the general term $\matrixelement{I}{\quantity(T^{+}_{I})^{n}}{\Psi^{+}}$ goes to zero as $n$ goes to infinity.
  * Note (2) As $\tau\rightarrow 0$ and $N\rightarrow\infty$ with $N\tau=t$, the operator $T^{N}$ converges to $e^{-t(H-E\mathds{1})}$. We then have $G^{E}_{ij}\rightarrow\intop\ilimits@_{0}^{\infty}dt\matrixelement{i}{e^{-t(H-E\mathds{1})}}{j}$, which is the Laplace transform of the time-dependent Green’s function $\matrixelement{i}{e^{-t(H-E\mathds{1})}}{j}$.
  * Gutzwiller (1963) M. C. Gutzwiller, Phys. Rev. Lett. 10, 159 (1963).
  * Huron  _et al._ (1973) B. Huron, J. P. Malrieu, and P. Rancurel, J. Chem. Phys. 58, 5745 (1973).
  * Harrison (1991) R. J. Harrison, J. Chem. Phys. 94, 5021 (1991).
  * Giner  _et al._ (2013) E. Giner, A. Scemama, and M. Caffarel, Can. J. Chem. 91, 879 (2013).
  * Holmes  _et al._ (2016) A. A. Holmes, N. M. Tubman, and C. J. Umrigar, J. Chem. Theory Comput. 12, 3674 (2016).
  * Schriber and Evangelista (2016) J. B. Schriber and F. A. Evangelista, J. Chem. Phys. 144, 161106 (2016).
  * Tubman  _et al._ (2020) N. M. Tubman, C. D. Freeman, D. S. Levine, D. Hait, M. Head-Gordon, and K. B. Whaley, J. Chem. Theory Comput. 16, 2139 (2020).
  * Booth  _et al._ (2009) G. H. Booth, A. J. W. Thom, and A. Alavi, J. Chem. Phys. 131, 054106 (2009).
  * Cleland  _et al._ (2010) D. Cleland, G. H. Booth, and A. Alavi, J. Chem. Phys. 132, 041103 (2010).
  * Zhang and Krakauer (2003) S. Zhang and H. Krakauer, Phys. Rev. Lett. 90, 136401 (2003).
  * White and Martin (1999) S. R. White and R. L. Martin, J. Chem. Phys. 110, 4127 (1999).
  * Chan and Sharma (2011) G. K.-L. Chan and S. Sharma, Annu. Rev. Phys. Chem. 62, 465 (2011).