g.bib

@@ -1,7 +1,7 @@
 %% This BibTeX bibliography file was created using BibDesk.
 %% http://bibdesk.sourceforge.net/
 
-%% Created for Pierre-Francois Loos at 2022-12-16 13:28:28 +0100 
+%% Created for Pierre-Francois Loos at 2022-10-05 14:42:12 +0200 
 
 
 %% Saved with string encoding Unicode (UTF-8) 
@@ -11,7 +11,7 @@
 @article{Ceperley_1983,
 	abstract = {Random walks with branching have been used to calculate exact properties of the ground state of quantum many-body systems. In this paper, a more general Green's function identity is derived which relates the potential energy, a trial wavefunction, and a trial density matrix to the rules of a branched random walk. It is shown that an efficient algorithm requires a good trial wavefunction, a good trial density matrix, and a good sampling of this density matrix. An accurate density matrix is constructed for Coulomb systems using the path integral formula. The random walks from this new algorithm diffuse through phase space an order of magnitude faster than the previous Green's Function Monte Carlo method. In contrast to the simple diffusion Monte Carlo algorithm, it is an exact method. Representative results are presented for several molecules.},
 	author = {D Ceperley},
-	date-modified = {2022-12-16 13:14:50 +0100},
+	date-modified = {2022-10-05 14:41:32 +0200},
 	doi = {https://doi.org/10.1016/0021-9991(83)90161-4},
 	issn = {0021-9991},
 	journal = {J. Comput. Phys.},
@@ -58,7 +58,7 @@
 
 @article{Moskowitz_1986,
 	author = {Moskowitz,Jules W. and Schmidt,K. E.},
-	date-modified = {2022-12-16 13:14:34 +0100},
+	date-modified = {2022-10-05 14:41:57 +0200},
 	doi = {10.1063/1.451046},
 	journal = {J. Chem. Phys.},
 	number = {5},
@@ -68,6 +68,9 @@
 	year = {1986},
 	bdsk-url-1 = {https://doi.org/10.1063/1.451046}}
 
+@misc{note,
+	note = {As $\tau \rightarrow 0$ and $N \rightarrow \infty$ with $N\tau=t$, the operator $T^N$ converges to $e^{-t(H-E \Id)}$. We then have $G^E_{ij} \rightarrow \int_0^{\infty} dt \mel{i}{e^{-t(H-E \Id)}}{j}$, which is the Laplace transform of the time-dependent Green's function $\mel{i}{e^{-t(H-E \Id)}}{j}$.}}
+
 @article{Willow_2012,
 	author = {Willow,Soohaeng Yoo and Kim,Kwang S. and Hirata,So},
 	date-modified = {2022-09-30 16:10:46 +0200},
@@ -1695,10 +1698,10 @@
 
 @article{Assaraf_2000,
 	author = {Assaraf, Roland and Caffarel, Michel and Khelif, Anatole},
-	date-modified = {2022-12-16 13:13:55 +0100},
+	date-modified = {2022-09-15 10:17:18 +0200},
 	doi = {10.1103/physreve.61.4566},
 	issn = {1095-3787},
-	journal = {Phys. Rev. E},
+	journal = {Phsys. Rev. E},
 	month = {Apr},
 	number = {4},
 	pages = {4566--4575},

g.tex

@@ -1,9 +1,6 @@
-%\documentclass[aps,prb,reprint]{revtex4-1}
+\documentclass[aps,prb,reprint,showkeys,superscriptaddress]{revtex4-1}
-%\documentclass[aps,prb,reprint,showkeys,superscriptaddress]{revtex4-1}
-\documentclass[reprint,superscriptaddress,nobibnotes,amsmath,amssymb,twocolumn,prb]{revtex4-1}
 %\usepackage{subcaption}
-%\usepackage{bm,graphicx,tabularx,array,booktabs,dcolumn,xcolor,microtype,multirow,amscd,amsmath,amssymb,amsfonts,physics,siunitx}
+\usepackage{bm,graphicx,tabularx,array,booktabs,dcolumn,xcolor,microtype,multirow,amscd,amsmath,amssymb,amsfonts,physics,siunitx}
-\usepackage{bm,graphicx,tabularx,array,dcolumn,xcolor,multirow,amsmath,amssymb,amsfonts,physics,siunitx}
 \usepackage[version=4]{mhchem}
 \usepackage[utf8]{inputenc}
 \usepackage[T1]{fontenc}
@@ -11,6 +8,12 @@
 \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}}
 
@@ -26,6 +29,7 @@ breaklinks=true
 
 \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.}}
@@ -117,17 +121,22 @@ breaklinks=true
 \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}
-       \affiliation{\LCPQ}
+       \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[Corresponding author: ]{caffarel@irsamc.ups-tlse.fr}
+        \email{caffarel@irsamc.ups-tlse.fr}
         \affiliation{\LCPQ}
 \begin{abstract}
 
@@ -139,7 +148,7 @@ Here, we extend this idea to the general case of a walker trapped within domains
 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 (discrete) finite linear spaces, it can be generalized without fundamental difficulties to continuous configuration spaces.
+Although this work presents the method for finite linear spaces, it can be generalized without fundamental difficulties to continuous configuration spaces.
 \end{abstract}
 
 \maketitle
@@ -154,7 +163,7 @@ DMC can be used either for systems defined in a continuous configuration space (
 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 well-established numerical approaches able to extract the largest or smallest eigenvalues of a matrix (see, \eg, Ref.~\onlinecite{Golub_2012}). 
+In essence, DMC is based on \textit{stochastic} power methods, a family of well established 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} 
@@ -176,21 +185,22 @@ It is shown how to define an effective stochastic dynamics describing walkers mo
 The equations of the effective dynamics are derived and a numerical application for a model (one-dimensional) problem is presented. 
 In particular, it shows that the statistical convergence of the energy can be greatly enhanced when domains associated with large average trapping times are considered.
 
-It should be noted that the use of domains in quantum Monte Carlo is not new. Domains have been introduced within the context of 
+It should be noted that the use of domains in quantum Monte Carlo is not new. Domains have been introduced within the context of the 
-Green's function Monte Carlo (GFMC) pioneered by Kalos \cite{Kalos_1962,Kalos_1970} 
+Green's function Monte Carlo (GFMC) method pioneered by Kalos\cite{Kalos_1962,Kalos_1970} 
-and later developed and applied by Kalos and others. \cite{Kalos_1974,Ceperley_1979,Ceperley_1983,Moskowitz_1986} 
+and later developped and applied by Kalos and others.\cite{Kalos_1974,Ceperley_1979,Ceperley_1983,Moskowitz_1986} 
-In GFMC, an approximate Green's function that can be sampled is required for the stochastic propagation of the wave function. 
+In GFMC, an approximate Green's function that can be sampled is needed to propagate stochastically the wavefunction. 
-In the so-called domain GFMC version of GFMC introduced in Ref.~\onlinecite{Kalos_1970} and \onlinecite{Kalos_1974} 
+In the so-called Domain's GFMC version of GFMC introduced in Ref.\onlinecite{Kalos_1970} and \onlinecite{Kalos_1974} 
 the sampling is realized by using the restriction of the Green's function 
-to a small domain consisting of the cartesian product of small spheres around each particle, the potential being considered constant within the domain. Fundamentally, the method presented in this work is closely related to the domain GFMC, although the way we present the formalism in terms of walkers
+to a small domain consisting of the Cartesian product of small spheres around each particle, the potential being considered as constant within the domain.
+Fundamentally, the method presented in this work is closely related to the Domain's GFMC, although the way we present the formalism in terms of walkers
 trapped within domains
-and derive the equations that may appear different. 
+and derive the equations may appear as different. 
-However, we show here how to use 
+However, a key difference here is that we show how to use 
-domains of arbitrary size, a new feature that greatly enhances the efficiency of the simulations when domains are suitably chosen, as illustrated in our numerical 
+domains of arbitrary size, a point which can greatly enhance the efficiency of the simulations when domains are suitably chosen, as illustrated in our numerical 
 application.
 
 
-Finally, from a general perspective, it is interesting to emphasize that the present method illustrates how suitable combinations of stochastic and deterministic techniques lead to a more efficient and valuable method.
+Finally, from a general perspective, it is interesting to mention that the method proposed here is an illustration of how valuable and efficient can be the combination of stochastic and deterministic techniques.
 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 hybrid algorithms for evaluating the second-order perturbation energy in selected configuration interaction methods, \cite{Garniron_2017b,Sharma_2017} the approach of Willow \textit{et al.} for computing stochastically second-order many-body perturbation energies, \cite{Willow_2012} or the zero variance Monte Carlo scheme for evaluating two-electron integrals in quantum chemistry. \cite{Caffarel_2019}
 
@@ -295,7 +305,7 @@ To derive a probabilistic expression for the Green's matrix, we introduce a guid
 \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}$.
 
-Now, the key idea to take advantage of probabilistic techniques is to rewrite the matrix elements of $\bar{T}$ as those of a stochastic matrix multiplied by some residual weights (here, not necessarily positive), namely
+Next, the matrix elements of $\bar{T}$ are expressed as those of a stochastic matrix multiplied by some residual weights, namely,
 \be
 \label{eq:defTij}
 	\bar{T}_{ij}= p_{i \to j} w_{ij}.
@@ -305,42 +315,7 @@ Here, we recall that a stochastic matrix is defined as a matrix with positive en
 \label{eq:sumup}
 	\sum_j p_{i \to j}=1.
 \ee
-Using this representation for $\bar{T}_{ij}$ the similarity-transformed Green's matrix components can be rewritten as
+To build the transition probability density, the following operator is introduced
-\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 state, 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}$ occurring 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.
-
-The rewriting of $\bar{T}_{ij}$ as a product of a stochastic matrix times some general real weight 
-does not introduce any constraint on the choice of the stochastic matrix, so that, in theory, any stochastic matrix could be used.  However, 
-in practice, it is highly desirable that the magnitude of the fluctuations of the weight during the Monte Carlo simulation be as small as
-possible. A natural solution is to choose a stochastic matrix as close as possible to $\bar{T}_{ij}$. This is done as follows.
-
-Let us introduce the  following operator
 \be
 \label{eq:T+}
 	T^+= \Id - \tau \qty( H^+ - \EL^+ \Id ), 
@@ -403,14 +378,54 @@ Note that one can eschew this condition via a simple generalization of the trans
 \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}$, defined from the Hamiltonian $H^{+}$ via $T^+$ and the operator $T$ associated with the exact Hamiltonian $H$. 
