\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}{\bf #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}{\mathbf{H}}
\newcommand{\bV}{\mathbf{V}}
\newcommand{\bh}{\mathbf{h}}
\newcommand{\bQ}{\mathbf{Q}}
\newcommand{\bSig}{\mathbf{\Sigma}}
\newcommand{\br}{\mathbf{r}}
\newcommand{\bp}{\mathbf{p}}
\newcommand{\cP}{\mathcal{P}}
\newcommand{\cS}{\mathcal{S}}
\newcommand{\cT}{\mathcal{T}}
\newcommand{\cC}{\mathcal{C}}
\newcommand{\PT}{\mathcal{PT}}

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

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


\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}{\rightarrow}

% 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 \titou{state} obeys a Poisson law and that the \titou{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, although it is usually not the most efficient variant of DMC. \titou{Why?}
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_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 a previous work,\cite{assaraf_99,caffarel_00} 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 \titou{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 cnosidered.

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.
Finally, note that some general equations for arbitrary domains in continuous space have also been proposed in Ref.~\onlinecite{Assaraf_1999B}.

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_1999} 
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.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\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 = \mathds{1} -\tau (H-E\mathds{1}),
\ee 
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. Starting from some initial vector, $\ket{\Psi_0}$, we have
\be
\lim_{N \to \infty} T^N \ket{\Psi_0} = \ket{\Phi_0}
\ee
where $\ket{\Phi_0}$ is the ground-state wave function. 
The equality is up to a global phase factor playing no role in physical quantum averages. 
This result is true for any $\ket{\Psi_0}$ provided that $\braket{\Phi_0}{\Psi_0} \ne 0$ and for $\tau$ sufficiently small.
At large but finite $N$, the vector $T^N \ket{\Psi_0}$ differs from $\ket{\Phi_0}$ only by an exponentially small correction, making easy the extrapolation of the finite-N results to $N=\infty$.\\

Ground-state properties may be obtained at large $N$. For example, in the important case of the energy one can use the formula
\be
E_0 = \lim_{N\rightarrow \infty } \frac{\langle \Psi_T|H T^N|\Psi_0 \rangle}
                               {\langle \Psi_T|T^N|\Psi_0 \rangle}
\label{E0}
\ee
where $|\Psi_T\rangle$ is some trial vector (some approximation of the true ground-state) on which $T^N|\Psi_0 \rangle$ is projected out.

To proceed further we introduce the time-dependent Green's matrix $G^{(N)}$ defined as
\be
G^{(N)}_{ij}=\langle j|T^N |i\rangle.
\ee
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 convenient notation, $i_k$, for the $N-1$ indices of the intermediate states in the $N$-th product of $T$, $G^{(N)}$ can be written in 
the expanded form
\be
 G^{(N)}_{i_0 i_N} = \sum_{i_1} \sum_{i_2} ... \sum_{i_{N-1}} \prod_{k=0}^{N-1} T_{i_{k} i_{k+1}}.
\label{cn}
\ee
Here, each index $i_k$ runs over all basis vectors.\\

In quantum physics, Eq.(\ref{cn}) is referred to as
the path-integral representation of the Green's matrix (function).
The series of states $|i_0\rangle,...,|i_N\rangle$ is interpreted as a "path" in the Hilbert space starting at vector $|i_0\rangle$ and ending at 
vector $|i_N\rangle$ where $k$ plays the role of 
a time index. To each path is associated the weight $\prod_{k=0}^{N-1} T_{i_{k} i_{k+1}}$ 
and the path integral expression of $G$ is written in the more suggestive form
\be
 G^{(N)}_{i_0 i_N}= \sum_{\rm all \; paths\; |i_1\rangle,...,|i_{N-1}\rangle} \prod_{k=0}^{N-1} T_{i_{k} i_{k+1}}
\label{G}
\ee
This expression allows a simple and vivid interpretation of the solution: In the $N \rightarrow \infty$-limit 
the $i$-th component of the ground-state (obtained as $\lim_{N\rightarrow \infty} G^{(N)}_{i i_0})$ is the weighted sum over all possible paths arriving 
at vector $|i\rangle$. This result is independent of the initial
vector $|i_0\rangle$, apart from some irrelevant global phase factor.
When the size of the linear space is small the explicit calculation of the full sums involving $M^N$ terms (here, $M$ is the size of the Hilbert space)
can be performed. We are in the realm of what can be called 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 used.

\subsection{Probabilistic framework}
\label{sec:proba}
In order to derive 
a probabilistic expression for the Green's matrix we introduce a so-called guiding vector, $|\Psi^+\rangle$,
having strictly positive components, $\Psi^+_i > 0$, and apply a similarity transformation to the operators $G^{(N)}$ and $T$
\be
{\bar T}_{ij}= \frac{\Psi^+_j}{\Psi^+_i} T_{ij}
\label{defT}
\ee
and
\be
{\bar G}^{(N)}_{ij}= \frac{\Psi^+_j}{\Psi^+_i} G^{(N)}_{ij}
\ee
Note that under the similarity transformation the path integral expression, Eq.(\ref{G}), relating $G^{(N)}$ and $T$ remains 
unchanged for the similarity-transformed operators, ${\bar G}^{(N)}$ and ${\bar T}$.

Next, the matrix elements of ${\bar T}$ are expressed as those of a stochastic matrix (or transition probability 
matrix) multiplied by some residual weight, namely
\be
{\bar T_{ij}}= p(i \rightarrow j) w_{ij} 
\label{defTij}
\ee
Here, we recall that a stochastic matrix is defined as a matrix with positive entries and obeying
\be
\sum_j p(i \rightarrow j)=1.
\label{sumup}
\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
T^+=\mathds{1} - \tau [ H^+-{\rm diag}(E_L^+)]
\ee
where 
$H^+$ is the matrix obtained from $H$ by imposing the off-diagonal elements to be negative
\be
H^+_{ii}=H_{ii} \;\;\;{\rm and} \;\;\;  H^+_{ij}=-|H_{ij}| \;\;\; {\rm for} \;\;\;i \ne j.
\ee
Here, ${\rm diag}(E_L^+)$ is the diagonal matrix whose diagonal elements are defined as
\be
E^+_{Li}= \frac{\sum_j H^+_{ij}\Psi^+_j}{\Psi^+_i}.
\ee
The vector $|E^+_L\rangle$ is known as the local energy vector associated with $|\Psi^+\rangle$.\\

Actually, the operator $H^+-diag(E^+_L)$ in the definition of the operator $T^+$ has been chosen to admit by construction $|\Psi^+ \rangle$ as ground-state with zero eigenvalue
\be
[H^+ - {\rm diag}(E_L^+)]|\Psi^+\rangle=0,
\label{defTplus}
\ee
leading to the relation
\be
T^+ |\Psi^+\rangle=|\Psi^+\rangle.
\label{relT+}
\ee

The stochastic matrix is now defined as
\be
p(i \rightarrow j) = \frac{\Psi^+_j}{\Psi^+_i} T^+_{ij}.
\label{pij}
\ee
The diagonal matrix elements of the stochastic matrix write
\be
p(i \rightarrow i) = 1 -\tau (H^+_{ii}- E^+_{Li})
\ee
while for $i \ne j$
\be
p(i \rightarrow j) = \tau \frac{\Psi^+_{j}}{\Psi^+_{i}} |H_{ij}|  \ge 0  \;\;\; i \ne j.
\ee
As seen, the off-diagonal terms, $p(i \rightarrow j)$ are positive while
the diagonal ones, $p(i \rightarrow i)$, can be made positive if $\tau$ is chosen
sufficiently small. More precisely, the condition writes
\be
\tau \leq \frac{1}{{\rm Max}_i| H^+_{ii}-E^+_{Li}|}
\label{cond}
\ee
The sum-over-states condition, Eq.(\ref{sumup}), follows from the fact that $|\Psi^+\rangle$ is eigenvector of $T^+$, Eq.(\ref{relT+})
\be
\sum_j p(i \rightarrow j)= \frac{1}{\Psi^+_{i}} \langle i |T^+| \Psi^ +\rangle =1.
\ee
We have then verified that $p(i \rightarrow j)$ is indeed a stochastic matrix.\\

At first sight, the condition defining the maximum value of $\tau$ allowed, Eq.(\ref{cond}), may appear as rather tight 
since for very large matrices it may impose an extremely small value for the time step. However, in practice during the simulation only a (tiny) 
fraction of the linear space is sampled, and the maximum value of $|H^+_{ii} -E^+_{Li}|$ for the sampled states turns out to be not too large, so reasonable values of $\tau$ 
can be used without violating the positivity of the transition probability matrix.
Note that we can even escape from this condition by slightly generalizing the transition probability
matrix used as follows
\be
p(i \rightarrow j) = \frac{ \frac{\Psi^+_{j}}{\Psi^+_{i}} |\langle i | T^+ | j\rangle| } { \sum_j \frac{\Psi^+_{j}}{\Psi^+_{i}} |\langle i | T^+ | j\rangle|}
= \frac{ \Psi^+_{j} |\langle i | T^+ | j\rangle| }{\sum_j \Psi^+_{j} |\langle i | T^+ | j\rangle|}
\ee
This new transition probability matrix with positive entries reduces to Eq.(\ref{pij}) when $T^+_{ij}$ is positive.\\

Now, using Eqs.(\ref{defT},\ref{defTij},\ref{pij}) the residual weight $w_{ij}$ writes
\be
w_{ij}=\frac{T_{ij}}{T^+_{ij}}.
\ee
Using these notations the Green's matrix components can be rewritten as 
\be
{\bar G}^{(N)}_{i i_0}=\sum_{i_1,..., i_{N-1}} \Big[ \prod_{k=0}^{N-1} p(i_{k} \rightarrow i_{k+1})\Big] \prod_{k=0}^{N-1} w_{i_{k}i_{k+1}} 
\ee
where $i$ is identified to $i_N$.\\

The product $\prod_{k=0}^{N-1} p(i_{k} \rightarrow i_{k+1})$ is the probability, denoted ${\rm Prob}_{i_0 \rightarrow i_N}(i_1,...,i_{N-1})$, 
for the path starting at $|i_0\rangle$ and ending at $|i_N\rangle$ to occur.
Using the fact that $p(i \rightarrow j) \ge 0$ and Eq.(\ref{sumup}) we verify that ${\rm Prob}_{i_0 \rightarrow i_N}$ is positive and obeys
\be
\sum_{i_1,..., i_{N-1}} {\rm Prob}_{i_0 \rightarrow i_N}(i_1,...,i_{N-1})=1
\ee
as it should be. 
The probabilistic average associated with this probability for the path, denoted here as, $ \Big \langle  ... \Big \rangle$ is then defined as
\be
\Big \langle F \Big \rangle =  \sum_{i_1,..., i_{N-1}} F(i_0,...,i_N) {\rm Prob}_{i_0 \rightarrow i_N}(i_1,...,i_{N-1}),
\label{average}
\ee
where $F$ is an arbitrary function.
Finally, the path-integral expressed as a probabilistic average writes
\be
{\bar G}^{(N)}_{ii_0}= \Big \langle \prod_{k=0}^{N-1} w_{i_{k}i_{k+1}}  \Big \rangle
\label{cn_stoch}
\ee
To calculate the probabilistic average, Eq.(\ref{average}),
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 \rightarrow i_{k+1})$, thus realizing the path of probability 
${\rm Prob}(i_0 \rightarrow i_n)$. 
The energy, Eq.(\ref{E0}) is given as
\be
E_0 = \lim_{N \rightarrow \infty } \frac{ \Big \langle \prod_{k=0}^{N-1} w_{i_{k}i_{k+1}} {(H\Psi_T)}_{i_N} \Big \rangle}
                                        { \Big \langle \prod_{k=0}^{N-1} w_{i_{k}i_{k+1}} {\Psi_T}_{i_N}      \Big \rangle}
