used for the nucleons, Greek ones refer to the Cartesian components of the
operators):

$\displaystyle A_{ij}^{[\tau]}$ $\displaystyle=v_{2}\left(r_{ij}\right)\;,$
$\displaystyle A_{i\alpha,j\beta}^{[\sigma]}$
$\displaystyle=v_{3}\left(r_{ij}\right)\delta_{\alpha\beta}+v_{5}\left(r_{ij}\right)\left(3\,\hat{r}_{ij}^{\alpha}\,\hat{r}_{ij}^{\beta}-\delta_{\alpha\beta}\right)\;,$
(3.127) $\displaystyle A_{i\alpha,j\beta}^{[\sigma\tau]}$
$\displaystyle=v_{4}\left(r_{ij}\right)\delta_{\alpha\beta}+v_{6}\left(r_{ij}\right)\left(3\,\hat{r}_{ij}^{\alpha}\,\hat{r}_{ij}^{\beta}-\delta_{\alpha\beta}\right)\;,$

that come from the decomposition of the operators in Cartesian coordinates:

$\displaystyle\bm{\sigma}_{i}\cdot\bm{\sigma}_{j}$
$\displaystyle=\sum_{\alpha\beta}\sigma_{i\alpha}\,\sigma_{j\beta}\,\delta_{\alpha\beta}\;,$
(3.128) $\displaystyle S_{ij}$
$\displaystyle=\sum_{\alpha\beta}\left(3\,\sigma_{i\alpha}\,\hat{r}_{ij}^{\alpha}\,\sigma_{j\beta}\,\hat{r}_{ij}^{\beta}-\sigma_{i\alpha}\,\sigma_{j\beta}\,\delta_{\alpha\beta}\right)\;.$
(3.129)

Being real and symmetric, the $A$ matrices have real eigenvalues and real
orthogonal eigenstates

$\displaystyle\sum_{j\beta}A_{i\alpha,j\beta}^{[\sigma]}\,\psi_{n,j\beta}^{[\sigma]}$
$\displaystyle=\lambda_{n}^{[\sigma]}\,\psi_{n,i\alpha}^{[\sigma]}\;,$
$\displaystyle\sum_{j\beta}A_{i\alpha,j\beta}^{[\sigma\tau]}\,\psi_{n,j\beta}^{[\sigma\tau]}$
$\displaystyle=\lambda_{n}^{[\sigma\tau]}\,\psi_{n,i\alpha}^{[\sigma\tau]}\;,$
(3.130) $\displaystyle\sum_{j}A_{ij}^{[\tau]}\,\psi_{n,j}^{[\tau]}$
$\displaystyle=\lambda_{n}^{[\tau]}\,\psi_{n,i}^{[\tau]}\;.$

Let us expand $\sigma_{i\alpha}$ on the complete set of eigenvectors
$\psi_{n,i\alpha}^{[\sigma]}$ of the matrix $A_{i\alpha,j\beta}^{[\sigma]}$ :

$\displaystyle\sigma_{i\alpha}=\sum_{n}c_{n}^{[\sigma]}\,\psi_{n,i\alpha}^{[\sigma]}=\sum_{n}\left(\sum_{j\beta}\psi_{n,j\beta}^{[\sigma]}\,\sigma_{j\beta}\right)\psi_{n,i\alpha}^{[\sigma]}\;,$
(3.131)

where we have used the orthogonality condition

$\displaystyle\sum_{i\alpha}\psi_{n,i\alpha}^{[\mathcal{O}]}\,\psi_{m,i\alpha}^{[\mathcal{O}]}=\delta_{nm}\;.$
(3.132)

Using Eq. (3.131) we can recast the first term of Eq. (3.126) in the following
form:

$\displaystyle\frac{1}{2}\sum_{i\alpha,j\beta}\sigma_{i\alpha}\,A_{i\alpha,j\beta}^{[\sigma]}\,\sigma_{j\beta}=$
$\displaystyle\,=\frac{1}{2}\sum_{i\alpha,j\beta}\left\\{\left[\sum_{n}\left(\sum_{k\gamma}\sigma_{k\gamma}\,\psi_{n,k\gamma}^{[\sigma]}\right)\psi_{n,i\alpha}^{[\sigma]}\right]A_{i\alpha,j\beta}^{[\sigma]}\left[\sum_{m}\left(\sum_{k\gamma}\sigma_{k\gamma}\,\psi_{m,k\gamma}^{[\sigma]}\right)\psi_{m,j\beta}^{[\sigma]}\right]\right\\}\;,$
$\displaystyle\,=\frac{1}{2}\sum_{i\alpha}\left\\{\left[\sum_{n}\left(\sum_{k\gamma}\sigma_{k\gamma}\,\psi_{n,k\gamma}^{[\sigma]}\right)\psi_{n,i\alpha}^{[\sigma]}\right]\left[\sum_{m}\lambda_{m}^{[\sigma]}\left(\sum_{k\gamma}\sigma_{k\gamma}\,\psi_{m,k\gamma}^{[\sigma]}\right)\psi_{m,i\alpha}^{[\sigma]}\right]\right\\}\;,$
$\displaystyle\,=\frac{1}{2}\sum_{n}\left(\sum_{k\gamma}\sigma_{k\gamma}\,\psi_{n,k\gamma}^{[\sigma]}\right)^{2}\\!\lambda_{n}^{[\sigma]}\;.$
(3.133)

Similar derivation can be written for the terms
$\tau_{i\gamma}\,\sigma_{i\alpha}$ and $\tau_{i\gamma}$ and we can define a
new set of operators expressed in terms of the eigenvectors of the matrices
$A$:

$\displaystyle\mathcal{O}_{n}^{[\sigma]}$
$\displaystyle=\sum_{j\beta}\sigma_{j\beta}\,\psi_{n,j\beta}^{[\sigma]}\;,$
$\displaystyle\mathcal{O}_{n,\alpha}^{[\sigma\tau]}$
$\displaystyle=\sum_{j\beta}\tau_{j\alpha}\,\sigma_{j\beta}\,\psi_{n,j\beta}^{[\sigma\tau]}\;,$
(3.134) $\displaystyle\mathcal{O}_{n,\alpha}^{[\tau]}$
$\displaystyle=\sum_{j}\tau_{j\alpha}\,\psi_{n,j}^{[\tau]}\;.$

The spin dependent part of the $NN$ interaction can be thus expressed as
follows:

$\displaystyle\\!\\!\\!\\!V_{NN}^{SD}=\frac{1}{2}\sum_{n=1}^{3A}\lambda_{n}^{[\sigma]}\\!\left(\mathcal{O}_{n}^{[\sigma]}\right)^{2}\\!+\frac{1}{2}\sum_{n=1}^{3A}\sum_{\alpha=1}^{3}\lambda_{n}^{[\sigma\tau]}\\!\left(\mathcal{O}_{n\alpha}^{[\sigma\tau]}\right)^{2}\\!+\frac{1}{2}\sum_{n=1}^{A}\sum_{\alpha=1}^{3}\lambda_{n}^{[\tau]}\\!\left(\mathcal{O}_{n\alpha}^{[\tau]}\right)^{2}\,.\\!$
(3.135)

$V_{NN}^{SD}$ is now written in a suitable form for the application of the
Hubbard-Stratonovich transformation of Eq. (3.105). The propagator
$\operatorname{e}^{-V_{NN}^{SD}\,d\tau}$ can be finally recast as:

$\displaystyle\operatorname{e}^{-\frac{1}{2}\sum_{n}\lambda_{n}(\mathcal{O}_{n})^{2}d\tau}$
$\displaystyle=\prod_{n}\operatorname{e}^{-\frac{1}{2}\lambda_{n}(\mathcal{O}_{n})^{2}d\tau}\,+\operatorname{o}\left(d\tau^{2}\right)\;,$
$\displaystyle\simeq\prod_{n}\frac{1}{\sqrt{2\pi}}\int\\!dx_{n}\operatorname{e}^{\frac{-x_{n}^{2}}{2}+\sqrt{-\lambda_{n}d\tau}\,x_{n}\mathcal{O}_{n}}\;,$
(3.136)

where we have used the compact notation $\mathcal{O}_{n}$ to denote the $3A$
$\mathcal{O}_{n}^{[\sigma]}$, the $9A$ $O_{n,\alpha}^{[\sigma\tau]}$ and the
$3A$ $\mathcal{O}_{n,\alpha}^{[\tau]}$ operators including the summation over
the coordinate index $\alpha$ where needed. The first step of the above
equation comes to the fact that in general the operators $\mathcal{O}_{n}$ do
not commute and so, due to Eq. (3.17), the equality is correct only at order
$d\tau^{2}$.

The standard DMC imaginary time propagation of Eq. (3.9) needs to be extended
to the spin-isospin space, as done in the GFMC algorithm via the projection of
Eqs. (3.82) and (3.83). In the AFDMC method, spin-isospin coordinates
$\\{S\\}$ are added to the spacial coordinates $\\{R\\}$, defining a set of
walkers which represents the single-particle wave function to be evolved in
imaginary time:

$\displaystyle\psi(R,S,\tau+d\tau)=\int
dR^{\prime}dS^{\prime}\,G(R,S,R^{\prime},S^{\prime},d\tau)\,\psi(R^{\prime},S^{\prime},\tau)\;.$
(3.137)

Including the integration over the Hubbard-Stratonovich auxiliary fields, the
Auxiliary Field DMC Green’s function reads (recall Eqs. (3.18) and (3.19)):

$\displaystyle G(R,S,R^{\prime},S^{\prime},d\tau)$ $\displaystyle=\langle
R,S|\operatorname{e}^{-(T+V-E_{T})d\tau}|R^{\prime},S^{\prime}\rangle\;,$
$\displaystyle\simeq\left(\frac{1}{4\pi
Dd\tau}\right)^{\frac{3\mathcal{N}}{2}}\\!\operatorname{e}^{-\frac{(R-R^{\prime})^{2}}{4Dd\tau}}\operatorname{e}^{-\left(\frac{V_{NN}^{SI}(R)+V_{NN}^{SI}(R^{\prime})}{2}-E_{T}\right)d\tau}\times$
$\displaystyle\quad\,\times\langle
S|\prod_{n=1}^{15A}\frac{1}{\sqrt{2\pi}}\int\\!dx_{n}\operatorname{e}^{\frac{-x_{n}^{2}}{2}+\sqrt{-\lambda_{n}d\tau}\,x_{n}\mathcal{O}_{n}}|S^{\prime}\rangle\;,$
(3.138)

Each operator $\mathcal{O}_{n}$ involves the sum over the particle index $j$,
as shown in Eq. (3.134). However, in the $A$-body tensor product space, each
$j$ sub-operator is a one-body operator acting on a different single particle
spin-isospin states, as in Eq. (3.123). Therefore the $j$-dependent terms
commute and we can represent the exponential of the sum over $j$ as a tensor
product of exponentials. The result is that the propagation of a spin-isospin
state $|S^{\prime}\rangle$ turns into a product of independent rotations, one
for each spin-isospin state. Considering just a spin wave function we have for
example:

$\displaystyle\operatorname{e}^{\sqrt{-\lambda_{n}d\tau}\,x_{n}\mathcal{O}_{n}^{[\sigma]}}|S^{\prime}\rangle=$
$\displaystyle\,=\operatorname{e}^{\sqrt{-\lambda_{n}d\tau}\,x_{n}\sum_{\beta}\sigma_{1\beta}\,\psi_{n,1\beta}^{[\sigma]}}\left(\begin{array}[]{c}a_{1\uparrow}\\\
a_{1\downarrow}\end{array}\right)_{1}\\!\otimes\cdots\otimes\operatorname{e}^{\sqrt{-\lambda_{n}d\tau}\,x_{n}\sum_{\beta}\sigma_{A\beta}\,\psi_{n,A\beta}^{[\sigma]}}\left(\begin{array}[]{c}a_{A\uparrow}\\\
a_{A\downarrow}\end{array}\right)_{A}\;,$ (3.143)
$\displaystyle\,=\left(\begin{array}[]{c}\widetilde{a}_{1\uparrow}\\\
\widetilde{a}_{1\downarrow}\end{array}\right)_{1}\\!\otimes\cdots\otimes\left(\begin{array}[]{c}\widetilde{a}_{A\uparrow}\\\
\widetilde{a}_{A\downarrow}\end{array}\right)_{A}\;.$ (3.148)

We can thus propagate spin-isospin dependent wave functions remaining inside
the space of single particle states.

The new coefficients $\widetilde{a}_{j\uparrow}$ and
$\widetilde{a}_{j\downarrow}$ are calculated at each time step for each
$\mathcal{O}_{n}$ operator. For neutron systems, i.e. for two-component
spinors for which only the operator $\mathcal{O}_{n}^{[\sigma]}$ is active,
there exists an explicit expression for these coefficients. Consider the
Landau relations

$\displaystyle\operatorname{e}^{i\,\vec{b}\cdot\vec{\sigma}}$
$\displaystyle=\cos(|\vec{b}|)+i\sin(|\vec{b}|)\frac{\vec{b}\cdot\vec{\sigma}}{|\vec{b}|}\;,$
(3.149) $\displaystyle\operatorname{e}^{\vec{b}\cdot\vec{\sigma}}$
$\displaystyle=\cosh(|\vec{b}|)+\sinh(|\vec{b}|)\frac{\vec{b}\cdot\vec{\sigma}}{|\vec{b}|}\;,$
(3.150)

and identify the $\vec{b}$ vector with

$\displaystyle\vec{b}=\sqrt{|\lambda_{n}|d\tau}\,x_{n}\vec{\psi}_{n,j}^{\,[\sigma]}\quad\quad\quad
b_{\beta}=\sqrt{|\lambda_{n}|d\tau}\,x_{n}\psi_{n,j\beta}^{\,[\sigma]}\;.$
(3.151)

The following expressions for the coefficients of the rotated spinors can be
then written, distinguishing the case $\lambda_{n}<0$ (Eq. (3.154)) and the
case $\lambda_{n}>0$ (Eq. (3.157)):

$\displaystyle\begin{array}[]{l}\widetilde{a}_{j\uparrow}\\!=\\!\\!\Biggl{[}\cosh(|\vec{b}|)+\sinh(|\vec{b}|)\operatorname{sgn}\\!\left(x_{n}\right)\frac{\psi_{n,jz}^{[\sigma]}}{|\vec{\psi}_{n,j}^{\,[\sigma]}|}\Biggr{]}a_{j\uparrow}\\!+\sinh(|\vec{b}|)\operatorname{sgn}\\!\left(x_{n}\right)\\!\\!\Biggl{[}\frac{\psi_{n,jx}^{[\sigma]}-i\,\psi_{n,jy}^{[\sigma]}}{|\vec{\psi}_{n,j}^{\,[\sigma]}|}\Biggr{]}a_{j\downarrow}\;,\\!\\!\\\\[13.00005pt]
\widetilde{a}_{j\downarrow}\\!=\\!\\!\Biggl{[}\cosh(|\vec{b}|)-\sinh(|\vec{b}|)\operatorname{sgn}\\!\left(x_{n}\right)\frac{\psi_{n,jz}^{[\sigma]}}{|\vec{\psi}_{n,j}^{\,[\sigma]}|}\Biggr{]}a_{j\downarrow}\\!+\sinh(|\vec{b}|)\operatorname{sgn}\\!\left(x_{n}\right)\\!\\!\Biggl{[}\frac{\psi_{n,jx}^{[\sigma]}+i\,\psi_{n,jy}^{[\sigma]}}{|\vec{\psi}_{n,j}^{\,[\sigma]}|}\Biggr{]}a_{j\uparrow}\;,\\!\\!\end{array}$
(3.154)
$\displaystyle\begin{array}[]{l}\widetilde{a}_{j\uparrow}\\!=\\!\\!\Biggl{[}\cos(|\vec{b}|)+i\,\sin(|\vec{b}|)\operatorname{sgn}\\!\left(x_{n}\right)\frac{\psi_{n,jz}^{[\sigma]}}{|\vec{\psi}_{n,j}^{\,[\sigma]}|}\Biggr{]}a_{j\uparrow}\\!+\sin(|\vec{b}|)\operatorname{sgn}\\!\left(x_{n}\right)\\!\\!\Biggl{[}\frac{i\,\psi_{n,jx}^{[\sigma]}+\psi_{n,jy}^{[\sigma]}}{|\vec{\psi}_{n,j}^{\,[\sigma]}|}\Biggr{]}a_{j\downarrow}\;,\\!\\!\\\\[13.00005pt]
\widetilde{a}_{j\downarrow}\\!=\\!\\!\Biggl{[}\cos(|\vec{b}|)+i\,\sin(|\vec{b}|)\operatorname{sgn}\\!\left(x_{n}\right)\frac{\psi_{n,jz}^{[\sigma]}}{|\vec{\psi}_{n,j}^{\,[\sigma]}|}\Biggr{]}a_{j\downarrow}\\!+\sin(|\vec{b}|)\operatorname{sgn}\\!\left(x_{n}\right)\\!\\!\Biggl{[}\frac{i\,\psi_{n,jx}^{[\sigma]}-\psi_{n,jy}^{[\sigma]}}{|\vec{\psi}_{n,j}^{\,[\sigma]}|}\Biggr{]}a_{j\uparrow}\;.\\!\\!\end{array}$
(3.157)

When we are dealing with the full 4-dimension single particle spinors, the
four coefficients $\widetilde{a}$ do not have an explicit expression. The
exponential of the $\mathcal{O}_{n}$ operators acting on the spinors is
calculated via a diagonalization procedure. Consider the general $4\times 4$
rotation matrix $B_{j}$ and its eigenvectors $\Psi_{m,j}\neq 0$ and
eigenvalues $\mu_{m,j}$:

$\displaystyle
B_{j}\,\Psi_{m,j}=\mu_{m,j}\,\Psi_{m,j}\quad\Rightarrow\quad\Psi_{m,j}^{-1}\,B_{j}\,\Psi_{m,j}=\mu_{m,j}\quad\quad
m=1,\ldots,4\;.$ (3.158)

Using the formal notation $\vec{\Psi}_{j}$ and $\vec{\mu}_{j}$ to denote the
$4\times 4$ matrix of eigenvectors and the 4-dimension vector of eigenvalues,
it is possible to write the action of $\operatorname{e}^{B_{j}}$ on a
4-dimensional single particle spinor $|S^{\prime}\rangle_{j}$ as follows:

$\displaystyle\operatorname{e}^{B_{j}}|S^{\prime}\rangle_{j}$
$\displaystyle=\vec{\Psi}_{j}\vec{\Psi}_{j}^{-1}\operatorname{e}^{B_{j}}\vec{\Psi}_{j}\vec{\Psi}_{j}^{-1}|S^{\prime}\rangle_{j}\;,$
$\displaystyle=\vec{\Psi}_{j}\operatorname{e}^{\vec{\Psi}_{j}^{-1}B_{j}\,\vec{\Psi}_{j}}\vec{\Psi}_{j}^{-1}|S^{\prime}\rangle_{j}\;,$
$\displaystyle=\vec{\Psi}_{j}\operatorname{e}^{\text{diag}\left(\vec{\mu}_{j}\right)}\vec{\Psi}_{j}^{-1}|S^{\prime}\rangle_{j}\;,$
$\displaystyle=\vec{\Psi}_{j}\,\text{diag}\left(\operatorname{e}^{\vec{\mu}_{j}}\right)\vec{\Psi}_{j}^{-1}|S^{\prime}\rangle_{j}\;.$
(3.159)