+Now, we need to make the connection between the transition probability matrix, $p_{i \to j}$, defined from the approximate Hamiltonian $H^{+}$ via $T^+$ and the operator $T$ associated with the exact Hamiltonian $H$. 
 This is done thanks to Eq.~\eqref{eq:defTij} that connects $p_{i \to j}$ and $T_{ij}$ through the weight
 \be
 	w_{ij}=\frac{T_{ij}}{T^+_{ij}},
 \ee
 derived from Eqs.~\eqref{eq:defT} and \eqref{eq:pij}.
+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}$. 
 
-To calculate the probabilistic averages, an artificial (mathematical) ``particle'' called a walker (or psi-particle) is introduced.
+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 state, 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}$ occurring 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}$. 
 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.
@@ -438,8 +453,8 @@ We shall not insist here on these practical details that are discussed, for exam
 \label{sec:single_domains}
 %=======================================%
 
-During the simulation, walkers move from state to state with the possibility of being trapped a certain number of times in the same state before 
+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 parts of the dynamics, thus leading to a reduction of the statistical 
+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 $ n \ge 1 $) and then exits to a different state $j$ (with $j \neq i$) is
@@ -494,7 +509,7 @@ Let us write an arbitrary path of length $N$ as
 	\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 path belongs to the set of paths that can be represented as follows