\ee
Note that, instead of using a single walker, it is possible to introduce a population of independent walkers and to calculate the averages over the population.
In addition, thanks to the ergodic property of the stochastic matrix (see, Refs \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 can be 
found, 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{DMC with domains}
\label{sec:DMC_domains}
\subsection{Domains consisting of a single state}
\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\cite{assaraf_99,Assaraf_1999B,caffarel_00} and applied to the SU(N) one-dimensional Hubbard model.

Let us consider a given state $|i\rangle$. The probability that the walker remains exactly $n$ times on $|i\rangle$ ($n$ from                                      
1 to $\infty$) and then exits to a different state $j$ is
\be
{\cal P}(i \rightarrow j, n) = [p(i \rightarrow i)]^{n-1}  p(i \rightarrow j) \;\;\;\; j \ne i.
\ee
Using the relation $\sum_{n=1}^{\infty} p^{n-1}(i \rightarrow i)=\frac{1}{1-p(i \rightarrow i)}$ and the normalization 
of the $p(i \rightarrow j)$, Eq.(\ref{sumup}), we verify that
the probability is normalized to one 
\be
\sum_{j \ne i} \sum_{n=1}^{\infty} {\cal P}(i \rightarrow j,n) = 1.
\ee

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} {\cal P}(i \rightarrow j,n) = [p(i \rightarrow i)]^{n-1}(1-p(i \rightarrow i)).
\ee
This probability defines a Poisson law
with an  average number $\bar{n}_i$ of trapping events given by
\be
\bar{n}_i= \sum_{n=1}^{\infty} n P_i(n) = \frac{1}{1 -p(i \rightarrow i)}.
\ee
Introducing the continuous time $t_i=n_i\tau$ the average trapping time is given by
\be
\bar{t_i}= \frac{1}{H^+_{ii}-E^+_{Li}}.
\ee
Taking the limit $\tau \rightarrow 0$ the Poisson probability takes the usual form
\be
P_{i}(t) = \frac{1}{\bar{t}_i} e^{-\frac{t}{\bar{t}_i}}
\ee
The time-averaged contribution of the on-state weight can 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 the effective dynamics can be in found in Refs. (\onlinecite{assaraf_99},\onlinecite{caffarel_00}).
\subsection{General domains}
\label{sec:general_domains}
Let us now extend the results of the preceding section to a general domain. For that, 
let us associate to each state $|i\rangle$ a set of states, called the domain of $|i\rangle$ and
denoted ${\cal D}_i$, consisting of the state $|i\rangle$ plus a certain number of states. No particular constraints on the type of domains 
are imposed, for example domains associated with different states can be identical, or may have or not 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
|i_0 \rangle \rightarrow |i_1 \rangle \rightarrow ... \rightarrow |i_N \rangle 
\ee
where the successive states are drawn using the transition probability matrix, $p(i \rightarrow j)$. This series can be rewritten 
\be
(|I_0\rangle,n_0) \rightarrow  (|I_1 \rangle,n_1) \rightarrow... \rightarrow (|I_p\rangle,n_p) 
\label{eff_series}
\ee
where $|I_0\rangle=|i_0\rangle$ is the initial state, 
$n_0$ the number of times the walker remains within the domain of $|i_0\rangle$ ($n_0=1$ to $N+1$), $|I_1\rangle$ is the first exit state, 
that is not belonging to 
${\cal D}_{i_0}$, $n_1$ is the number of times the walker remains within ${\cal D}_{i_1}$ ($n_1=1$ to $N+1-n_0$), $|I_2\rangle$ the second exit state, and so on.
Here, the integer $p$ goes from 0 to $N$ and indicates the number of exit events occurring along the path. The two extreme cases, $p=0$ and $p=N$, 
correspond to the cases where the walker remains for ever within the initial domain, and to the case where the walker leaves the current domain at each step, 
respectively. 
In what follows, we shall systematically write the integers representing the exit states in capital letter.