Each component of the rotated spinor $|\widetilde{S}^{\prime}\rangle_{j}$ is
thus derived from the eigenvectors and eigenvalues of the rotation matrix
$B_{j}$, which is built starting from the $\mathcal{O}_{n}$ operators. Moving
from neutrons to nucleons, i.e. adding the isospin degrees of freedom to the
system, the computational time spent to rotate each single particle spin-
isospin state during the propagation is increased by the time for the
diagonalization of the $4\times 4$ Hubbard-Stratonovich rotation matrices.
However, the total time for the propagation of the wave function as $A$
becomes large, is dominated by the diagonalization of the potential matrices.
Since the cost of this operation goes as the cube of the number of matrix rows
(columns), the AFDMC computational time is proportional to $A^{3}$, which is
much slower than the scaling factor $A!$ of GFMC.

In addition to the diagonalization of the AV6 potential matrices and the
spinor rotation matrices, we have to deal with the evaluation of the integral
over the auxiliary fields $x_{n}$. The easiest way, in the spirit of Monte
Carlo, is to sample the auxiliary fields from the Gaussian of Eq. (3.138),
which is interpreted as a probability distribution. The sampled values are
then used to determine the action of the operators on the spin-isospin part of
the wave function as described above. The integral over all the spin-isospin
rotations induced by the auxiliary fields eventually recovers the action of
the quadratic spin-isospin operators on a trial wave function containing all
the possible good spin-isospin states, as the GFMC one.

In this scheme, the integration over the auxiliary fields is performed jointly
with the integration over the coordinates. This generally leads to a large
variance. The integral of Eq. (3.138) should be indeed evaluated for each
sampled position and not simply estimated “on the fly”. A more refined
algorithm, in which for each sampled configuration the integral over $x_{n}$
is calculated by sampling more than one auxiliary variable, has been tested.
The energy values at convergence are the same for both approaches. However, in
the latter case the variance is much reduced, although the computational time
for each move is increased due to the iteration over the newly sampled
auxiliary points.

As done in the DMC method, see § 3.1.1, we can introduce an importance
function to guide the diffusion in the coordinate space also in the AFDMC
algorithm. The drift term (3.32) is added to the $R-R^{\prime}$ Gaussian
distribution of Eq. (3.138) and the branching weight $\widetilde{\omega}_{i}$
is given by the local energy as in Eq. (3.33). The idea of the importance
sampling can be applied to guide the rotation of the spin-isospin states in
the Hubbard-Stratonovich transformation. This can be done by properly shifting
the Gaussian over the auxiliary fields of Eq. (3.138) by means of a drift term
$\bar{x}_{n}$:

$\displaystyle\operatorname{e}^{-\frac{x_{n}^{2}}{2}+\sqrt{-\lambda_{n}d\tau}\,x_{n}\mathcal{O}_{n}}=\operatorname{e}^{-\frac{(x_{n}-\bar{x}_{n})^{2}}{2}}\operatorname{e}^{\sqrt{-\lambda_{n}d\tau}\,x_{n}\mathcal{O}_{n}}\operatorname{e}^{-\bar{x}_{n}\left(x_{n}-\frac{\bar{x}_{n}}{2}\right)}\;,$
(3.160)

where

$\displaystyle\bar{x}_{n}=\operatorname{Re}\left[\sqrt{-\lambda_{n}d\tau}\langle\mathcal{O}_{n}\rangle_{m}\right]\;,$
(3.161)

and $\langle\mathcal{O}_{n}\rangle_{m}$ is the mixed expectation value of
$\mathcal{O}_{n}$ (Eq. (3.24)) calculated on the old spin-isospin
configurations. The mixed estimator is introduced in order to guide the
rotations, by maximizing the overlap between the walker and the trial
function, which is not generally picked around $x_{n}=0$.

The last factor of Eq. (3.160) can be interpreted as an additional weight term
that has to be included in the total weight. By combining diffusion, rotation
and all the additional factors we can derive two different algorithms.

  * _v1_

In the first one, the ratio between the importance functions in the new and
old configurations (see Eq. (3.28)) is kept explicit. However the drifted
Gaussian $\widetilde{G}_{0}(R,R^{\prime},d\tau)$ of Eq. (3.31) is used for the
diffusion in the coordinate space and the drifted Gaussian of Eq. (3.160) for
the sampling of auxiliary fields. The weight for the branching process
$\omega_{i}$ defined in Eq. (3.21) takes then an overall factor

$\displaystyle\frac{\langle\psi_{I}|RS\rangle}{\langle\psi_{I}|R^{\prime}S^{\prime}\rangle}\operatorname{e}^{-\frac{d(R^{\prime})\left[d(R^{\prime})+2(R-R^{\prime})\right]}{4Dd\tau}}\prod_{n}\operatorname{e}^{-\bar{x}_{n}\left(x_{n}-\frac{\bar{x}_{n}}{2}\right)}\;,$
(3.162)

due to the counter terms coming from the coordinate drift
$d(R)=\bm{v}_{d}(R)Dd\tau$ added in the original $G_{0}(R,R^{\prime},d\tau)$
and from the auxiliary field shift $\bar{x}_{n}$.

  * _v2_

The second algorithm corresponds the local energy scheme described in § 3.1.1.
Again the coordinates are diffused via the drifted Gaussian
$\widetilde{G}_{0}(R,R^{\prime},d\tau)$ of Eq. (3.31) and the auxiliary fields
are sampled from the shifted Gaussian of Eq. (3.160). The branching weight
$\widetilde{w}_{i}$ is instead given by the local energy as in Eq. (3.33). The
counter terms related to $\bar{x}_{n}$ are automatically included in the
weight because the local energy
$E_{L}(R,S)=\frac{H\psi_{I}(R,S)}{\psi_{I}(R,S)}$ takes now contributions from
all the spin-isospin operators of the full potential $V_{NN}$. Actually, the
term $\operatorname{e}^{-\bar{x}_{n}x_{n}}$ vanishes during the auxiliary
field integration because $x_{n}$ can take positive and negative values. The
term $\frac{\bar{x}_{n}^{2}}{2}$ is nothing but the
$-\frac{1}{2}\lambda_{n}\langle\mathcal{O}_{n}\rangle_{m}^{2}d\tau$
contribution already included in the weight via $E_{L}(R,S)$.

Given the same choice for the drift term, that depends, for example, on the
constraint applied to deal with the sign problem, the two algorithms are
equivalent and should sample the same Green’s function.

In both versions, the steps that constitute the AFDMC algorithm are almost the
same of the DMC one, reported in § 3.1. The starting point is the initial
distribution of walkers, step 1. In step 2 the diffusion of the coordinates is
performed including the drift factor. Now also the spin-isospin degrees of
freedom are propagated, by means of the Hubbard-Stratonovich rotations and the
integral over the auxiliary fields. As in step 3, a weight is assigned to each
walker, choosing one of the two equivalent solutions proposed above (explicit
$\psi_{I}$ ratio or local energy). Both propagation and weight depend on the
prescription adopted in order to keep under control the sign problem. Usually
the fixed phase approximation (see § 3.1.2) is applied with the evaluation of
local operators. The branching process follows then the DMC version described
in step 4 and the procedure is iterated in the same way with the computation
of expectation values at convergence.

####  3.2.2 Propagator for neutrons: spin-orbit terms

In the previous section we have seen how to deal in an efficient way with a
propagator containing the first six components of the Argonne two-body
potential. Next terms in Eq. (2.16) are the spin-orbit contributions for
$p=7,8$. Although an attempt to treat the spin-orbit terms for nucleon systems
has been reported by Armani in his Ph.D. thesis [224] (together with a
possible $\bm{L}_{ij}^{2}$ inclusion for $p=9$), at present the
$\bm{L}_{ij}\cdot\bm{S}_{ij}$ operator is consistently employed in the AFDMC
algorithm only for neutron systems. No other terms of the $NN$ interaction are
included in the full propagator, neither for nucleons nor for neutrons,
although a perturbative treatment of the remaining terms of AV18 is also
possible [136]. The full derivation of the neutron spin-orbit propagator is
reported in Ref. [37]. Here we want just to sketch the idea behind the
treatment of this non local term for which the corresponding Green’s function
is not trivial to be derived.

Consider the spin-orbit potential for neutrons:

$\displaystyle
v_{ij}^{LS}=v_{LS}(r_{ij})\,\bm{L}_{ij}\cdot\bm{S}_{ij}=v_{LS}(r_{ij})\left(\bm{L}\cdot\bm{S}\right)_{ij}\;,$
(3.163)

where

$\displaystyle v_{LS}(r_{ij})=v_{7}(r_{ij})+v_{8}(r_{ij})\;,$ (3.164)

and $\bm{L}_{ij}$ and $\bm{S}_{ij}$ are defined respectively by Eqs. (2.19)
and (2.20). As reported in Ref. [228], one way to evaluate the propagator for
$\bm{L}\cdot\bm{S}$ is to consider the expansion at first order in $d\tau$

$\displaystyle\operatorname{e}^{-v_{LS}(r_{ij})\left(\bm{L}\cdot\bm{S}\right)_{ij}d\tau}\simeq\left[1-v_{LS}(r_{ij})\left(\bm{L}\cdot\bm{S}\right)_{ij}d\tau\right]\;,$
(3.165)

acting on the free propagator $G_{0}(R,R^{\prime},d\tau)$ of Eq. (3.11). The
derivatives terms of the above expression give

$\displaystyle\left(\bm{\nabla}_{i}-\bm{\nabla}_{j}\right)G_{0}(R,R^{\prime},d\tau)=-\frac{1}{2Dd\tau}\left(\Delta\bm{r}_{i}-\Delta\bm{r}_{j}\right)G_{0}(R,R^{\prime},d\tau)\;,$
(3.166)

where $\Delta\bm{r}_{i}=\bm{r}_{i}-\bm{r}^{\prime}_{i}$ . We can then write:

$\displaystyle\left(\bm{L}\cdot\bm{S}\right)_{ij}G_{0}(R,R^{\prime},d\tau)=$
$\displaystyle\quad\quad=-\frac{1}{4i}\frac{1}{2Dd\tau}\left(\bm{r}_{i}-\bm{r}_{j}\right)\times\left(\Delta\bm{r}_{i}-\Delta\bm{r}_{j}\right)\cdot\left(\bm{\sigma}_{i}+\bm{\sigma}_{j}\right)G_{0}(R,R^{\prime},d\tau)\;,$
$\displaystyle\quad\quad=-\frac{1}{4i}\frac{1}{2Dd\tau}\left(\bm{\Sigma}_{ij}\times\bm{r}_{ij}\right)\cdot\left(\Delta\bm{r}_{i}-\Delta\bm{r}_{j}\right)G_{0}(R,R^{\prime},d\tau)\;,$
(3.167)

where $\bm{\Sigma}_{ij}=\bm{\sigma}_{i}+\bm{\sigma}_{j}$ and
$\bm{r}_{ij}=\bm{r}_{i}-\bm{r}_{j}$, and the relation
$\bm{a}\cdot\left(\bm{b}\times\bm{c}\right)=\bm{c}\cdot\left(\bm{a}\times\bm{b}\right)$
has been used.

By inserting the last expression in Eq. (3.165) and re-exponentiating,
including also the omitted sum over particle indices $i$ and $j$, the
following propagator is obtained:

$\displaystyle\mathcal{P}_{LS}\simeq\operatorname{e}^{\sum_{i\neq
j}\frac{v_{LS}(r_{ij})}{8iD}\left(\bm{\Sigma}_{ij}\times\bm{r}_{ij}\right)\cdot\left(\Delta\bm{r}_{i}-\Delta\bm{r}_{j}\right)}\;.$
(3.168)

The effect of $\mathcal{P}_{LS}$ can be studied starting from the formal
solution

$\displaystyle\psi(R,S,\tau+d\tau)\stackrel{{\scriptstyle LS}}{{\simeq}}\int
dR^{\prime}dS^{\prime}\,G_{0}(R,R^{\prime},d\tau)\,\mathcal{P}_{LS}\,\psi(R^{\prime},S^{\prime},\tau)\;,$
(3.169)

and expanding the propagator to the second order and the wave function
$\psi(R^{\prime},S^{\prime},\tau)$ to the first order in $R-R^{\prime}$. It is
possible to show (see Ref. [37] for the details) that the spin-orbit
contribution of the propagator takes a simple form, but two- and three-body
extra corrections appear. However, in the case of neutrons these additional
terms contain quadratic spin operators and so they can be handled by the
Hubbard-Stratonovich transformation and the rotations over new auxiliary
fields.

####  3.2.3 Propagator for neutrons: three-body terms

As reported in § 2.1.2, the Illinois (Urbana IX) TNI can be written as a sum
of four different terms:

$\displaystyle
V_{ijk}=A_{2\pi}^{P}\,\mathcal{O}^{2\pi,P}_{ijk}+A_{2\pi}^{S}\,\mathcal{O}^{2\pi,S}_{ijk}+A_{3\pi}\,\mathcal{O}^{3\pi}_{ijk}+A_{R}\,\mathcal{O}^{R}_{ijk}\;.$
(3.170)

For neutron systems, being $\bm{\tau}_{i}\cdot\bm{\tau}_{j}=1$, the operator
structure simplify in such a way that $V_{ijk}$ can be recast as a sum of two-
body terms only [37, 33]. We can therefore handle also the TNI in the AFDMC
propagator by means of the Hubbard-Stratonovich transformation. Let analyze
how each term of the above relation can be conveniently rewritten for neutron
systems.

  * •

$\mathcal{O}^{2\pi,P}_{ijk}$  _term_. The $P$-wave 2$\pi$ exchange term (and
also the 3$\pi$ exchange one) of Eq. (2.26) includes the OPE operator
$X_{ij}$, defined in Eq. (2.7). $X_{ij}$ involves the
$\bm{\sigma}_{i}\cdot\bm{\sigma}_{j}$ and the $S_{ij}$ operators that can be
decomposed via Eqs. (3.128) and (3.129) in order to define a $3A\times 3A$
matrix $X_{i\alpha,j\beta}$ analogous to the $A_{i\alpha,j\beta}^{[\sigma]}$
of Eq. (3.127), where $v_{3}(r_{ij})\\!\rightarrow Y_{\pi}(r_{ij})$ and
$v_{5}(r_{ij})\\!\rightarrow T_{\pi}(r_{ij})$. The OPE operator can be thus
expressed as

$\displaystyle
X_{ij}=\sigma_{i\alpha}\,X_{i\alpha,j\beta}\,\sigma_{j\beta}\;,$ (3.171)

where the matrix $X_{i\alpha,j\beta}$ is real with zero diagonal and has the
same symmetries of $A_{i\alpha,j\beta}^{[\sigma]}$. The commutator over the
$\bm{\tau}_{i}$ operators vanishes, while the anticommutator gives simply a
factor 2. Recalling that $X_{ij}=X_{ji}$ we can derive the following relation:

$\displaystyle\sum_{i<j<k}\mathcal{O}^{2\pi,P}_{ijk}$
$\displaystyle=\frac{1}{3!}\sum_{i\neq j\neq
k}\sum_{cyclic}2\phantom{\frac{1}{4}}\\!\\!\\!\\!\Bigl{\\{}X_{ij},X_{jk}\Bigr{\\}}\;,$
$\displaystyle=2\sum_{i\neq j\neq k}X_{ik}X_{kj}\;,$
$\displaystyle=2\sum_{i\neq
j}\sum_{\alpha\beta}\sigma_{i\alpha}\Biggl{(}\sum_{k\gamma}X_{i\alpha,k\gamma}\,X_{k\gamma,j\beta}\Biggr{)}\sigma_{j\beta}\;,$
$\displaystyle=2\sum_{i\neq
j}\sum_{\alpha\beta}\sigma_{i\alpha}\,X^{2}_{i\alpha,j\beta}\,\sigma_{j\beta}\;.$
(3.172)

  * •

$\mathcal{O}^{2\pi,S}_{ijk}$  _term_. In the $S$-wave TPE term the isospin
operators do not contribute and we are left with

$\displaystyle\sum_{i<j<k}\mathcal{O}_{ijk}^{2\pi,S}$
$\displaystyle=\frac{1}{3!}\sum_{i\neq j\neq
k}\sum_{cyclic}Z_{\pi}(r_{ij})Z_{\pi}(r_{jk})\,\bm{\sigma}_{i}\cdot\hat{\bm{r}}_{ij}\,\bm{\sigma}_{k}\cdot\hat{\bm{r}}_{kj}\;,$
$\displaystyle=\frac{1}{2}\sum_{i\neq
j}\sum_{\alpha\beta}\sigma_{i\alpha}\Biggl{[}\sum_{k}Z_{\pi}(r_{ik})\,\hat{r}_{ik}^{\alpha}\,Z_{\pi}(r_{jk})\,\hat{r}_{jk}^{\beta}\Biggr{]}\sigma_{j\beta}\;,$
$\displaystyle=\frac{1}{2}\sum_{i\neq
j}\sum_{\alpha\beta}\sigma_{i\alpha}\,Z_{i\alpha,j\beta}\,\sigma_{j\beta}\;.$
(3.173)

  * •

$\mathcal{O}^{3\pi}_{ijk}$  _term_. The 3$\pi$ exchange term, even with the
isospin reduction for neutrons, still keeps a very complicated operator
structure. As reported in Ref. [33], this factor can be conveniently written
as a sum of a spin independent and a spin dependent components

$\displaystyle\sum_{i<j<k}\mathcal{O}_{ijk}^{3\pi}=V_{c}^{3\pi}+V_{\sigma}^{3\pi}\;,$
(3.174)

with

$\displaystyle V_{c}^{3\pi}$ $\displaystyle=\frac{400}{18}\sum_{i\neq
j}X_{i\alpha,j\beta}^{2}\,X_{i\alpha,j\beta}\;,$ (3.175) $\displaystyle
V_{\sigma}^{3\pi}$ $\displaystyle=\frac{200}{54}\sum_{i\neq
j}\sum_{\alpha\beta}\sigma_{i\alpha}\Biggl{(}\sum_{\gamma\delta\mu\nu}X_{i\gamma,j\mu}^{2}\,X_{i\delta,j\nu}\,\varepsilon_{\alpha\gamma\delta}\,\varepsilon_{\beta\mu\nu}\Biggr{)}\sigma_{j\beta}\;,$
$\displaystyle=\frac{200}{54}\sum_{i\neq
j}\sum_{\alpha\beta}\sigma_{i\alpha}\,W_{i\alpha,j\beta}\,\sigma_{j\beta}\;,$
(3.176)

where $\varepsilon_{\alpha\beta\gamma}$ is the full antisymmetric tensor.

  * •

$\mathcal{O}^{R}_{ijk}$  _term_. The last spin independent term can be recast
as a two body operator as follows

$\displaystyle\sum_{i<j<k}\mathcal{O}^{R}_{ijk}=G_{0}^{R}+\frac{1}{2}\sum_{i}\left(G_{i}^{R}\right)^{2}\;,$
(3.177)

with

$\displaystyle G_{0}^{R}$ $\displaystyle=-\sum_{i<j}T_{\pi}^{4}(r_{ij})\;,$
(3.178) $\displaystyle G_{i}^{R}$ $\displaystyle=\sum_{k\neq
i}T_{\pi}^{2}(r_{ik})\;.$ (3.179)

Finally, for neutron systems we can still write the spin dependent part of the
$NN$ potential in the form of Eq. (3.126), with the inclusion of TNI
contributions:

$\displaystyle V_{NN}^{SD}$ $\displaystyle=\frac{1}{2}\sum_{i\neq
j}\sum_{\alpha\beta}\sigma_{i\alpha}\,A_{i\alpha,j\beta}^{[\sigma]}\,\sigma_{j\beta}\;,$
(3.180)

where now

$\displaystyle A_{i\alpha,j\beta}^{[\sigma]}\longrightarrow
A_{i\alpha,j\beta}^{[\sigma]}+2A_{2\pi}^{P}\,X^{2}_{i\alpha,j\beta}+\frac{1}{2}A_{2\pi}^{S}\,Z_{i\alpha,j\beta}+\frac{200}{54}\,A_{3\pi}\,W_{i\alpha,j\beta}\;.$
(3.181)

The central term of the two-body potential of Eq. (3.125) keeps also
contributions from the TNI 3$\pi$ exchange term and from the phenomenological
term, and it reads now:

$\displaystyle V_{NN}^{SI}\longrightarrow
V_{NN}^{SI}+A_{3\pi}V_{c}^{3\pi}+A_{R}\left[G_{0}^{R}+\frac{1}{2}\sum_{i}\left(G_{i}^{R}\right)^{2}\right]\;.$
(3.182)

####  3.2.4 Wave functions

In this section the trial wave functions used in AFDMC calculations for
nuclear and hypernuclear systems will be presented, distinguishing between the
case of finite and infinite systems. Restoring the convention of Chapter 2,
which is commonly used in the literature for hypernuclear systems, $A$ will
refer to the total number of baryons, $\mathcal{N}_{N}$ nucleons plus
$\mathcal{N}_{\Lambda}$ lambda particles. Latin indices will be used for the
nucleons, Greek $\lambda$, $\mu$ and $\nu$ indices for the lambda particles.
Finally, the first letters of the Greek alphabet
($\alpha,\beta,\gamma,\delta,\ldots$) used as indices will refer to the
Cartesian components of the operators.

##### Non strange finite and infinite systems

As already sketched in § 3.2, the AFDMC wave function is written in the single
particle state representation. The trial wave function for nuclear systems,
which is used both as projection and importance function
$|\psi_{I}\rangle\equiv|\psi_{T}\rangle$, is assumed of the form [31, 38]

$\displaystyle\psi_{T}^{N}(R_{N},S_{N})=\prod_{i<j}f_{c}^{NN}(r_{ij})\,\Phi_{N}(R_{N},S_{N})\;,$
(3.183)

where $R_{N}=\\{\bm{r}_{1},\ldots,\bm{r}_{\mathcal{N}_{N}}\\}$ are the
Cartesian coordinates and $S_{N}=\\{s_{1},\ldots,s_{\mathcal{N}_{N}}\\}$ the
spin-isospin coordinates, represented as complex 4- or 2-component vectors:

nucleons: $\displaystyle\quad s_{i}=\left(\begin{array}[]{c}a_{i}\\\ b_{i}\\\
c_{i}\\\
d_{i}\end{array}\right)_{i}\\!=a_{i}|p\uparrow\rangle_{i}+b_{i}|p\downarrow\rangle_{i}+c_{i}|n\uparrow\rangle_{i}+d_{i}|n\downarrow\rangle_{i}\;,$
(3.188) neutrons: $\displaystyle\quad s_{i}=\left(\begin{array}[]{c}a_{i}\\\
b_{i}\end{array}\right)_{i}\\!=a_{i}|n\uparrow\rangle_{i}+b_{i}|n\downarrow\rangle_{i}\;,$
(3.191)

with
$\left\\{|p\uparrow\rangle,|p\downarrow\rangle,|n\uparrow\rangle,|n\downarrow\rangle\right\\}$
the proton-neutron-up-down basis.

The function $f_{c}^{NN}(r)$ is a symmetric and spin independent Jastrow
correlation function, solution of the Schrödinger-like equation for
$f_{c}^{NN}(r<d)$

$\displaystyle-\frac{\hbar^{2}}{2\mu_{NN}}\nabla^{2}f_{c}^{NN}(r)+\eta\,v_{c}^{NN}(r)f_{c}^{NN}(r)=\xi
f_{c}^{NN}(r)\;,$ (3.192)

where $v_{c}^{NN}(r)$ is the spin independent part of the two-body $NN$
interaction, $\mu_{NN}=m_{N}/2$ the reduced mass of the nucleon pair and
$\eta$ and the healing distance $d$ are variational parameters. For distances
$r\geq d$ we impose $f_{c}^{NN}(r)=1$. The role of the Jastrow function is to
include the short-range correlations in the trial wave function. In the AFDMC
algorithm the effect is simply a reduction of the overlap between pairs of
particles, with the reduction of the energy variance. Since there is no change
in the phase of the wave function, the $f_{c}^{NN}$ function does not
influence the computed energy value in the long imaginary time projection.

The antisymmetric part $\Phi_{N}(R_{N},S_{N})$ of the trial wave function
depends on the system to be studied (finite or infinite). As already seen, it
is generally built starting from single particle states
$\varphi_{\epsilon}^{N}(\bm{r}_{i},s_{i})$, where $\epsilon$ is the set of
quantum numbers describing the state and $\bm{r}_{i}$, $s_{i}$ the single
particle nucleon coordinates. The antisymmetry property is then realized by
taking the Slater determinant of the $\varphi_{\epsilon}^{N}$:

$\displaystyle\Phi_{N}(R_{N},S_{N})=\mathcal{A}\Bigg{[}\prod_{i=1}^{\mathcal{N}_{N}}\varphi_{\epsilon}^{N}(\bm{r}_{i},s_{i})\Bigg{]}=\det\Bigl{\\{}\varphi_{\epsilon}^{N}(\bm{r}_{i},s_{i})\Bigr{\\}}\;.$
(3.193)

For nuclei and neutron drops [31] a good set of quantum number is given by
$\epsilon=\\{n,j,m_{j}\\}$. The single particle states are written as:

$\displaystyle\varphi_{\epsilon}^{N}(\bm{r}_{i},s_{i})=R_{n,j}^{N}(r_{i})\Bigl{[}Y_{l,m_{l}}^{N}(\Omega)\,\chi_{s,m_{s}}^{N}(s_{i})\Bigr{]}_{j,m_{j}}\;,$
(3.194)

where $R_{n,j}^{N}$ is a radial function, $Y_{l,m_{l}}^{N}$ the spherical
harmonics depending on the solid angle $\Omega$ and $\chi_{s,m_{s}}^{N}$ the
spinors in the proton-neutron-up-down basis. The angular functions are coupled
to the spinors using the Clebsh-Gordan coefficients to have orbitals in the
$\\{n,j,m_{j}\\}$ basis according to the usual shell model classification of
the nuclear single particle spectrum. For finite systems, in order to make the
wave function translationally invariant, the single particle orbitals have to
be defined with respect to the center of mass (CM) of the system. We have
thus:

$\displaystyle\varphi_{\epsilon}^{N}(\bm{r}_{i},s_{i})\longrightarrow\varphi_{\epsilon}^{N}(\bm{r}_{i}-\bm{r}_{CM},s_{i})\quad\quad\text{with}\quad\bm{r}_{CM}=\frac{1}{\mathcal{N}_{N}}\sum_{i=1}^{\mathcal{N}_{N}}\bm{r}_{i}\;.$
(3.195)

In order to deal with new shifted coordinates, we need to correct all the
first and second derivatives of trial wave function with respect to
$\bm{r}_{i}$. The derivation of such corrections is reported in Appendix A.
The choice of the radial functions $R_{n,j}^{N}$ depends on the system studied
and, typically, solutions of the self-consistent Hartree-Fock problem with
Skyrme interactions are adopted. For nuclei the Skyrme effective interactions
of Ref. [229] are commonly used. For neutron drops, the Skyrme SKM force of
Ref. [230] has been considered.

An additional aspect to take care when dealing with finite systems, is the
symmetry of the wave function. Because the AFDMC projects out the lowest
energy state not orthogonal to the starting trial wave function, it is
possible to study a state with given symmetry imposing to the trial wave
function the total angular momentum $J$ experimentally observed. This can be
achieved by taking a sum over a different set of determinants

$\displaystyle\det\Bigl{\\{}\varphi_{\epsilon}^{N}(\bm{r}_{i},s_{i})\Bigr{\\}}\longrightarrow\left[\sum_{\kappa}c_{\kappa}\,\text{det}_{\kappa}\Bigl{\\{}\varphi_{\epsilon_{\kappa}}^{N}(\bm{r}_{i},s_{i})\Bigr{\\}}\right]_{J,M_{J}}\;,$
(3.196)

where the $c_{\kappa}$ coefficients are determined in order to have the
eigenstate of total angular momentum $J=j_{1}+\ldots+j_{\mathcal{N}_{N}}$.

For nuclear and neutron matter [38], the antisymmetric part of the wave
function is given by the ground state of the Fermi gas, built from a set of
plane waves. The infinite uniform system at a given density is simulated with
$\mathcal{N}_{N}$ nucleons in a cubic box of volume $L^{3}$ replicated into
the space by means of periodic boundary conditions (PBC):

$\displaystyle\varphi_{\epsilon}^{N}(\bm{r}_{1}+L\hat{\bm{r}},\bm{r}_{2},\ldots,s_{i})=\varphi_{\epsilon}^{N}(\bm{r}_{1},\bm{r}_{2},\ldots,s_{i})\;.$
(3.197)

Working in a discrete space, the momentum vectors are quantized and can be
expressed as

$\displaystyle\bm{k}_{\epsilon}=\frac{2\pi}{L}\left(n_{x},n_{y},n_{z}\right)_{\epsilon}\;,$
(3.198)

where $\epsilon$ labels the quantum state and $n_{x}$, $n_{y}$ and $n_{z}$ are
integer numbers labelling the momentum shell. The single particle states are
then given by

$\displaystyle\varphi_{\epsilon}^{N}(\bm{r}_{i},s_{i})=\operatorname{e}^{-i\bm{k}_{\epsilon}\cdot\bm{r}_{i}}\chi_{s,m_{s}}^{N}(s_{i})\;.$
(3.199)

In order to meet the requirement of homogeneity and isotropy, the shell
structure of the system must be closed. The total number of Fermions in a
particular spin-isospin configuration that can be correctly simulated in a box
corresponds to the closure of one of the
$\left(n_{x},n_{y},n_{z}\right)_{\epsilon}$ shells. The list of the first
closure numbers is

$\displaystyle\mathcal{N}_{c}=1,7,19,27,33,57,81,93\ldots\;.$ (3.200)

Given a closure number $\mathcal{N}_{c}^{N}$, we can thus simulate an infinite
system by means of a periodic box with $2\,\mathcal{N}_{c}^{N}$ neutrons (up
and down spin) or $4\,\mathcal{N}_{c}^{N}$ nucleons (up and down spin and
isospin). Although the use of PBC should reduce the finite-size effects, in
general there are still sizable errors in the kinetic energy arising from
shell effects in filling the plane wave orbitals, even at the closed shell
filling in momentum space. However, in the thermodynamical limit
$\mathcal{N}_{c}^{N}\rightarrow\infty$, exact results should be obtained. For
symmetric nuclear matter (SNM), 28, 76 and 108 nucleons have been used [35],
resulting in comparable results for the energy per particle at a given
density. In the case of pure neutron matter (PNM), finite-size effects are
much more evident [38] and the thermodynamical limit is not reached
monotonically. Typically, PNM is simulated using 66 neutrons, which was found
to give the closest kinetic energy compared to the Fermi gas in the range of
$\mathcal{N}_{c}^{N}$ corresponding to feasible computational times.

As reported in Ref. [231], twist-averaged boundary conditions (TABC) can be
imposed on the trial wave function to reduce the finite-size effects. One can
allow particles to pick up a phase $\theta$ when they wrap around the periodic
boundaries:

$\displaystyle\varphi_{\epsilon}^{N}(\bm{r}_{1}+L\hat{\bm{r}},\bm{r}_{2},\ldots,s_{i})=\operatorname{e}^{i\theta}\varphi_{\epsilon}^{N}(\bm{r}_{1},\bm{r}_{2},\ldots,s_{i})\;.$
(3.201)

The boundary condition $\theta=0$ corresponds to the PBC, $\theta\neq 0$ to
the TABC. It has been shown that if the twist phase is integrated over, the
finite size effects are substantially reduced. TABC has been used in PNM
calculations [38], showing a small discrepancy in the energy per particle for
38, 45, 66 and 80 neutrons at fixed density. A remarkable result is that the
PNM energy for 66 neutrons using PBC is very close to the extrapolated result
obtained employing the TABC, validating then the standard AFDMC calculation
for 66 particles. Compare to PBC, employing the TABC results in a more
computational time and they have not been used in this work.

##### Strange finite and infinite systems

The $\Lambda$ hyperon, having isospin zero, does not participate to the
isospin doublet of nucleons. Referring to hypernuclear systems, we can
therefore treat the additional strange baryons as distinguishable particles
writing a trial wave function of the form

$\displaystyle\psi_{T}(R,S)=\prod_{\lambda i}f_{c}^{\Lambda N}(r_{\lambda
i})\,\psi_{T}^{N}(R_{N},S_{N})\,\psi_{T}^{\Lambda}(R_{\Lambda},S_{\Lambda})\;,$
(3.202)

where
$R=\\{\bm{r}_{1},\ldots,\bm{r}_{\mathcal{N}_{N}},\bm{r}_{1},\ldots,\bm{r}_{\mathcal{N}_{\Lambda}}\\}$
and
$S=\\{s_{1},\ldots,s_{\mathcal{N}_{N}},s_{1},\ldots,s_{\mathcal{N}_{\Lambda}}\\}$
refer to the coordinates of all the baryons and $\psi_{T}^{N}(R_{N},S_{N})$ is
the nucleon single particle wave function of Eq. (3.183).
$\psi_{T}^{\Lambda}(R_{\Lambda},S_{\Lambda})$ is the lambda single particle
wave function that takes the same structure of the nucleon one:

$\displaystyle\psi_{T}^{\Lambda}(R_{\Lambda},S_{\Lambda})=\prod_{\lambda<\mu}f_{c}^{\Lambda\Lambda}(r_{\lambda\mu})\,\Phi_{\Lambda}(R_{\Lambda},S_{\Lambda})\;,$
(3.203)

with

$\displaystyle\Phi_{\Lambda}(R_{\Lambda},S_{\Lambda})=\mathcal{A}\Bigg{[}\prod_{\lambda=1}^{\mathcal{N}_{\Lambda}}\varphi_{\epsilon}^{\Lambda}(\bm{r}_{\lambda},s_{\lambda})\Bigg{]}=\det\Bigl{\\{}\varphi_{\epsilon}^{\Lambda}(\bm{r}_{\lambda},s_{\lambda})\Bigr{\\}}\;.$
(3.204)

$R_{\Lambda}=\\{\bm{r}_{1},\ldots,\bm{r}_{\mathcal{N}_{\Lambda}}\\}$ are the
hyperon Cartesian coordinates and
$S_{\Lambda}=\\{s_{1},\ldots,s_{\mathcal{N}_{\Lambda}}\\}$ the hyperon spin
coordinates, represented by the 2-dimension spinor in the lambda-up-down
basis:

$\displaystyle s_{\lambda}=\left(\begin{array}[]{c}u_{\lambda}\\\
v_{\lambda}\end{array}\right)_{\lambda}\\!=u_{\lambda}|\Lambda\uparrow\rangle_{\lambda}+v_{\lambda}|\Lambda\downarrow\rangle_{\lambda}\;.$
(3.207)

The $\Lambda\Lambda$ Jastrow correlation function $f_{c}^{\Lambda\Lambda}(r)$
is calculated by means of Eq. (3.192) for the hyperon-hyperon pair using the
central channel of the $\Lambda\Lambda$ potential of Eq. (2.48). Eq. (3.192)
is also used to calculate the hyperon-nucleon correlation function
$f_{c}^{\Lambda N}(r)$ of the hypernuclear wave function (3.202) by
considering the pure central term of the $\Lambda N$ potential of Eq. (2.35)
and using the reduced mass

$\displaystyle\mu_{\Lambda N}=\frac{m_{\Lambda}\,m_{N}}{m_{\Lambda}+m_{N}}\;.$
(3.208)

For $\Lambda$ hypernuclei (and $\Lambda$ neutron drops) the hyperon single
particle states take the same structure as the nuclear case, and they read:

$\displaystyle\varphi_{\epsilon}^{\Lambda}(\bm{r}_{\lambda},s_{\lambda})=R_{n,j}^{\Lambda}(r_{\lambda})\Bigl{[}Y_{l,m_{l}}^{\Lambda}(\Omega)\,\chi_{s,m_{s}}^{\Lambda}(s_{\lambda})\Bigr{]}_{j,m_{j}}\;.$
(3.209)

Although the AFDMC code for hypernuclei is set up for an arbitrary number of
hyperons, we focused on single and double $\Lambda$ hypernuclei. Having just
two hyperons to deal with, only one radial function $R_{n,j}^{\Lambda}$ is
needed. Being the mass difference between the neutron and the $\Lambda$
particle small, we used the neutron $1s_{1/2}$ radial function also for the
hyperon.

Dealing with finite systems, the coordinates of all the baryons must be
related to the CM, that now is given by the coordinates of particles with
different mass. Nucleon and hyperon single particle orbitals are thus defined
as:

$\displaystyle\begin{array}[]{rcl}\varphi_{\epsilon}^{N}(\bm{r}_{i},s_{i})\\!\\!\\!&\longrightarrow&\\!\\!\varphi_{\epsilon}^{N}(\bm{r}_{i}-\bm{r}_{CM},s_{i})\\\\[5.0pt]
\varphi_{\epsilon}^{\Lambda}(\bm{r}_{\lambda},s_{\lambda})\\!\\!\\!&\longrightarrow&\\!\\!\varphi_{\epsilon}^{\Lambda}(\bm{r}_{\lambda}-\bm{r}_{CM},s_{\lambda})\end{array}$
(3.212)

where

$\displaystyle\bm{r}_{CM}=\frac{1}{M}\left(m_{N}\sum_{i=1}^{\mathcal{N}_{N}}\bm{r}_{i}+m_{\Lambda}\sum_{\lambda=1}^{\mathcal{N}_{\Lambda}}\bm{r}_{\lambda}\right)\quad\text{with}\quad
M=\mathcal{N}_{N}\,m_{N}+\mathcal{N}_{\Lambda}\,m_{\Lambda}\;.$ (3.213)

As in the nuclear case, the use of relative coordinates introduces corrections
in the calculation of the derivatives of the trial wave function. For
hypernuclei such corrections, and in general the evaluation of derivatives,
are more complicated than for nuclei. This is because we have to deal with two
set of spacial coordinates ($R_{N}$ and $R_{\Lambda}$) and the Jastrow
function $f_{c}^{\Lambda N}$ depends on both. The derivatives of the trial
wave function including CM corrections are reported in Appendix A.

For $\Lambda$ neutron (nuclear) matter the antisymmetric part of the hyperon
wave function is given by the ground state of the Fermi gas, as for nucleons.
We are thus dealing with two Slater determinants of plane waves with
$\bm{k}_{\epsilon}$ vectors quantized in the same $L^{3}$ box (see Eq.
(3.198)). The dimension of the simulation box, and thus the quantization of
the $\bm{k}_{\epsilon}$ vectors, is given by the total numeric baryon density