+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),
@@ -540,11 +555,7 @@ Since the sums are restricted to states belonging to the domain, it is convenien
 \label{eq:pi}
 	P_I = \sum_{\ket{i} \in \cD_I} \dyad{i}{i},
 \ee
-as well as the projection of $T^+$ over $\cD_I$
+as well as the projection of $T^+$ over $\cD_I$, \ie, $T^+_I = P_I T^+ P_I$, which governs the dynamics of the walkers trapped in this domain.
-\be
-T^+_I = P_I T^+ P_I, 
-\ee
-which governs the dynamics of the walkers trapped in this domain.
 %see Eq.~\eqref{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
@@ -555,14 +566,14 @@ where the operator $F^+_I =  P_I T^+ (1-P_I)$, corresponding to the last move co
  \be
 	(F^+_I)_{ij} = 
 	\begin{cases}
-		T^+_{ij},	&	\qif* \ket{i} \in \cD_{I} \;\;{\rm and}\;\; \ket{j} \notin \cD_{I}, 
+		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, exiting to some state, it is possible to obtain
+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:
 \be
 \label{eq:PiN}
@@ -573,9 +584,7 @@ The normalization of this probability can be verified using the fact that
 \label{eq:relation}
 	\qty(T^+_{I})^{n-1} F^+_I = \qty(T^+_{I})^{n-1} T^+ - \qty(T^+_I)^n,
 \ee
-leading to\footnote{
+leading to
-The property results from the fact that the series is a telescoping series and that the general term
-$\mel{ I }{ \qty(T^+_{I})^{n} }{ \PsiG }$ goes to zero as $n$ goes to infinity.}
 \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 }
