\documentclass[aps,prb,reprint,showkeys,superscriptaddress]{revtex4-1}
%\usepackage{subcaption}
\usepackage{bm,graphicx,tabularx,array,booktabs,dcolumn,xcolor,microtype,multirow,amscd,amsmath,amssymb,amsfonts,physics,siunitx}
\usepackage[version=4]{mhchem}
\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\usepackage{txfonts,dsfont}
\usepackage{xspace}

\usepackage[normalem]{ulem}
\newcommand{\roland}[1]{\textcolor{cyan}{\bf #1}}
\newcommand{\manu}[1]{\textcolor{blue}{\bf #1}}
\newcommand{\vijay}[1]{\textcolor{green}{\bf #1}}
\newcommand{\titou}[1]{\textcolor{red}{#1}}
\newcommand{\toto}[1]{\textcolor{purple}{\bf #1}}
\newcommand{\mimi}[1]{\textcolor{orange}{\bf #1}}
\newcommand{\be}{\begin{equation}}
\newcommand{\ee}{\end{equation}}

\usepackage[
	colorlinks=true,
    citecolor=blue,
    linkcolor=blue,
    filecolor=blue,      
    urlcolor=blue,	
    breaklinks=true
	]{hyperref}
\urlstyle{same}

\newcommand{\ctab}{\multicolumn{1}{c}{---}}
\newcommand{\latin}[1]{#1}
%\newcommand{\latin}[1]{\textit{#1}}
\newcommand{\ie}{\latin{i.e.}}
\newcommand{\eg}{\latin{e.g.}}
\newcommand{\etal}{\textit{et al.}}

\newcommand{\mc}{\multicolumn}
\newcommand{\fnm}{\footnotemark}
\newcommand{\fnt}{\footnotetext}
\newcommand{\mcc}[1]{\multicolumn{1}{c}{#1}}
\newcommand{\mr}{\multirow}

% operators
\newcommand{\bH}{\boldsymbol{H}}
\newcommand{\bV}{\boldsymbol{V}}
\newcommand{\bh}{\boldsymbol{h}}
\newcommand{\bQ}{\boldsymbol{Q}}
\newcommand{\br}{\boldsymbol{r}}
\newcommand{\bp}{\boldsymbol{p}}
\newcommand{\cP}{\mathcal{P}}
\newcommand{\cS}{\mathcal{S}}
\newcommand{\cT}{\mathcal{T}}
\newcommand{\cC}{\mathcal{C}}
\newcommand{\cD}{\mathcal{D}}
\newcommand{\cE}{\mathcal{E}}
\newcommand{\cI}{\mathcal{I}}

\newcommand{\EPT}{E_{\PT}}
\newcommand{\laPT}{\lambda_{\PT}}

\newcommand{\EEP}{E_\text{EP}}
\newcommand{\laEP}{\lambda_\text{EP}}

\newcommand{\PsiT}{\Psi_\text{T}}
\newcommand{\PsiG}{\Psi^{+}}
\newcommand{\EL}{E_\text{L}}
\newcommand{\Id}{\mathds{1}}

\newcommand{\Ne}{N} % Number of electrons
\newcommand{\Nn}{M} % Number of nuclei
\newcommand{\hI}{\Hat{I}}
\newcommand{\hH}{\Hat{H}}
\newcommand{\hS}{\Hat{S}}
\newcommand{\hT}{\Hat{T}}
\newcommand{\hW}{\Hat{W}}
\newcommand{\hV}{\Hat{V}}
\newcommand{\hc}[2]{\Hat{c}_{#1}^{#2}}
\newcommand{\hn}[1]{\Hat{n}_{#1}}
\newcommand{\n}[1]{n_{#1}}
\newcommand{\Dv}{\Delta v}

\newcommand{\ra}{\to}

% Center tabularx columns
\newcolumntype{Y}{>{\centering\arraybackslash}X}

% HF rotation angles
\newcommand{\ta}{\theta^{\,\alpha}}
\newcommand{\tb}{\theta^{\,\beta}}
\newcommand{\ts}{\theta^{\sigma}}

% Some constants
\renewcommand{\i}{\mathrm{i}} % Imaginary unit
\newcommand{\e}{\mathrm{e}} % Euler number
\newcommand{\rc}{r_{\text{c}}}
\newcommand{\lc}{\lambda_{\text{c}}}
\newcommand{\lp}{\lambda_{\text{p}}}
\newcommand{\lep}{\lambda_{\text{EP}}}

% Some energies
\newcommand{\Emp}{E_{\text{MP}}}

% Blackboard bold
\newcommand{\bbR}{\mathbb{R}}
\newcommand{\bbC}{\mathbb{C}}

% Making life easier
\newcommand{\Lup}{\mathcal{L}^{\uparrow}}
\newcommand{\Ldown}{\mathcal{L}^{\downarrow}}
\newcommand{\Lsi}{\mathcal{L}^{\sigma}}
\newcommand{\Rup}{\mathcal{R}^{\uparrow}}
\newcommand{\Rdown}{\mathcal{R}^{\downarrow}}
\newcommand{\Rsi}{\mathcal{R}^{\sigma}}
\newcommand{\vhf}{\Hat{v}_{\text{HF}}}
\newcommand{\whf}{\Psi_{\text{HF}}}

\newcommand{\LCPQ}{Laboratoire de Chimie et Physique Quantiques (UMR 5626), Universit\'e de Toulouse, CNRS, UPS, France}
\newcommand{\LCT}{Laboratoire de Chimie Th\'eorique, Sorbonne-Universit\'e, Paris, France}

\begin{document}	
\title{Diffusion Monte Carlo using Domains in Configuration Space}
\author{Roland Assaraf}
        \email{assaraf@lct.jussieu.fr}
        \affiliation{\LCT}
\author{Emmanuel Giner}
        \email{giner@lct.jussieu.fr}
        \affiliation{\LCT}
\author{Vijay Gopal Chilkuri}
       \email{vijay.gopal.c@gmail.com}
       \affiliation{\LCT}
\author{Pierre-Fran\c{c}ois Loos}
        \email{loos@irsamc.ups-tlse.fr}
        \affiliation{\LCPQ}
\author{Anthony Scemama}
        \email{scemama@irsamc.ups-tlse.fr}
        \affiliation{\LCPQ}
\author{Michel Caffarel}
        \email{caffarel@irsamc.ups-tlse.fr}
        \affiliation{\LCPQ}
\begin{abstract}

\noindent
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 [\href{https://doi.org/10.1103/PhysRevB.60.2299}{Assaraf et al.~Phys.~Rev.~B \textbf{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.
\end{abstract}

\maketitle

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{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,\cite{Foulkes_2001,Kolorenc_2011} quantum liquids,\cite{Holzmann_2006} nuclear physics,\cite{Carlson_2015,Carlson_2007} theoretical chemistry,\cite{Austin_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 \textit{stochastic} power methods, a family of old and widely employed numerical approaches able to extract the largest or smallest eigenvalues of a matrix (see, \eg, Ref.~\onlinecite{Golub_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, \eg, the Lancz\`os algorithm (based on Krylov subspaces) \cite{Golub_2012} or Davidson's method where a diagonal preconditioning is performed. \cite{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, \cite{Caffarel_1988} DMC with branching, \cite{Reynolds_1982} reptation Monte Carlo, \cite{Baroni_1999} stochastic reconfiguration Monte Carlo, \cite{Sorella_1998,Assaraf_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 the 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 \cite{Assaraf_1999A}) 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,\cite{Assaraf_1999B,Caffarel_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 within 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. 
In their pioneering work, \cite{Kalos_1974} Kalos and collaborators introduced the so-called domain Green's function Monte Carlo approach in continuous space that they applied to a system of bosons with hard-sphere interaction. 
The domain used was the Cartesian product of small spheres around each particle, the Hamiltonian being approximated by the kinetic part only within the domain.
Some time later, Kalos proposed to extend these ideas to more general domains such as rectangular and/or cylindrical domains. \cite{Kalos_2000} In both works, the size of the domains is infinitely small in the limit of a vanishing time-step. 
Here, the domains are of arbitrary size, thus greatly increasing the efficiency of the approach.
Note also that some general equations for arbitrary domains in continuous space have also been proposed by one of us in Ref.~\onlinecite{Assaraf_1999B}.

Finally, from a general perspective, it is interesting to mention that the method proposed here is an example 
of how combining stochastic and deterministic techniques can be valuable.
%Here, a new stochastic process is built by taking advantage of a partial exact (deterministic) summation of local quantities.
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 {\it et al.},\cite{Petruzielo_2012} 
two different methods for evaluating stochastically the second-order perturbation energy of selected CI methods\cite{Garniron_2017b,Sharma_2017}, 
the stochastic approach of Willow {\it et al.} for the second-order many-body correction to the Hartree-Fock energy,\cite{Willow_2012} 
or the zero variance Monte Carlo scheme for evaluating the two-electron integrals of quantum chemistry.\cite{Caffarel_2019}

The paper is organized as follows. 
Section \ref{sec:DMC} presents the basic equations and notations of DMC. 
First, the path integral representation of the Green's function is given in Subsec.~\ref{sec:path}. 
The probabilistic framework allowing the Monte Carlo calculation of the Green's function is presented in Subsec.~\ref{sec:proba}. 
Section \ref{sec:DMC_domains} is devoted to the use of domains in DMC. 
First, we recall in Subsec.~\ref{sec:single_domains} the case of a domain consisting of a single state. \cite{Assaraf_1999B} 
The general case is then treated in Subsec.~\ref{sec:general_domains}. 
In Subsec.~\ref{sec:Green}, both the time- and energy-dependent Green's function using domains is derived. 
Section \ref{sec:numerical} presents the application of the approach to the one-dimensional Hubbard model. 
Finally, in Sec.\ref{sec:conclu}, some conclusions and perspectives are given.
Atomic units are used throughout.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Diffusion Monte Carlo}
\label{sec:DMC}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%=======================================%
\subsection{Path-integral representation}
\label{sec:path}
%=======================================%

As previously mentioned, DMC is a stochastic implementation of the power method defined by the following operator:
\be
	T = \Id -\tau (H-E\Id),
\ee 
where $\Id$ 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 $\braket{\Phi_0}{\Psi_0} \ne 0$ and for $\tau$ sufficiently small, we have
\be
\label{eq:limTN}
	\lim_{N \to \infty} T^N \ket{\Psi_0} = \ket{\Phi_0},
\ee
where $\ket{\Phi_0}$ is the ground-state wave function, \ie, $H \ket{\Phi_0} = E_0 \ket{\Phi_0}$.
The equality in Eq.~\eqref{eq:limTN} 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{\PsiT}$, as follows
\be
\label{eq:E0}
	E_0 = \lim_{N \to \infty } \frac{\mel{\Psi_0}{T^N}{H\Psi_T}}{\mel{\Psi_0}{T^N}{\Psi_T}}.
\ee
$\ket{\PsiT}$ 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
\be
	G^{(N)}_{ij}=\mel{i}{T^N}{j}.
\ee
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 \le k \le N-1}$, in $T^N$, $G^{(N)}$ can be written in the following expanded form
\be
\label{eq:cn}
	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}},
\ee
where $T_{ij} =\mel{i}{T}{j}$.
Here, each index $i_k$ runs over all basis vectors.

In quantum physics, Eq.~\eqref{eq:cn} 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:
\be
\label{eq:G}
	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}}.
\ee

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.\ref{fig1}.
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 the ``deterministic'' power methods, such as the Lancz\`os or Davidson approaches. 
If not, probabilistic techniques for generating only the paths contributing significantly to the sums are to be preferred.
This is the central theme of the present work.

%%% FIG 0 %%%
\begin{figure}
\includegraphics[width=\columnwidth]{fig1}
\caption{Path integral representation of the exact coefficient $c_i=\braket{i}{\Psi_0}$ of the ground-state wavefunction obtained as an infinite sum of paths starting 
from $\ket{i_0}$ and ending at $\ket{i}$, Eq.(\ref{eq:G}). Each path carries a weight $\prod_k T_{i_{k} i_{k+1}}$ computed along the path.
The result is independent of the choice of the initial state  $\ket{i_0}$, provided that $\braket{i_0}{\Psi_0}$ is non-zero.
Here, only four paths of infinite length have been represented. 
}
\label{fig1}
\end{figure}


%=======================================%
\subsection{Probabilistic framework}
\label{sec:proba}
%=======================================%

To derive a probabilistic expression for the Green's matrix, we introduce a guiding wave function, $\ket{\PsiG}$, having strictly positive components, \ie, $\PsiG_i > 0$, in order to perform a similarity transformation to the operators $G^{(N)}$ and $T$, as follows:
\begin{align}
\label{eq:defT}
	\bar{T}_{ij} & = \frac{\PsiG_j}{\PsiG_i} T_{ij},
	&
	\bar{G}^{(N)}_{ij}& = \frac{\PsiG_j}{\PsiG_i} G^{(N)}_{ij}.
\end{align}
Note that, thanks to the properties of similarity transformations, the path integral expression relating $G^{(N)}$ and $T$  [see Eq.~\eqref{eq:G}] remains unchanged for $\bar{G}^{(N)}$ and $\bar{T}$.

Next, the matrix elements of $\bar{T}$ are expressed as those of a stochastic matrix multiplied by some residual weight, namely,
\be
\label{eq:defTij}
	\bar{T}_{ij}= p_{i \to j} w_{ij}.
\ee
Here, we recall that a stochastic matrix is defined as a matrix with positive entries that obeys
\be
\label{eq:sumup}
	\sum_j p_{i \to j}=1.
\ee
To build the transition probability density, the following operator is introduced
%As known, there is a natural way of associating a stochastic matrix to a matrix having a positive ground-state vector (here, a positive vector is defined here as 
%a vector with all components positive). 
\be
\label{eq:T+}
	T^+= \Id - \tau \qty( H^+ - \EL^+ \Id ), 
\ee
where 
$H^+$ is the matrix obtained from $H$ by imposing the off-diagonal elements to be negative
\be
    H^+_{ij}=
    \begin{cases}
    	\phantom{-}H_{ij},	&	\text{if $i=j$},	
    	\\
    	-\abs{H_{ij}},		&	\text{if $i\neq j$}.
    \end{cases}
\ee
Here, $\EL^+ \Id$ is the diagonal matrix whose diagonal elements are defined as
\be
	(\EL^+)_{i}= \frac{\sum_j H^+_{ij}\PsiG_j}{\PsiG_i}.
\ee
The vector $\EL^+$ is known as the local energy vector associated with $\PsiG$.
By construction, the operator $H^+ - \EL^+ \Id$ in the definition of the operator $T^+$ [see Eq.~\eqref{eq:T+}] has been chosen to admit $\ket{\PsiG}$ as a ground-state wave function with zero eigenvalue, \ie, $\qty(H^+ - E_L^+ \Id) \ket{\PsiG} = 0$, leading to the relation
\be
\label{eq:relT+}
	T^+ \ket{\PsiG} = \ket{\PsiG}.
\ee

We are now in the position to define the stochastic matrix as
\be
\label{eq:pij}
	p_{i \to j} 
	= \frac{\PsiG_j}{\PsiG_i} T^+_{ij}
	= 
	\begin{cases}
		1 - \tau \qty[ H^+_{ii}- (\EL^+)_{i} ],					&	\qif* i=j,	
		\\
		\tau \frac{\PsiG_{j}}{\PsiG_{i}} \abs{H_{ij}}  \ge 0, 	&	\qif* i\neq j.
	\end{cases}
\ee
As readily seen in Eq.~\eqref{eq:pij}, 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 
\be
\label{eq:cond}
	\tau \leq \frac{1}{\max_i\abs{H^+_{ii}-(\EL^+)_{i}}}.
\ee
The sum-over-states condition [see Eq.~\eqref{eq:sumup}],
\be
	\sum_j p_{i \to j}= \frac{\mel{i}{T^+}{\PsiG}}{\PsiG_{i}} = 1,
\ee
follows from the fact that $\ket{\PsiG}$ is eigenvector of $T^+$, as evidenced by Eq.~\eqref{eq:relT+}.
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.~\eqref{eq:cond}] 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}-(\EL^+)_{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:
\be
	p_{i \to j} 
	= \frac{ \frac{\PsiG_{j}}{\PsiG_{i}} \abs{\mel{i}{T^+}{j}} } 
	{ \sum_j \frac{\PsiG_{j}}{\PsiG_{i}} \abs{\mel{i}{T^+}{j}} }
	= \frac{ \PsiG_{j} \abs*{T^+_{ij}} }
	{ \sum_j \PsiG_{j} \abs*{T^+_{ij}} }.
\ee
This new transition probability matrix with positive entries reduces to Eq.~\eqref{eq:pij} when $T^+_{ij}$ is positive as $\sum_j \PsiG_{j} T^+_{ij} = 1$.

Now, we need to make the connection between the transition probability matrix, $p_{i \to j}$, just defined from the operator $T^+$ corresponding 
to the approximate Hamiltonian 
$H^{+}$ and the operator $T$ associated with the exact Hamiltonian $H$. 
This is done thanks to the relation, Eq.(\ref{eq:defTij}), connecting $p_{i \to j}$ and $T_{ij}$ through a weight. 
Using Eqs.~\eqref{eq:defT} and \eqref{eq:pij}, the weights read
\be
	w_{ij}=\frac{T_{ij}}{T^+_{ij}}.
\ee
Using these notations the similarity-transformed Green's matrix components can be rewritten as 
\be
\label{eq:GN_simple}
	\bar{G}^{(N)}_{i_0 i_N} = 
	\sum_{i_1,\ldots,i_{N-1}} \qty( \prod_{k=0}^{N-1} p_{i_{k} \to i_{k+1}} ) \prod_{k=0}^{N-1} w_{i_{k}i_{k+1}},
\ee
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.~\eqref{eq:E0}. Taking $\ket{\Psi_0}=\ket{i_0}$ as initial condition, we have
\be
        E_0 = \lim_{N \to \infty }
        \frac{ \sum_{i_N} G^{(N)}_{i_0 i_N} {(H\PsiT)}_{i_N} }
        { \sum_{i_N} {G}^{(N)}_{i_0 i_N} {\PsiT}_{i_N} }.
\ee
which can be rewritten probabilistically as 
\be
        E_0 = \lim_{N \to \infty }
        \frac{ \expval{ \prod_{k=0}^{N-1} w_{i_{k}i_{k+1}} \frac{ {(H\PsiT)}_{i_N} }{ \PsiG_{i_N} }}}
        { \expval{ \prod_{k=0}^{N-1} w_{i_{k}i_{k+1}} \frac{ {\PsiT}_{i_N} } {\PsiG_{i_N}} }}.
\ee
where $\expval{...}$ is the probabilistic average defined over the set of paths $\ket{i_1},\ldots,\ket{i_N}$ occuring with probability
\be
 \text{Prob}_{i_0}(i_1,\ldots,i_{N}) = \prod_{k=0}^{N-1} p_{i_{k} \to i_{k+1}}
\ee
Using Eq.~\eqref{eq:sumup} and the fact that $p_{i \to j} \ge 0$, one can easily verify that $\text{Prob}_{i_0}$ is positive and obeys
\be
	\sum_{i_1,\ldots,i_{N}} \text{Prob}_{i_0}(i_1,\ldots,i_{N}) = 1,
\ee
as it should. 
%For a given path $i_1,\ldots,i_{N-1}$, the probabilistic average associated with this probability is then defined as
%\be
%\label{eq:average}
%	\expval{F} =  \sum_{i_1,\ldots,i_{N-1}} F(i_0,\ldots,i_N) \text{Prob}_{i_0 \to i_N}(i_1,\ldots,i_{N-1}),
%\ee
%where $F$ is an arbitrary function.
%Finally, the path-integral expressed as a probabilistic average reads
%\be
%\label{eq:cn_stoch}
%	\bar{G}^{(N)}_{i_0 i_N}= \expval{ \prod_{k=0}^{N-1} w_{i_{k}i_{k+1}}}.
%\ee
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}$. 
%In this framework, the energy defined in Eq.~\eqref{eq:E0} is given by
%\be
%	E_0 = \lim_{N \to \infty } 
%	\frac{ \expval{ \prod_{k=0}^{N-1} w_{i_{k}i_{k+1}} {(H\PsiT)}_{i_N}} }
%	{ \expval{ \prod_{k=0}^{N-1} w_{i_{k}i_{k+1}} {\PsiT}_{i_N} }}.
%\ee
%A schematic algorithm is presented in Fig.\ref{scheme1B}.
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 ergodic property of the stochastic matrix (see, for example, Ref.~\onlinecite{Caffarel_1988}), a unique path of infinite length from which sub-paths of length $N$ can be extracted may also be used. 
We shall not here insist on these practical details that are discussed, for example, in Refs.~\onlinecite{Foulkes_2001,Kolorenc_2011}.

%{\it Spawner representation} In this representation, we no longer consider moving particles but occupied or non-occupied states $|i\rangle$.
%To each state is associated the (positive or negative) quantity $c_i$.
%At iteration $n$ each state $|i\rangle$ with weight $c^n_i$
%"spawnes" (or deposits) the weight $T_{ij} c_i^{(n)}$ on a number of its connected sites $j$
%(that is, $T_{ij} \ne 0$). At the end of iteration the new weight on site $|i_n\rangle$  is given as the sum
%of the contributions spawn on this site by all states of the previous iteration
%\be
%c_i^{(n+1)}=\sum_{i_n} T_{ij} c^{(n)}_{i_n}
%\ee
%In the numerical applications to follow, we shall use the walker representation.

%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Diffusion Monte Carlo with domains}
\label{sec:DMC_domains}
%%%%%%%%%%%%%%%%%%%%%%%%%%%

%=======================================%
\subsection{Single-state domains}
\label{sec:single_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.\cite{Assaraf_1999A,Assaraf_1999B,Caffarel_2000} 

Considering a given state $\ket{i}$, the probability that a walker remains exactly $n$ times in $\ket{i}$ (with $1 \le n < \infty$) and then exits to a different state $j$ (with $j \neq i$) is
\be
	\cP_{i \to j}(n) = \qty(p_{i \to i})^{n-1}  p_{i \to j}.
\ee
Using the relation 
\be
	\sum_{n=1}^{\infty} (p_{i \to i})^{n-1}=\frac{1}{1-p_{i \to i}}
\ee
and the normalization of the $p_{i \to j}$'s [see Eq.~\eqref{eq:sumup}], one can check that the probability is properly normalized, \ie,
\be
	\sum_{j \ne i} \sum_{n=1}^{\infty} \cP_{i \to j}(n) = 1.
\ee

Naturally, the probability of being trapped during $n$ steps is obtained by summing over all possible exit states
\be
P_i(n)=\sum_{j \ne i} \cP_{i \to j}(n) = \qty(p_{i \to i})^{n-1} \qty( 1 - p_{i \to i} ),
\ee
and this defines a Poisson law with an average number of trapping events
\be
	\bar{n}_i= \sum_{n=1}^{\infty} n P_i(n) = \frac{1}{1 -p_{i \to i}}.
\ee
Introducing the continuous time $t_i = n \tau$,  the average trapping time is thus given by
\be
	\bar{t_i}= \qty[ H^+_{ii}-(\EL^+)_{i} ]^{-1},
\ee
and, in the limit $\tau \to 0$, the Poisson probability takes the usual form
\be
	P_{i}(t) = \frac{1}{\bar{t}_i} \exp(-\frac{t}{\bar{t}_i}).
\ee
The time-averaged contribution of the on-state weight can then be easily calculated to be
\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}}
\ee
Details of the implementation of this effective dynamics can be in found in Refs.~\onlinecite{Assaraf_1999B} and \onlinecite{Caffarel_2000}.

%=======================================%
\subsection{Multi-state domains}
\label{sec:general_domains}
%=======================================%

Let us now extend the results of Sec.~\ref{sec:single_domains} 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 $\cD_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 
\be
	\ket{i_0} \to \ket{i_1} \to \cdots \to \ket{i_N},
\ee
where the successive states are drawn using the transition probability matrix, $p_{i \to j}$. 
This series can be recast
\be
\label{eq:eff_series}
	(\ket*{I_0},n_0) \to  (\ket*{I_1},n_1) \to \cdots \to (\ket*{I_p},n_p),
\ee
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 \le n_0 \le N+1$), $\ket{I_1}$ is the first exit state that does not belong to $\cD_{i_0}$, $n_1$ is the number of times the walker remains in $\cD_{i_1}$ (with $1 \le n_1 \le N+1-n_0$), $\ket{I_2}$ is the second exit state, and so on.
Here, the integer $p$ (with $0 \le p \le 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 write the integers representing the exit states in capital letter, while small letters will be used for 
denoting the elementary states $\ket{i_k}$ generated with $p_{i \to j}$. Making this distinction is important since the effective 
stochastic dynamics used in practical Monte Carlo calculations will only involve exit states, the contribution of the elementary states, $\ket{i_k}$, being 
exactly integrated out.
Figure \ref{fig2} exemplifies how a path can be decomposed as proposed in Eq.(\ref{eq:eff_series}).
To make things as clear as possible, let us explicit in detail how the path drawn in Figure \ref{fig2} evolves in time.
The walker realizing the path starts at $\ket{i_0}$ within the domain $\cD_{i_0}$. It then makes two steps to arrive at $\ket{i_1}$, then $\ket{i_2}$ and, finally, leaves the domain 
at $\ket{i_3}$. The state $\ket{i_3}$ is the first exit state and is denoted as $\ket{I_1}(=\ket{i_3})$ following our convention of denoting exit states 
with capital letters. The trapping time in $\cD_{i_0}$ is $n_0=3$ since three states of the domain have been visited 
(namely, $\ket{i_0}$,$\ket{i_1}$,and $\ket{i_2}$).
During the next steps the domains $\cD_{I_1}$, $\cD_{I_2}$, and $\cD_{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. This work takes advantage of 
the fact that each possible path can be decomposed in this way.

%Generalizing what has been done for domains consisting of only one single state, the general idea here is to integrate out exactly the stochastic dynamics over the 
%set of all paths having the same representation, Eq.(\ref{eff_series}). As a consequence, an effective Monte Carlo dynamics including only exit states 
%averages for renormalized quantities will be defined.\\

%%% FIG 0B %%%
\begin{figure}
\includegraphics[width=\columnwidth,angle=0]{fig2}
\caption{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 with small black solid circles and the exit states, $\ket{I_k}$, with larger black solid squares. 
By convention the initial state, $\ket{i_0}$ is denoted using a capital letter, $\ket{I_0}$ since it will be the first state of the effective dynamics 
involving only exit states. 
See text, for some comments on the time evolution of the path.}
\label{fig2}
\end{figure}


Now, generalizing what has been done previously for a single-state domain, let us define the probability of remaining $n$ times in the domain of $\ket{I_0}$ and to exit at $\ket{I} \notin \cD_{I_0}$. This probability is given by
\be
\label{eq:eq1C}
	\cP_{I_0 \to I}(n) 
	= \sum_{\ket{i_1} \in \cD_{I_0}} \cdots \sum_{\ket{i_{n-1}} \in \cD_{I_0}} 
	p_{I_0 \to i_1} \ldots p_{i_{n-2} \to i_{n-1}} p_{i_{n-1} \to I}.
\ee
Since the sums are restricted to states belonging to the domain it is convenient to introduce a projector for over each domain
\be
\label{eq:pi}
	P_I = \sum_{\ket{i} \in \cD_I} \dyad{i}{i}.
\ee
and, also,
the projection of the operator $T^+$ over the domain, \ie,
\be
	T^+_I= P_I T^+ P_I.
\ee
Th operator $T^+_I$  governs the dynamics of the walkers trapped within the domain $\cD_{I}$,
see Eq.(\ref{eq:pij}) where $T^+$ is now restricted to the domain.
Using Eqs.~\eqref{eq:pij} and \eqref{eq:eq1C}, the probability can be rewritten as
\be
\label{eq:eq3C}
	\cP_{I_0 \to I}(n) = \frac{1}{\PsiG_{I_0}} \mel{I_0}{\qty(T^+_{I_0})^{n-1} F^+_{I_0}}{I} \PsiG_{I},
\ee
where the operator 
\be
\label{eq:Fi}
	F^+_I =  P_I T^+ (1-P_I),
\ee
corresponding to the last move connecting the inside and outside regions of the domain, has the following matrix elements:
 \be
	(F^+_I)_{ij} = 
	\begin{cases}
		T^+_{ij},	&	\qif* \ket{i} \in \cD_{I} \land \ket{j} \notin \cD_{I}, 
		\\
		0,			&	\text{otherwise}.
	\end{cases}
\ee
Physically, $F$ may be seen as a flux operator through the boundary of $\cD_{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 $\cD_{I}$, just by summing over all possible exit states. Thus, we get
\be
\label{eq:PiN}
	P_{I}(n) = \frac{1}{\PsiG_{I}} \mel{ I }{ \qty(T^+_{I})^{n-1} F^+_{I} }{ \PsiG }.
\ee
The normalization to one of this probability can be verified 
by using 
the fact that 
\be
\label{eq:relation}
	\qty(T^+_{I})^{n-1} F^+_I = \qty(T^+_{I})^{n-1} T^+ - \qty(T^+_I)^n,
\ee
leading to
\be
	\sum_{n=0}^{\infty} P_{I}(n) 
	= \frac{1}{\PsiG_{I}}  \sum_{n=1}^{\infty} \qty[ \mel{ I }{ \qty(T^+_{I})^{n-1} }{ \PsiG }
	-  \mel{ I }{ \qty(T^+_{I})^{n} }{ \PsiG } ] = 1.
\ee
The average trapping time defined as $t_{I}={\bar n}_{I} \tau$ where $ {\bar n}_{I}=\sum_n n P_{I}(n)$  is calculated to be
\be 
	t_{I}=\frac{1}{\PsiG_I}  \mel{I}{ { \qty[ P_I \qty( H^+ - \EL^+ \Id ) P_I ] }^{-1} }{ \PsiG }. 
\ee
In practice, the various quantities restricted to the domain will be computed by diagonalizing the matrix $(H^+-\EL^+ \Id)$ in $\cD_{I}$. 
Note that it is possible only if the dimension of the domains is not too large (say, less than a few thousands).

%=======================================%
\subsection{Time-dependent Green's matrix using domains}
\label{sec:Green}
%=======================================%

%--------------------------------------------%
%\subsubsection{Time-dependent Green's matrix}
%\label{sec:time}
%--------------------------------------------%
In this section we generalize the path-integral expression of the Green's matrix, Eq.~\eqref{eq:G}, to the case where domains are used.
To do so, we introduce the Green's matrix associated with each domain as follows:
\be
	G^{(N),\cD}_{ij}= \mel{ i }{ T_i^N }{ j }.
\ee
%Starting from Eq.~\eqref{eq:cn}
%\be
%G^{(N)}_{i_0 i_N}= \sum_{i_1,...,i_{N-1}} \prod_{k=0}^{N-1} \langle i_k|  T |i_{k+1} \rangle.
%\ee
Starting from Eq.~\eqref{eq:cn} and using the representation of the paths in terms of exit states and trapping times [see Eq.~\eqref{eq:eff_series}], we get
\be
	G^{(N)}_{i_0 i_N} = \sum_{p=0}^N 
	\sum_{\cC_p} \sum_{(i_1,...,i_{N-1}) \in \cC_p}
	\prod_{k=0}^{N-1} \mel{ i_k }{ T }{ i_{k+1} } ,
\ee
where $\cC_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 \cD_{k-1}$ (with $1 \le n_k \le N+1$ and $\sum_{k=0}^p n_k= N+1$).
It follows that 
\begin{multline}
\label{eq:Gt}
	G^{(N)}_{i_0 i_N}= G^{(N),\cD}_{i_0 i_N} +
	\sum_{p=1}^{N}
	\sum_{\ket{I_1} \notin \cD_{I_0}, \ldots , \ket{I_p} \notin \cD_{I_{p-1}}  }
	\sum_{n_0 \ge 1} \cdots \sum_{n_p \ge 1}
	\delta_{\sum_{k=0}^p n_k,N+1} 
	\\
	\times 
	\qty[ \prod_{k=0}^{p-1} \mel{ I_k }{ \qty(T_{I_k})^{n_k-1} F_{I_k} }{ I_{k+1} } ]
	G^{(n_p-1),\cD}_{I_p i_N},
\end{multline}
where $\delta_{i,j}$ is a Kronecker delta. Note that the first contribution, $G^{(N),\cD}_{i_0 i_N}$, corresponds to the case
$p=0$ and collects all contributions to the Green's matrix coming from paths remaining for ever within the domain of $\ket{i_0}$ (no exit event). 
This contribution is here isolated from the $p$-sum since, as a domain Green's matrix, it will be calculated exactly and will not be suject 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 representation of the paths (if not, it would be 
associated to the next contribution coresponding to $p+1$ exit events).

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, Eq.~\eqref{eq:G}, 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 \le k \le 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} \mel{ I_k }{ F_{I_k} }{ I_{k+1} }$, with $F=T$, thus recovering Eq.~\eqref{eq:G}.

In order to compute $G^{(N)}_{i_0 i_N}$ by resorting to  Monte Carlo techniques, let us reformulate Eq.~\eqref{eq:Gt}
using the transition probability $\cP_{I \to J}(n)$ introduced above,
Eq.~\eqref{eq:eq3C}. We first rewrite Eq.~\eqref{eq:Gt} under the form
\begin{multline}
	{G}^{(N)}_{i_0 i_N}={G}^{(N),\cD}_{i_0 i_N} + {\PsiG_{i_0}}
	\sum_{p=1}^{N}
	\sum_{\ket{I_1} \notin \cD_{I_0}, \ldots , \ket{I_p} \notin \cD_{I_{p-1}}}
	\sum_{n_0 \ge 1} \cdots \sum_{n_p \ge 1}
	\\
	\times
	\delta_{\sum_{k=0}^p n_k,N+1} \qty{ \prod_{k=0}^{p-1} \qty[ \frac{\PsiG_{I_{k+1}}}{\PsiG_{I_k}} \mel{ I_k }{ \qty(T_{I_k})^{n_k-1} F_{I_k} }{ I_{k+1} 	} ] }
	\frac{{G}^{(n_p-1),\cD}_{I_p i_N}}{\PsiG_{I_p}}.
\end{multline}
Introducing the weights 
\be
	W_{I_k I_{k+1}} = \frac{\mel{ I_k }{ \qty(T_{I_k})^{n_k-1} F_{I_k} }{ I_{k+1} }}{\mel{ I_k }{ \qty(T^{+}_{I_k})^{n_k-1} F^+_{I_k} }{ I_{k+1} }},
\ee
and using the effective transition probability, Eq.~\eqref{eq:eq3C}, we get
\begin{multline}
\label{eq:Gbart}
	{G}^{(N)}_{i_0 i_N} = {G}^{(N),\cD}_{i_0 i_N} 
\\
+ {\PsiG_{i_0}}\sum_{p=1}^{N} \sum_{{(I,n)}_{p,N}} 
	{
	\qty( \prod_{k=0}^{p-1} W_{I_k I_{k+1}} ) \qty( \prod_{k=0}^{p-1}\cP_{I_k \to I_{k+1}}(n_k)  )
	 \frac{1}{\PsiG_{I_p}} {G}^{(n_p-1), \cD}_{I_p i_N} },
\end{multline}
where, for clarity, $\sum_{{(I,n)}_{p,N}}$ is used as a short-hand notation for the sum,
$ \sum_{\ket{I_1} \notin \cD_{I_0}, \ldots , \ket{I_p} \notin \cD_{I_{p-1}}}
        \sum_{n_0 \ge 1} \cdots \sum_{n_p \ge 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, $\cP_{I \to J}(n)$. For example, in the case of the energy, we start from
\be
        E_0 = \lim_{N \to \infty }
        \frac{ \sum_{i_N} {G}^{(N)}_{i_0 i_N} {(H\PsiT)}_{i_N} }
        { \sum_{i_N} {G}^{(N)}_{i_0 i_N}  {\PsiT}_{i_N} }
\ee
which can be rewritten probabilistically as
\be
        E_0 = \lim_{N \to \infty }
        \frac{ {G}^{(N),\cD}_{i_0 i_N} + {\PsiG_{i_0}}  \sum_{p=1}^{N}  \expval{ \qty( \prod_{k=0}^{p-1} W_{I_k I_{k+1}} ) {\cal H}_{n_p,I_p} }_p}
             { {G}^{(N),\cD}_{i_0 i_N} + {\PsiG_{i_0}} \sum_{p=1}^{N}  \expval{ \qty( \prod_{k=0}^{p-1} W_{I_k I_{k+1}} ) {\cal S}_{n_p,I_p} }_p}
\ee
where $\expval{...}_p$ is the probabilistic average defined over the set of paths $p$ exit events of probability
$\prod_{k=0}^{p-1} \cP_{I_k \to I_{k+1}}(n_k) $ 
and $({\cal H}_{n_p,I_p},{\cal S}_{n_p,I_p})$ two quantities taking into account the contribution  of the trial wave function at the end of the path and defined as follows
\be
{\cal H}_{n_p,I_p}=   \frac{1}{\PsiG_{I_p}} \sum_{i_N} {G}^{(n_p-1),\cD}_{I_p i_N} (H \Psi_T)_{i_N} 
\ee
and
\be
{\cal S}_{n_p,I_p}=   \frac{1}{\PsiG_{I_p}} \sum_{i_N} {G}^{(n_p-1), \cD}_{I_p i_N} (\Psi_T)_{i_N} 
\ee
In practice, the Monte Carlo algorithm is a simple generalization of that used in standard diffusion Monte Carlo calculations.
Stochastic paths starting at $\ket{I_0}$ are generated using 
the probability $\cP_{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 
products of weights, $ \prod_{k=0}^{p-1} W_{I_k I_{k+1}} $ times the end-point contributions, ${{(\cal H}/{\cal 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.~\eqref{eq:PiN}]
and, then,
the exit state, $\ket{I_{k+1}}$, is drawn using the conditional probability $\frac{\cP_{I_k \to I_{k+1}}(n_k)}{P_{I_k}(n_k)}$.
%See fig.\ref{scheme1C}.
%\titou{i) Choose some initial vector $\ket{I_0}$\\
%ii) Generate a stochastic path by running over $k$ (starting at $k=0$) as follows.\\
%$\;\;\;\bullet$ Draw $n_k$ using the probability $P_{I_k}(n)$ [see Eq.~\eqref{eq:PiN}]\\
%$\;\;\;\bullet$ Draw the exit state, $\ket{I_{k+1}}$, using the conditional probability $$\frac{\cP_{I_k \to I_{k+1}}(n_k)}{P_{I_k}(n_k)}$$\\
%iii) Terminate the path when $\sum_k n_k=N$ is greater than some target value $N_\text{max}$ and compute $F(I_0,n_0;\ldots;I_N,n_N)$\\
%iv) Go to step ii) until some maximum number of paths is reached.\\
%\\
%%At the end of the simulation, an estimate of the average for a few values of $N$ greater but close to $N_\text{max}$ is obtained. 
%At large $N_\text{max}$ where the convergence of the average as a function of $p$ is reached, such values can be averaged.}
%
%%==============================================%
\subsection{Domain Green's Function Monte Carlo}
\label{sec:energy}
%==============================================%

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.


%--------------------------------------------%
\subsubsection{The pedestrian way} 
%--------------------------------------------%
Let us define the energy-dependent Green's matrix
\be
	G^E_{ij}= \tau \sum_{N=0}^{\infty} \mel{ i }{ T^N }{ j} = \mel{i}{ \qty( H-E \Id )^{-1} }{j}.
\ee
The denomimation ``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.\cite{note}
The remarkable property is that, thanks to the summation over $N$ up to infinity, the constrained multiple sums appearing in Eq.~\eqref{eq:Gt} can be factorized in terms of a product of unconstrained single sums, as follows
\begin{multline}
	\sum_{N=1}^\infty \sum_{p=1}^N \sum_{n_0 \ge 1} \cdots \sum_{n_p \ge 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).
\end{multline}
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.~\eqref{eq:Gt} and summing over the $n_k$-variables, we get
\begin{multline}
\label{eq:eqfond}
  {G}^{E}_{i_0 i_N}
	= {G}^{E,\cD}_{i_0 i_N}
	+ \sum_{p=1}^{\infty} \sum_{I_1 \notin \cD_0, \hdots , I_p \notin \cD_{p-1}} \\
	\qty[ \prod_{k=0}^{p-1} \mel{ I_k }{ {\qty[ P_k \qty( H-E \Id ) P_k ] }^{-1} (-H)(\Id-P_k) }{ I_{k+1} } ]
       {G}^{E,\cD}_{I_p i_N}
\end{multline}
where, ${G}^{E,\cD}$ is the energy-dependent domain's Green matrix defined as ${G}^{E,\cD}_{ij} = \tau \sum_{N=0}^{\infty} \mel{ i }{ T^N_i }{ j}$.

As a didactical example, Appendix \ref{app:A} reports the exact derivation of this formula in the case of a two-state system.

%----------------------------%
\subsubsection{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
\be
	\qty(H-E\Id)^{-1} = \qty(H_0-E\Id)^{-1} + \qty(H_0-E\Id)^{-1} (H_0-H) \qty(H-E\Id)^{-1},
\ee
where $H_0$ is some arbitrary reference Hamiltonian, we have the Dyson equation
\be
\label{eq:GE}
	G^E_{ij} = G^E_{0,ij} + \sum_{kl} G^{E}_{0,ik} (H_0-H)_{kl} G^E_{lj},
\ee
with $G^E_{0,ij} = \mel{i}{ \qty( H_0-E \Id )^{-1} }{j}$.
Let us choose $H_0$ such that $\mel{ i }{ H_0 }{ j } = \mel{ i }{ P_i H P_i }{ j }$ for all $i$ and $j$.
Then, the Dyson equation \eqref{eq:GE} becomes
\begin{multline}
{G}^{E}_{ij}
	= {G}^{E,\cD}_{ij}
	\\
	+ \sum_k  \mel{ i }{ {\qty[P_i \qty(H-E \Id)P_i ]}^{-1} (H_0-H) }{ k } {G}^{E}_{kj}.
\end{multline}
Using the following identity
\be
\begin{split}
{\qty[	P_i \qty(H-E \Id) P_i ]}^{-1} (H_0-H) 
	& = { \qty[ P_i \qty(H-E \Id)P_i ]}^{-1} (P_i H P_i - H) 
	\\
	& = {\qty[ P_i \qty(H-E \Id) P_i]}^{-1} (-H) ( \Id -P_i),
\end{split}
\ee
the Dyson equation may be written under the form
\begin{multline}
 G^E_{ij} =  {G}^{E,\cD}_{ij}
	\\
	+ \sum_{k \notin \cD_i}  \mel{ i }{ {\qty[P_i \qty(H-E \Id) P_i]}^{-1}  (-H)( \Id -P_i) }{ k }  {G}^{E}_{kj}
\end{multline}
which is identical to Eq.~\eqref{eq:eqfond} when $G^E_{ij}$ is expanded iteratively.
Let us use as effective transition probability density
\be
	P_{I \to J} = \frac{1} {\PsiG(I)} \mel{ I }{ {\qty[ P_I \qty(H^+ - \EL^+ \Id) P_I]}^{-1} (-H^+) (\Id -P_I) }{ J } \PsiG(J)
\ee
and the weight
\be
	W^E_{IJ} =
	\frac{ \mel{ I }{ {\qty[ P_I \qty(H-E \Id) P_I]}^{-1} (-H)( \Id -P_I) }{ J} }
	{\mel{ I }{ {\qty[ P_I \qty(H^+ - \EL^+ \Id) P_I]}^{-1} (-H^+)( \Id -P_I) }{ J} }
\ee
Using Eqs.~\eqref{eq:eq1C}, \eqref{eq:eq3C} and \eqref{eq:relation}, one can easily verify that $P_{I \to J} \ge 0$ and $\sum_J P_{I \to J}=1$. 
Finally, the probabilistic expression reads
\begin{multline}
\label{eq:final_E}
G^E_{i_0 i_N}= {G}^{E,\cD}_{i_0 i_N}
	+ \sum_{p=1}^{\infty} \PsiG_{i_0} \expval{ \qty( \prod_{k=0}^{p-1} W^E_{I_k I_{k+1}} ) 
\frac{{G}^{E,\cD}_{I_p i_N}} { \PsiG_{I_p}} }
\end{multline}
% HERE

%----------------------------%
\subsubsection{Energy estimator} 
%----------------------------%
To calculate the energy, we introduce the following estimator 
\be
	\cE(E) = \frac{ \mel{ I_0 }{ \qty(H-E \Id)^{-1} }{ H\PsiT } } {\mel{ I_0 }{ \qty(H-E \Id)^{-1} }{ \PsiT } }.
\ee
which is the counterpart of the quantity ${\cal E}_N =\frac{ \mel{ I_0 }{ T^N }{ H\PsiT } } {\mel{ I_0 }{T^N}{ \PsiT } }$ 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.(\ref{eq:E0}).
Here, the situation is not as simple and we need a way to extract the energy from $\cE(E)$.

Using the spectral decomposition of $H$, we have
\be
\label{eq:calE}
        \cE(E) = \frac{ \sum_i \frac{E_i c_i}{E_i-E}}{\sum_i \frac{c_i}{E_i-E}},
\ee
where $c_i = \braket{ I_0 }{ \Phi_i } \braket{ \Phi_i}{  \PsiT }$ and $\Phi_i$ are the eigenstates of $H$, \ie, $H \ket{\Phi_i} = E_i \ket{\Phi_i}$.
The important observation is that for all eigenstates we have $\cE(E_i)= E_i$. Thus, to get the ground-state energy we propose to search for the solution 
of the non-linear equation $\cE(E)= E$ in the vicinity of $E_0$.

In practical Monte Carlo calculations the DMC energy will be obtained by computing a finite number of components $H_p$ and $S_p$ defined as follows
\be
\cE^\text{DMC}(E,p_{max})= \frac{ H_0+ \sum_{p=1}^{p_{max}} H^\text{DMC}_p }
                         {S_1 +\sum_{p=1}^{p_{max}} S^\text{DMC}_p }
\ee
For $ p\ge 1$, Eq~\eqref{eq:final_E} gives
\begin{align}
        H^\text{DMC}_p & = \PsiG_{i_0}\expval{ \qty(\prod_{k=0}^{p-1} W^E_{I_k I_{k+1}}) {\cal H}_{I_p} }
        \\
        S^\text{DMC}_p & = \PsiG_{i_0} \expval{ \qty(\prod_{k=0}^{p-1} W^E_{I_k I_{k+1}}) {\cal S}_{I_p} }.
\end{align}
where
\be
{\cal H}_{I_p}=   \frac{1}{\PsiG_{I_p}} \sum_{i_N} {G}^{(E,\cD}_{I_p i_N} (H \Psi_T)_{i_N}
\ee
and
\be
{\cal S}_{I_p}=   \frac{1}{\PsiG_{I_p}} \sum_{i_N} {G}^{E, \cD}_{I_p i_N} (\Psi_T)_{i_N}
\ee
For $p=0$, the two components are exactly evaluated as
\begin{align}
        H_0 & = \mel{ I_0 }{ {\qty[ P_{I_0} \qty(H-E \Id) P_{I_0} ]}^{-1} }{ H\PsiT } 
        \\
        S_0 & = \mel{ I_0 }{ {\qty[ P_{I_0} \qty(H-E \Id) P_{I_0} ]}^{-1} }{ \PsiT }.
\end{align}
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.
To avoid the calculation by Monte Carlo of a quantity, such as $H_0$ or $S_0$, which can be evaluated analytically is computationally 
very appealing (the statistical error is removed). We thus propose to extend 
the exact calculation of $H_p$ and $S_p$ to higher values of $p$,
up to the point where the exponential increase 
of the number of intermediate states involved in the explicit sums makes the calculation unfeasible. 
%\begin{multline}
%H_p=   \sum_{I_1 \notin \cD_0, \hdots , I_p \notin \cD_{p-1}}
%        \qty[ \prod_{k=0}^{p-1} \mel{ I_k }{ {\qty[ P_k \qty( H-E \Id ) P_k ] }^{-1} (-H)(\Id-P_k) }{ I_{k+1} } ]
%\\
%\times \mel{ I_p }{ {\qty[ P_p ( H-E \Id) P_p ] }^{-1} }{ I_N }
%\end{multline}
%with a similar formula for $S_p$. In practice, $\cE^\text{DMC}(E,p_{ex},p_{max})$

Finally, $\cE^\text{DMC}(E,p_{ex},p_{max})$
is evaluated in practice as follows
\be
\cE^\text{DMC}(E,p_{ex},p_{max})= \frac{   \sum_{p=0}^{p_{ex}-1} H_p +\sum_{p=p_{ex}}^{p_{max}} H^\text{DMC}_p }
                         {                 \sum_{p=0}^{p_{ex}-1} S_p +\sum_{p=p_{ex}}^{p_{max}} S^\text{DMC}_p }
\ee
where $p_{ex}$ is the number of components exactly computed. 
Let us emphasize that, since the magnitude of $H_p$ and $S_p$ decreases as a function of $p$, to remove the statistical error 
for the first most important contributions can lead 
to important gains, as illustrated in the numerical application to follow.

It is easy to check that, in the vicinity of the exact energy $E_0$, $\cE(E)$ is a linear function of $E - E_0$.
Therefore, in practice, we compute the value of $\cE(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 $\cE(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.~\eqref{eq:calE} 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. 


%Let us describe the 3 functions used here for the fit.
%i) Linear fit\\
%We write
%\be
%\cE(E)= a_0 + a_1 E
%\ee
%and search the best value of $(a_0,a_1)$ by fitting the data.
%Then, $\cE(E)=E$ leads to
%$$
%E_0= \frac{a_0}{1-a_1}
%$$
%ii) Quadratic fit\\
%At the quadratic level we write
%$$
%\cE(E)= a_0 + a_1 E + a_2 E^2
%$$
%leading to 
%$$
%E_0 = \frac{1 - a_1 \pm \sqrt{ (a_1 -1)^2 - 4 a_0 a_2}}{2 a_2}
%$$
%iii) Two-component fit\\
%We take the advantage that the exact expression of $\cE(E)$ is known as an infinite series, Eq.(\ref{calE}).
%If we limit ourselves to the first two-component, we write
%\be
%\cE(E) = \frac{ \frac{E_0 c_0}{E_0-E} + \frac{E_1 c_1}{E_1-E}}{\frac{c_0}{E_0-E} + \frac{c_1}{E_1-E} }
%\ee
%Here, the variational parameters used for the fit are ($c_0,E_0,c_1,E_1)$.\\


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Numerical application to the Hubbard model}
\label{sec:numerical}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%======================================
\subsection{Hamiltonian and trial wave function}
%======================================

Let us consider the one-dimensional Hubbard Hamiltonian for a chain of $N$ sites 
\be
	H= -t \sum_{\expval{ i j } \sigma} \hat{a}^+_{i\sigma} \hat{a}_{j\sigma} 
	+ U \sum_i \hat{n}_{i\uparrow} \hat{n}_{i\downarrow},
\ee
where $\expval{ i j }$ 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
\be
	\ket{n} = \ket{n_{1 \uparrow},\ldots,n_{N \uparrow},n_{1 \downarrow},\ldots,n_{N \downarrow}},
\ee
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 \ge 0$). 
It is then possible to equate the guiding and trial vectors, \ie, $\ket{c^+} = \ket{c_\text{T}}$.
As trial wave function, we shall employ a generalization of the Gutzwiller wave function \cite{Gutzwiller_1963}
\be
	\braket{ n }{ c_\text{T} }= e^{-\alpha n_D(n)-\beta n_A(n)},
\ee
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
\be
	n_A(n)= \sum_{\expval{ i j }} n_{i\uparrow} n _{j\downarrow} \delta_{n_{i\uparrow},1} \delta_{n_{j\downarrow},1}.
\ee
The parameters $\alpha$ and $\beta$ are optimized by minimizing the variational energy of the trial wave function, \ie, 
\be
	E_\text{T}(\alpha,\beta) = \frac{\mel{ c_\text{T} }{H }{c_\text{T}}} {\braket{ c_\text{T} }{ c_\text{T} }}.
\ee

%======================================
\subsection{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 kinetic energy. 
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 diagonalizations to perform, we shall consider only one non-trivial domain called here the main domain and denoted as $\cD$. 
This domain will be chosen common to all states belonging to it, that is,
\be
	\cD_i= \cD \qq{for} \ket{i} \in \cD.
\ee
For the other states, we choose a single-state domain as
\be
	\cD_i= \qty{ \ket{i} } \qq{for} \ket{i} \notin \cD.
\ee
To define $\cD$, let us introduce the following set of states
\be
	{\cS_{ij} = \qty{ \ket{n} : n_D(n)=i \land n_A(n)=j }}.
\ee
which means that $\cD$ contains 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)$ will not be varied and taken to be the maximum possible for a given $m$, $n_A^\text{max}(m)= \max_{N}(N-1-2m,0)$. 
Using these definitions, the main domain is taken as the union of some elementary domains
\be
	\cD = \bigcup_{n_D=0}^{n_D^\text{max}}\cD(n_D,n_A^\text{min}(n_D)),
\ee
where the elementary domain is defined as
\be
	\cD(n_D,n_A^\text{min}(n_D))=\bigcup_{ j = n_A^\text{min}(n_D)}^{n_A^\text{max}(n_D)} \cS_{n_D j}.
\ee
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 
\begin{align*}
	\cD(0,3) & = \cS_{03}, 
	&
	\cD(0,2) & = \cS_{03} \cup \cS_{02}, 
	\\
	\cD(0,1) & = \cS_{03} \cup \cS_{02} \cup \cS_{01},
	&
	\cD(1,1) & = \cS_{11},
	\\
	\cD(1,0) & = \cS_{11} \cup \cS_{10},
	&
	\cD(2,0) & = \cS_{20}, 
\end{align*}
where
\begin{align*}
	\cS_{03} & = \qty{ \ket{\uparrow,\downarrow,\uparrow,\downarrow }, \ket{\downarrow,\uparrow,\downarrow,\uparrow } }, 
	\qq{(the two N\'eel states)}
\\
	\cS_{02} & = \qty{ \ket{\uparrow, \downarrow, \downarrow, \uparrow }, \ket{\downarrow, \uparrow, \uparrow, \downarrow } },
\\
	\cS_{01} & = \qty{ \ket{\uparrow, \uparrow, \downarrow, \downarrow }, \ket{\downarrow, \downarrow, \uparrow, \uparrow } },
\\
	\cS_{11} & = \qty{ \ket{\uparrow \downarrow, \uparrow ,\downarrow, 0 }, \ket{\uparrow \downarrow, 0, \uparrow,\uparrow } + \ldots },
\\
	\cS_{10} & = \qty{ \ket{\uparrow \downarrow, \uparrow, 0, \downarrow }, \ket{\uparrow \downarrow, 0, \uparrow, \downarrow } + \ldots },
\\
	\cS_{20} & = \qty{ \ket{\uparrow \downarrow, \uparrow \downarrow, 0 ,0 } + \ldots }.
\end{align*}
For the three last cases, the dots indicate the remaining states obtained by permuting the position of the pairs.

%======================================
\subsection{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 of $E_\text{T}=-0.495361\ldots$. 
In what follows $\ket{I_0}$ will be systematically chosen as one of the two N\'eel states, \eg, $\ket{I_0} = \ket{\uparrow,\downarrow, \uparrow,\ldots}$. 

Figure \ref{fig3} 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_{ex}=0$). 
At the scale of the figure, error bars are too small to be seen.
When $E$ is far from the exact value of $E_0=-0.768068\ldots$, 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. 
Note the oscillations of the curves as a function of $p$ 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, $\cE^\text{DMC}(E,p_{ex}=1,p_{max})$ as a function of $E$ is presented in Fig.~\ref{fig4}. 
Here, $p_{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:
\be
\label{eq:2comp}
	\cE(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} },
\ee
which considers only the first two terms of the exact expression of $\cE(E)$ [see Eq.~\eqref{eq:calE}].
%If we limit ourselves to the first two-component, we write
%\be
%\cE(E) = \frac{ \frac{E_0 c_0}{E_0-E} + \frac{E_1 c_1}{E_1-E}}{\frac{c_0}{E_0-E} + \frac{c_1}{E_1-E} }
%\ee
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 $\cE(E)=E$ is $-0.76807(5)$ in full agreement with the exact value of $-0.768068\ldots$.

%%% FIG 1 %%%
\begin{figure}
\includegraphics[width=\columnwidth]{fig3}
\caption{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.}
\label{fig3}
\end{figure}

%%% FIG 2 %%%
\begin{figure}
\includegraphics[width=\columnwidth]{fig4}
\caption{One-dimensional Hubbard model with $N=4$ and $U=12$. 
$\cE(E)$ as a function of $E$.
The horizontal and vertical lines are at $\cE(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.~\eqref{eq:2comp}.
Error bars are smaller than the symbol size.}
\label{fig4}
\end{figure}

Table \ref{tab1} 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 using a minimal 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_{conv}$, defined as the smallest value of $p$ for which the energy is converged with six decimal places. 
The smaller $p_{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 $\cD(0,3) \cup \cD(1,0)$ has a trapping time for the N\'eel state of $6.2$, while the domain $\cD(0,3) \cup \cD(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, $\cE(E=-1)=-0.75272390\ldots$, can be found at the last row of Table \ref{tab1} 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, nicely illustrating the critical importance of the domains employed.

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

As explained above, it is very advantageous to calculate exactly as many $(H_p,S_p)$ as possible in order to avoid the sttistical error on the heaviest 
components. 
Table \ref{tab2} shows the results both for the case of a single-state main domain and for the domain having the largest average trapping time, namely $\cD(0,1) \cup \cD(1,1)$ (see Table \ref{tab1}). 
Table \ref{tab2} reports the statistical fluctuations of the energy for the simulation of Table \ref{tab1}. 
Results show that it is indeed interesting to compute exactly as many components as possible. 
For the single-state domain, a factor 2 of reduction of the statistical error is obtained when passing from no analytical computation, $p_{ex}=0$, to the case where eight components for $H_p$ and $S_p$ are exactly computed. 
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 \ref{tab3} reports the energies converged as a function of $p$ with 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.~\ref{fig4}. 
As seen on the figure the behavior of $\cE$ as a function of $E$ is very close to the linearity. 
The extrapolated values obtained from the five values of the energy with the three fitting functions are reported. 
Using the linear fitting function leads to an energy of $-0.7680282(5)$ to compare 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 with the quadratic fit. 
It is also true when using the two-component function. 
Our final value is in full agreement with the exact value with about six decimal places.

Table \ref{tab4} shows the evolution of the average trapping times and extrapolated energies as a function of $U$ when using $\cD(0,1) \cup \cD(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 \ge 12$) the increase of $\bar{t}_{I_0}$ becomes particularly steep.

Finally, in Table \ref{tab4} we report the results obtained for larger systems at $U=12$ for a size of the chain 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 (not 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 domains chosen 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 called for.
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 for 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 size increases. 
To get lower error bars we need to use better trial wave functions, 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 II %%%
\begin{table}
\caption{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 Table \ref{tab1}. 
Results are presented when a single-state domain is used for all states and when $\cD(0,1) \cup \cD(1,0)$ is chosen as the main domain.}
\label{tab2}
\begin{ruledtabular}
\begin{tabular}{lcc}
$p_{ex}$ & single-state & $\cD(0,1) \cup \cD(1,0)$ \\
\hline
$0$ & $4.3 \times 10^{-5}$ &$ 347 \times 10^{-8}$ \\
$1$ & $4.0  \times10^{-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  \times10^{-5}$ &$ 6.0  \times 10^{-8}$\\
$6$ & $2.3  \times10^{-5}$ &$ 0.7  \times 10^{-8}$\\
$7$ & $2.2 \times 10^{-5}$ &$ 0.6  \times 10^{-8}$\\
$8$ & $2.2  \times10^{-5}$ &$  0.05 \times 10^{-8}$\\
\end{tabular}
\end{ruledtabular}
\end{table}

%%% TABLE III %%%
\begin{table}
\caption{One-dimensional Hubbard model with $N=4$, $U=12$, $\alpha=1.292$, $\beta=0.552$, and $p_{ex}=4$. 
The main domain is $\cD(0,1) \cup \cD(1,0)$. 
The simulation is performed with 20 independent blocks and $10^6$ stochastic paths starting from the N\'eel state. 
The fits are performed with the five values of $E$.}
\label{tab3}
\begin{ruledtabular}
\begin{tabular}{lc}
$E$ & $E^\text{DMC}$ \\
\hline
$-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$\\
\end{tabular}
\end{ruledtabular}
\end{table}

%%% TABLE IV %%%
\begin{table}
\caption{One-dimensional Hubbard model with $N=4$.
The main domain is $\cD(0,1) \cup \cD(1,0)$.}
\label{tab4}
\begin{ruledtabular}
\begin{tabular}{rcccr}
$U$ & $(\alpha,\beta)$ & $E_\text{T}$ & $E_{ex}$  & $\bar{t}_{I_0}$ \\
\hline
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  \\
\end{tabular}
\end{ruledtabular}
\end{table}

%%% TABLE V %%%
\begin{table*}
\caption{One-dimensional Hubbard model with $U=12$. 
The fits to extrapolate the DMC energies are done using the two-component function}
\label{tab5}
\begin{ruledtabular}
\begin{tabular}{crcrccccr}
$N$ & Hilbert space size & Domain & Domain size & $(\alpha,\beta)$ &$\bar{t}_{I_0}$ & $E_\text{T}$ & $E_{ex}$ & $E^\text{DMC}$\\
\hline
4  & 36     & $\cD(0,1) \cup \cD(1,0)$ & 30 &(1.292,\,0.552)& 108.7  & $-0.495361$ & $-0.768068$ & $-0.7680676(5)$\\
6  & 400    & $\cD(0,1) \cup \cD(1,0)$ &200 &(1.124,\,0.689)&57.8 & $-0.633297$ & $-1.215395$&  $-1.215389(9)$\\
8  & 4 900   & $\cD(0,1) \cup \cD(1,0)$ & 1 190 &(0.984,\,0.788)& 42.8 & $-0.750995 $& $-1.66395$&  $-1.6637(2)$\\
10 & 63 504  & $\cD(0,5) \cup \cD(1,4)$ & 2 682 &(0.856,\,0.869)& 31.0  & $-0.855958$  & $-2.113089$& $-2.1120(7)$\\
12 & 853 776 & $\cD(0,8) \cup \cD(1,7)$ & 1 674 &(0.739,\,0.938)& 16.7 & $-0.952127$ & $-2.562529$& $-2.560(6)$\\
\end{tabular}
\end{ruledtabular}
\end{table*}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Summary and perspectives}
\label{sec:conclu}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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 domain Green's functions 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, \cite{Huron_1973,Harrison_1991,Giner_2013,Holmes_2016,Schriber_2016,Tubman_2020}, FCIQMC, \cite{Booth_2009,Cleland_2010}, AFQMC, \cite{Zhang_2003} or DMRG. \cite{White_1999,Chan_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.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\acknowledgments{
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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Particular case of a $2\times2$ matrix}
\label{app:A}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

For the simplest case of a system containing only two states, $\ket{1}$ and $\ket{2}$, the fundamental equation given in Eq.~\eqref{eq:eqfond} simplifies to
\be
	\cI
	= \mel{ I_0 }{ \qty(H-E \Id)^{-1} }{ \Psi } 
	=  \mel{ I_0 }{ { \qty[ P_0 \qty(H-E \Id) P_0 ]}^{-1} }{ \Psi }
	+ \sum_{p=1}^{\infty} \cI_p,
\ee
with
\begin{multline}
	\cI_p= \sum_{I_1 \notin \cD_0, \ldots , I_p \notin \cD_{p-1}}
	\qty[ \prod_{k=0}^{p-1} \mel{ I_k }{ { \qty[P_k \qty(H-E \Id) P_k ]}^{-1} (-H)(1-P_k) }{ I_{k+1} } ]
	\\
	\times \mel{ I_p }{ { \qty[ P_p \qty(H-E \Id) P_p ]}^{-1} }{ \Psi }.
\end{multline}
To treat simultaneously the two possibilities for the final state, \ie, $\ket{I_N} = \ket{1}$ or $\ket{2}$, Eq.~\eqref{eq:eqfond} has been slightly generalized to the case of a general vector for the final state written as
\be
	\ket{\Psi} = \Psi_1 \ket{1} +  \Psi_2 \ket{2}.
\ee

Let us choose a single-state domain for both states, namely $\cD_1 = \qty{ \ket{1} }$  and $\cD_2 = \qty{ \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
\begin{align}
\label{eq:defA1}
	A_1 & = \mel{ 1 }{ {\qty[ P_1 \qty(H-E \Id) P_1 ]}^{-1} (-H)(1-P_1) }{ 2 },
	\\
\label{eq:defA2}
	A_2 & = \mel{ 2 }{ {\qty[ P_2 \qty(H-E \Id) P_2]}^{-1} (-H) (1-P_2) }{ 1 }, 
\end{align}
and
\begin{align}
	C_1 & = \mel{ 1 }{ {\qty[ P_1 \qty(H-E \Id) P_1]}^{-1}  }{ \Psi },
	\\
	C_2 & = \mel{ 2 }{ {\qty[ P_2 \qty(H-E \Id) P_2]}^{-1} }{ \Psi }.
\end{align}
Without loss of generality, let us choose, for example, $\ket{I_0} = \ket{1}$.
Then, it is easy to show that
%\begin{align}
%	\cI_1 & = A_1 C_2, 
%	& 
%	\cI_2 & = A_1 A_2 C_1, 
%	&
%	\cI_3 & = A_1 A_2 A_1 C_2, 
%	&
%	\qq{etc.}
%\end{align}
%that is
\begin{align}
	\cI_{2k+1} & = \frac{C_2}{A_2} (A_1 A_2)^{2k+1},
	& 
	\cI_{2k} & = C_1 (A_1 A_2)^{2k},
\end{align}
which yields
\be
	\sum_{p=1}^{\infty}
	\cI_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},
\ee
and 
\be
	\mel{ 1 }{ \qty(H-E \Id)^{-1} }{ \Psi } 
	= C_1 +  A_1 \frac{C_2 + C_1 A_2}{1-A_1 A_2}.
\ee
For a $2\times2$ matrix of the form
\be
	H = 
	\begin{pmatrix}
	H_{11}	&	H_{12}	\\
	H_{12}	&	H_{22}	\\
	\end{pmatrix},
\ee
it is easy to evaluate the $A_i$'s.
Using Eqs.~\eqref{eq:defA1} and \eqref{eq:defA2}, one gets, for $i = 1$ or $2$,
\begin{align}
	A_i & = -\frac{H_{12}}{H_{ii}-E},
	&
	C_i & = \frac{1}{H_{ii}-E} \Psi_i.
\end{align}
which finally yields
\be
	\mel{ 1 }{ \qty(H-E \Id)^{-1} }{ \Psi} 
%	= \qty(H_{11}-E)^{-1} \Psi_1 
%	- H_{12} \frac{ \qty( \Psi_2 - \frac{ H_{12}\Psi_1}{H_{11}-E} )} 
%                                          { (H_{11}-E)(H_{22}-E) - H^2_{12}}
%	\\
	= \frac{ H_{22}-E}{\Delta}  \Psi_1 - \frac{H_{12}}{\Delta} \Psi_2,
\label{final}
\ee
where $\Delta$ is the determinant of $H$. 

On the other hand, the quantity $\mel{ 1 }{ \qty(H-E \Id)^{-1} }{ \Psi}$ of the L.H.S of Eq.(\ref{final}) can be directly calculated using the inverse of the $2\times2$ matrix 
\be
{ \qty(H-E \Id)^{-1} }{ |\Psi\rangle}=
              \frac{1}{H-E \Id}      \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}
\ee
As seen, the first component of the vector ${ \qty(H-E \Id)^{-1} }{ |\Psi\rangle}$ is identical to the one given in 
Eq.~\eqref{final}, thus confirming independently the validity of this equation.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\bibliography{g}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\end{document}