$\displaystyle\rho_{b}=\frac{\mathcal{N}_{b}}{L^{3}}=\frac{\mathcal{N}_{N}+\mathcal{N}_{\Lambda}}{L^{3}}=\rho_{N}+\rho_{\Lambda}\;,$
(3.214)

and the number of nucleons and lambda particles. The hyperon single particle
states correspond then to

$\displaystyle\varphi_{\epsilon}^{\Lambda}(\bm{r}_{\lambda},s_{\lambda})=\operatorname{e}^{-i\bm{k}_{\epsilon}\cdot\bm{r}_{\lambda}}\chi_{s,m_{s}}^{\Lambda}(s_{\lambda})\;,$
(3.215)

where the $\bm{k}_{\epsilon}$ structure derived from $\rho_{b}$ is used also
for the the nuclear part. The requirements of homogeneity and isotropy
discussed in the previous section still apply and so the lambda plane waves
have to close their own momentum shell structure. Given the list of closure
numbers (3.200), we can add $2\,\mathcal{N}_{c}^{\Lambda}$ hyperons (up and
down spin) to the $2\,\mathcal{N}_{c}^{N}$ neutrons or
$4\,\mathcal{N}_{c}^{N}$ nucleons in the periodic box.

The wave functions described so far are appropriate only if we consider
nucleons and hyperons as distinct particles. In this way, it is not possible
to include the $\Lambda N$ exchange term of Eq. (2.35) directly in the
propagator, because it mixes hyperon and nucleon states. The complete
treatment of this factor would require a drastic change in the AFDMC code
and/or a different kind of nuclear-hypernuclear interactions, as briefly
discussed in Appendix B. A perturbative analysis of the $v_{0}(r_{\lambda
i})\,\varepsilon\,\mathcal{P}_{x}$ term is however possible and it is reported
in the next section.

####  3.2.5 Propagator for hypernuclear systems

Consider a many-body system composed by nucleons and hyperons, interacting via
the full Hamiltonian (2.1) and described by the trial wave function (3.202).
Suppose to switch off all the spin-isospin interactions in all the channels
and keep only the central terms:

$\displaystyle
H=T_{N}+T_{\Lambda}+V_{NN}^{c}+V_{\Lambda\Lambda}^{c}+V_{\Lambda N}^{c}\;,$
(3.216)

where also the central contributions from the three-body interactions are
included. Neglecting the spin-isospin structure of the trial wave function we
can follow the idea of the standard DMC described in § 3.1 and write the
analog of Eq. (3.18):

$\displaystyle\psi(R_{N},R_{\Lambda},\tau+d\tau)\simeq\int
dR^{\prime}_{N}\,dR^{\prime}_{\Lambda}\langle
R_{N},R_{\Lambda}|\operatorname{e}^{-\left(V_{NN}^{c}+V_{\Lambda\Lambda}^{c}+V_{\Lambda
N}^{c}\right)\frac{d\tau}{2}}\operatorname{e}^{-T_{N}d\tau}\operatorname{e}^{-T_{\Lambda}d\tau}\times$
$\displaystyle\hskip
14.22636pt\times\operatorname{e}^{-\left(V_{NN}^{c}+V_{\Lambda\Lambda}^{c}+V_{\Lambda
N}^{c}\right)\frac{d\tau}{2}}\operatorname{e}^{E_{T}d\tau}|R^{\prime}_{N},R^{\prime}_{\Lambda}\rangle\,\psi(R^{\prime}_{N},R^{\prime}_{\Lambda},\tau)\;,$
$\displaystyle\simeq\int
dR^{\prime}_{N}\,dR^{\prime}_{\Lambda}\underbrace{\langle
R_{N}|\operatorname{e}^{-T_{N}d\tau}|R^{\prime}_{N}\rangle}_{G_{0}^{N}(R_{N},R^{\prime}_{N},d\tau)}\underbrace{\langle
R_{\Lambda}|\operatorname{e}^{-T_{\Lambda}d\tau}|R^{\prime}_{\Lambda}\rangle}_{G_{0}^{\Lambda}(R_{\Lambda},R^{\prime}_{\Lambda},d\tau)}\times$
$\displaystyle\hskip
14.22636pt\times\underbrace{\phantom{\langle}\\!\\!\operatorname{e}^{-\left(\widetilde{V}_{NN}^{c}(R_{N},R^{\prime}_{N})+\widetilde{V}_{\Lambda\Lambda}^{c}(R_{\Lambda},R^{\prime}_{\Lambda})+\widetilde{V}_{\Lambda
N}^{c}(R_{N},R_{\Lambda},R^{\prime}_{N},R^{\prime}_{\Lambda})-E_{T}\right)d\tau}}_{G_{V}(R_{N},R_{\Lambda},R^{\prime}_{N},R^{\prime}_{\Lambda},d\tau)}\psi(R^{\prime}_{N},R^{\prime}_{\Lambda},\tau)\;,$
$\displaystyle\simeq\left(\frac{1}{4\pi
D_{N}d\tau}\right)^{\frac{3\mathcal{N}_{N}}{2}}\\!\\!\left(\frac{1}{4\pi
D_{\Lambda}d\tau}\right)^{\frac{3\mathcal{N}_{\Lambda}}{2}}\\!\\!\int
dR^{\prime}_{N}\,dR^{\prime}_{\Lambda}\operatorname{e}^{-\frac{(R_{N}-R^{\prime}_{N})^{2}}{4D_{N}d\tau}}\operatorname{e}^{-\frac{(R_{\Lambda}-R^{\prime}_{\Lambda})^{2}}{4D_{\Lambda}d\tau}}\times$
$\displaystyle\hskip
14.22636pt\times\operatorname{e}^{-\left(\widetilde{V}_{NN}^{c}(R_{N},R^{\prime}_{N})+\widetilde{V}_{\Lambda\Lambda}^{c}(R_{\Lambda},R^{\prime}_{\Lambda})+\widetilde{V}_{\Lambda
N}^{c}(R_{N},R_{\Lambda},R^{\prime}_{N},R^{\prime}_{\Lambda})-E_{T}\right)d\tau}\psi(R^{\prime}_{N},R^{\prime}_{\Lambda},\tau)\;,$
(3.217)

where

$\displaystyle\widetilde{V}_{NN}^{c}(R_{N},R^{\prime}_{N})$
$\displaystyle=\frac{1}{2}\Bigl{[}V_{NN}^{c}(R_{N})+V_{NN}^{c}(R^{\prime}_{N})\Bigr{]}\;,$
$\displaystyle\widetilde{V}_{\Lambda\Lambda}^{c}(R_{\Lambda},R^{\prime}_{\Lambda})$
$\displaystyle=\frac{1}{2}\Bigl{[}V_{\Lambda\Lambda}^{c}(R_{\Lambda})+V_{\Lambda\Lambda}^{c}(R^{\prime}_{\Lambda})\Bigr{]}\;,$
(3.218) $\displaystyle\widetilde{V}_{\Lambda
N}^{c}(R_{N},R_{\Lambda},R^{\prime}_{N},R^{\prime}_{\Lambda})$
$\displaystyle=\frac{1}{2}\Bigl{[}V_{\Lambda
N}^{c}(R_{N},R_{\Lambda})+V_{\Lambda
N}^{c}(R^{\prime}_{N},R^{\prime}_{\Lambda})\Bigr{]}\;,$

and $D_{N}=\hbar^{2}/2m_{N}$ and $D_{\Lambda}=\hbar^{2}/2m_{\Lambda}$ are the
diffusion constants of the Brownian motion of nucleons and lambda particles.

The evolution in imaginary time is thus performed in the same way of the
standard DMC algorithm. A set of walkers, which now contains nucleon and
hyperon coordinates, is diffused according to
$G_{0}^{N}(R_{N},R^{\prime}_{N},d\tau)$ and
$G_{0}^{\Lambda}(R_{\Lambda},R^{\prime}_{\Lambda},d\tau)$. A weight
$\omega_{i}=G_{V}(R_{N},R_{\Lambda},R^{\prime}_{N},R^{\prime}_{\Lambda},d\tau)$
is assigned to each waker and it is used for the estimator contributions and
the branching process. We can also apply the importance function technique,
the result of which is the inclusion of a drift term in the diffusion of each
type of baryon and the use of the local energy for the branching weight. The
drift velocities take the same form of Eq. (3.32), but now the derivatives are
calculated with respect to nucleon or hyperon coordinates, including all the
possible CM (for finite systems) and Jastrow corrections, as reported in
Appendix A.

Reintroduce now the spin-isospin structure in the wave function and consider
then the spin-isospin dependent interactions. For the nuclear part, we can
still deal with AV6 like potentials for nucleon systems by means of the
Hubbard-Stratonovich transformation, as discussed in § 3.2.1. In the case of
pure neutron systems, we can also include spin-orbit and three-body
contributions as reported in § 3.2.2 and § 3.2.3. In the next we will discuss
how to deal with the spin-isospin dependent part of the hypernuclear
potentials, in both two- and three-body channels.

##### Two-body terms

Consider the full two-body $\Lambda N$ interaction described in the previous
chapter:

$\displaystyle V_{\Lambda N}$ $\displaystyle=\sum_{\lambda i}\left(v_{\lambda
i}+v_{\lambda i}^{CSB}\right)\;,$ $\displaystyle=\sum_{\lambda
i}v_{0}(r_{\lambda i})(1-\varepsilon)+\sum_{\lambda i}v_{0}(r_{\lambda
i})\,\varepsilon\,\mathcal{P}_{x}+\sum_{\lambda
i}\frac{1}{4}v_{\sigma}T^{2}_{\pi}(r_{\lambda
i})\,\bm{\sigma}_{\lambda}\cdot\bm{\sigma}_{i}$
$\displaystyle\quad+\sum_{\lambda i}C_{\tau}\,T_{\pi}^{2}\left(r_{\lambda
i}\right)\tau_{i}^{z}\;,$ $\displaystyle=\sum_{\lambda i}v_{0}(r_{\lambda
i})(1-\varepsilon)+\sum_{\lambda i}B_{\lambda
i}^{[\mathcal{P}_{x}]}\,\mathcal{P}_{x}+\sum_{\lambda
i}\sum_{\alpha}\sigma_{\lambda\alpha}\,B_{\lambda
i}^{[\sigma]}\,\sigma_{i\alpha}+\sum_{i}B_{i}^{[\tau]}\,\tau_{i}^{z}\;,$
(3.219)

where

$\displaystyle B_{\lambda i}^{[\mathcal{P}_{x}]}=v_{0}(r_{\lambda
i})\,\varepsilon\quad\quad B_{\lambda
i}^{[\sigma]}=\frac{1}{4}v_{\sigma}T^{2}_{\pi}(r_{\lambda i})\quad\quad
B_{i}^{[\tau]}=\sum_{\lambda}C_{\tau}\,T_{\pi}^{2}\left(r_{\lambda
i}\right)\;.$ (3.220)

The first term of Eq. (3.219) is simply a spin independent factor and can be
included in the $V_{\Lambda N}^{c}$ contribution of Eq. (3.217). The remaining
terms involve operators acting on nucleons and hyperons and need to be
discussed separately.

  * •

$\bm{\sigma}_{\lambda}\cdot\bm{\sigma}_{i}$  _term_. The quadratic spin-spin
term of the $\Lambda N$ interaction is written in same form of the nucleon-
nucleon one of Eq. (3.126). However, in general the matrix $B_{\lambda
i}^{[\sigma]}$ is not a square matrix ($\dim B_{\lambda
i}^{[\sigma]}=\mathcal{N}_{\Lambda}\times\mathcal{N}_{N}$) and so we can not
follow the derivation of § 3.2.1. Recalling that we are working with single
particle wave functions and that each spin-isospin operator is the
representation in the $A$-body tensor product space of a one-body operator as
in Eq. (3.123), we can write

$\displaystyle\\!\\!\\!\\!\sum_{\alpha}\sigma_{\lambda\alpha}\otimes\sigma_{i\alpha}=\frac{1}{2}\sum_{\alpha}\left[\left(\sigma_{\lambda\alpha}\oplus\sigma_{i\alpha}\right)^{2}-\left(\sigma_{\lambda\alpha}\otimes\mathbb{I}_{i\alpha}\right)^{2}-\left(\mathbb{I}_{\lambda\alpha}\otimes\sigma_{i\alpha}\right)^{2}\right]\;.\\!$
(3.221)

The square of the Pauli matrices of the last two terms gives the identity with
respect to the single particle state $\lambda$ or $i$, so that they can be
simply written as a spin independent contribution

$\displaystyle\sum_{\alpha}\sigma_{\lambda\alpha}\otimes\sigma_{i\alpha}=-3+\frac{1}{2}\sum_{\alpha}\left(\mathcal{O}_{\lambda
i,\alpha}^{[\sigma_{\Lambda N}]}\right)^{2}\;,$ (3.222)

where we have defined a new spin-spin operator

$\displaystyle\mathcal{O}_{\lambda i,\alpha}^{[\sigma_{\Lambda
N}]}=\sigma_{\lambda\alpha}\oplus\sigma_{i\alpha}\;.$ (3.223)

We can now write the the $\bm{\sigma}_{\lambda}\cdot\bm{\sigma}_{i}$ term as
follows

$\displaystyle V_{\Lambda N}^{\sigma\sigma}$ $\displaystyle=\sum_{\lambda
i}\sum_{\alpha}\sigma_{\lambda\alpha}\,B_{\lambda
i}^{[\sigma]}\,\sigma_{i\alpha}\;,$ $\displaystyle=-3\sum_{\lambda
i}B_{\lambda i}^{[\sigma]}+\frac{1}{2}\sum_{\lambda i}\sum_{\alpha}B_{\lambda
i}^{[\sigma]}\left(\mathcal{O}_{\lambda i,\alpha}^{[\sigma_{\Lambda
N}]}\right)^{2}$ (3.224)

The first term is a central factor that can be included in $V_{\Lambda
N}^{c}$. The second term is written in the same way of the spin-isospin
dependent part of the nuclear interaction of Eq. (3.135). With a little abuse
of notation
$n=\\{\lambda,i\\}={1,\ldots,\mathcal{N}_{N}\,\mathcal{N}_{\Lambda}}$, the
spin dependent part of the propagator for the
$\bm{\sigma}_{\lambda}\cdot\bm{\sigma}_{i}$ takes a suitable form for the
application of the Hubbard-Stratonovich transformation:

$\displaystyle\operatorname{e}^{-\sum_{\lambda
i}\sum_{\alpha}\sigma_{\lambda\alpha}\,B_{\lambda
i}^{[\sigma]}\,\sigma_{i\alpha}\,d\tau}$
$\displaystyle=\operatorname{e}^{-3\sum_{n}B_{n}^{[\sigma]}-\frac{1}{2}\sum_{n\alpha}B_{n}^{[\sigma]}\left(\mathcal{O}_{n,\alpha}^{[\sigma_{\Lambda
N}]}\right)^{2}d\tau}\;,$ $\displaystyle=\operatorname{e}^{-V_{\Lambda
N}^{c\,[\sigma]}}\prod_{n\alpha}\operatorname{e}^{-\frac{1}{2}B_{n}^{[\sigma]}\left(\mathcal{O}_{n\alpha}^{[\sigma_{\Lambda
N}]}\right)^{2}d\tau}\,+\operatorname{o}\\!\left(d\tau^{2}\right)\;,$
$\displaystyle\simeq\operatorname{e}^{-V_{\Lambda
N}^{c\,[\sigma]}}\prod_{n\alpha}\frac{1}{\sqrt{2\pi}}\int\\!dx_{n\alpha}\operatorname{e}^{\frac{-x_{n\alpha}^{2}}{2}+\sqrt{-B_{n}^{[\sigma]}d\tau}\,x_{n\alpha}\mathcal{O}_{n\alpha}^{[\sigma_{\Lambda
N}]}}\;.$ (3.225)

Recalling Eq. (3.138), we can write the hyperon spin dependent part of the
AFDMC propagator for hypernuclear systems as

$\displaystyle\langle
S_{N}S_{\Lambda}|\prod_{n\alpha=1}^{3\mathcal{N}_{N}\mathcal{N}_{\Lambda}}\\!\frac{1}{\sqrt{2\pi}}\int\\!dx_{n\alpha}\operatorname{e}^{\frac{-x_{n\alpha}^{2}}{2}+\sqrt{-B_{n}^{[\sigma]}d\tau}\,x_{n\alpha}\mathcal{O}_{n\alpha}^{[\sigma_{\Lambda
N}]}}|S^{\prime}_{N}S^{\prime}_{\Lambda}\rangle\;.$ (3.226)

By the definition of Eq. (3.223), it comes out that the action of the operator
$\mathcal{O}_{n\alpha}^{[\sigma_{\Lambda N}]}$ on the spinor
$|S^{\prime}_{N},S^{\prime}_{\Lambda}\rangle$ factorizes in a
$\sigma_{i\alpha}$ rotation for the nucleon spinor $|S_{N}\rangle$ and a
$\sigma_{\lambda\alpha}$ rotation for the $\Lambda$ spinor
$|S_{\Lambda}\rangle$, coupled by the same coefficient
$\sqrt{-B_{n}^{[\sigma]}d\tau}\,x_{n\alpha}$.

  * •

$\tau_{i}^{z}$  _term_. As already seen in § 3.2, the single particle wave
function is closed with respect to the application of a propagator containing
linear spin-isospin operators. The action of the CSB potential corresponds to
the propagator

$\displaystyle\operatorname{e}^{-\sum_{i}B_{i}^{[\tau]}\,\tau_{i}^{z}\,d\tau}=\prod_{i}\operatorname{e}^{-B_{i}^{[\tau]}\,\tau_{i}^{z}\,d\tau}\,+\operatorname{o}\left(d\tau^{2}\right)\;,$
(3.227)

that, acting on the trial wave function, simply produces a rotation of the
nucleon spinors, as in Eq. (3.148). Being the CSB term already linear in
$\tau_{i}^{z}$, there is no need for Hubbard-Stratonovich transformation. The
$\tau_{i}^{z}$ rotations can be applied after the integration over auxiliary
fields on the spinors modified by the Hubbard-Stratonovich rotations. In the
$\psi_{I}$ ratio AFDMC algorithm (_v1_) there are no additional terms in the
weight coming from the CSB rotations. If we use the local energy version of
the algorithm (_v2_), we need to subtract the CSB contribution to $E_{L}(R)$
(Eq. (3.34)) because there are no counter terms coming from the importance
sampling on auxiliary fields. Note that in neutron systems, $\tau_{i}^{z}$
gives simply a factor $-1$, so the CSB becomes a positive central contribution
($C_{\tau}<0$) to be added in $V_{\Lambda N}^{c}$.

  * •

$\mathcal{P}_{x}$  _term_. As discussed in § 3.2.4, the structure of our trial
wave function for hypernuclear systems prevents the straightforward inclusion
of the $\Lambda N$ space exchange operator in the AFDMC propagator. We can try
to treat this contribution perturbatively: $\mathcal{P}_{x}$ is not included
in the propagator but its effect is calculated as a standard estimator on the
propagated wave function. The action of $\mathcal{P}_{x}$ is to exchange the
position of one nucleon and one hyperon, modifying thus the CM of the whole
system due to the mass difference between the baryons. To compute this
potential contribution we have thus to sum over all the hyperon-nucleon pairs.
For each exchanged pair, all the positions are referred to the new CM and the
wave function is evaluated and accumulated. Then, particles are moved back to
the original positions and a new pair is processed. At the end of the sum the
contribution $\sum_{\lambda i}\mathcal{P}_{x}\,\psi$ is obtained. As reported
in Refs. [191, 188, 195, 189], the $\Lambda N$ space exchange operator induces
strong correlations and thus a perturbative approach might not be appropriate.
A possible non perturbative extension of the AFDMC code for the space exchange
operator is outlined in Appendix B.