@@ -585,8 +594,8 @@ The average trapping time defined as ${\bar t}_{I}={\bar n}_{I} \tau$ where $ {\
 \be 
 	{\bar 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 inverting the matrix $(H^+-\EL^+ \Id)$ in $\cD_{I}$. 
+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 thousand).
+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}
@@ -597,40 +606,23 @@ Note that it is possible only if the dimension of the domains is not too large (
 %\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.
+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 need to introduce the Green's matrix restricted to each domain as follows:
+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
-where $T_i$ is the projection of the operator $T$ over the domain $\cD_i$ of $\ket{i}$
+%Starting from Eq.~\eqref{eq:cn}
-\be
+%\be
-T_i = P_i T^+ P_i,
+%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
+%\ee
-and $P_i$ is the projector, $P_i=\sum_{\ket{k} \in \cD_i} \dyad{k}{k}$, Eq.~\eqref{eq:pi}.
+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
-
-Using the representation of the paths in terms of exit states and trapping times [see Eq.~\eqref{eq:eff_series}], 
-the set of paths can be partitioned according to the number $p$ of exit states.
-As already noted above, the number of exit states ranges from $p=0$ (no exit of the initial domain) to $p=N$ (exit of the current domain at each time step).
-We shall denote $\cC_p$ the set of paths with $p$ exit states.
-The time-dependent Green's matrix, Eq.~\eqref{eq:cn}, is then rewritten as
 \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} } .
+	\prod_{k=0}^{N-1} \mel{ i_k }{ T }{ i_{k+1} } ,
 \ee
-$\cC_p$ can be partitioned further by considering the subset of paths, denoted by $\cC_p^{I,n}$, 
+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$).
-having $(\ket{I_k}, n_k)$, for $0 \le k \le p$, as exit states and trapping times.
+It follows that 
-We recall that, by definition of the exit states, $\ket{I_k} \notin \cD_{I_{k-1}}$. The total time must be conserved, so the relation $\sum_{k=0}^p n_k= N+1$ 
-must be fulfilled.
-Now, the contribution of $\cC_p^{I,n}$ to the path integral is obtained by summing over all elementary paths $(i_1,...,i_{N-1})$ 
-of this set, thus giving
-\begin{multline}
-\sum_{(i_1,...,i_{N-1}) \in \cC_p^{I,n}} \prod_{k=0}^{N-1} \mel{ i_k }{ T }{ i_{k+1} }
-\\
-= \qty[ \prod_{k=0}^{p-1} \mel{ I_k }{ \qty(T_{I_k})^{n_k-1} F_{I_k} }{ I_{k+1} } ]
-         \mel{ I_k }{ \qty(T_{I_p})^{n_p} }{ i_N }
-\end{multline}
-Finally, collecting all contributions the Green's matrix is written as
 \begin{multline}
 \label{eq:Gt}
 	G^{(N)}_{i_0 i_N}= G^{(N),\cD}_{i_0 i_N} +