%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.\\

Let us define the probability of being $n$ times within the domain of $|I_0\rangle$ and, then, to exit at $|I\rangle \notin {\cal D}_{I_0}$. 
It is given by
$$
{\cal P}(I_0 \rightarrow I,n)= \sum_{|i_1\rangle \in {\cal D}_{I_0}} ... \sum_{|i_{n-1}\rangle \in {\cal D}_{I_0}} 
$$
\be
p(I_0 \rightarrow i_1) ... p(i_{n-2} \rightarrow i_{n-1}) p(i_{n-1} \rightarrow I)
\label{eq1C}
\ee
To proceed we need to introduce the projector associated with each domain
\be
P_I=  \sum_{|k\rangle \in {\cal D}_I}  |k\rangle \langle k|
\label{pi}
\ee
and to define the restriction of the operator $T^+$ to the domain
\be
T^+_I= P_I T^+ P_I.
\ee
$T^+_I$ is the operator governing the dynamics of the walkers moving within ${\cal D}_{I}$.
Using Eqs.(\ref{eq1C}) and (\ref{pij}), the probability can be rewritten as
\be
{\cal P}(I_0 \rightarrow I,n)=
\frac{1}{\Psi^+_{I_0}} \langle I_0 | {T^+_{I_0}}^{n-1} F^+_{I_0}|I\rangle \Psi^+_{I}
\label{eq3C}
\ee
where the operator $F$, corresponding to the last move connecting the inside and outside regions of the 
domain, is given by
\be
F^+_I =  P_I T^+ (1-P_I),
\label{Fi}
\ee
that is, $(F^+_I)_{ij}= T^+_{ij}$ when $(|i\rangle \in {\cal D}_{I}, |j\rangle \notin {\cal D}_{I})$, and zero
otherwise.
Physically, $F$ may be seen as a flux operator through the boundary of ${\cal D}_{I}$.