In the two body hypernuclear sector a $\Lambda\Lambda$ interaction is also
employed, as reported in § 2.2.3. The potential described in Eq. (2.48) can be
recast as

$\displaystyle V_{\Lambda\Lambda}$
$\displaystyle=\sum_{\lambda<\mu}\sum_{k=1}^{3}\left(v_{0}^{(k)}+v_{\sigma}^{(k)}\,{\bm{\sigma}}_{\lambda}\cdot{\bm{\sigma}}_{\mu}\right)\operatorname{e}^{-\mu^{(k)}r_{\lambda\mu}^{2}}\;,$
$\displaystyle=\sum_{\lambda<\mu}\sum_{k=1}^{3}v_{0}^{(k)}\operatorname{e}^{-\mu^{(k)}r_{\lambda\mu}^{2}}+\frac{1}{2}\sum_{\lambda\neq\mu}\sum_{\alpha}\sigma_{\lambda\alpha}\,C_{\lambda\mu}^{[\sigma]}\,\sigma_{\mu\alpha}\;,$
(3.228)

where

$\displaystyle
C_{\lambda\mu}^{[\sigma]}=\sum_{k=1}^{3}v_{\sigma}^{(k)}\operatorname{e}^{-\mu^{(k)}r_{\lambda\mu}^{2}}\;.$
(3.229)

The first term of $V_{\Lambda\Lambda}$ is a pure central factor to be included
in $V_{\Lambda\Lambda}^{c}$, while the second part has exactly the same form
of the isospin component of Eq. (3.126). We can thus diagonalize the $C$
matrix and define a new operator $\mathcal{O}_{n,\alpha}^{[\sigma_{\Lambda}]}$
starting from the eigenvectors $\psi_{n,\lambda}^{[\sigma_{\Lambda}]}$:

$\displaystyle\mathcal{O}_{n,\alpha}^{[\sigma_{\Lambda}]}$
$\displaystyle=\sum_{\lambda}\sigma_{\lambda\alpha}\,\psi_{n,\lambda}^{[\sigma_{\Lambda}]}\;.$
(3.230)

In this way, the spin dependent part of the hyperon-hyperon interaction
becomes

$\displaystyle
V_{\Lambda\Lambda}^{SD}=\frac{1}{2}\sum_{n=1}^{\mathcal{N}_{\Lambda}}\sum_{\alpha=1}^{3}\lambda_{n}^{[\sigma_{\Lambda}]}\\!\left(\mathcal{O}_{n\alpha}^{[\sigma_{\Lambda}]}\right)^{2}\;,$
(3.231)

and we can apply the Hubbard-Stratonovich transformation to linearize the
square dependence of $\mathcal{O}_{n\alpha}^{[\sigma_{\Lambda}]}$ introducing
the integration over $3\,\mathcal{N}_{\Lambda}$ new auxiliary fields and the
relative $|S^{\prime}_{\Lambda}\rangle$ rotations.

At the end, using the diagonalization of the potential matrices and the
derivation reported in this section, the spin-isospin dependent part of the
nuclear and hypernuclear two-body potentials (but spin-orbit term for
simplicity) reads:

$\displaystyle V_{NN}^{SD}+V_{\Lambda N}^{SD}$
$\displaystyle=\frac{1}{2}\sum_{n=1}^{3\mathcal{N}_{N}}\lambda_{n}^{[\sigma]}\\!\left(\mathcal{O}_{n}^{[\sigma]}\right)^{2}\\!+\frac{1}{2}\sum_{n=1}^{3\mathcal{N}_{N}}\sum_{\alpha=1}^{3}\lambda_{n}^{[\sigma\tau]}\\!\left(\mathcal{O}_{n\alpha}^{[\sigma\tau]}\right)^{2}\\!+\frac{1}{2}\sum_{n=1}^{\mathcal{N}_{N}}\sum_{\alpha=1}^{3}\lambda_{n}^{[\tau]}\\!\left(\mathcal{O}_{n\alpha}^{[\tau]}\right)^{2}$
$\displaystyle\,+\frac{1}{2}\sum_{n=1}^{\mathcal{N}_{\Lambda}}\sum_{\alpha=1}^{3}\lambda_{n}^{[\sigma_{\Lambda}]}\\!\left(\mathcal{O}_{n\alpha}^{[\sigma_{\Lambda}]}\right)^{2}\\!+\frac{1}{2}\sum_{n=1}^{\mathcal{N}_{N}\mathcal{N}_{\Lambda}}\\!\sum_{\alpha=1}^{3}B_{n}^{[\sigma]}\\!\left(\mathcal{O}_{n\alpha}^{[\sigma_{\Lambda
N}]}\right)^{2}\\!+\sum_{i=1}^{\mathcal{N}_{N}}B_{i}^{[\tau]}\,\tau_{i}^{z}\;.$
(3.232)

Using a compact notation, the AFDMC propagator for hypernuclear systems of Eq.
(3.217) with the inclusion of spin-isospin degrees of freedom becomes:

$\displaystyle G(R,S,R^{\prime},S^{\prime},d\tau)=\langle
R,S|\operatorname{e}^{-(T_{N}+T_{\Lambda}+V_{NN}+V_{\Lambda\Lambda}+V_{\Lambda
N}-E_{T})d\tau}|R^{\prime},S^{\prime}\rangle\;,$ $\displaystyle\hskip
3.99994pt\simeq\left(\frac{1}{4\pi
D_{N}d\tau}\right)^{\frac{3\mathcal{N}_{N}}{2}}\\!\\!\left(\frac{1}{4\pi
D_{\Lambda}d\tau}\right)^{\frac{3\mathcal{N}_{\Lambda}}{2}}\\!\\!\operatorname{e}^{-\frac{(R_{N}-R^{\prime}_{N})^{2}}{4D_{N}d\tau}}\operatorname{e}^{-\frac{(R_{\Lambda}-R^{\prime}_{\Lambda})^{2}}{4D_{\Lambda}d\tau}}\times$
$\displaystyle\hskip
17.07182pt\times\operatorname{e}^{-\left(\widetilde{V}_{NN}^{c}(R_{N},R^{\prime}_{N})+\widetilde{V}_{\Lambda\Lambda}^{c}(R_{\Lambda},R^{\prime}_{\Lambda})+\widetilde{V}_{\Lambda
N}^{c}(R_{N},R_{\Lambda},R^{\prime}_{N},R^{\prime}_{\Lambda})-E_{T}\right)d\tau}$
$\displaystyle\hskip 17.07182pt\times\langle
S_{N},S_{\Lambda}|\prod_{i=1}^{\mathcal{N}_{N}}\operatorname{e}^{-B_{i}^{[\tau]}\,\tau_{i}^{z}\,d\tau}\prod_{n=1}^{\mathcal{M}}\frac{1}{\sqrt{2\pi}}\int\\!dx_{n}\operatorname{e}^{\frac{-x_{n}^{2}}{2}+\sqrt{-\lambda_{n}d\tau}\,x_{n}\mathcal{O}_{n}}|S^{\prime}_{N},S^{\prime}_{\Lambda}\rangle\;,$
(3.233)

where $|R,S\rangle\equiv|R_{N},R_{\Lambda},S_{N},S_{\Lambda}\rangle$ is the
state containing all the coordinates of the baryons and
$\widetilde{V}_{NN}^{c}$, $\widetilde{V}_{\Lambda\Lambda}^{c}$ and
$\widetilde{V}_{\Lambda N}^{c}$ defined in Eqs. (3.218) contain all the
possible central factors. Formally,
$\mathcal{M}=15\,\mathcal{N}_{N}+3\,\mathcal{N}_{\Lambda}+3\,\mathcal{N}_{N}\mathcal{N}_{\Lambda}$
and $\mathcal{O}_{n}$ stays for the various operators of Eq. (3.232), which
have a different action on the spinors
$|S^{\prime}_{N},S^{\prime}_{\Lambda}\rangle$. The
$\mathcal{O}_{n}^{[\sigma]}$, $\mathcal{O}_{n\alpha}^{[\sigma\tau]}$ and
$\mathcal{O}_{n\alpha}^{[\tau]}$ act on the nucleon spinor
$|S^{\prime}_{N}\rangle$. The $\mathcal{O}_{n\alpha}^{[\sigma_{\Lambda}]}$
rotates the lambda spinor $|S^{\prime}_{\Lambda}\rangle$.
$\mathcal{O}_{n\alpha}^{[\sigma_{\Lambda N}]}$ acts on both baryon spinors
with a separate rotation for nucleons and hyperons coupled by the same
coefficient $(-B_{n}^{[\sigma]}d\tau)^{1/2}\,x_{n}$ (recall Eq. (3.223)). The
algorithm follows then the nuclear version (§ 3.2.1) with the sampling of the
nucleon and hyperon coordinates and of the auxiliary fields, one for each
linearized operator. The application of the propagator of Eq. (3.233) has the
effect to rotate the spinors of the baryons. The weight for each walker is
then calculate starting from the central part of the interaction with possible
counter terms coming from the importance sampling on spacial coordinates and
on auxiliary fields (algorithm  _v1_), or by means of the local energy
(algorithm  _v2_). Fixed phase approximation, branching process and
expectation values are the same discussed in § 3.1.

##### Three-body terms

We have already shown in § 3.2.3 that for neutron systems the three-body
nucleon force can be recast as a sum of two-body terms only. In the case of
the three-body $\Lambda NN$ interaction it is possible to verify that the same
reduction applies both for nucleon and neutron systems. Let consider the full
potential of Eqs. (2.43) and (2.47)

$\displaystyle V_{\Lambda NN}=\sum_{\lambda,i<j}\left(v_{\lambda
ij}^{2\pi,P}+v_{\lambda ij}^{2\pi,S}+v_{\lambda ij}^{D}\right)\;,$ (3.234)

and assume the notations:

$\displaystyle T_{\lambda i}=T_{\pi}(r_{\lambda i})\quad\quad Y_{\lambda
i}=Y_{\pi}(r_{\lambda i})\quad\quad Q_{\lambda i}=Y_{\lambda i}-T_{\lambda
i}\;.$ (3.235)

By expanding the operators over the Cartesian components as done in Eqs.
(3.128) and (3.129), it is possible to derive the following relations:

$\displaystyle V_{\Lambda NN}^{2\pi,S}$ $\displaystyle=\frac{1}{2}\sum_{i\neq
j}\sum_{\alpha\beta\gamma}\tau_{i\gamma}\,\sigma_{i\alpha}\left(-\frac{C_{P}}{3}\sum_{\lambda}\sum_{\delta}\Theta_{\lambda
i}^{\alpha\delta}\,\Theta_{\lambda
j}^{\beta\delta}\right)\sigma_{j\beta}\,\tau_{j\gamma}\;,$ (3.236)
$\displaystyle V_{\Lambda NN}^{2\pi,P}$ $\displaystyle=\frac{1}{2}\sum_{i\neq
j}\sum_{\alpha\beta\gamma}\tau_{i\gamma}\,\sigma_{i\alpha}\,\Xi_{i\alpha,j\beta}\,\sigma_{j\beta}\,\tau_{j\gamma}\;,$
(3.237) $\displaystyle V_{\Lambda NN}^{D}$
$\displaystyle=W_{D}\\!\sum_{\lambda,i<j}T^{2}_{\lambda i}\,T^{2}_{\lambda
j}+\frac{1}{2}\sum_{\lambda i}\sum_{\alpha}\sigma_{\lambda\alpha}\,D_{\lambda
i}^{[\sigma]}\,\sigma_{i\alpha}\;,$ (3.238)

where

$\displaystyle\Theta_{\lambda i}^{\alpha\beta}$ $\displaystyle=Q_{\lambda
i}\,\delta^{\alpha\beta}+3\,T_{\lambda i}\hat{r}_{\lambda
i}^{\alpha}\,\hat{r}_{\lambda i}^{\beta}\;,$ (3.239)
$\displaystyle\Xi_{i\alpha,j\beta}$
$\displaystyle=\frac{1}{9}C_{S}\,\mu_{\pi}^{2}\sum_{\lambda}\,Q_{i\lambda}\,Q_{\lambda
j}\,|r_{i\lambda}||r_{j\lambda}|\,\hat{r}_{i\lambda}^{\alpha}\,\hat{r}_{j\lambda}^{\beta}\;,$
(3.240) $\displaystyle D_{\lambda i}^{[\sigma]}$
$\displaystyle=\frac{1}{3}W_{D}\\!\sum_{j,j\neq i}T^{2}_{\lambda
i}\,T^{2}_{\lambda j}\;.$ (3.241)

By combining then Eq. (3.239) and Eq. (3.240), the TPE term of the three-body
hyperon-nucleon interaction can be recast as:

$\displaystyle V_{\Lambda NN}^{2\pi}=\frac{1}{2}\sum_{i\neq
j}\sum_{\alpha\beta\gamma}\tau_{i\gamma}\,\sigma_{i\alpha}\,D_{i\alpha,j\beta}^{[\sigma\tau]}\,\sigma_{j\beta}\,\tau_{j\gamma}\;,$
(3.242)

where

$\displaystyle D_{i\alpha,j\beta}^{[\sigma\tau]}$
$\displaystyle=\sum_{\lambda}\Bigg{\\{}\\!-\frac{1}{3}C_{P}Q_{\lambda
i}Q_{\lambda j}\delta_{\alpha\beta}-C_{P}Q_{\lambda i}T_{\lambda
j}\,\hat{r}_{j\lambda}^{\,\alpha}\,\hat{r}_{j\lambda}^{\,\beta}-C_{P}Q_{\lambda
j}T_{\lambda i}\,\hat{r}_{i\lambda}^{\,\alpha}\,\hat{r}_{i\lambda}^{\,\beta}$
$\displaystyle\quad+\left[-3\,C_{P}T_{\lambda i}T_{\lambda
j}\left({\sum_{\delta}}\,\hat{r}_{i\lambda}^{\,\delta}\,\hat{r}_{j\lambda}^{\,\delta}\right)+\frac{1}{9}C_{S}\mu_{\pi}^{2}Q_{\lambda
i}Q_{\lambda
j}\,|r_{i\lambda}||r_{j\lambda}|\right]\,\hat{r}_{i\lambda}^{\,\alpha}\,\hat{r}_{j\lambda}^{\,\beta}\Bigg{\\}}\;.$
(3.243)

Finally, the $\Lambda NN$ interaction takes the following form:

$\displaystyle V_{\Lambda NN}$
$\displaystyle=W_{D}\\!\sum_{\lambda,i<j}T^{2}_{\lambda i}\,T^{2}_{\lambda
j}+\frac{1}{2}\sum_{\lambda i}\sum_{\alpha}\sigma_{\lambda\alpha}\,D_{\lambda
i}^{[\sigma]}\,\sigma_{i\alpha}$ $\displaystyle\quad\,+\frac{1}{2}\sum_{i\neq
j}\sum_{\alpha\beta\gamma}\tau_{i\gamma}\,\sigma_{i\alpha}\,D_{i\alpha,j\beta}^{[\sigma\tau]}\,\sigma_{j\beta}\,\tau_{j\gamma}\;.$
(3.244)

The first term is a pure central factor that can be included in $V_{\Lambda
N}^{c}$. The second factor is analogous to the
$\bm{\sigma}_{\lambda}\cdot\bm{\sigma}_{i}$ term (3.219) of the two body
hyperon-nucleon interaction. The last term acts only on the spin-isospin of
the two nucleons $i$ and $j$ and has the same structure of the nuclear
$\bm{\sigma}\cdot\bm{\tau}$ contribution described by the matrix
$A_{i\alpha,j\beta}^{[\sigma\tau]}$. The three-body hyperon-nucleon
interaction is then written as a sum of two-body operators only, of the same
form of the ones already discussed for the $NN$ and $\Lambda N$ potentials. We
can therefore include also these contributions in the AFDMC propagator of Eq.
(3.233) by simply redefining the following matrices:

$\displaystyle B_{\lambda i}^{[\sigma]}$ $\displaystyle\longrightarrow
B_{\lambda i}^{[\sigma]}+D_{\lambda i}^{[\sigma]}\;,$ (3.245) $\displaystyle
A_{i\alpha,j\beta}^{[\sigma\tau]}$ $\displaystyle\longrightarrow
A_{i\alpha,j\beta}^{[\sigma\tau]}+D_{i\alpha,j\beta}^{[\sigma\tau]}\;.$
(3.246)

The algorithm follows then the steps already discussed in the previous
section. Note that in the case of pure neutron systems, the last term of Eq.
(3.244) simply reduces to a $\bm{\sigma}_{i}\cdot\bm{\sigma}_{j}$ contribution
that is included in the propagator by redefining the nuclear matrix
$A_{i\alpha,j\beta}^{[\sigma]}$.

With the AFDMC method extended to the hypernuclear sector, we can study finite
and infinite lambda-nucleon and lambda-neutron systems. In the first case we
can treat Hamiltonians that include the full hyperon-nucleon, hyperon-nucleon-
nucleon and hyperon-hyperon interaction of Chapter 2, but we are limited to
the Argonne V6 like potentials for the nuclear sector. However it has been
shown that this approach gives good results for finite nuclei [31] and nuclear
matter [35, 36]. In the latter case, instead, we can also add the nucleon
spin-orbit contribution, so AV8 like potentials, and the three-neutron force.
The neutron version of the AFDMC code has been extensively and successfully
applied to study the energy differences of oxygen [30] and calcium [32]
isotopes, the properties of neutron drops [33, 20, 34] and the properties of
neutron matter in connection with astrophysical observables [37, 38, 39, 40].
Very recently, the AFDMC algorithm has been also used to perform calculations
for neutron matter using chiral effective field theory interactions [132].

## Chapter 4 Results: finite systems

This chapter reports on the analysis of finite systems, nuclei and
hypernuclei. For single $\Lambda$ hypernuclei a direct comparison of energy
calculations with experimental results is given for the $\Lambda$ separation
energy, defined as:

$\displaystyle
B_{\Lambda}\left(\,{}^{A}_{\Lambda}\text{Z}\,\right)=E\left(\,{}^{A-1}\text{Z}\,\right)-E\left(\,{}^{A}_{\Lambda}\text{Z}\,\right)\;,$
(4.1)

where, using the notation of the previous chapters, ${}^{A}_{\Lambda}\text{Z}$
refers to the hypernucleus and ${}^{A-1}\text{Z}$ to the corresponding
nucleus. $E$ is the binding energy of the system, i.e. the expectation value
of the Hamiltonian on the ground state wave function

$\displaystyle
E(\kappa)=\frac{\langle\psi_{0,\kappa}|H_{\kappa}|\psi_{0,\kappa}\rangle}{\langle\psi_{0,\kappa}|\psi_{0,\kappa}\rangle}\;,\quad\quad\kappa=\text{nuc},\text{hyp}\;,$
(4.2)

that we can compute by means of the AFDMC method. In the case of double
$\Lambda$ hypernuclei, the interesting experimental observables we can have
access are the double $\Lambda$ separation energy