@@ -726,7 +718,7 @@ First, the integer $n_k$ is drawn using the probability $P_{I_k}(n)$ [see Eq.~\e
 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 making the explicit integration over the $n_k$'s. 
+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.
 
 
@@ -739,8 +731,8 @@ Let us define the energy-dependent Green's matrix
 \ee
 The denomination ``energy-dependent'' is chosen here since 
 this quantity is the discrete version of the Laplace transform of the time-dependent Green's function in a continuous space,
-usually known under this name.\footnote{As $\tau \to 0$ and $N \to \infty$ with $N\tau=t$, the operator $T^N$ converges to $e^{-t(H-E \Id)}$. We then have $G^E_{ij} \to \int_0^{\infty} dt \mel{i}{e^{-t(H-E \Id)}}{j}$, which is the Laplace transform of the time-dependent Green's function $\mel{i}{e^{-t(H-E \Id)}}{j}$.}
+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 sums, as follows
+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 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)
 	\\
@@ -843,10 +835,8 @@ In practical Monte Carlo calculations, the DMC energy is obtained by computing a
 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} },
-\label{eq:defHp}
         \\
         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} },
-\label{eq:defSp}
 \end{align}
 where
 \begin{align}
@@ -863,7 +853,7 @@ For $p=0$, the two components are exactly evaluated as
 \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.
 Avoiding the stochastic calculation of quantities, such as $H_0$ or $S_0$, that can be evaluated analytically is computationally very appealing as the statistical error associated with these quantities is completely removed. 
-We thus suggest extending the exact calculation of $H_p$ and $S_p$ to higher $p$ values, up to the point where the exponential increase of the number of intermediate states involved in the explicit sums makes the calculation unfeasible.
+We thus suggest to extend the exact calculation of $H_p$ and $S_p$ to higher $p$ values, up to the point where the exponential increase of the number of intermediate states involved in the explicit sums makes the calculation unfeasible.
 %\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} } ]
@@ -876,7 +866,6 @@ Finally, $\cE^\text{DMC}(E,p_\text{ex},p_\text{max})$ is evaluated in practice a
 \be
 \cE^\text{DMC}(E,p_\text{ex},p_\text{max})= \frac{   \sum_{p=0}^{p_\text{ex}-1} H_p +\sum_{p=p_\text{ex}}^{p_\text{max}} H^\text{DMC}_p }
                          {                 \sum_{p=0}^{p_\text{ex}-1} S_p +\sum_{p=p_\text{ex}}^{p_\text{max}} S^\text{DMC}_p },
-\label{eq:E_pex}
 \ee
 where $p_\text{ex}$ is the number of components computed exactly. 
 Let us emphasize that, since the magnitude of $H_p$ and $S_p$ decreases as a function of $p$, most of the statistical error is removed by computing the dominant terms analytically.
@@ -886,7 +875,7 @@ It is easy to check that, in the vicinity of the exact energy $E_0$, $\cE(E)$ is
 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 in statistical fluctuations.
+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. 
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -903,7 +892,7 @@ Let us consider the one-dimensional Hubbard Hamiltonian for a chain of $N$ sites
 	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 a 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.
+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
@@ -913,7 +902,7 @@ where $n_{i \sigma}=0$ or $1$ is the number of electrons of spin $\sigma$ on sit
 
 For the one-dimensional Hubbard model with open boundary conditions, the components of the ground-state vector have the same sign (say, $c_i > 0$). 
 It is then possible to equate the guiding and trial vectors, \ie, $\ket{c^+} = \ket{c_\text{T}}$.
-As a trial wave function, we shall employ a generalization of the Gutzwiller wave function \cite{Gutzwiller_1963}
+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
@@ -923,7 +912,7 @@ where $n_D(n)$ is the number of doubly occupied sites for the configuration $\ke
 \ee
 The parameters $\alpha$ and $\beta$ are optimized by minimizing the variational energy of the trial wave function, \ie, 
 \be
-	E_\text{v}(\alpha,\beta) = \frac{\mel{ c_\text{T} }{H }{c_\text{T}}} {\braket{ c_\text{T} }{ c_\text{T} }}.
+	E_\text{T}(\alpha,\beta) = \frac{\mel{ c_\text{T} }{H }{c_\text{T}}} {\braket{ c_\text{T} }{ c_\text{T} }}.
 \ee
 
 %======================================