Now, the probability of being trapped $n$ times within ${\cal D}_{I}$ is given by
\be
P_{I}(n)=
\frac{1}{\Psi^+_{I}} \langle I | {T^+_{I}}^{n-1} F^+_{I}|\Psi^+ \rangle.
\label{PiN}
\ee
Using the fact that 
\be
{T^+_I}^{n-1} F^+_I= {T^+_I}^{n-1} T^+ - {T^+_I}^n 
\label{relation}
\ee
we have
\be
\sum_{n=0}^{\infty} P_{I}(n) = \frac{1}{\Psi^+_{I}}  \sum_{n=1}^{\infty} \Big( \langle I | {T^+_{I}}^{n-1} |\Psi^+\rangle 
-  \langle I | {T^+_{I}}^{n} |\Psi^+\rangle \Big) = 1
\ee
and the average trapping time
\be 
t_{I}={\bar n}_{I} \tau=  \frac{1}{\Psi^+_{I}}  \langle I | P_{I} \frac{1}{H^+ -E_L^+}  P_{I} | \Psi^+\rangle 
\ee
In practice, the various quantities restricted to the domain are computed by diagonalizing the matrix $(H^+-E_L^+)$ in ${\cal D}_{I}$. Note that 
it is possible only if the dimension of the domains is not too large (say, less than a few thousands).
\subsection{Expressing the Green's matrix using domains}
\label{sec:Green}
\subsubsection{Time-dependent Green's matrix}
\label{time}
In this section we generalize the path-integral expression of the Green's matrix, Eqs.(\ref{G}) and (\ref{cn_stoch}), to the case where domains are used.
For that we introduce the Green's matrix associated with each domain 
\be
G^{(N),{\cal D}}_{IJ}= \langle J| T_I^N| I\rangle.
\ee
Starting from Eq.(\ref{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
and using the representation of the paths in terms of exit states and trapping times we write
\be
G^{(N)}_{I_0 I_N} = \sum_{p=0}^N 
\sum_{{\cal C}_p} \sum_{(i_1,...,i_{N-1}) \in \;{\cal C}_p}
\prod_{k=0}^{N-1} \langle i_k|T |i_{k+1} \rangle 
\ee
where ${\cal C}_p$ is the set of paths with $p$ exit states, $|I_k\rangle$, and trapping times $n_k$ with the
constraints that $|I_k\rangle \notin {\cal D}_{k-1}$, $1 \le n_k \le N+1$ and $\sum_{k=0}^p n_k= N+1$.
We then have
$$
G^{(N)}_{I_0 I_N}= G^{(N),{\cal D}}_{I_0 I_N} +
$$
\be
\sum_{p=1}^{N}
\sum_{|I_1\rangle \notin {\cal D}_{I_0}, \hdots , |I_p\rangle \notin {\cal D}_{I_{p-1}}  }
\sum_{n_0 \ge 1} ... \sum_{n_p \ge 1}
\ee
\be
\delta(\sum_{k=0}^p n_k=N+1) \Big[ \prod_{k=0}^{p-1} \langle I_k|T^{n_k-1}_{I_k} F_{I_k} |I_{k+1} \rangle \Big]
G^{(n_p-1),{\cal D}}_{I_p I_N}
\label{Gt}
\ee
This expression is the path-integral representation of the Green's matrix using only the variables $(|I_k\rangle,n_k)$ of the effective dynamics defined over the set
of domains. The standard formula derived above, Eq.(\ref{G}) may be considered as the particular case where the domain associated with each state is empty, 
In that case, $p=N$ and $n_k=1$ for $k=0$ to $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} \langle I_k|F_{I_k}|I_{k+1} \rangle $
where $F=T$.

To express the fundamental equation for $G$ under the form of a probabilistic average, we write the importance-sampled version of the equation
$$
{\bar G}^{(N)}_{I_0 I_N}={\bar G}^{(N),{\cal D}}_{I_0 I_N} +
$$
\be
\sum_{p=1}^{N}
\sum_{|I_1\rangle \notin {\cal D}_{I_0}, \hdots , |I_p\rangle \notin {\cal D}_{I_{p-1}}} 
\sum_{n_0 \ge 1} ... \sum_{n_p \ge 1}
\ee
\be
\delta(\sum_k n_k=N+1) \Big[ \prod_{k=0}^{p-1} [\frac{\Psi^+_{I_{k+1}}}{\Psi^+_{I_k}} \langle I_k|  T^{n_k-1}_{I_k} F_{I_k} |I_{k+1} \rangle \Big]
{\bar G}^{(n_p-1),{\cal D}}_{I_p I_N}.
\label{Gbart}
\ee
Introducing the weight 
\be
W_{I_k I_{k+1}}=\frac{\langle I_k|T^{n_k-1}_{I_k} F_{I_k} |I_{k+1}\rangle}{\langle I_k|T^{+\;n_k-1}_{I_k} F^+_{I_k} |I_{k+1} \rangle}
\ee
and using the effective transition probability, Eq.(\ref{eq3C}), we get
\be
{\bar G}^{(N)}_{I_0 I_N}={\bar G}^{(N),{\cal D}}_{I_0 I_N}+ \sum_{p=1}^{N}
\bigg \langle
\Big( \prod_{k=0}^{p-1} W_{I_k I_{k+1}} \Big)
{\bar G}^{(n_p-1), {\cal D}}_{I_p I_N}
\bigg \rangle
\label{Gbart}
\ee
where the average is defined as
$$
\bigg \langle F \bigg \rangle
= \sum_{|I_1\rangle \notin {\cal D}_{I_0}, \hdots , |I_p\rangle \notin {\cal D}_{I_{p-1}}} 
\sum_{n_0 \ge 1} ... \sum_{n_p \ge 1}
\delta(\sum_k n_k=N+1) 
$$
\be
\prod_{k=0}^{N-1}{\cal P}(I_k \rightarrow I_{k+1},n_k-1)  F(I_0,n_0;...;I_N,n_N)
\ee
In practice, a schematic DMC algorithm to compute the average is as follows.\\
i) Choose some initial vector $|I_0\rangle$\\
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)$, Eq.(\ref{PiN})\\
$\;\;\;\bullet$ Draw the exit state, $|I_{k+1}\rangle$, using the conditional probability $$\frac{{\cal P}(I_k \rightarrow 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_{\rm max}$ and compute $F(I_0,n_0;...;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_{max}$ is obtained. At large $N_{max}$ where the 
convergence of the average as a function of $p$ is reached, such values can be averaged.
\subsubsection{Integrating out the trapping times : The Domain Green's Function Monte Carlo approach}
\label{energy}
Now, let us 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 within the time-evolution of all 
stochastic paths trapped within 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$. The second, more direct and elegant, is based on the Dyson equation.\\
\\
{\it $\bullet$ The pedestrian way}. Let us define the quantity\\
$$
G^E_{ij}= \tau \sum_{N=0}^{\infty} \langle i|T^N|j\rangle.
$$
By summing over $N$ we obtain
\be
G^E_{ij}= \langle i | \frac{1}{H-E}|j\rangle.
\ee
This quantity, which no longer depends on the time-step, is referred to as the energy-dependent Green's matrix. Note that in the continuum this quantity is 
essentially the Laplace transform of the time-dependent Green's function. Here, we then use the same denomination. The remarkable property 
is that, thanks to the summation over $N$ up to the 
infinity the constrained multiple sums appearing in Eq.(\ref{Gt}) can be factorized in terms of a product of unconstrained single sums as follows
$$
\sum_{N=1}^\infty \sum_{p=1}^N \sum_{n_0 \ge 1} ...\sum_{n_p \ge 1} \delta(n_0+...+n_p=N+1)
$$
$$
= \sum_{p=1}^{\infty} \sum_{n_0=1}^{\infty} ...\sum_{n_p=1}^{\infty}.
$$
It is then a trivial matter to integrate out exactly the $n_k$ variables, leading to
$$
\langle I_0|\frac{1}{H-E}|I_N\rangle =  \langle I_0|P_0\frac{1}{H-E} P_0|I_N\rangle
+ \sum_{p=1}^{\infty}
\sum_{I_1 \notin {\cal D}_0, \hdots , I_p \notin {\cal D}_{p-1}}
$$
\be
\Big[ \prod_{k=0}^{p-1} \langle I_k| P_k \frac{1}{H-E} P_k (-H)(1-P_k)|I_{k+1} \rangle \Big]
\langle I_p| P_p \frac{1} {H-E} P_p|I_N\rangle
\label{eqfond}
\ee
\noindent
As an illustration, Appendix \ref{A} reports the exact derivation of this formula in the case of a two-state system.\\
\\
{\it $\bullet$ 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
\frac{1}{H-E} = \frac{1}{H_0-E}
+ \frac{1}{H_0-E} (H_0-H)\frac{1}{H-E}
\ee
where $H_0$ is some arbitrary reference Hamiltonian, we have
the Dyson equation
$$
G^E_{ij}=  G^E_{0,ij}  + \sum_{k,l} G^{E}_{0,ik} (H_0-H)_{kl} G^E_{lj}
$$
Let us choose as $H_0$
$$
\langle i |H_0|j\rangle= \langle i|P_i H P_i|j\rangle  \;\;\; {\rm for \; all \; i \;j}.
$$
The Dyson equation becomes
$$
\langle i| \frac{1}{H-E}|j\rangle 
= \langle i| P_i \frac{1}{H-E} P_i|j\rangle
$$
\be
+ \sum_k  \langle i| P_i \frac{1}{H-E} P_i(H_0-H)|k\rangle \langle k|\frac{1}{H-E}|j\rangle
\ee
Now, we have
$$
P_i \frac{1}{H-E} P_i(H_0-H) = P_i \frac{1}{H-E} P_i (P_i H P_i - H) 
$$
$$
= P_i \frac{1}{H-E} P_i (-H) (1-P_i)
$$
and the Dyson equation may be written under the form
$$
\langle i| \frac{1}{H-E}|j\rangle 
= \langle i| P_i \frac{1}{H-E} P_i|j\rangle
$$
$$
+ \sum_{k \notin {\cal D}_i}  \langle i| P_i \frac{1}{H-E} P_i (-H)(1-P_i)|k\rangle \langle k|\frac{1}{H-E}|j\rangle
$$
which is identical to Eq.(\ref{eqfond}) when $G$ is expanded iteratively.\\
\\
Let us use as effective transition probability density
\be
P(I \rightarrow J) = \frac{1} {\Psi^+(I)} \langle I| P_I \frac{1}{H^+-E^+_L} P_I (-H^+) (1-P_I)|J\rangle \Psi^+(J)
\ee
and the weight
\be
W^E_{IJ} =
\frac{\langle I|\frac{1}{H-E} P_I (-H)(1-P_I) |J\rangle }{\langle I|\frac{1}{H^+-E^+_L} P_I (-H^+)(1-P_I) |J\rangle}
\ee
Using Eqs.(\ref{eq1C},\ref{eq3C},\ref{relation}) we verify that $P(I \rightarrow J) \ge 0$ and $\sum_J P(I \rightarrow J)=1$. 
Finally, the probabilistic expression writes
$$
\langle I_0| \frac{1}{H-E}|I_N\rangle 
= \langle I_0| P_{I_0} \frac{1}{H-E} P_{I_0}|I_N\rangle
$$
\be
+ \sum_{p=0}^{\infty}  \Bigg \langle \Big( \prod_{k=0}^{p-1} W^E_{I_k I_{k+1}} \Big) \langle I_p| P_{I_p} \frac{1}{H-E} P_{I_p}|I_N\rangle \Bigg \rangle
\label{final_E}
\ee
{\it Energy estimator.} To calculate the energy we introduce the following quantity
\be
{\cal E}(E) = \frac{ \langle I_0|\frac{1}{H-E}|H\Psi_T\rangle} {\langle I_0|\frac{1}{H-E}|\Psi_T\rangle}.
\label{calE}
\ee
and search for the solution $E=E_0$ of
\be
{\cal E}(E)= E
\ee
Using the spectral decomposition of $H$ we have
\be
{\cal E}(E) = \frac{ \sum_i \frac{E_i c_i}{E_i-E}}{\sum_i \frac{c_i}{E_i-E}}
\label{calE}
\ee
with
\be
c_i = \langle  I_0| \Phi_0\rangle \langle  \Phi_0| \Psi_T\rangle
\ee
It is easy to check that in the vicinity of $E=E_0$, ${\cal E}(E)$ is a linear function of $E-E_0$.
So, in practice we compute a few values of ${\cal E}(E^{k})$ and fit the data using some function of $E$ close to the linearity
to extrapolate the exact value of $E_0$. Let us describe the 3 functions used here for the fit.\\

i) Linear fit\\
We write
\be
{\cal E}(E)= a_0 + a_1 E
\ee
and search the best value of $(a_0,a_1)$ by fitting the data.
Then, ${\cal E}(E)=E$ leads to
$$
E_0= \frac{a_0}{1-a_1}
$$
ii) Quadratic fit\\
At the quadratic level we write
$$
{\cal E}(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 ${\cal E}(E)$ is known as an infinite series, Eq.(\ref{calE}).
If we limit ourselves to the first two-component, we write
\be
{\cal E}(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)$.\\

In order to have a precise extrapolation of the energy, it is interesting to compute the ratio 
${\cal E}(E)$ for values of $E$ as close as possible to the exact energy. However, in that case 
the numerators and denominators computed diverge. This is reflected by the fact that we need to compute 
more and more $p$-components with an important increase of statistical fluctuations. So, in practice 
a tradoff has to be found between the possible bias in the extrapolation and the amount of simulation time 
required. 

\section{Numerical application to the Hubbard model}
\label{sec:numerical}
Let us consider the one-dimensional Hubbard Hamiltonian for a chain of $N$ sites 
\be
\hat{H}= -t \sum_{\langle i j\rangle \sigma} \hat{a}^+_{i\sigma} \hat{a}_{j\sigma} 
+ U \sum_i \hat{n}_{i\uparrow} \hat{n}_{i\downarrow}
\ee
where $\langle i j\rangle$ denotes the summation over two neighboring sites,
$\hat{a}_{i\sigma} (\hat{a}_{i\sigma})$ is the fermionic creation (annihilation) operator of
an electron of spin $\sigma$ ($=\uparrow$ or $\downarrow$) at 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 (OBC)
at half-filling, that is, $N_{\uparrow}=N_{\downarrow}=\frac{N}{2}$.
In the site-representation, a general vector of the Hilbert space will be written as
\be
|n\rangle = |n_{1 \uparrow},...,n_{N \uparrow},n_{1 \downarrow},...,n_{N \downarrow}\rangle
\ee
where $n_{i \sigma}=0,1$ is the number of electrons of spin $\sigma$ at site $i$.

For the 1D Hubbard model and OBC the components of the ground-state vector have the same sign (say, $c_i \ge 0$). 
It is then possible to identify the guiding and trial vectors, that is, $|c^+\rangle=|c_T\rangle$.
As trial wave function we shall employ a generalization of the Gutzwiller wave function\cite{Gutzwiller_1963} written under the form
\be
\langle n|c_T\rangle= e^{-\alpha n_D(n)-\beta n_A(n)}
\ee
where $n_D(n)$ is the number of doubly occupied sites for the configuration $|n\rangle$ and $n_A(n)$ 
the number of nearest-neighbor antiparallel pairs defined as
\be
n_A(n)= \sum_{\langle i j\rangle} n_{i\uparrow} n _{j\downarrow} \delta(n_{i\uparrow}-1) \delta(n_{j\downarrow}-1)
\ee
where the kronecker function $\delta(k)$ is equal to 1 when $k=0$ and 0 otherwise.
The parameters $\alpha, \beta$ are optimized by minimizing the variational energy,
$E_v(\alpha,\beta)=\frac{\langle c_T|H|c_T\rangle} {\langle c_T|c_T\rangle}$.\\

{\it Domains}. 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. 
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 (electrons move 
more easily). Therefore, we build domains associated with small $n_D$ and large $n_A$ in a hierarchical way as described below. 
For simplicity and decreasing the number of diagonalizations to perform,
we shall consider only one non-trivial domain called here the main domain and denoted as ${\cal D}$. 
This domain will be chosen common to all states belonging to it, that is 
\be
{\cal D}_i= {\cal D} \;\; {\rm for } \; |i\rangle \in {\cal D}
\ee
For the other states we choose a single-state domain
\be
{\cal D}_i= \{ |i\rangle \} \;\; {\rm for} \; |i\rangle \notin {\cal D}
\ee
To define ${\cal D}$, let us introduce the following set of states
\be
{\cal S}_{ij} = \{ |n\rangle; n_D(n)=i {\rm \; and\;} n_A(n)=j \}.
\ee
${\cal D}$ is defined as the set of states having up to $n_D^{\rm 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^{\rm min}(m)$ and $n_A^{\rm max}(m)$. 
Here, $n_A^{\rm max}(m)$ will not be varied and taken to be the maximum possible for a given $m$, 
$n_A^{\rm max}(m)= {\rm Max}(N-1-2m,0)$. 
Using these definitions, the main domain is taken as the union of some elementary domains
\be
{\cal D} = \cup_{n_D=0}^{n_D^{\rm max}}{\cal D}(n_D,n_A^{\rm min}(n_D))
\ee
where the elementary domain is defined as
\be
{\cal D}(n_D,n_A^{\rm min}(n_D))=\cup_{ n_A^{\rm min}(n_D) \leq j \leq n_A^{\rm max}(n_D)} {\cal S}_{n_D j}
\ee
The two quantities defining the main domain are thus $n_D^{\rm max}$ and
$n_A^{\rm min}(m)$.
To give an illustrative example, let us consider the 4-site case. There are 6 possible elementary domains 
$$
{\cal D}(0,3)= {\cal S}_{03} 
$$
$$
{\cal D}(0,2)= {\cal S}_{03} \cup {\cal S}_{02}
$$
$$
{\cal D}(0,1)= {\cal S}_{03} \cup {\cal S}_{02} \cup {\cal S}_{01}
$$
$$
{\cal D}(1,1)= {\cal S}_{11}
$$
$$
{\cal D}(1,0)= {\cal S}_{11} \cup {\cal S}_{10}
$$
$$
{\cal D}(2,0)= {\cal S}_{20} 
$$
where
$$
{\cal S}_{03} = \{\; |\uparrow,\downarrow,\uparrow,\downarrow \rangle, |\downarrow,\uparrow,\downarrow,\uparrow \rangle \} 
\;\;({\rm the\; two\; N\acute eel\; states})
$$
$$
{\cal S}_{02} = \{\; |\uparrow, \downarrow, \downarrow, \uparrow \rangle, |\downarrow, \uparrow, \uparrow, \downarrow \rangle \}
$$
$$
{\cal S}_{01} = \{\; |\uparrow, \uparrow, \downarrow, \downarrow \rangle, |\downarrow, \downarrow, \uparrow, \uparrow \rangle \}
$$
$$
{\cal S}_{11} = \{\; |\uparrow \downarrow, \uparrow ,\downarrow, 0 \rangle, |\uparrow \downarrow, 0, \uparrow,\uparrow \rangle + ... 
\}
$$
$$
{\cal S}_{10} = \{\; |\uparrow \downarrow, \uparrow, 0, \downarrow \rangle, |\uparrow \downarrow, 0, \uparrow, \downarrow \rangle + ... \}
$$
$$
{\cal S}_{20} = \{\; |\uparrow \downarrow, \uparrow \downarrow, 0 ,0 \rangle + ... \}
$$
For the three last cases, the dots indicate the remaining states obtained by permuting the position of the pairs.\\
\\
Let us now present our QMC 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.
\\

Following Eqs.(\ref{final_E},\ref{calE}), the practical formula used for computing the QMC energy is written as
\be
{\cal E}_{QMC}(E,p_{ex},p_{max})= \frac{H_0 +...+ H_{p_{ex}} + \sum_{p=p_{ex}+1}^{p_{max}} H^{QMC}_p } 
                         {S_0 +...+ S_{p_{ex}} + \sum_{p=p_{ex}+1}^{p_{max}} S^{QMC}_p }
\label{calE}
\ee
where $p_{ex}+1$ is the number of pairs, ($H_p$, $S_p$), computed analytically. For $p_{ex} < p \le p_{max}$
the Monte Carlo estimates are written as
\be
H^{QMC}_p= \Bigg \langle \Big( \prod_{k=0}^{p-1} W^E_{I_k I_{k+1}} \Big) \langle I_p| P_{I_p} \frac{1}{H-E} P_{I_p}|H\Psi_T\rangle \Bigg \rangle
\ee
and
\be
S^{QMC}_p= \Bigg \langle \Big(\prod_{k=0}^{p-1} W^E_{I_k I_{k+1}} \Big) \langle I_p| P_{I_p} \frac{1}{H-E} P_{I_p}|\Psi_T\rangle \Bigg \rangle.
\ee

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...$. 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_v=-0.495361...$. In what follows 
$|I_0\rangle$ will be systematically chosen as one of the two N\'eel states, {\it e.g.} $|I_0\rangle =|\uparrow,\downarrow, \uparrow,...\rangle$. 

Figure \ref{fig1} 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.9$, and
$E=-0.8$. Only $H_0$ is computed analytically ($p_{ex}=0$). At the scale of the figure error bars are too small to be perceptible.
When $E$ is far from the exact value of $E_0=-0.768...$ the convergence is very rapid and only a few terms of the $p$-expansion are necessary. 
In constrast, 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, ${\cal E}_{QMC}(E,p_{ex}=1,p_{max})$ as a function of $E$ is presented in figure \ref{fig2}. 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 smaller $E$ the curve is extrapolated using 
the two-component expression. The estimate of the energy obtained from ${\cal E}(E)=E$ is $-0.76807(5)$ in full agreement with the exact value of $-0.768068...$.

\begin{figure}[h!]
\includegraphics[width=\columnwidth]{fig1}
\caption{1D-Hubbard model, $N=4$, $U=12$. $H_p$ as a function of $p$ for $E=-1.6,-1.2,-1.,-0.9,-0.8$. $H_0$ is
computed analytically and $H_p$ (p > 0) by Monte Carlo. Error bars are smaller than the symbol size.}
\label{fig1}
\end{figure}


\begin{figure}[h!]
\includegraphics[width=\columnwidth]{fig2}
\caption{1D-Hubbard model, $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$ is the exact energy of -0.768068.... The dotted line is the two-component extrapolation.
Error bars are smaller than the symbol size.}
\label{fig2}
\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 column 2, the average trapping time 
for the state $|I_0\rangle$ 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 ${\cal D}(0,3) \cup {\cal D}(1,0)$ has a trapping time for the N\'eel state of 6.2, while
the domain ${\cal D}(0,3) \cup {\cal 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, ${\cal E}(E=-1)=-0.75272390...$, can be found at the last row of the Table 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.

\begin{table}[h!]
\centering
\caption{$N$=4, $U$=12, $E$=-1, $\alpha=1.292$, $\beta=0.552$,$p_{ex}=4$. Simulation 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_{\rm conv}$ is a measure of the convergence of ${\cal E}_{QMC}(p)$ as a function of $p$, see text.}
\label{tab1}
\begin{ruledtabular}
\begin{tabular}{lcccl}
Domain & Size & $\bar{t}_{I_0}$ & $p_{\rm conv}$ & $\;\;\;\;\;\;{\cal E}_{QMC}$ \\
\hline
Single & 1 & 0.026 & 88 &$\;\;\;\;$-0.75276(3)\\
${\cal D}(0,3)$ & 2 & 2.1  & 110 &$\;\;\;\;$-0.75276(3)\\
${\cal D}(0,2)$ & 4 & 2.1 & 106  &$\;\;\;\;$-0.75275(2)\\
${\cal D}(0,1)$ & 6 & 2.1&  82   &$\;\;\;\;$-0.75274(3)\\
${\cal D}(0,3)$ $\cup$ ${\cal D}(1,1)$ &14 &4.0& 60 & $\;\;\;\;$-0.75270(2)\\
${\cal D}(0,3)$ $\cup$ ${\cal D}(1,0)$ &26 &6.2& 45 & $\;\;\;\;$-0.752730(7) \\
${\cal D}(0,2)$ $\cup$ ${\cal D}(1,1)$ &16 &10.1 & 36 &$\;\;\;\;$-0.75269(1)\\
${\cal D}(0,2)$ $\cup$ ${\cal D}(1,0)$ &28 &34.7  & 14&$\;\;\;\;$-0.7527240(6)\\
${\cal D}(0,1)$ $\cup$ ${\cal D}(1,1)$ &18 & 10.1 & 28 &$\;\;\;\;$-0.75269(1)\\
${\cal D}(0,1)$ $\cup$ ${\cal D}(1,0)$ &30 & 108.7 & 11&$\;\;\;\;$-0.75272400(5) \\
${\cal D}(0,3)$ $\cup$ ${\cal D}(1,1)$ $\cup$ ${\cal D}$(2,0) &20 & 4.1 & 47 &$\;\;\;\;$-0.75271(2)\\
${\cal D}(0,3)$ $\cup$ ${\cal D}(1,0)$ $\cup$ ${\cal D}$(2,0) &32 & 6.5 & 39 &$\;\;\;\;$-0.752725(3)\\
${\cal D}(0,2)$ $\cup$ ${\cal D}(1,1)$ $\cup$ ${\cal D}$(2,0) &22 & 10.8 & 30 &$\;\;\;\;$-0.75270(1)\\
${\cal D}(0,2)$ $\cup$ ${\cal D}(1,0)$ $\cup$ ${\cal D}$(2,0) &34 & 52.5 & 13&$\;\;\;\;$-0.7527236(2)\\
${\cal D}(0,1)$ $\cup$ ${\cal D}(1,1)$ $\cup$ ${\cal D}$(2,0) & 24 & 10.8 & 26&$\;\;\;\;$-0.75270(1)\\
${\cal D}(0,1)$ $\cup$ ${\cal D}(1,0)$ $\cup$ ${\cal D}$(2,0) & 36 & $\infty$&1&$\;\;\;\;$-0.75272390\\
\end{tabular}
\end{ruledtabular}
\end{table}

As a general rule, it is always good to avoid the Monte Carlo calculation of a quantity which is computable analytically. Here, we apply 
this idea to the case of the energy, Eq.(\ref{calE}), where the first $p$-components can be evaluated exactly up to some maximal value of $p$, $p_{ex}$. 
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 ${\cal D}(0,1) \cup {\cal D}(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 on Fig.\ref{fig2}. As seen on the figure the behavior of ${\cal E}$ 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... 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 ${\cal D}(0,1) \cup {\cal 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 \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 QMC 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 QMC 
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. 
laptop. 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.

\begin{table}[h!]
\caption{$N=4$, $U=12$, and $E=-1$. Dependence of the statistical error on the energy with the number of $p$-components calculated 
analytically. Same simulation as for Table \ref{tab1}. Results are presented when a single-state domain 
is used for all states and when 
${\cal D}(0,1) \cup {\cal D}(1,0)$ is chosen as main domain.}
\label{tab2}
\begin{ruledtabular}
\begin{tabular}{lcc}
$p_{ex}$ & single-state & ${\cal D}(0,1) \cup {\cal D}(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}


\begin{table}[h!]
\caption{$N=4$, $U=12$, $\alpha=1.292$, $\beta=0.552$. Main domain = ${\cal D}(0,1) \cup {\cal D}(1,0)$. Simulation with 20 independent blocks and $10^6$ paths.
$p_{ex}=4$. The various fits are done with the five values of $E$}
\label{tab3}
\begin{ruledtabular}
\begin{tabular}{lc}
$E$ & $E_{QMC}$ \\
\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...\\
\end{tabular}
\end{ruledtabular}
\end{table}

\begin{table}[h!]
\caption{$N$=4, Domain ${\cal D}(0,1) \cup {\cal D}(1,0)$}
\label{tab4}
\begin{ruledtabular}
\begin{tabular}{cccccc}
$U$ & $\alpha,\beta$ & $E_{var}$ & $E_{ex}$  & $\bar{t}_{I_0}$ \\
\hline
8  & 0.908,\;0.520 & -0.770342... &-1.117172... & 33.5\\
10 & 1.116,\;0.539 & -0.604162... &-0.911497... & 63.3\\
12 & 1.292,\;0.552 & -0.495361... &-0.768068... & 108.7\\
14 & 1.438,\;0.563 & -0.419163... &-0.662871... & 171.7  \\
20 & 1.786,\;0.582 & -0.286044... &-0.468619... & 504.5  \\
50 & 2.690,\;0.609 & -0.110013... &-0.188984... & 8040.2  \\
200 & 4.070,\;0.624& -0.026940... &-0.047315... & 523836.0  \\
\end{tabular}
\end{ruledtabular}
\end{table}

\begin{table*}[h!]
\caption{$U=12$. The fits to extrapolate the QMC energies are done using the two-component function}
\label{tab5}
\begin{ruledtabular}
\begin{tabular}{crcrccccc}
$N$ & Size Hilbert space & Domain & Domain size & $\alpha,\beta$ &$\bar{t}_{I_0}$ & $E_{var}$ & $E_{ex}$ & $E_{QMC}$\\
\hline
4  & 36     & ${\cal D}(0,1) \cup {\cal D}(1,0)$ & 30 &1.292, \; 0.552& 108.7  & -0.495361 & -0.768068 & -0.7680676(5)\\\
6  & 400    & ${\cal D}(0,1) \cup {\cal D}(1,0)$ &200 & 1.124,\;  0.689 &57.8 & -0.633297 & -1.215395&  -1.215389(9)\\
8  & 4 900   & ${\cal D}(0,1) \cup {\cal D}(1,0)$ & 1 190 & 0.984,\;  0.788 & 42.8 & -0.750995 & -1.66395&  -1.6637(2)\\
10 & 63 504  & ${\cal D}(0,5) \cup {\cal D}(1,4)$ & 2 682 & 0.856,\;  0.869& 31.0  & -0.855958  & -2.113089& -2.1120(7)\\
12 & 853 776 & ${\cal D}(0,8) \cup {\cal D}(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 and the corresponding general equations have been derived.
In this way a new effective stochastic dynamics connecting only the domains and not the individual states is defined.
A key property of the effective dynamics is that it does not depend on the shape of the domains used for each state, so rather general 
domains overlapping or not can be used. To obtain the effective transition probability giving the probability of going from 
one domain to another and the renormalized estimators,
the Green's functions limited to the domains are to be computed analytically, that is, in practice, matrices of size the number 
of states of the sampled domains need to be inverted. This is the main computationally intensive step of the method. 
The efficiency of the method is directly related to 
the importance of the average time spent by the stochastic trajectories within 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 application.
So a tradeoff has to be found in the complexity of the domains maximizing the average trapping time and the cost of computing 
the local Green's functions associated with such domains. In practice, there is no general rule to construct efficient domains.  
For each system at hand, we need to determine on physical grounds which regions of the 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 1D-Hubbard model 
we exploit the physics of the large $U$ regime which is known to approach the Heisenberg limit where double occupations have a small weight. 
This simple example has been chosen to illustrate the various aspects of the approach. However, our goal is of course to
tackle much larger systems, like those treated by state-of-the-art methods, such as selected CI, (for example, \cite{Huron_1973,Harrison_1991,Giner_2013,Holmes_2016,
Schriber_2016,Tubman_2020}, QMCFCI,\cite{Booth_2009,Cleland_2010}, AFQMC\cite{Zhang_2003}, or DMRG\cite{White_1999,Chan_2011} approaches).
Here, we have mainly focused on the theoretical aspects of the approach. To reach larger systems will require some 
more elaborate implementation of the method in order to keep under control the cost of the simulation.
Doing this was out of the scope of the present work and will be presented in a forthcoming work.

\section*{Acknowledgement}
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 the $2\times2$ matrix}
\label{A}
For the simplest case of a two-state system the fundamental equation (\ref{eqfond}) writes
$$
{\cal I}=
\langle I_0|\frac{1}{H-E}|\Psi\rangle =  \langle I_0|P_0\frac{1}{H-E} P_0|\Psi\rangle
+ \sum_{p=1}^{\infty} {\cal I}_p$$
with
$$
 {\cal I}_p=
\sum_{I_1 \notin {\cal D}_0, \hdots , I_p \notin {\cal D}_{p-1}}
$$
\be
\Big[ \prod_{k=0}^{p-1} \langle I_k| P_k \frac{1}{H-E} P_k (-H)(1-P_k)|I_{k+1} \rangle \Big]
\langle I_p| P_p \frac{1} {H-E} P_p|\Psi\rangle
\label{eqfond}
\ee
To treat simultaneously the two possible cases for the final state, $|I_N\rangle =|1\rangle$ or $|2\rangle$,
the equation has been slightly generalized to the case of a general vector for the final state written as
\be
|\Psi\rangle = \Psi_1 |1\rangle +  \Psi_2 |2\rangle
\ee
where
$|1\rangle$ and $|2\rangle$ denote the two states. Let us choose a single-state domain for both states, namely ${\cal D}_1=\{|1\rangle \}$  and ${\cal D}_2=\{ |2\rangle \}$. Note that the single exit state for each state is the other state. Thus there are only two  possible deterministic "alternating" paths, namely either 
$|1\rangle \rightarrow |2\rangle \rightarrow  |1\rangle,...$ or $|2\rangle \rightarrow |1\rangle \rightarrow  |2\rangle,...$
We introduce the following quantities
\be
A_1= \langle 1| P_1 \frac{1}{H-E} P_1 (-H)(1-P_1)|2\rangle \;\;\;
A_2= \langle 2| P_2 \frac{1}{H-E} P_2 (-H) (1-P_2)|1\rangle 
\label{defA}
\ee
and
\be
C_1= \langle 1| P_1 \frac{1}{H-E} P_1  |\Psi\rangle  \;\;\;
C_2= \langle 2| P_2 \frac{1}{H-E} P_2  |\Psi\rangle 
\ee
Let us choose, for example, $|I_0\rangle =|1\rangle$, we then have
\be
{\cal I}_1= A_1 C_2  \;\; {\cal I}_2 = A_1 A_2 C_1 \;\;  {\cal I}_3 = A_1 A_2 A_1 C_2 \;\; etc.
\ee
that is
\be
{\cal I}_{2k+1} = \frac{C_2}{A_2} (A_1 A_2)^{2k+1}
\ee
and
\be
{\cal I}_{2k} = C_1 (A_1 A_2)^{2k}
\ee
Then
\be
\sum_{p=1}^{\infty}
{\cal I}_p
= \frac{C_2}{A_2} \big[ \sum_{p=1}^{\infty} (A_1 A_2)^p \big] + {C_1}  \big[ \sum_{p=1}^{\infty} (A_1 A_2)^p \big]
\ee
which gives
\be
\sum_{p=1}^{\infty}
{\cal I}_p
= A_1 \frac{C_2 + C_1 A_2}{1-A_1 A_2}
\ee
and thus
\be
\langle 1 | \frac{1} {H-E} |\Psi\rangle 
=
C_1 +  A_1 \frac{C_2 + C_1 A_2}{1-A_1 A_2}
\ee
For a two by two matrix it is easy to evaluate the $A_i$'s, Eq.\ref{defA}, we
have
$$
A_i = -\frac{H_{12}}{H_{ii}-E} \;\;\; \;\;\; C_i= \frac{1}{H_{ii}-E} \Psi_i \;\;\; i=1,2
$$
Then
$$
\langle 1 | \frac{1} {H-E} |\Psi\rangle 
= \frac{1}{H_{11}-E} \Psi_1  -H_{12} \frac{ \big( \Psi_2 - \frac{ H_{12}\Psi_1}{H_{11}-E} \big)} 
                                          { (H_{11}-E)(H_{22}-E) - H^2_{12}}
$$
$$
= \frac{ (H_{22}-E) \Psi_1}{\Delta} -\frac{H_{12} \Psi_2}{\Delta}
$$
where $\Delta$ is the determinant of the matrix $H$. It is easy to check that the RHS is equal to the LHS, thus validating the fundamental equation 
for this particular case.

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