$\displaystyle B_{\Lambda\Lambda}\left(\,{}^{\leavevmode\nobreak\
\,A}_{\Lambda\Lambda}\text{Z}\,\right)=E\left(\,{}^{A-2}\text{Z}\,\right)-E\left(\,{}^{\leavevmode\nobreak\
\,A}_{\Lambda\Lambda}\text{Z}\,\right)\;,$ (4.3)

and the incremental $\Lambda\Lambda$ energy

$\displaystyle\Delta B_{\Lambda\Lambda}\left(\,{}^{\leavevmode\nobreak\
\,A}_{\Lambda\Lambda}\text{Z}\,\right)=B_{\Lambda\Lambda}\left(\,{}^{\leavevmode\nobreak\
\,A}_{\Lambda\Lambda}\text{Z}\,\right)-2B_{\Lambda}\left(\,{}^{A-1}_{\quad\;\Lambda}\text{Z}\,\right)\;.$
(4.4)

The calculation of these quantities proceeds thus with the computation of the
binding energies for both strange and non strange systems. Moreover it is
interesting to compare other observables among the systems with strangeness
$0$, $-1$ and $-2$, such as the single particle densities. By looking at the
densities in the original nucleus and in the one modified by the addition of
the lambda particles, information about the hyperon-nucleon interaction can be
deduced.

As reported in Ref. [31], the ground state energies of 4He, 8He, 16O and 40Ca
have been computed using the AV6’ potential (§ 2.1.1). Due to the limitations
in the potential used, the results cannot reproduce the experimental energies
and all the nuclei result less bound than expected. However, given the same
simplified interaction, the published AFDMC energies are close to the GFMC
results, where available. AFDMC has also been used to compute the energy
differences between oxygen [30] and calcium [32] isotopes, by studying the
external neutrons with respect to a nuclear core obtained from Hartree-Fock
calculations using Skyrme forces. In this case the results are close to the
experimental ones.

The idea behind the AFDMC analysis of $\Lambda$ hypernuclei follows in some
sense the one assumed in the study of oxygen and calcium isotopes by the
analysis of energy differences. The two-body nucleon interaction employed is
limited to the first six operators of AV18. However, if we use the same
potential for the nucleus and the core of the corresponding hypernucleus, and
take the difference between the binding energies of the two systems, the
uncertainties in the $NN$ interaction largely cancel out. We shall see that
this assumption, already used in other works [148, 181], is indeed consistent
with our results, thereby confirming that the specific choice of the nucleon
Hamiltonian does not significantly affect the results on $B_{\Lambda}$. On the
grounds of this observation, we can focus on the interaction between hyperons
and nucleons, performing QMC simulations with microscopic interactions in a
wide mass range.

###  4.1 Nuclei

Let us start with the AFDMC study of finite nuclei. In the previous chapter,
we have seen that two versions of the AFDMC algorithm, that should give the
same results, are available (_v1_ and _v2_). Before including the strange
degrees of freedom, we decided to test the stability and accuracy of the two
algorithms, within the fixed phase approximation, by performing some test
simulations on 4He. The result of $-27.13(10)$ MeV for the AV6’ potential
reported in Ref. [31], was obtained employing the algorithm  _v2_ using single
particle Skyrme orbitals and a particular choice of the parameters for the
solution of the Jastrow correlation equation (see § 3.2.4). In order to check
the AFDMC projection process, we tried to modify the starting trial wave
function:

  * •

we changed the healing distance $d$ and the quencher parameter $\eta$ for the
Jastrow function $f_{c}^{NN}$;

  * •

we used a different set of radial functions, labelled as HF-B1 [232], coming
from Hartree-Fock calculations for the effective interaction B1 of Ref. [233].
The B1 is a phenomenological two-body nucleon-nucleon potential fitted in
order to reproduce the binding energies and densities of various light nuclei
and of nuclear matter in the HF approximation.

Although a central correlation function should not affect the computed energy
value, in the version _v2_ of the algorithm an unpleasant dependence on
$f_{c}^{NN}$ was found, and in particular as a function of the quencher
parameter $\eta$. This dependence is active for both the AV4’ and the AV6’
potentials, regardless of the choice of the single particle orbitals. The time
step extrapolation ($d\tau\rightarrow 0$) of the energy does not solve the
issue. Energy differences are still more than 1 MeV among different setups for
the trial wave function. By varying the parameter $\eta$ from zero (no Jastrow
at all) to one (full central channel of the $NN$ potential used for the
solution of Eq. (3.192)), the energies increase almost linearly. For example,
in the case of AV4’ for the Skyrme orbital functions, the energy of 4He goes
from $-31.3(2)$ MeV for $\eta=0$, to $-27.2(2)$ MeV for $\eta=1$. Same effect
is found for the HF-B1 orbitals with energies going from $-32.5(2)$ MeV to
$-28.4(2)$ MeV. The inclusion of the pure central Jastrow introduces thus
strong biases in the evaluation of the total energy. Moreover, there is also a
dependence on the choice of the single particle orbitals, as shown from the
results for AV4’. Same conclusions follow for the AV6’ potential.

On the grounds of these observations we moved from the AFDMC local energy
scheme to the importance function ratio scheme (version _v1_), with no
importance sampling on auxiliary fields. In this case the bias introduced by
the Jastrow correlation function is still present but reduced to $0.3\div 0.4$
MeV for AV4’ and $0.1\div 0.2$ MeV for AV6’. In spite of the improvement with
respect to the previous case, we decided to remove this source of uncertainty
from the trial wave function and proceed with the test of the _v1_ algorithm
with no Jastrow. It has to be mentioned that a new sampling procedure, for
both coordinates and auxiliary fields, capable to reduce the dependence on
central correlations is being studied.

As shown in Figs. 4.1 and 4.2, the _v1_ extrapolated energies obtained using
different single particle orbitals are consistent within the Monte Carlo
statistical errors, both for the AV4’ and the AV6’ potentials.

Figure 4.1: Binding energy of 4He as a function of the Monte Carlo imaginary
time step. Results are obtained using the AV4’ $NN$ potential. Red dots are
the AFDMC results for the Skyrme radial orbitals. Blue triangles the ones for
the HF-B1 orbitals. For comparison, the GFMC result of Ref. [198], corrected
by the Coulomb contribution (see text for details), is reported with the green
band. Figure 4.2: Binding energy of 4He as a function of the Monte Carlo
imaginary time step. Results are obtained using the AV6’ $NN$ potential. As in
Fig. 4.1, red dots refers to the AFDMC results for the Skyrme radial functions
and blue triangles for the HF-B1 orbitals. The green arrow points to the GFMC
result.

For $\mathcal{N}_{N}=4$ we can compare the AFDMC results with the GFMC ones.
In our calculations the Coulomb interaction is not included. A precise VMC
estimate, that should be representative also for the GFMC estimate, of the
Coulomb expectation value for 4He is $0.785(2)$ MeV [234]. The AFDMC values of
$-32.67(8)$ MeV (Skyrme) and $-32.7(1)$ MeV (HF-B1) for 4He with AV4’ are thus
very close to the GFMC $-32.11(2)$ MeV of Ref. [198] for the same potential
once the Coulomb contribution is subtracted. Our results are still $\sim
0.1\div 0.2$ MeV above the GFMC one, most likely due to the removal of the
sign problem constraint applied at the end of the GFMC runs (release node
procedure [235]).

Although AFDMC and GFMC energies for 4He described by the AV4’ potential are
consistent, a clear problem appears using the AV6’ interaction (Fig. 4.2).
With the two sets of radial functions, the energies are $-19.59(8)$ MeV
(Skyrme) and $-19.53(13)$ MeV (HF-B1) and thus the AFDMC actually projects out
the same ground state. However, the GFMC estimate is $-26.15(2)$ MeV minus the
Coulomb contribution. This large difference in the energies cannot be
attributed to the GFMC release node procedure. The difference in using AV4’
and AV6’ is the inclusion of the tensor term $S_{ij}$ of Eq. (2.8). The
Hamiltonian moves then from real to complex and this might result in a phase
problem during the imaginary time propagation. There might be some issues with
the fixed phase approximation or with the too poor trial wave function (or
both), which does not include operatorial correlations. This is still an
unsolved question but many ideas are being tested. According to the lack of
control on the AFDMC simulations for the AV6’ potential, from now on we will
limit the study to AV4’. As we shall see, this choice does not affect the
result on energy differences as the hyperon separation energy, which is the
main observable of this study for finite systems.

In order to complete the check of the accuracy of the algorithm _v1_ for 4He,
we performed simulations using the Minnesota potential of Ref. [236]. This
two-nucleon interaction has the same operator structure of AV4’ but much
softer cores. Our AFDMC result for the energy is $-30.69(7)$ MeV. It has to be
compared with the $-29.937$ MeV ($-30.722$ MeV with the Coulomb subtraction)
obtained with the Stochastic Variational Method (SVM) [237], that has been
proven to give consistent results with the GFMC algorithm for 4He [21]. The
agreement of the results is remarkable.

Moreover, we tested the consistency of the _v1_ algorithm for the AV4’
potential by studying the deuteron, tritium and oxygen nuclei.

  * •

The AFDMC binding energy for 2H is $-2.22(5)$ MeV, in agreement with the
experimental $-2.225$ MeV. The result is significant because, although the
Argonne V4’ was exactly fitted in order to reproduce the deuteron energy, our
starting trial wave function is just a Slater determinant of single particle
orbitals, with no correlations.

  * •

The result for 3H is $-8.74(4)$ MeV, close to the GFMC $-8.99(1)$ MeV of Ref.
[198]. As for 4He, the small difference in the energies is probably due to the
release node procedure in GFMC. Without the Coulomb contribution, we obtained
the same energy $-8.75(4)$ MeV also for 3He. In AV4’ there are no charge
symmetry breaking terms. Therefore, this result can be seen as a consistency
test on the correct treatment of the spin-isospin operators acting on the wave
function during the Hubbard-Stratonovich rotations.

  * •

For 16O we found the energy values of $-176.8(5)$ MeV for the Skyrme orbitals
and $-174.3(8)$ MeV for the HF-B1 radial functions. The energy difference is
of order 1% even for a medium mass nucleus. The projection mechanism is
working accurately regardless the starting trial function. GFMC results are
limited to 12 nucleons [17, 18, 19], so we cannot compare the two methods for
$\mathcal{N}_{N}=16$. The binding energy cannot be compared with the
experimental data due to the poor employed Hamiltonian. However the AFDMC
results are consistent with the overbinding predicted by the available GFMC
energies for AV4’ [198] and the nucleus results stable under alpha particle
break down, as expected.

On the grounds of the results of these consistency checks, in the present work
we adopt the version _v1_ of the AFDMC algorithm employing the nuclear
potential AV4’ for both nuclei and hypernuclei.

System | $E_{\text{AFDMC}}$ | $E_{\text{GFMC}}$ | $E_{\text{exp}}$ | $E_{\text{AFDMC}}/\mathcal{N}_{N}$ | $E_{\text{exp}}/\mathcal{N}_{N}$
---|---|---|---|---|---
2H | -2.22(5) | — | -2.225 | -1.11 | -1.11
3H | -8.74(4) | -8.99(1) | -8.482 | -2.91 | -2.83
3He | -8.75(4) | — | -7.718 | -2.92 | -2.57
4He | -32.67(8) | -32.90(3) | -28.296 | -8.17 | -7.07
5He | -27.96(13) | -31.26(4) | -27.406 | -5.59 | -5.48
6He | -29.87(14) | -33.00(5) | -29.271 | -4.98 | -4.88
12C | -77.31(25)* | — | -92.162 | -6.44 | -7.68
15O | -144.9(4) | — | -111.955 | -9.66 | -7.46
16O | -176.8(5) | — | -127.619 | -11.05 | -7.98
17O | -177.0(6) | — | -131.762 | -10.41 | -7.75
40Ca | -597(3) | — | -342.052 | -14.93 | -8.55
48Ca | -645(3) | — | -416.001 | -13.44 | -8.67
90Zr | -1457(6) | — | -783.899 | -16.19 | -8.71
Table 4.1: Binding energies (in MeV) for different nuclei. AFDMC and GFMC
results are obtained using the the AV4’ $NN$ potential. The GFMC data are from
Ref. [198] corrected by the Coulomb contribution (see text for details). In
the fourth column the experimental results are from Ref. [238]. Errors are
less than 0.1 KeV. In the last two columns the calculated and experimental
binding energies per particle. For the note * on 12C see the text.

As reported in Tab. 4.1, the resulting absolute binding energies using AV4’
are not comparable with experimental ones, as expected, due to the lack of
information about the nucleon interaction in the Hamiltonian. With the
increase of the number of particles, the simulated nuclei become more an more
bound until the limit case of 90Zr, for which the estimated binding energy is
almost twice the experimental one. Looking at the results for helium isotopes,
we can see that for $\mathcal{N}_{N}=3$ and $4$ the energies are compatible
with GFMC calculations, once the Coulomb contribution is removed. For 5He and
6He instead, we obtained discrepancies between the results for the two
methods. However this is an expected result. When moving to open shell
systems, as 5He and 6He with one or two neutrons out of the first $s$ shell,
the structure of the wave function becomes more complicated and results are
more dependent on the employed $\psi_{T}$. For example, in the case of 6He, in
order to have total angular momentum zero, the two external neutrons can
occupy the $1p_{3/2}$ or the $1p_{1/2}$ orbitals of the nuclear shell model
classification. By using just one of the two $p$ shells, one gets the
unphysical result $E(^{5}\text{He})<E(^{6}\text{He})$. The reported binding
energy has been instead obtained by considering the linear combination of the
Slater determinants giving $J=0$

$\displaystyle\Phi(R_{N},S_{N})=(1-c)\det\Bigl{\\{}\varphi_{\epsilon}^{N}(\bm{r}_{i},s_{i})\Bigr{\\}}_{1p_{3/2}}\\!\\!+c\,\det\Bigl{\\{}\varphi_{\epsilon}^{N}(\bm{r}_{i},s_{i})\Bigr{\\}}_{1p_{1/2}}\;,$
(4.5)

and minimizing the energy with respect to the mixing parameter $c$, as shown
in Fig. 4.3. However the final result is still far from the GFMC data. This is
a clear indication that better wave functions are needed for open shell
systems. A confirmation of that is the non physical result obtained for the
12C nucleus (marked in Tab. 4.1 with *), which is even less bound than
expected, although the employed AV4’ potential, resulting thus unstable under
$\alpha$ break down. In the case of $\mathcal{N}_{N}=12$ indeed, the 8
additional neutrons and protons to the alpha core have just been placed in the
$1p_{3/2}$ shell without any linear combination of the other possible setups
giving zero total angular momentum. This result will be useful in the hyperon
separation energy estimate anyway. In fact, we shall see in the next section
that regardless the total binding energies, by using the same nucleon
potential to describe nuclei and the core of hypernuclei, the obtained hyperon
separation energy is in any case realistic.

Last comment on a technical detail regarding the computation of AFDMC
observables. As shown in Fig. 4.1 and 4.2, the extrapolation of the energy
values in the limit $d\tau\rightarrow 0$ is linear. This is consistent with
the application of the Trotter-Suzuki formula of Eq. (3.17) in the Hubbard-
Stratonovich transformation (3.136), that is thus correct at order
$\sqrt{d\tau}^{\,2}$. Focusing on the AV4’ case, for 4He the time step
extrapolation is almost flat. The differences between the final results and
the energies computed at large $d\tau$ are less than $0.5\%$ and almost within
the statistical errors of the Monte Carlo run. The situation dramatically
changes with the increase of the particle number. For $\mathcal{N}_{N}=16$
this difference is around $2\%$. For 40 and 48 particles, large time step
values and the extrapolated ones are, respectively, $6\%$ and $8.5\%$
different. Therefore, the binding energies must always be carefully studied by
varying the time step of the AFDMC run. The same behavior has been found for
observables other than the total energy (single particle densities and radii).
Each reported result in this chapter has been thus obtained by means of a
computationally expensive procedure of imaginary time extrapolation.

Figure 4.3: 6He binding energy as a function of the mixing parameter $c$ of
Eq. (4.5). The arrows point to the results for the pure $1p_{3/2}$
($-27.65(8)$ MeV) and $1p_{1/2}$ ($-25.98(8)$ MeV) configurations used for the
two external neutrons. The green line is the GFMC result of Ref. [198]
corrected by the VMC Coulomb expectation contribution $0.776(2)$ MeV [234].

###  4.2 Single $\Lambda$ hypernuclei

When a single $\Lambda$ particle is added to a core nucleus, the wave function
of Eq. (3.202) is given by

$\displaystyle\psi_{T}(R,S)=\prod_{i}f_{c}^{\Lambda N}(r_{\Lambda
i})\,\psi_{T}^{N}(R_{N},S_{N})\,\varphi_{\epsilon}^{\Lambda}(\bm{r}_{\Lambda},s_{\Lambda})\;.$
(4.6)

The structure of the nucleon trial wave function is the same of Eq. (3.183),
used in the AFDMC calculations for nuclei. The hyperon Slater determinant is
simply replaced by the single particle state $\varphi_{\epsilon}^{\Lambda}$ of
Eq. (3.209), assumed to be the neutron $1s_{1/2}$ radial function, as already
described in the previous chapter. In order to be consistent with the
calculations for nuclei, we neglected the Jastrow $\Lambda N$ correlation
function which was found to produce a similar but smaller bias on the total
energy. As radial functions we used the same Skyrme set employed in the
calculations for the nuclei of Tab. 4.1.

The $\Lambda$ separation energies defined in Eq. (4.1), are calculated by
taking the difference between the nuclei binding energies presented in the
previous section, and the AFDMC energies for hypernuclei, given the same
nucleon potential. By looking at energy differences, we studied the
contribution of the $\Lambda N$ and $\Lambda NN$ terms defined in Chapter 2.
By comparing AFDMC results with the expected hyperon separation energies,
information about the hyperon-nucleon interaction are deduced. Some
qualitative properties have been also obtained by studying the nucleon and
hyperon single particle densities and the root mean square radii.

####  4.2.1 Hyperon separation energies

We begin the study of $\Lambda$ hypernuclei with the analysis of closed shell
hypernuclei, in particular ${}^{5}_{\Lambda}$He and
${}^{17}_{\leavevmode\nobreak\ \Lambda}$O. We have seen in the previous
section that the AFDMC algorithm is most accurate in describing closed shell
nuclei. Results for 4He and 16O with the AV4’ potential are indeed consistent
and under control. This give us the possibility to realistically describe the
hyperon separation energy for such systems and deduce some general properties
of the employed hyperon-nucleon force.

The step zero of this study was the inclusion in the Hamiltonian of the $NN$
AV4’ interaction and the two-body $\Lambda N$ charge symmetric potential of
Eq. (2.35). The employed parameters $\bar{v}$ and $v_{\sigma}$ are reported in
Tab 2.1. The exchange parameter $\varepsilon$ has been initially set to zero
due to the impossibility of including the space exchange operator directly in
the AFDMC propagator (see § 3.2.4). As reported in Tab. 4.2, the AV4’
$\Lambda$ separation energy for ${}^{5}_{\Lambda}$He is more than twice the
expected value. For the heavier ${}^{17}_{\leavevmode\nobreak\ \Lambda}$O the
discrepancy is even larger. Actually, this is an expected result. As firstly
pointed out by Dalitz [148], $\Lambda N$ potentials, parameterized to account
for the low-energy $\Lambda N$ scattering data and the binding energy of the
$A=3,4$ hypernuclei, overbind ${}^{5}_{\Lambda}$He by $2\div 3$ MeV. That is,
the calculated $A=5$ $\Lambda$ separation energy is about a factor of 2 too
large. This fact is usually reported as _$A=5$ anomaly_. With only a $\Lambda
N$ potential fitted to $\Lambda p$ scattering, the heavier hypernuclei result
then strongly overbound.