@@ -939,7 +928,7 @@ The contribution of the other states vanishes as $U$ increases with a rate incre
 In addition, for a given number of double occupations, configurations with large values of $n_A(n)$ are favored due to their high electronic mobility.
 Therefore, we build domains associated with small $n_D$ and large $n_A$ in a hierarchical way as described below. 
 
-For simplicity and reducing the number of matrix inversions to perform, we shall consider only one non-trivial domain called here the main domain and denoted as $\cD$. 
+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.
@@ -1005,21 +994,20 @@ In what follows, we shall restrict ourselves to the case of the Green's function
 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 
+The two variational parameters of the trial vector have been optimized and fixed at the values of $\alpha=1.292$ and $\beta=0.552$ with a variational energy $E_\text{T}=-0.495361\ldots$. 
-variational energy $E_\text{v}=-0.495361\ldots$. 
 In what follows, $\ket{I_0}$ is 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^\text{DMC}_p$, Eq.~\eqref{eq:defHp}, as a function of $p$ for different values of the reference energy $E$.
+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_\text{ex}=0$). 
-At the scale of the figure, the error bars are too small to be seen.
+At the scale of the figure, error bars are too small to be seen.
 
 When $E$ is far from the exact value, the convergence is very rapid and only a few terms of the $p$-expansion are necessary. 
 In contrast, when $E$ approaches the exact energy, a slower convergence is observed, as expected from the divergence of the matrix elements of the Green's matrix at $E=E_0$ where the expansion does not converge at all. 
 The oscillations of the curves as a function of $p$ are due to a parity effect specific to this system. 
 In practice, it is not too much of a problem since a smoothly convergent behavior is nevertheless observed for the even- and odd-parity curves.
-The ratio, $\cE^\text{DMC}(E,p_\text{ex}=1,p_\text{max})$, Eq.~\eqref{eq:E_pex}, as a function of $E$ is presented in Fig.~\ref{fig4}. 
+The ratio, $\cE^\text{DMC}(E,p_\text{ex}=1,p_\text{max})$ as a function of $E$ is presented in Fig.~\ref{fig4}. 
 Here, $p_\text{max}$ is taken sufficiently large so that the convergence at large $p$ is reached. 
 The values of $E$ are $-0.780$, $-0.790$, $-0.785$, $-0.780$, and $-0.775$. 
 For small $E$, the curve is extrapolated using the so-called two-component expression:
@@ -1032,7 +1020,7 @@ which considers only the first two terms of the exact expression of $\cE(E)$ [se
 %\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 need to be determined are $c_0$, $\epsilon_0$, $c_1$, and $\epsilon_1$.
+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 %%%
@@ -1060,10 +1048,9 @@ Error bars are smaller than the symbol size.}
 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 one uses a single-state domain for all states and the last one for the maximal domain containing the full linear space. 
+The first line of the table gives the results when one uses a single-state domain for all states, and the last one for the maximal domain containing the full linear space. 
 The size of the various domains is given in the second column, the average trapping time for the state $\ket{I_0}$ in the third column, and an estimate of the speed of convergence of the $p$-expansion for the energy in the fourth column. 
-To quantify the rate of convergence, we report the quantity, $p_\text{conv}$, defined as the smallest value of $p$ for which the energy is roughly stabilized 
+To quantify the rate of convergence, we report the quantity, $p_\text{conv}$, defined as the smallest value of $p$ for which the energy is converged to six decimal places. 
-to five decimal places. 
 The smaller $p_\text{conv}$, the better the convergence is.
 Although this is a rough estimate, it is sufficient here for our purpose.
 As clearly seen, the speed of convergence is directly related to the magnitude of $\bar{t}_{I_0}$. 
@@ -1085,13 +1072,12 @@ Using this domain leads to a reduction in the statistical error as large as abou
 \caption{One-dimensional Hubbard model with $N=4$, $U=12$, $E=-1$, $\alpha=1.292$, $\beta=0.552$, and $p_\text{ex}=4$. 
 The simulation is performed with 20 independent blocks and $10^5$ stochastic paths starting from the N\'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}$ as a function of $p$ 
+$p_\text{conv}$ is a measure of the convergence of $\cE^\text{DMC}(p)$ as a function of $p$. 
-(smallest value of $p$ for which the energy is roughly stabilized to five decimal places)
 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}(E=-1)$} \\