$NN$ potential |  ${}^{5}_{\Lambda}$He |  ${}^{17}_{\leavevmode\nobreak\ \Lambda}$O
---|---|---
$V_{\Lambda N}$ |  $V_{\Lambda N}$+$V_{\Lambda NN}$ |  $V_{\Lambda N}$ |  $V_{\Lambda N}$+$V_{\Lambda NN}$
Argonne V4’ |  7.1(1) | 5.1(1) |  43(1) | 19(1)
Argonne V6’ |  6.3(1) | 5.2(1) |  34(1) | 21(1)
Minnesota |  7.4(1) | 5.2(1) |  50(1) | 17(2)
Expt. | 3.12(2) | 13.0(4)
Table 4.2: $\Lambda$ separation energies (in MeV) for ${}^{5}_{\Lambda}$He and
${}^{17}_{\leavevmode\nobreak\ \Lambda}$O obtained using different nucleon
potentials (AV4’, AV6’, Minnesota) and different hyperon-nucleon interaction
(two-body alone and two-body plus three-body, set of parameters (I)) [41]. In
the last line the experimental $B_{\Lambda}$ for ${}^{5}_{\Lambda}$He is from
Ref. [77]. Since no experimental data for ${}^{17}_{\leavevmode\nobreak\
\Lambda}$O exists, the reference separation energy is the semiempirical value
reported in Ref. [192].

As suggested by the same Dalitz [148] and successively by Bodmer and Usmani
[181], the inclusion of a $\Lambda$-nucleon-nucleon potential may solve the
overbinding problem. This is indeed the case, as reported for instance in
Refs. [192, 184, 189]. Therefore, in our AFDMC calculations we included the
three-body $\Lambda NN$ interaction developed by Bodmer, Usmani and Carlson
and described in § 2.2.2. Among the available parametrizations coming from
different VMC studies of light hypernuclei, the set of parameters for the
$\Lambda NN$ potential has been originally taken from Ref. [185], being the
choice that made the variational $B_{\Lambda}$ for ${}_{\Lambda}^{5}$He and
${}^{17}_{\leavevmode\nobreak\ \Lambda}$O compatible with the expected
results. It reads:

$\hypertarget{par_I}{(\text{I})}\phantom{I}\quad\left\\{\begin{array}[]{rcll}C_{P}&\\!=\\!&0.60&\\!\text{MeV}\\\
C_{S}&\\!=\\!&0.00&\\!\text{MeV}\\\
W_{D}&\\!=\\!&0.015&\\!\text{MeV}\end{array}\right.$

The inclusion of the $\Lambda NN$ force reduces the overbinding and thus the
hyperon separation energies, as reported in Tab. 4.2. Although the results are
still not compatible with the experimental ones, the gain in energy due to the
inclusion of the three-body hypernuclear force is considerable.

It has to be pointed out that this result might in principle depend on the
particular choice of the $NN$ interaction used to describe both nucleus and
hypernucleus. One of the main mechanisms that might generate this dependence
might be due to the different environment experienced by the hyperon in the
hypernucleus because of the different nucleon densities and correlations
generated by each $NN$ potential. To discuss this possible dependence, we
performed calculations with different $NN$ interactions having very different
saturation properties. As it can be seen from Tab. 4.2, for
${}^{5}_{\Lambda}$He the extrapolated $B_{\Lambda}$ values with the two-body
$\Lambda N$ interaction alone are about 10% off and well outside statistical
errors. In contrast, the inclusion of the three-body $\Lambda NN$ force gives
a similar $\Lambda$ binding energy independently to the choice of the $NN$
force. On the grounds of this observation, we feel confident that the use of
AV4’, for which AFDMC calculations for nuclei are under control, will in any
case return realistic estimates of $B_{\Lambda}$ for larger masses when
including the $\Lambda NN$ interaction. We checked this assumption performing
simulations in ${}^{17}_{\leavevmode\nobreak\ \Lambda}$O, where the
discrepancy between the $\Lambda$ separation energy computed using the
different $NN$ interactions and the full $\Lambda N$+$\Lambda NN$ force is
less than few per cent (last column of Tab.4.2). The various $NN$ forces
considered here are quite different. The AV6’ includes a tensor force, while
AV4’ and Minnesota have a simpler structure with a similar operator structure
but very different intermediate- and short-range correlations. The fact that
the inclusion of the $\Lambda NN$ force does not depend too much on the
nuclear Hamiltonian is quite remarkable, because the different $NN$ forces
produce a quite different saturation point for the nuclear matter EoS,
suggesting that our results are pretty robust.

Figure 4.4: $\Lambda$ separation energy as a function of $A$ for closed shell
hypernuclei, adapted from Ref. [41]. Solid green dots (dashed curve) are the
available $B_{\Lambda}$ experimental or semiempirical values. Empty red dots
(upper banded curve) refer to the AFDMC results for the two-body $\Lambda N$
interaction alone. Empty blue diamonds (middle banded curve) are the results
with the inclusion also of the three-body hyperon-nucleon force in the
parametrization (I).

For ${}^{5}_{\Lambda}$He the hyperon separation energy with the inclusion of
the $\Lambda NN$ force with the set of parameters (I) reduces of a factor
$\sim 1.4$. For ${}^{17}_{\leavevmode\nobreak\ \Lambda}$O the variation is
around $40\div 50\%$. In order to check the effect of the three-body force
with increasing the particle number, we performed simulations for the next
heavier closed or semi-closed shell $\Lambda$ hypernuclei,
${}^{41}_{\leavevmode\nobreak\ \Lambda}$Ca and ${}^{91}_{\leavevmode\nobreak\
\Lambda}$Zr. The $\Lambda$ separation energies for all the studied closed
shell hypernuclei are shown in Fig. 4.4. While the results for lighter
hypernuclei might be inconclusive in terms of the physical consistency of the
$\Lambda NN$ contribution to the hyperon binding energy in AFDMC calculations,
the computations for ${}^{41}_{\leavevmode\nobreak\ \Lambda}$Ca and
${}^{91}_{\leavevmode\nobreak\ \Lambda}$Zr reveal a completely different
picture. The saturation binding energy provided by the $\Lambda N$ force alone
is completely unrealistic, while the inclusion of the $\Lambda NN$ force gives
results that are much closer to the experimental behavior. Therefore, the
$\Lambda$-nucleon-nucleon force gives a very important repulsive contribution
towards a realistic description of the saturation of medium-heavy hypernuclei
[41]. However, with the given parametrization, only a qualitative agreement
wiht the expected separation energies is reproduced. A refitting procedure for
the three-body hyperon-nucleon interaction might thus improve the quality of
the results.

As already discussed in § 2.2.2, the $C_{S}$ parameter can be estimated by
comparing the $S$-wave term of the $\Lambda NN$ force with the Tucson-
Melbourne component of the $NNN$ interaction. We take the suggested
$C_{S}=1.50$ MeV value [189], in order to reduce the number of fitting
parameters. This choice is justified because the $S$-wave component of the
three-body $\Lambda NN$ interaction is sub-leading. We indeed verified that a
change in the $C_{S}$ value yields a variation of the total energy within
statistical error bars, and definitely much smaller than the variation in
energy due to a change of the $W_{D}$ parameter.

Figure 4.5: $\Lambda$ separation energy for ${}^{5}_{\Lambda}$He as a function
of strengths $W_{D}$ and $C_{P}$ of the three-body $\Lambda NN$ interaction
[43]. The red grid represents the experimental $B_{\Lambda}=3.12(2)$ MeV [77].
The dashed yellow curve is the interception between the expected result and
the $B_{\Lambda}$ surface in the $W_{D}-C_{P}$ parameter space. Statistical
error bars on AFDMC results (solid black dots) are of the order of $0.10\div
0.15$ MeV. Figure 4.6: Projection of Fig. 4.5 on the $W_{D}-C_{P}$ plane [43].
Error bars come from a realistic conservative estimate of the uncertainty in
the determination of the parameters due to the statistical errors of the Monte
Carlo calculations. Blue and green dashed, long-dashed and dot-dashed lines
(lower curves) are the variational results of Ref. [189] for different
$\varepsilon$ and $\bar{v}$ (two-body $\Lambda N$ potential). The dashed box
corresponds to the parameter domain of Fig. 4.5. Black empty dots and the red
band (upper curve) are the projected interception describing the possible set
of parameters reproducing the experimental $B_{\Lambda}$.

In Fig. 4.5 we report the systematic study of the $\Lambda$ separation energy
of ${}_{\Lambda}^{5}$He as a function of both $W_{D}$ and $C_{P}$. Solid black
dots are the AFDMC results. The red grid represents the experimental
$B_{\Lambda}=3.12(2)$ MeV [77]. The dashed yellow curve follows the set of
parameters reproducing the expected $\Lambda$ separation energy. The same
curve is also reported in Fig. 4.6 (red banded curve with black empty dots and
error bars), that is a projection of Fig. 4.5 on the $W_{D}-C_{P}$ plane. The
dashed box represents the $W_{D}$ and $C_{P}$ domain of the previous picture.
For comparison, also the variational results of Ref. [189] are reported. Green
curves are the results for $\bar{v}=6.15$ MeV and $v_{\sigma}=0.24$ MeV, blue
ones for $\bar{v}=6.10$ MeV and $v_{\sigma}=0.24$ MeV. Dashed, long-dashed and
dot-dashed lines correspond respectively to $\varepsilon=0.1$, $0.2$ and
$0.3$. In our calculations we have not considered different combinations for
the parameters of the two-body $\Lambda N$ interaction, focusing on the three-
body part. We have thus kept fixed $\bar{v}$ and $v_{\sigma}$ to the same
values of the green curves of Fig. 4.6 which are the same reported in Tab.
2.1. Moreover, we have set $\varepsilon=0$ for all the hypernuclei studied due
to the impossibility of exactly including the $\mathcal{P}_{x}$ exchange
operator in the propagator. A perturbative analysis of the effect of the
$v_{0}(r)\varepsilon(\mathcal{P}_{x}-1)$ term on the hyperon separation energy
is reported in § 4.2.1.

As it can be seen from Fig. 4.5, $B_{\Lambda}$ significantly increases with
the increase in $C_{P}$, while it decreases with $W_{D}$. This result is
consistent with the attractive nature of $V_{\Lambda NN}^{2\pi,P}$ and the
repulsion effect induced by $V_{\Lambda NN}^{D}$. It is also in agreement with
all the variational estimates on ${}^{5}_{\Lambda}$He (see for instance Refs.
[184, 189]). Starting from the analysis of the results in the $W_{D}-C_{P}$
space for ${}_{\Lambda}^{5}$He, we performed simulations for the next closed
shell hypernucleus ${}^{17}_{\leavevmode\nobreak\ \Lambda}$O. Using the
parameters in the red band of Fig. 4.6 we identified a parametrization able to
reproduce the experimental $B_{\Lambda}$ for both ${}_{\Lambda}^{5}$He and
${}^{17}_{\leavevmode\nobreak\ \Lambda}$O at the same time within the AFDMC
framework:

$\hypertarget{par_II}{(\text{II})}\quad\left\\{\begin{array}[]{rcll}C_{P}&\\!=\\!&1.00&\\!\text{MeV}\\\
C_{S}&\\!=\\!&1.50&\\!\text{MeV}\\\
W_{D}&\\!=\\!&0.035&\\!\text{MeV}\end{array}\right.$

Given the set (II), the $\Lambda$ separation energy of the closed shell
hypernuclei reported in Fig. 4.4 has been re-calculated. We have seen that
$B_{\Lambda}$ is not sensitive neither to the details of the $NN$ interaction,
nor to the total binding energies of nuclei and hypernuclei, as verified by
the good results in Tab. 4.2 even for the problematic case of AV6’ (see §
4.1). On the grounds of this observation, we tried to simulate also open shell
hypernuclei, using the $\Lambda N$, $\Lambda NN$ set (I) and $\Lambda NN$ set
(II) potentials. The binding energies for these systems might not be accurate,
as in the case of the corresponding nuclei. The hyperon separation energy is
expected to be in any case realistic. All the results obtained so far in the
mass range $3\leq A\leq 91$ are summarized in Fig. 4.7 and Fig. 4.8.

Figure 4.7: $\Lambda$ separation energy as a function of $A$. Solid green dots
(dashed curve) are the available $B_{\Lambda}$ experimental or semiempirical
values. Empty red dots (upper banded curve) refer to the AFDMC results for the
two-body $\Lambda N$ interaction alone. Empty blue diamonds (middle banded
curve) and empty black triangles (lower banded curve) are the results with the
inclusion also of the three-body hyperon-nucleon force, respectively for the
parametrizations (I) and (II). Figure 4.8: $\Lambda$ separation energy as a
function of $A^{-2/3}$, adapted from Ref. [43]. The key is the same of Fig.
4.7.

We report $B_{\Lambda}$ as a function of $A$ and $A^{-2/3}$, which is an
approximation of the $A$ dependence of the kinetic term of the Hamiltonian.
Solid green dots are the available experimental data, empty symbols the AFDMC
results. The red curve is obtained using only the two-body hyperon-nucleon
interaction in addition to the nuclear AV4’ potential. The blue curve refers
to the results for the same systems when also the three-body $\Lambda NN$
interaction with the old set of parameters (I) is included. The black lower
curve shows the results obtained by including the three-body hyperon-nucleon
interaction described by the new parametrization (II). A detailed comparison
between numerical and experimental results for the hyperon-separation energy
is given in Tab. 4.3.

System | $E$ | $B_{\Lambda}$ | Expt. $B_{\Lambda}$
---|---|---|---
${}^{3}_{\Lambda}$H | -1.00(14) | -1.22(15) | 0.13(5) [77]
${}^{4}_{\Lambda}$H | -9.69(8) | 0.95(9) | 2.04(4) [77]
${}^{4}_{\Lambda}$He | -9.97(8) | 1.22(9) | 2.39(3) [77]
${}^{5}_{\Lambda}$He | -35.89(12) | 3.22(14) | 3.12(2) [77]
${}^{6}_{\Lambda}$He | -32.72(15) | 4.76(20) | 4.25(10) [77]
${}^{7}_{\Lambda}$He | -35.82(15) | 5.95(25) | 5.68(28) [86]
${}^{13}_{\leavevmode\nobreak\ \Lambda}$C | -88.5(26)* | 11.2(4) | 11.69(12) [78]
${}^{16}_{\leavevmode\nobreak\ \Lambda}$O | -157.5(6) | 12.6(7) | 12.50(35) [80]
${}^{17}_{\leavevmode\nobreak\ \Lambda}$O | -189.2(4) | 12.4(6) | 13.0(4) [192]
${}^{18}_{\leavevmode\nobreak\ \Lambda}$O | -189.7(6) | 12.7(9) | —
${}^{41}_{\leavevmode\nobreak\ \Lambda}$Ca | -616(3) | 19(4) | 19.24(0) [181]
${}^{49}_{\leavevmode\nobreak\ \Lambda}$Ca | -665(4) | 20(5) | —
${}^{91}_{\leavevmode\nobreak\ \Lambda}$Zr | -1478(7) | 21(9) | 23.33(0) [181]
Table 4.3: Binding energies and $\Lambda$ separation energies (in MeV)
obtained using the two-body plus three-body hyperon-nucleon interaction with
the set of parameters (II) [43]. The results already include the CSB
contribution. The effect is evident only for light systems, as discussed in
the next section. In the last column, the expected $B_{\Lambda}$ values. Since
no experimental data for $A=17,41,91$ exists, the reference separation
energies are semiempirical values.

From Fig. 4.7 and Fig. 4.8 we can see that the new parametrization for the
three-body hyperon-nucleon interaction correctly reproduces the experimental
saturation property of the $\Lambda$ separation energy. All the separation
energies for $A\geq 5$ are compatible or very close to the expected results,
where available, as reported in Tab. 4.3. Since for
${}^{18}_{\leavevmode\nobreak\ \Lambda}$O and ${}^{49}_{\leavevmode\nobreak\
\Lambda}$Ca no experimental data have been found, the values of $12.7(9)$ MeV
and $20(5)$ MeV are AFDMC predictions, that follows the general trend of the
experimental curve. Although for $A\geq 41$ the Monte Carlo statistical error
bars become rather large, the extrapolation of the $\Lambda$ binding energy
for $A\rightarrow\infty$ points to the correct region for the expected value
$D_{\Lambda}\sim 30$ MeV of $s_{\Lambda}$ states in nuclear matter.

We can find the same problems discussed in the case of nuclei (§ 4.1) in the
analysis of the total hypernuclear binding energies for $A\geq 5$. For
instance, the binding energy of ${}^{13}_{\leavevmode\nobreak\ \Lambda}$C is
non physical, as for the energy of the core nucleus 12C. However, the energy
difference is consistent with the expected result. Moreover, for the core wave
function of ${}^{7}_{\Lambda}$He we have used the same mixing parameter
adopted in the description of 6He (see Eq. (4.5) and Fig. 4.3), in order to
have at least the correct ordering in the hypernuclear energy spectrum.
However, the same hyperon separation energy can be found by just using the
$1p_{3/2}$ shell for the outer neutrons for both strange and non strange
nucleus. Our working hypothesis regarding the computation of the hyperon
separation energy is thus correct, at least for medium-heavy hypernuclei.

For $A<5$ our results are more than 1 MeV off from experimental data. For
${}^{3}_{\Lambda}$H, the $\Lambda$ separation energy is even negative, meaning
that the hypernucleus is less bound than the corresponding nucleus 2H. We can
ascribe this discrepancy to the lack of accuracy of our wave function for few-
body systems. Since the $\Lambda$ hyperon does not suffer from Pauli blocking
by the other nucleons, it can penetrate into the nuclear interior and form
deeply bound hypernuclear states. For heavy systems the $\Lambda$ particle can
be seen as an impurity that does not drastically alter the nuclear wave
function. Therefore, the trial wave function of Eq. (4.6) with the single
particle state $\varphi_{\epsilon}^{\Lambda}$ described by the $1s_{1/2}$
neutron orbital, is accurate enough as starting point for the imaginary time
propagation. For very light hypernuclei, for which the first nucleonic $s$
shell is not closed, this might not be the case. In order to have a correct
projection onto the ground state, the single particle orbitals of both
nucleons and lambda might need to be changed when the hyperon is added to the
nucleus. Moreover, in very light hypernuclei, the neglected nucleon-nucleon
and hyperon-nucleon correlations, might result in non negligible contributions
to the $\Lambda$ binding energy. A study of these systems within a few-body
method or a different projection algorithm like the GFMC, might solve this
issue.

##### Effect of the charge symmetry breaking term

The effect of the CSB potential has been studied for the $A=4$ mirror
hypernuclei. As reported in Tab. 4.4, without the CSB term there is no
difference in the $\Lambda$ binding energy of ${}^{4}_{\Lambda}$H and
${}^{4}_{\Lambda}$He. When CSB is active, a splitting appears due to the
different behavior of the $\Lambda p$ and $\Lambda n$ channels. The strength
of the difference $\Delta B_{\Lambda}^{CSB}$ is independent on the parameters
of the three-body $\Lambda NN$ interaction and it is compatible with the
experimental result [77], although the $\Lambda$ separation energies are not
accurate.

Parameters | System | $B_{\Lambda}^{sym}$ | $B_{\Lambda}^{CSB}$ |  $\Delta B_{\Lambda}^{CSB}$
---|---|---|---|---
Set (I) |  ${}^{4}_{\Lambda}$H | 1.97(11) | 1.89(9) |  0.24(12)
${}^{4}_{\Lambda}$He | 2.02(10) | 2.13(8)
Set (II) |  ${}^{4}_{\Lambda}$H | 1.07(8) | 0.95(9) |  0.27(13)
${}^{4}_{\Lambda}$He | 1.07(9) | 1.22(9)
Expt. [77] |  ${}^{4}_{\Lambda}$H | — | 2.04(4) |  0.35(5)0
${}^{4}_{\Lambda}$He | — | 2.39(3)
Table 4.4: $\Lambda$ separation energies (in MeV) for the $A=4$ mirror $\Lambda$ hypernuclei with (fourth column) and without (third column) the inclusion of the charge symmetry breaking term [43]. In the last column the difference in the separation energy induced by the CSB interaction. First and second rows refer to different set of parameters for the $\Lambda NN$ interaction, while the last row is the experimental result. System | $p$ | $n$ | $\Delta_{np}$ |  $\Delta B_{\Lambda}$
---|---|---|---|---
${}^{4}_{\Lambda}$H | 1 | 2 | $+1$ |  $-0.12(8)$
${}^{4}_{\Lambda}$He | 2 | 1 | $-1$ |  $+0.15(9)$
${}^{5}_{\Lambda}$He | 2 | 2 |  $0$ |  $+0.02(9)$
${}^{6}_{\Lambda}$He | 2 | 3 | $+1$ |  $-0.06(8)$
${}^{7}_{\Lambda}$He | 2 | 4 | $+2$ |  $-0.18(8)$
${}^{16}_{\leavevmode\nobreak\ \Lambda}$O | 8 | 7 | $-1$ |  $+0.27(35)$
${}^{17}_{\leavevmode\nobreak\ \Lambda}$O | 8 | 8 |  $0$ |  $+0.15(35)$
${}^{18}_{\leavevmode\nobreak\ \Lambda}$O | 8 | 9 | $+1$ |  $-0.74(49)$
Table 4.5: Difference (in MeV) in the hyperon separation energies induced by
the CSB term (Eq. (2.42)) for different hypernuclei [43]. The fourth column
reports the difference between the number of neutrons and protons. Results are
obtained with the full two- plus three-body (set (II)) hyperon-nucleon
interaction. In order to reduce the errors, $\Delta B_{\Lambda}$ has been
calculated by taking the difference between total hypernuclear binding
energies, instead of the hyperon separation energies.

The same CSB potential of Eq. (2.42) has been included in the study of
hypernuclei for $A>4$. In Tab. 4.5 the difference in the hyperon separation
energies $\Delta B_{\Lambda}=B_{\Lambda}^{CSB}-B_{\Lambda}^{sym}$ is reported
for different hypernuclei up to $A=18$. The fourth column shows the difference
between the number of neutrons and protons
$\Delta_{np}=\mathcal{N}_{n}-\mathcal{N}_{p}$. For the symmetric hypernuclei
${}^{5}_{\Lambda}$He and ${}^{17}_{\leavevmode\nobreak\ \Lambda}$O the CSB
interaction has no effect, being this difference zero. In the systems with
neutron excess ($\Delta_{np}>0$), the effect of the CSB consists in decreasing
the hyperon separation energy compared to the charge symmetric case. When
$\Delta_{np}$ becomes negative, $\Delta B_{\Lambda}>0$ due to the attraction
induced by the CSB potential in the $\Lambda p$ channel, that produces more
bound hypernuclei. Being $\Delta_{np}$ small, these effects are in any case
rather small and they become almost negligible compared to the statistical
errors on $B_{\Lambda}$ when the number of baryons becomes large enough
($A>16$). However, in the case of $\Lambda$ neutron matter, the CSB term might
have a relevant effect for large enough $\Lambda$ fraction.

##### Effect of the hyperon-nucleon space-exchange term

As already mentioned in the previous chapter, the inclusion of the $\Lambda N$
space exchange operator of Eq. (2.35) in the AFDMC propagator is not yet
possible. In § 3.2.5 we presented a possible perturbative approach for the
treatment of such term. In Tab. 4.6 we report the results of this analysis.

All the results for ${}^{41}_{\leavevmode\nobreak\ \Lambda}$Ca are consistent
within the statistical errors. On the contrary, for lighter systems the
$\Lambda$ separation energy seems rather sensitive to the value of the
exchange parameter $\varepsilon$. Considering larger values for $\varepsilon$,
$B_{\Lambda}$ generally increases. This trend is opposite to what is found for
instance in Ref. [188]. We recall that only the computation of the Hamiltonian
expectation value by means of Eq. (3.24) gives exact results. For other
operators, like the space exchange $\mathcal{P}_{x}$, the pure estimators have
to be calculated with the extrapolation method via the two relations (3.25) or
(3.26). The variational estimate $\langle\mathcal{P}_{x}\rangle_{v}$ is thus
needed. In the mentioned reference, the importance of space exchange
correlations for variational estimates is discussed. Being these correlations
neglected in this work, our perturbative treatment of the $\mathcal{P}_{x}$
contribution might not be accurate. Moreover, the evidence of the importance
of space exchange correlations might invalid the perturbative approach itself.
An effective but more consistent treatment of this term could consist in a
slight change in the strength of the central $\Lambda N$ potential. However,
due to the very limited information about the space exchange parameter and its
effect on single $\Lambda$ hypernuclei heavier than ${}^{5}_{\Lambda}$He, this
approach has not been considered in the present work. Recent calculations of
many hadron systems within an EFT treatment at NLO for the full $SU(3)$
hadronic spectrum confirmed indeed that exchange terms are sub-leading [170].

System | $\varepsilon=0.0$ | $\varepsilon=0.1$ |  $\varepsilon=0.3$
---|---|---|---
${}^{5}_{\Lambda}$He | 3.22(14) | 3.89(15) | 4.67(25)
${}^{17}_{\leavevmode\nobreak\ \Lambda}$O | 12.4(6) | 12.9(9) | 14.0(9)
${}^{41}_{\leavevmode\nobreak\ \Lambda}$Ca | 19(4) | 21(5) | 25(7)
Table 4.6: Variation of the $\Lambda$ separation energy as a consequence of
the exchange potential $v_{0}(r)\varepsilon(\mathcal{P}_{x}-1)$ in the
$\Lambda N$ interaction of Eq. (2.35). The contribution of $\mathcal{P}_{x}$
is treated perturbatively for different value of the parameter $\varepsilon$.
The interaction used is the full AV4’+$\Lambda N$+$\Lambda NN$ set (II).
Results are expressed in MeV.

####  4.2.2 Single particle densities and radii

Single particle densities can be easily computed in Monte Carlo calculations
by considering the expectation value of the density operator

$\displaystyle\hat{\rho}_{\kappa}(r)=\sum_{i}\delta(r-r_{i})\quad\quad\kappa=N,\Lambda\;,$
(4.7)

where $i$ is the single particle index running over nucleons for
$\rho_{N}=\langle\hat{\rho}_{N}\rangle$ or hyperons for
$\rho_{\Lambda}=\langle\hat{\rho}_{\Lambda}\rangle$. The normalization is
given by

$\displaystyle\int dr4\pi r^{2}\rho_{\kappa}(r)=1\;.$ (4.8)

Root mean square radii $\langle r_{\kappa}^{2}\rangle^{1/2}$ are simply
calculated starting from the Cartesian coordinates of nucleons and hyperons. A
consistency check between AFDMC densities and radii is then taken by verifying
the relation

$\displaystyle\langle r_{\kappa}^{2}\rangle=\int dr4\pi
r^{4}\rho_{\kappa}(r)\;.$ (4.9)

Before reporting the results we recall that also for densities and radii the
AFDMC calculation can only lead to mixed estimators. The pure estimators are
thus approximated by using Eq. (3.25) or Eq. (3.26). The two relations should
lead to consistent results. This is the case for the nucleon and hyperon
radii. In computing the densities instead, the low statistics for
$r\rightarrow 0$ generates differences in the two approaches. For nucleons
these discrepancies are almost within the statistical errors. For hyperons,
the much reduced statistics (1 over $A-1$ for single $\Lambda$ hypernuclei)
and the fact that typically the $\Lambda$ density is not peaked in $r=0$,
create some uncertainties in the region for small $r$, in particular for the
first estimator. We therefore chose to adopt the pure estimator of Eq. (3.26)
to have at least a positive definite estimate. Finally, it has to be pointed
out that the pure extrapolated results are sensitive to the quality of the
variational wave function and the accuracy of the projection sampling
technique. Although we successfully tested the AFDMC propagation, we are
limited in the choice of the VMC wave function. In order to be consistent with
the mixed estimators coming from AFDMC calculations, we considered the same
trial wave functions also for the variational runs. This might introduce some
biases in the evaluation of pure estimators. Therefore, the results presented
in the following have to be considered as a qualitative study on the general
effect of the hypernuclear forces on the nucleon and hyperon distributions.

In Fig. 4.9 we report the results for the single particle densities for 4He
and ${}^{5}_{\Lambda}$He. The green curves are the densities of nucleons in
the nucleus, while the red and blue curves are, respectively, the density of
nucleons and of the lambda particle in the hypernucleus. In the left panel the
results are obtained using AV4’ for the nuclear part and the two-body $\Lambda
N$ interaction alone for the hypernuclear component. In the right panel the
densities are calculated with the full two- plus three-body (set (II))
hyperon-nucleon interaction.

Figure 4.9: Single particle densities for nucleons in 4He [green, upper banded
curve] and for nucleons [red, middle banded curve] and the lambda particle
[blue, lower banded curve] in ${}^{5}_{\Lambda}$He [43]. In the left panel the
results for the two-body $\Lambda N$ interaction alone. In the right panel the
results with the inclusion also of the three-body hyperon-nucleon force in the
parametrization (II). The AV4’ potential has been used for the nuclear core.
Figure 4.10: Single particle densities for the $\Lambda$ particle in different
hypernuclei [43]. Top panel reports the results for the two-body $\Lambda N$
interaction alone. Bottom panel shows the results when the three-body hyperon-
nucleon interaction with the set of parameters (II) is also included. The
nuclear core is described by the AV4’ potential.

The addition of the $\Lambda$ particle to the nuclear core of 4He has the
effect to slightly reduce the nucleon density in the center. The $\Lambda$
particle tries to localize close to $r=0$, enlarging therefore the nucleon
distribution. When the three-body $\Lambda NN$ interaction is turned on (right
panel of Fig. 4.9), the repulsion moves the nucleons to large distances but
the main effect is that the hyperon is pushed away from the center of the
system. As can be seen from Fig. 4.10, this effect is much more evident for
large $A$. When the hypernucleus is described by the $\Lambda N$ interaction
alone, the $\Lambda$ particle is localized near the center, in the range $r<2$
fm (left panel of Fig. 4.10). The inclusion of the three-body $\Lambda NN$
potential forces the hyperon to move from the center, in a region that roughly
correspond to the skin of nucleons (see Tab. 4.7). Although these densities
are strictly dependent to the nuclear interaction, by using the AV6’ potential
we found the same qualitative effects on the $\Lambda$ particle, confirming
the importance of the three-body hyperon-nucleon interaction and its repulsive
nature. Due to the limitations discussed above and the use of too simplified
interactions for the nucleon-nucleon force, the comparison with the available
VMC density profiles [187, 192] is difficult.

System | nucleus | hypernucleus
---|---|---
$r_{N}^{\text{exp}}$ | $r_{N}$ |  $r_{N}$ |  $r_{\Lambda}$
2H - ${}^{3}_{\Lambda}$H |  2.142 | 1.48(8) |  1.9(1) | 2.00(16)
3H - ${}^{4}_{\Lambda}$H |  1.759 | 1.5(1) |  1.77(9) | 2.12(15)
3He - ${}^{4}_{\Lambda}$He |  1.966 | 1.5(1) |  1.77(9) | 2.10(14)
4He - ${}^{5}_{\Lambda}$He |  1.676 | 1.57(9) |  1.58(7) | 2.2(2)
5He - ${}^{6}_{\Lambda}$He |  — | 2.02(16) |  2.16(17) | 2.43(17)
6He - ${}^{7}_{\Lambda}$He |  2.065 | 2.3(2) |  2.4(2) | 2.5(2)
15O - ${}^{16}_{\leavevmode\nobreak\ \Lambda}$O |  — | 2.20(12) |  2.3(1) | 3.2(3)
16O - ${}^{17}_{\leavevmode\nobreak\ \Lambda}$O |  2.699 | 2.16(12) |  2.23(11) | 3.3(3)
17O - ${}^{18}_{\leavevmode\nobreak\ \Lambda}$O |  2.693 | 2.26(13) |  2.32(14) | 3.3(3)
40Ca - ${}^{41}_{\leavevmode\nobreak\ \Lambda}$Ca |  3.478 | 2.8(2) |  2.8(2) | 4.2(5)
48Ca - ${}^{49}_{\leavevmode\nobreak\ \Lambda}$Ca |  3.477 | 3.1(2) |  3.1(2) | 4.3(5)
Table 4.7: Nucleon and hyperon root mean square radii (in fm) for nuclei and
corresponding $\Lambda$ hypernuclei. The employed nucleon-nucleon potential is
AV4’. For the strange sector we used the full two- plus three-body hyperon-
nucleon force in the parametrization (II). The experimental nuclear charge
radii are from Ref. [239]. Errors are on the fourth digit.

In Tab. 4.7 we report the nucleon and hyperon root mean square radii for
nuclei and hypernuclei. The experimental nuclear charge radii are reported as
a reference. AFDMC $r_{N}$, that do not distinguish among protons and
neutrons, are typically smaller than the corresponding experimental results.
This can be understood as a consequence of the employed AV4’ $NN$ interaction
that overbinds nuclei. The main qualitative information is that the hyperon
radii are systematically larger than the nucleon ones, as expected by looking
at the single particle densities. Starting from $A=5$, the nucleon radii in
the nucleus and the corresponding hypernucleus do not change, although the
differences in the nucleon densities for $r\rightarrow 0$. This is due to the
small contribution to the integral (4.9) given by the density for $r$ close to
zero. For the hypernuclei with $A<5$, AFDMC calculations predict larger
$r_{N}$ when the hyperon is added to the core nucleus. This is inconsistent
with the results of Ref. [184], where a shrinking of the core nuclei due to
the presence of the $\Lambda$ particle in $A\leq 5$ hypernuclei is found. We
need to emphasize once more that the results presented in this section are
most likely strictly connected to the employed nucleon-nucleon potential. For
instance, the shrinkage of hypernuclei has been investigated experimentally by
$\gamma$-ray spectroscopy [52, 240]. In the experiment of Ref. [240], by
looking at the electric quadrupole transition probability from the excited
$5/2^{+}$ state to the ground state in ${}^{7}_{\Lambda}$Li, a $19\%$
shrinkage of the intercluster distance was inferred, assuming the two-body
cluster structure core+deuteron. Therefore, the AFDMC study of densities and
radii, differently from the analysis of $\Lambda$ separation energies, cannot
lead to accurate results at this level. It has to be considered as a first
explorative attempt to get hypernuclear structure information from Diffusion
Monte Carlo simulations.

###  4.3 Double $\Lambda$ hypernuclei

In the single particle wave function representation, two $\Lambda$ particles
with antiparallel spin can be added to a core nucleus filling the first
hyperon $s$ shell, assumed to be the neutron $1s_{1/2}$ Skyrme radial function
as in the case of single $\Lambda$ hypernuclei. The complete hypernuclear wave
function is given by Eq. (3.202), where the nucleon trial wave function is the
same used in the AFDMC calculations for nuclei and in this case also the
hyperon Slater determinant is employed. Although the effect on the total
energy introduced by a $\Lambda\Lambda$ correlation function is found to be
negligible, for consistency with the calculations for nuclei and single
$\Lambda$ hypernuclei we neglected the central Jastrow correlations.

The double $\Lambda$ separation energy and the incremental $\Lambda\Lambda$
energy of Eqs. (4.3) and (4.4) are calculated starting from the energy of the
nucleus and the corresponding single and double $\Lambda$ hypernuclei
described by the same $NN$ AV4’ potential. Due to the difficulties in treating
open shell nuclei and the limited amount of data about double $\Lambda$
hypernuclei, we performed the AFDMC study for just the lightest
$\Lambda\Lambda$ hypernucleus for which energy experimental information are
available, ${}^{\;\;\,6}_{\Lambda\Lambda}$He.

####  4.3.1 Hyperon separation energies

In Tab. 4.8 we report the total binding energies for 4He, ${}^{5}_{\Lambda}$He
and ${}^{\;\;\,6}_{\Lambda\Lambda}$He in the second column, the single or
double hyperon separation energies in the third and the incremental binding
energy in the last column. The value of $B_{\Lambda\Lambda}$ confirms the weak
attractive nature of the $\Lambda\Lambda$ interaction [173, 150, 151, 152].
Starting from 4He and adding two hyperons with $B_{\Lambda}=3.22(14)$ MeV, the
energy of ${}^{\;\;\,6}_{\Lambda\Lambda}$He would be $1.0\div 1.5$ MeV less
than the actual AFDMC result. Therefore the $\Lambda\Lambda$ potential of Eq.
(2.48) induces a net attraction between hyperons, at least at this density.

System | $E$ | $B_{\Lambda(\Lambda)}$ |  $\Delta B_{\Lambda\Lambda}$
---|---|---|---
4He | -32.67(8) | — | —
${}^{5}_{\Lambda}$He | -35.89(12) | 3.22(14) | —
${}^{\;\;\,6}_{\Lambda\Lambda}$He | -40.6(3) | 7.9(3) | 1.5(4)
${}^{\;\;\,6}_{\Lambda\Lambda}$He | Expt. [91] | $7.25\pm 0.19^{+0.18}_{-0.11}$ |  $1.01\pm 0.20^{+0.18}_{-0.11}$
Table 4.8: Comparison between 4He and the corresponding single and double
$\Lambda$ hypernuclei [43]. In the second column the total binding energies
are reported. The third column shows the single or double $\Lambda$ separation
energies. In the last column the incremental binding energy $\Delta
B_{\Lambda\Lambda}$ is reported. All the results are obtained using the
complete two- plus three-body (set (II)) hyperon-nucleon interaction with the
addition of the $\Lambda\Lambda$ force of Eq. (2.48). The results are
expressed in MeV.

Our $B_{\Lambda\Lambda}$ and $\Delta B_{\Lambda\Lambda}$ are very close to the
expected results for which the potential has originally been fitted within the
cluster model. The latest data $B_{\Lambda\Lambda}=6.91(0.16)$ MeV and $\Delta
B_{\Lambda\Lambda}=0.67(0.17)$ MeV of Ref. [93] suggest a weaker attractive
force between the two hyperons. A refit of the interaction of the form
proposed in Eq. (2.48) would be required. It would be interesting to study
more double $\Lambda$ hypernuclei within the AFDMC framework with the $\Lambda
N$, $\Lambda NN$ and $\Lambda\Lambda$ interaction proposed. Some experimental
data are available in the range $A=7\div 13$, but there are uncertainties in
the identification of the produced double $\Lambda$ hypernuclei, reflecting in