+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)$\\
@@ -1113,12 +1099,11 @@ $\cD(0,1)\cup\cD(1,0)\cup\cD$(2,0) & 36 & $\infty$&1&$-0.75272390$\\
 \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 statistical error on the largest components. 
+As explained above, it is very advantageous to calculate exactly as many $(H_p,S_p)$ as possible in order to avoid the statistical error on the largest components. 
-Ideally, this should be done up to the value of $p$ for which the calculation of these quantities whose cost increases exponentially is possible.
 Table \ref{tab2} shows the results both for the case of a single-state 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 worth computing exactly as many components as possible. 
-For the single-state domain, the statistical error is reduced by a factor of two when passing from no analytical computation, $p_\text{ex}=0$, to the case where eight components for $H_p$ and $S_p$ are computed exactly. 
+For the single-state domain, the statistical error is reduced by a factor two when passing from no analytical computation, $p_\text{ex}=0$, to the case where eight components for $H_p$ and $S_p$ are computed exactly. 
 For the best domain, the impact is much more important with a huge reduction of about three orders of magnitude in the statistical error. 
 
 Table \ref{tab3} reports the energies convergence as a function of $p$ alongside their statistical error on the last digit for $E=
@@ -1126,10 +1111,9 @@ Table \ref{tab3} reports the energies convergence as a function of $p$ alongside
 The values are displayed in Fig.~\ref{fig4}. 
 As seen, the behavior of $\cE$ as a function of $E$ is almost perfectly linear. 
 The extrapolated values obtained from the five values of the energy with three different fitting functions are reported. 
-Fitting the data using a simple linear function leads to an energy of $-0.7680282(5)$ (to be compared with the exact value of $-0.768068\ldots$). 
+Using the linear fitting function leads to an energy of $-0.7680282(5)$ (to be compared with the exact value of $-0.768068\ldots$). 
 A small bias of about $4 \times 10^{-5}$ is observed. 
-This bias vanishes within the statistical error when resorting to more flexible fitting functions, such as a quadratic function of $E$ or the
+This bias vanishes within the statistical error when one relies on the quadratic and two-component fitting functions. 
-two-component representation given by Eq.~\eqref{eq:2comp}.
 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. 
@@ -1137,25 +1121,24 @@ We also report the variational and exact energies together with the values of th
 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{tab5}, we report the results obtained for larger systems at $U=12$ for a chain size ranging from $N=4$ (36 states) 
+Finally, in Table \ref{tab4}, we report the results obtained for larger systems at $U=12$ for a chain size ranging from $N=4$ (36 states) 
 to $N=12$ ($\sim 10^6$ states). 
 No careful construction of domains maximizing the average trapping time has been performed, we have merely chosen domains of reasonable size (no more than 2682) by taking not too large number of double occupations (only, $n_D=0,1$) and not too small number of nearest-neighbor antiparallel pairs. 
 As seen, as the number of sites increases, the average trapping time for the chosen domains decreases. 
 This point is, of course, undesirable since the efficiency of the approach may gradually deteriorate when considering large systems. 
 A more elaborate way of constructing domains is clearly desirable in this case.
-The exact energies extrapolated using the two-component function are also reported. 
+The exact DMC energies extrapolated using the two-component function are also reported. 
 Similarly to what has been done for $N=4$, the extrapolation is performed using about five values of the reference energy. 
-The extrapolated DMC energies are in full agreement with the exact value within the error bar. 
+The extrapolated DMC energies are in full agreement with the exact value within error bars. 
-However, an increase in statistical error is observed when the system size increases. 
+However, an increase of the statistical error is observed when the system size increases. 
-To get lower error bars, more accurate trial wave functions may be considered, better domains, and also larger simulation times. 
+To get lower error bars, a more accurate trial wave functions may be considered, better domains, and also larger simulation times. 
 Of course, it will also be particularly interesting to take advantage of the fully parallelizable character of the algorithm to get much lower error bars. 
 All these aspects will be considered in a forthcoming work.
 
 %%% TABLE 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 $p_\text{ex}$ calculated analytically in the expression 
+Dependence of the statistical error on the energy with the number of $p$-components calculated analytically. 
-of the energy, $\cE^\text{DMC}(E,p_\text{ex},p_\text{max})$, Eq.~\eqref{eq:E_pex}.
 The simulation is performed the same way as in 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}
@@ -1181,11 +1164,11 @@ $8$ & $2.2  \times10^{-5}$ &$  0.05 \times 10^{-8}$\\
 \caption{One-dimensional Hubbard model with $N=4$, $U=12$, $\alpha=1.292$, $\beta=0.552$, and $p_\text{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$ reported in this table.}
+The fits are performed with the five values of $E$.}
 \label{tab3}
 \begin{ruledtabular}
 \begin{tabular}{lc}
-$E$ & $\cE^\text{DMC}(E)$ \\
+$E$ & $E^\text{DMC}$ \\
 \hline
 $-0.8$  &$-0.7654686(2)$\\
 $-0.795$&$-0.7658622(2)$\\
@@ -1202,15 +1185,12 @@ $E_0$ exact & $-0.768068\ldots$\\
 
 %%% TABLE IV %%%
 \begin{table}
-\caption{One-dimensional Hubbard model with $N=4$. Average trapping time as a 
+\caption{One-dimensional Hubbard model with $N=4$.
-function of $U$. $(\alpha,\beta)$ and $E_\text{v}$ are the parameters and variational energy
-of the trial wave function, respectively. $E_\text{ex}$ is the exact energy obtained by diagonalization of the 
-Hamiltonian matrix.
 The main domain is $\cD(0,1) \cup \cD(1,0)$.}
 \label{tab4}
 \begin{ruledtabular}
 \begin{tabular}{rcccr}
-$U$ & $(\alpha,\beta)$ & $E_\text{v}$ & $E_\text{ex}$  & $\bar{t}_{I_0}$ \\
+$U$ & $(\alpha,\beta)$ & $E_\text{T}$ & $E_\text{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\\
@@ -1225,20 +1205,18 @@ $U$ & $(\alpha,\beta)$ & $E_\text{v}$ & $E_\text{ex}$  & $\bar{t}_{I_0}$ \\
 
 %%% TABLE V %%%
 \begin{table*}
-\caption{Ground-state energy of the one-dimensional Hubbard model for different sizes and $U=12$. 
+\caption{One-dimensional Hubbard model with $U=12$. 
-The parameters $(\alpha,\beta)$ of the trial wave function and the corresponding variational energy $E_\text{v}$ are reported.
+The extrapolated DMC energies are obtained using the two-component fitting function.}
-$\bar{t}_{I_0}$ is the average trapping time for the N\'eel state.
-The extrapolated DMC energies, $E_\text{DMC}$, are obtained using the two-component fitting function, Eq.~\eqref{eq:2comp}. $E_\text{DMC}$ values are in full agreement with the exact values $E_\text{ex}$ computed by diagonalization of the Hamiltonian matrix.}
 \label{tab5}
 \begin{ruledtabular}
 \begin{tabular}{crcrccccr}
-$N$ & Hilbert space size & Domain & Domain size & $(\alpha,\beta)$ &$\bar{t}_{I_0}$ & $E_\text{v}$ & $E_\text{DMC}$ & $ E_\text{ex}$\\
+$N$ & Hilbert space size & Domain & Domain size & $(\alpha,\beta)$ &$\bar{t}_{I_0}$ & $E_\text{T}$ & $E_\text{ex}$ & $E^\text{DMC}$\\
 \hline
-4  & 36     & $\cD(0,1) \cup \cD(1,0)$ & 30 &(1.292,\,0.552)& 108.7  & $-0.495361$    &    $-0.7680676(5)$  &  $-0.768068$  \\
+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.215389(9)$  &  $-1.215395$  \\
+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.6637(2)$    &  $-1.66395$  \\
+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.1120(7)$     &  $-2.113089$ \\
+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.560(6)$      &  $-2.562529$  \\
+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*}
@@ -1327,6 +1305,16 @@ and
 \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},
 	& 
@@ -1363,6 +1351,10 @@ Using Eqs.~\eqref{eq:defA1} and \eqref{eq:defA2}, one gets, for $i = 1$ or $2$,
 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{eq:final}
 \ee