UNIVERSITÀ DEGLI STUDI DI TRENTO
Facoltà di Scienze Matematiche, Fisiche e Naturali

Dipartimento di Fisica

Tesi di Dottorato di Ricerca in Fisica

Ph.D. Thesis in Physics

From Hypernuclei to Hypermatter:
a Quantum Monte Carlo Study of Strangeness
in Nuclear Structure and Nuclear Astrophysics

> _Tutta la materia di cui siamo fatti noi l’hanno costruita le stelle, tutti
> gli elementi dall’idrogeno all’uranio sono stati fatti nelle reazioni
> nucleari che avvengono nelle supernove, cioè queste stelle molto più grosse
> del Sole che alla fine della loro vita esplodono e sparpagliano nello spazio
> il risultato di tutte le reazioni nucleari avvenute al loro interno. Per cui
> noi siamo veramente figli delle stelle._

> _Il computer non è una macchina intelligente che aiuta le persone stupide,
> anzi, è una macchina stupida che funziona solo nelle mani delle persone
> intelligenti._

###### Contents

  1. Introduction
  2. 1 Strangeness in nuclear systems
    1. 1.1 Hyperons in finite nuclei
    2. 1.2 Hyperons in neutron stars
  3. 2 Hamiltonians
    1. 2.1 Interactions: nucleons
      1. 2.1.1 Two-body $NN$ potential
      2. 2.1.2 Three-body $NNN$ potential
    2. 2.2 Interactions: hyperons and nucleons
      1. 2.2.1 Two-body $\Lambda N$ potential
      2. 2.2.2 Three-body $\Lambda NN$ potential
      3. 2.2.3 Two-body $\Lambda\Lambda$ potential
  4. 3 Method
    1. 3.1 Diffusion Monte Carlo
      1. 3.1.1 Importance Sampling
      2. 3.1.2 Sign Problem
      3. 3.1.3 Spin-isospin degrees of freedom
    2. 3.2 Auxiliary Field Diffusion Monte Carlo
      1. 3.2.1 Propagator for nucleons: $\bm{\sigma}$, $\bm{\sigma}\cdot\bm{\tau}$ and $\bm{\tau}$ terms
      2. 3.2.2 Propagator for neutrons: spin-orbit terms
      3. 3.2.3 Propagator for neutrons: three-body terms
      4. 3.2.4 Wave functions
      5. 3.2.5 Propagator for hypernuclear systems
  5. 4 Results: finite systems
    1. 4.1 Nuclei
    2. 4.2 Single $\Lambda$ hypernuclei
      1. 4.2.1 Hyperon separation energies
      2. 4.2.2 Single particle densities and radii
    3. 4.3 Double $\Lambda$ hypernuclei
      1. 4.3.1 Hyperon separation energies
      2. 4.3.2 Single particle densities and radii
  6. 5 Results: infinite systems
    1. 5.1 Neutron matter
    2. 5.2 $\Lambda$ neutron matter
      1. 5.2.1 Test: fixed $\Lambda$ fraction
      2. 5.2.2 $\Lambda$ threshold density and the equation of state
      3. 5.2.3 Mass-radius relation and the maximum mass
  7. 6 Conclusions
  8. A AFDMC wave functions
    1. A.1 Derivatives of the wave function: CM corrections
    2. A.2 Derivatives of a Slater determinant
  9. B $\Lambda N$ space exchange potential
  10. C Acknowledgements

###### List of Figures

  1. _i_.1 Neutron star structure
  2. 1.1 Strangeness producing reactions
  3. 1.2 $\Lambda$ hypernuclei accessible via different experimental reactions
  4. 1.3 $\Lambda$ hypernuclear chart
  5. 1.4 Hyperon and nucleon chemical potentials
  6. 1.5 Neutron star mass-radius relation: Schulze 2011
  7. 1.6 Neutron star mass-radius relation: Massot 2012
  8. 1.7 Neutron star mass-radius relation: Miyatsu 2012
  9. 1.8 Neutron star mass-radius relation: Bednarek 2012
  10. 2.1 Two-pion exchange processes in the $NNN$ force
  11. 2.2 Three-pion exchange processes in the $NNN$ force
  12. 2.3 Short-range contribution in the $NNN$ force
  13. 2.4 Meson exchange processes in the $\Lambda N$ force
  14. 2.5 Two-pion exchange processes in the $\Lambda NN$ force
  15. 4.1 Binding energies: $E$ vs. $d\tau$ for 4He, Argonnne V4’
  16. 4.2 Binding energies: $E$ vs. $d\tau$ for 4He, Argonnne V6’
  17. 4.3 Binding energies: $E$ vs. mixing parameter for 6He
  18. 4.4 $\Lambda$ separation energy vs. $A$: closed shell hypernuclei
  19. 4.5 $\Lambda$ separation energy for ${}^{5}_{\Lambda}$He vs. $W_{D}-C_{P}$: 3D plot
  20. 4.6 $\Lambda$ separation energy for ${}^{5}_{\Lambda}$He vs. $W_{D}-C_{P}$: 2D plot
  21. 4.7 $\Lambda$ separation energy vs. $A$
  22. 4.8 $\Lambda$ separation energy vs. $A^{-2/3}$
  23. 4.9 Single particle densities: $N$ and $\Lambda$ in 4He and ${}^{5}_{\Lambda}$He
  24. 4.10 Single particle densities: $\Lambda$ in hypernuclei for $3\leq A\leq 91$
  25. 4.11 Single particle densities: $N$ and $\Lambda$ in 4He, ${}^{5}_{\Lambda}$He and ${}^{\;\;\,6}_{\Lambda\Lambda}$He
  26. 5.1 Energy per particle vs. baryon density at fixed $\Lambda$ fraction
  27. 5.2 Pair correlation functions at fixed $\Lambda$ fraction: $\rho_{b}=0.16\leavevmode\nobreak\ \text{fm}^{-3}$
  28. 5.3 Pair correlation functions at fixed $\Lambda$ fraction: $\rho_{b}=0.40\leavevmode\nobreak\ \text{fm}^{-3}$
  29. 5.4 YNM and PNM energy difference vs. $\Lambda$ fraction
  30. 5.5 Hyperon symmetry energy vs. baryon density
  31. 5.6 $x_{\Lambda}(\rho_{b})$ function and $\Lambda$ threshold density
  32. 5.7 YNM equation of state
  33. 5.8 YNM mass-radius relation
  34. 5.9 YNM mass-central density relation

###### List of Tables

  1. 1.1 Nucleon and hyperon properties
  2. 2.1 Parameters of the $\Lambda N$ and $\Lambda NN$ interaction
  3. 2.2 Parameters of the $\Lambda\Lambda$ interaction
  4. 4.1 Binding energies: nuclei, $2\leq A-1\leq 90$
  5. 4.2 $\Lambda$ separation energies: $\Lambda N+\Lambda NN$ set (I) for ${}^{5}_{\Lambda}$He and ${}^{17}_{\leavevmode\nobreak\ \Lambda}$O
  6. 4.3 $\Lambda$ separation energies: $\Lambda$ hypernuclei, $3\leq A\leq 91$
  7. 4.4 $\Lambda$ separation energies: $A=4$ mirror hypernuclei
  8. 4.5 $\Lambda$ separation energies: effect of the CSB potential
  9. 4.6 $\Lambda$ separation energies: effect of the $\Lambda N$ exchange potential
  10. 4.7 Nucleon and hyperon radii in hypernuclei for $3\leq A\leq 49$
  11. 4.8 $\Lambda$ separation energies: ${}^{\;\;\,6}_{\Lambda\Lambda}$He
  12. 4.9 Nucleon and hyperon radii for ${}^{\;\;\,6}_{\Lambda\Lambda}$He
  13. 5.1 Energy per particle: neutron matter
  14. 5.2 Energy per particle: $\Lambda$ neutron matter
  15. 5.3 Baryon number and $\Lambda$ fraction
  16. 5.4 Coefficients of the hyperon symmetry energy fit

Empty page

## Introduction

Neutron stars (NS) are among the densest objects in the Universe, with central
densities several times larger than the nuclear saturation density
$\rho_{0}=0.16\leavevmode\nobreak\ \text{fm}^{-3}$. As soon as the density
significantly exceeds this value, the structure and composition of the NS core
become uncertain. Moving from the surface towards the interior of the star,
the stellar matter undergoes a number of transitions, Fig.  _i_.1. From
electrons and neutron rich ions in the outer envelopes, the composition is
supposed to change to the $npe\mu$ matter in the outer core, a degenerate gas
of neutrons, protons, electrons and muons. At densities larger than $\sim
2\rho_{0}$ the $npe\mu$ assumption can be invalid due to the appearance of new
hadronic degrees of freedom or exotic phases.

Figure _i_.1: Schematic structure of a neutron star. Stellar parameters
strongly depend on the equation of state of the core. Figure taken from Ref.
[1]

In the pioneering work of 1960 [2], Ambartsumyan and Saakyan reported the
first theoretical evidence of hyperons in the core of a NS. Contrary to
terrestrial conditions, where hyperons are unstable and decay into nucleons
through the weak interaction, the equilibrium conditions in a NS can make the
inverse process happen. At densities of the order $2\div 3\rho_{0}$, the
nucleon chemical potential is large enough to make the conversion of nucleons
into hyperons energetically favorable. This conversion reduces the Fermi
pressure exerted by the baryons, and makes the equation of state (EoS) softer.
As a consequence, the maximum mass of the star is typically reduced.

Nowadays many different approaches of hyperonic matter are available, but
there is no general agreement among the predicted results for the EoS and the
maximum mass of a NS including hyperons. Some classes of methods extended to
the hyperonic sector predict that the appearance of hyperons at around $2\div
3\rho_{0}$ leads to a strong softening of EoS and consequently to a large
reduction of the maximum mass. Other approaches, instead, indicate much weaker
effects as a consequence of the presence of strange baryons in the core of the
star.

The situation has recently become even more controversial as a result of the
latest astrophysical observations. Until 2010, the value of $1.4\div
1.5M_{\odot}$ for the maximum mass of a NS, inferred from precise neutron star
mass determinations [4], was considered the canonical limit. First neutron
star matter calculations with the inclusion of hyperons seemed to better agree
with this value compared to the case of pure nucleonic EoS, that predicts
relatively large maximum masses ($1.8\div 2.4M_{\odot}$) [5]. The recent
measurements of unusually high masses of the millisecond pulsars PSR
J1614-2230 ($1.97(4)M_{\odot}$) [6] and PSR J1903+0327 ($2.01(4)M_{\odot}$)
[7], rule out almost all these results, making uncertain the appearance of
strange baryons in high-density matter. However, in the last three years new
models compatible with the recent observations have been proposed, but many
inconsistency still remain. The solution of this problem, known as _hyperon
puzzle_ , is far to be understood.

The difficulty of correctly describe the effect of strange baryons in the
nuclear medium, is that one needs a precise solution of a many-body problem
for a very dense system with strong and complicated interactions which are
often poorly known.

The determination of a realistic interaction among hyperons and nucleons
capable to reconcile the terrestrial measurements on hypernuclei and the NS
observations is still an unsolved question. The amount of data available for
nucleon-nucleon scattering and binding energies is enough to build
satisfactory models of nuclear forces, either purely phenomenological or built
on the basis of an effective field theory. Same approaches have been used to
derive potentials for the hyperon-nucleon and hyperon-hyperon interaction, but
the accuracy of these models is far from that of the non strange counterparts.
The main reason of this is the lack of experimental information due the
impossibility to collect hyperon-neutron and hyperon-hyperon scattering data.
This implies that interaction models must be fitted mostly on binding energies
(and possibly excitations) of hypernuclei. In the last years several
measurements of the energy of hypernuclei became available. These can be used
to validate or to constrain the hyperon-nucleon interactions within the
framework of many-body systems. The ultimate goal is then to constrain these
forces by reproducing at best the experimental energies of hypernuclei from
light systems made of few particles up to heavier systems.

The method used to accurately solve the many-body Schrödinger equation
represents the second part of the problem. Accurate calculations are indeed
limited to very few nucleons. The exact Faddeev-Yakubovsky equation approach
has been applied up to four particle systems [8]. Few nucleon systems can be
accurately described by means of techniques based on shell models calculations
like the No-Core Shell Model [9], on the Hyperspherical Harmonics approach
[10, 11, 12, 13, 14] or on QuantumMonte Carlo methods, like the Variational
Monte Carlo [15, 16] or Green Function Monte Carlo [17, 18, 19, 20]. These
methods have been proven to solve the nuclear Schrödinger equation in good
agreement with the Faddeev-Yakubovsky method [21]. For heavier nuclei,
Correlated Basis Function theory [22], Cluster Variational Monte Carlo [23,
24] and Coupled Cluster Expansion [25, 26] are typically adopted. In addition,
the class of method which includes the Brueckner-Goldstone [27] and the
Hartree-Fock [28] algorithms is widely used, also for nuclear matter
calculations. The drawback of these many-body methods is that they modify the
original Hamiltonian to a more manageable form, often introducing uncontrolled
approximations in the algorithm. In absence of an exact method for solving the
many-body Schrödinger equation for a large number of nucleons, the derivation
of model interactions and their applicability in different regimes is subject
to an unpleasant degree of arbitrariness.

In this work we address the problem of the hyperon-nucleon interaction from a
Quantum Monte Carlo point of view. We discuss the application of the Auxiliary
Field Diffusion Monte Carlo (AFDMC) algorithm to study a non relativistic
Hamiltonian based on a phenomenological hyperon-nucleon interaction with
explicit two- and three-body components. The method was originally developed
for nuclear systems [29] and it has been successfully applied to the study of
nuclei [30, 31, 32], neutron drops [33, 20, 34], nuclear matter [35, 36] and
neutron matter [37, 38, 39, 40]. We have extended this ab-initio algorithm in
order to include the lightest of the strange baryons, the $\Lambda$ particle.
By studying the ground state properties of single and double $\Lambda$
hypernuclei, information about the employed microscopic hyperon-nucleon
interaction are deduced.

The main outcome of the study on finite strange systems is that only the
inclusion of explicit $\Lambda NN$ terms provides the necessary repulsion to
realistically describe the separation energy of a $\Lambda$ hyperon in
hypernuclei of intermediate masses [41, 42, 43]. The analysis of single
particle densities confirms the importance of the inclusion of the $\Lambda
NN$ contribution. On the ground of this observation, the three-body hyperon-
nucleon interaction has been studied in detail. By refitting the coefficients
in the potential, it has been possible to reproduce at the same time the
available experimental data accessible with AFDMC calculations in a medium-
heavy mass range [43]. Other details of the hypernuclear force, like the
charge symmetry breaking contribution and the effect of a $\Lambda\Lambda$
interaction, have been successfully analyzed. The AFDMC study of $\Lambda$
hypernuclei results thus in a realistic phenomenological hyperon-nucleon
interaction accurate in describing the ground state physics of medium-heavy
mass hypernuclei.

The large repulsive contribution induced by the three-body $\Lambda NN$ term,
makes very clear the fact that the lack of an accurate Hamiltonian might be
responsible for the unrealistic predictions of the EoS, that would tend to
rule out the appearance of strange baryons in high-density matter. We
speculate that the application of the developed hyperon-nucleon interaction to
the study of the homogeneous medium would lead to a stiffer EoS for the
$\Lambda$ neutron matter. This fact might eventually reconcile the physically
expected onset of hyperons in the inner core of a NS with the observed masses
of order $2M_{\odot}$.

First steps in this direction have been taken. The study of $\Lambda$ neutron
matter at fixed $\Lambda$ fraction shows that the repulsive nature of the
three-body hyperon-nucleon interaction is still active and relevant at
densities larger than the saturation density. The density threshold for the
appearance of $\Lambda$ hyperons has then been derived and the EoS has been
computed. Very preliminary results suggest a rather stiff EoS even in the
presence of hyperons, implying a maximum mass above the observational limit.
The study of hypermatter is still work in progress.

The present work is organized as follows:

Chapter 1:

a general overview about strangeness in nuclear systems, from hypernuclei to
neutron stars, is reported with reference to the terrestrial experiments and
astronomical observations.

Chapter 2:

a description of nuclear and hypernuclear non-relativistic Hamiltonians is
presented, with particular attention to the hyperon-nucleon sector in the two-
and three-body channels.

Chapter 3:

the Auxiliary Field Diffusion Monte Carlo method is discussed in its original
form for nuclear systems and in the newly developed version with the inclusion
of strange degrees of freedom, both for finite and infinite systems.

Chapter 4:

the analysis and set up of a realistic hyperon-nucleon interaction are
reported in connection with the AFDMC results for the hyperon separation
energy. Qualitative information are also deduced from single particle
densities and root mean square radii for single and double $\Lambda$
hypernuclei.

Chapter 5:

using the interaction developed for finite strange systems, first Quantum
Monte Carlo calculations on $\Lambda$ neutron matter are presented and the
implications of the obtained results for the properties of neutron stars are
explored.

Chapter 6:

the achievements of this work are finally summarized and future perspective
are discussed.

Empty page

## Chapter 1 Strangeness in nuclear systems

Hyperons are baryons containing one or more strange quarks. They have masses
larger than nucleons and lifetimes characteristic of the weak decay. The
$\Lambda$ and $\Omega$ hyperons belong to an isospin singlet, the $\Sigma$s to
an isospin triplet and the $\Xi$ particles to an isospin doublet. In Tab. 1.1
we report the list of hyperons (excluding resonances and unnatural parity
states [44]), with their main properties. The isospin doublet of nucleons is
also shown for comparison.

Baryon | qqq | $S$ | $I$ |  $m$ [MeV] |  $\tau$ [$10^{-10}$ s] | Decay mode
---|---|---|---|---|---|---
$p$ | uud | $0$ | $\displaystyle\frac{1}{2}$ | $938.272\,05(2)$ |  $\sim 10^{32}$ y | many
$n$ | udd | $939.565\,38(2)$ | 808(1) s | $p\,e\,\bar{\nu}_{e}$
$\Lambda$ | uds | $-1$ | 0 | $1115.683(6)$ | $2.63(2)$ | $p\,\pi^{-},n\,\pi^{0}$
$\Sigma^{+}$ | uus | $-1$ | 1 | $1189.37(7)$ | $0.802(3)$ | $p\,\pi^{0},n\,\pi^{+}$
$\Sigma^{0}$ | uds | $1192.64(2)$ |  $7.4(7)$$\times 10^{-10}$ | $\Lambda\,\gamma$
$\Sigma^{-}$ | dds | $1197.45(3)$ | $1.48(1)$ | $n\,\pi^{-}$
$\Xi^{0}$ | uss | $-2$ | $\displaystyle\frac{1}{2}$ | $1314.9(2)$ | $2.90(9)$ | $\Lambda\,\pi^{0}$
$\Xi^{-}$ | dss | $1321.71(7)$ | $1.64(2)$ | $\Lambda\,\pi^{-}$
$\Omega^{-}$ | sss | $-3$ | 0 | $1672.5(3)$ | $0.82(1)$ | $\Lambda\,K^{-},\Xi^{0}\,\pi^{-},\Xi^{-}\,\pi^{0}$
Table 1.1: Nucleon and hyperon properties: quark components, strangeness,
isospin, mass, mean life and principal decay modes [44].

In the non strange nuclear sector many information are available for nucleon-
nucleon scattering. The Nijmegen $NN$ scattering database [45, 46] includes
1787 $pp$ and 2514 $np$ data in the range $0\div 350$ MeV. Due to the
instability of hyperons in the vacuum and the impossibility to collect
hyperon-neutron and hyperon-hyperon scattering data, the available information
in the strange nuclear sector are instead very limited. Although many events
have been reported both in the low and high energy regimes [47], the standard
set employed in the modern hyperon-nucleon interactions (see for example Ref.
[48]) comprises 35 selected $\Lambda p$ low energy scattering data [49] and
some $\Lambda N$ and $\Sigma N$ data at higher energies [50]. In addition
there are the recently measured $\Sigma^{+}p$ cross sections of the KEK-PS
E289 experiment [51], for a total of 52 $YN$ scattering data.

The very limited experimental possibilities of exploring hyperon-nucleon and
hyperon-hyperon interactions in elementary scattering experiments, makes the
detailed study of hypernuclei essential to understand the physics in the
strange sector. In the next, we will present a summary of the available
hypernuclei experimental data. These information are the key ingredient to
develop realistic hyperon-nucleon and hyperon-hyperon interactions, as
described in the next chapters. The theoretical evidence of the appearance of
hyperons in the core of a NS and the problem of the hyperon puzzle will then
be discussed, following the results of many-body calculations for the
available models of hypermatter.

###  1.1 Hyperons in finite nuclei

In high-energy nuclear reactions strange hadrons are produced abundantly, and
they are strongly involved in the reaction process. When hyperons are captured
by nuclei, hypernuclei are formed, which can live long enough in comparison
with nuclear reaction times. Extensive efforts have been devoted to the study
of hypernuclei. Among many strange nuclear systems, the single $\Lambda$
hypernucleus is the most investigated one [52].

The history of hypernuclear experimental research (see Refs. [53, 54, 52] for
a complete review) celebrates this year the sixtieth anniversary, since the
publication of the discovery of hypernuclei by Danysz and Pniewski in 1953
[55]. Their first event was an example of ${}^{3}_{\Lambda}$H decaying via

${}^{3}_{\Lambda}\text{H}\longrightarrow\,^{3}\text{He}+\pi^{-}\;,$ (1.1)

confirming that the bound particle was a $\Lambda$ hyperon. The event was
observed in an emulsion stack as a consequence of nuclear multifragmentation
induced by cosmic rays. This first evidence opened the study of light
$\Lambda$ hypernuclei ($A<16$) by emulsion experiments, by means of cosmic ray
observations at the beginning and then through proton and pion beams, although
the production rates were low and there was much background. In the early
70’s, the advent of kaon beam at CERN and later at Brookhaven National
Laboratory (BNL), opened the possibility of spectroscopic studies of
hypernuclei, including excited states, by means of the $(K^{-},\pi^{-})$
reaction (see Fig. 1.1). A third stage, which featured the use of the
$(\pi^{+},K^{+})$ reaction, began in the mid 1980’s at the Alternating
Gradient Synchrotron (AGS) of BNL first, and then at the proton synchrotron
(PS) of the High Energy Accelerator Organization (KEK) in Japan. Here, the
superconducting kaon spectrometer (SKS) played a key role in exploring
$\Lambda$ hypernuclear spectroscopy by the $(\pi^{+},K^{+})$ reaction.
$\gamma$-ray spectroscopy developed reaching unprecedented resolution through
the use of a germanium detector array, the Hyperball, and the high quality and
high intensity electron beams available at the Thomas Jefferson National
Accelerator Facility (JLab). This permitted the first successful
$(e,e^{\prime}K^{+})$ hypernuclear spectroscopy measurement (an historical
review of hypernuclear spectroscopy with electron beams can be found in Ref.
[56]. The detailed analysis of $\Lambda$ hypernuclei spectroscopy is reported
in Ref. [52]).

Figure 1.1: Schematic presentation of three strangeness producing reactions
used in the study of $\Lambda$ hypernuclei.

With the development of new facilities, like the japanese J-PARC (Proton
Accelerator Research Complex), other reaction channels for the production of
neutron rich $\Lambda$ hypernuclei became available. The candidates are the
single charge exchange (SCX) reactions $(K^{-},\pi^{0})$ and
$(\pi^{-},K^{0})$, and double charge exchange (DCX) reactions
$(\pi^{-},K^{+})$ and $(K^{-},\pi^{+})$. Fig. 1.2 nicely illustrates the
complementarity of the various production mechanisms and thus the need to
study hypernuclei with different reactions. Moreover, during the last 20 years
of research, great progress has been made in the investigation of
multifragmentation reactions associated with heavy ion collisions (see for
instance [57] and reference therein). This gives the opportunity to apply the
same reactions for the production of hypernuclei too [58, 59]. On the other
hand, it was noticed that the absorption of hyperons in spectator regions of
peripheral relativistic ion collisions is a promising way to produce
hypernuclei [60, 61]. Also, central collisions of relativistic heavy ions can
lead to the production of light hypernuclei [62]. Recent experiments have
confirmed observations of hypernuclei in such reactions, in both peripheral
[63, 64] and central collisions [65].

Figure 1.2: $\Lambda$ hypernuclei accessible by experiments for different
production channels. The boundaries at the neutron and proton rich side mark
the predicted drip lines by a nuclear mass formula extended to strange nuclei.
Figure taken from Ref. [66].

At the time of writing, many laboratories are making extensive efforts in the
study of $\Lambda$ hypernuclei. The status of the art together with future
prospects can be found in Refs. [67, 68, 69] for the J-PARC facility and in
Ref. [70] for the ALICE (A Large Ion Collider Experiment) experiment at the
LHC. Ref. [71] reports the status of the JLab’s Hall A program. In Ref. [72]
future prospects for the the PANDA (antiProton ANihilation at DArmstadt)
project at FAIR (Facility for Antiproton ad Ion Research) and the hypernuclear
experiments using the KAOS spectrometer at MAMI (Mainz Microtron) can be
found. Last results from the FINUDA (FIsica NUcleare a DA$\Phi$NE)
collaboration at DA$\Phi$NE, Italy, are reported in Ref. [73]. Recent interest
has been also focused on the $S=-2$ sector with the study of double $\Lambda$
hypernuclei [74] and the $S=-3$ sector with the search for $\Omega$
hypernuclei [75].

So far, there is no evidence for $\Lambda p$ and ${}^{3}_{\Lambda}$He bound
states. Only very recently the possible evidence of the three-body system
$\Lambda nn$ has been reported [76]. The first well established weakly bound
systems is ${}^{3}_{\Lambda}$H, with hyperon separation energy $B_{\Lambda}$
(the energy difference between the $A-1$ nucleus and the $A$ hypernucleus,
being $A$ the total number of baryons) of $0.13(5)$ MeV [77]. Besides the very
old experimental results [77, 78, 79], several measurements of single
$\Lambda$ hypernuclei became available in the last years trough the many
techniques described above [80, 81, 82, 83, 84, 85, 86, 73]. The update
determination of the lifetime of ${}_{\Lambda}^{3}$H and ${}_{\Lambda}^{4}$H
has been recently reported [87] and new proposals for the search of exotic
$\Lambda$ hypernuclei are constantly discussed (see for example the search for
${}_{\Lambda}^{9}$He [88]). One of the results of this investigation is the
compilation of the $\Lambda$ hypernuclear chart reported in Fig. 1.3. Although
the extensive experimental studies in the $S=-1$ strangeness sector, the
availability of information for hypernuclei is still far from the abundance of
data for the non strange sector.

Figure 1.3: $\Lambda$ hypernuclear chart presented at the XI International
Conference on Hypernuclear and Strange Particle Physics (HYP2012), October
2012, Spain. The figure has been updated from Ref. [52].

It is interesting to observe that with the increase of $A$, there is an
orderly increase of $B_{\Lambda}$ with the number of particles, of the order
of 1 MeV/nucleon (see Tab. 4.3 or the mentioned experimental references). Many
stable hypernuclei with unstable cores appears, as for example
${}^{6}_{\Lambda}$He, ${}^{8}_{\Lambda}$He, ${}^{7}_{\Lambda}$Be and
${}^{9}_{\Lambda}$Be. These evidences testify that the presence of a $\Lambda$
particle inside a nucleus has a glue like effect, increasing the binding
energy and stability of the system. This should be reflected by the attractive
behavior of the $\Lambda$-nucleon interaction, at least in the low density
regime of hypernuclei.

For $\Sigma$ hypernuclei, the situation is quite different. Up to now, only
one bound $\Sigma$ hypernucleus, ${}^{4}_{\Sigma}$He, was detected [89],
despite extensive searches. The analysis of experimental data suggests a
dominant $\Sigma$-nucleus repulsion inside the nuclear surface and a weak
attraction outside the nucleus. In the case of $\Xi$ hypernuclei, although
there is no definitive data for any $\Xi$ hypernucleus at present, several
experimental results suggest that $\Xi$-nucleus interactions are weakly
attractive [90]. No experimental indication exists for $\Omega$ hypernuclei.
It is a challenge to naturally explain the net attraction in $\Lambda$\- and
$\Xi$-nucleus potentials and at the same time the dominant repulsion in
$\Sigma$-nucleus potentials.

In addition to single hyperon nuclei, the binding energies of few double
$\Lambda$ hypernuclei (${}^{\;\;\,6}_{\Lambda\Lambda}$He [91, 92, 93],
${}^{\;10}_{\Lambda\Lambda}$Be, ${}^{\;12}_{\Lambda\Lambda}$Be and
${}^{\;12}_{\Lambda\Lambda}$Be [94, 92], ${}^{\;13}_{\Lambda\Lambda}$B [92])
have been measured. The indication is that of a weakly attractive
$\Lambda\Lambda$ interaction, which reinforces the glue like role of $\Lambda$
hyperons inside nuclei.

From the presented picture it is clear that experimental hypernuclear physics
has become a very active field of research. However there is still lack of
information, even in the most investigated sector of $\Lambda$ hypernuclei.
Due to the technical difficulties in performing scattering experiments
involving hyperons and nucleons, the present main goal is the extension of the
$\Lambda$ hypernuclear chart to the proton and neutron drip lines and for
heavier systems. Parallel studies on $\Sigma$, $\Xi$ and double $\Lambda$
hypernuclei have been and will we be funded in order to try to complete the
scheme. This will hopefully provide the necessary information for the
development of realistic hyperon-nucleon and hyperon-hyperon interactions.

###  1.2 Hyperons in neutron stars

The matter in the outer core of a NS is supposed to be composed by a
degenerate gas of neutrons, protons, electrons and muons, the $npe\mu$ matter,
under $\beta$ equilibrium. Given the energy density

$\displaystyle\mathcal{E}(\rho_{n},\rho_{p},\rho_{e},\rho_{\mu})=\mathcal{E}_{N}(\rho_{n},\rho_{p})+\mathcal{E}_{e}(\rho_{e})+\mathcal{E}_{\mu}(\rho_{\mu})\;,$
(1.2)

where $\mathcal{E}_{N}$ is the nucleon contribution, the equilibrium condition
at a given baryon density $\rho_{b}$ corresponds to the minimum of
$\mathcal{E}$ under the constraints

fixed baryon density: $\displaystyle\rho_{n}+\rho_{p}-\rho_{b}=0\;,$ (1.3a)
electrical neutrality: $\displaystyle\rho_{e}+\rho_{\mu}-\rho_{p}=0\;.$ (1.3b)

The result is the set of conditions

$\displaystyle\mu_{n}$ $\displaystyle=\mu_{p}+\mu_{e}\;,$ (1.4a)
$\displaystyle\mu_{\mu}$ $\displaystyle=\mu_{e}\;,$ (1.4b)

where $\mu_{j}=\partial\mathcal{E}/\partial\rho_{j}$ with $j=n,p,e,\mu$ are
the chemical potentials. These relations express the equilibrium with respect
to the weak interaction processes

$\displaystyle\begin{array}[]{rclrcl}n&\longrightarrow&p+e+\bar{\nu}_{e}\;,&p+e&\longrightarrow&n+\nu_{e}\;,\\\\[5.0pt]
n&\longrightarrow&p+\mu+\bar{\nu}_{\mu}\;,&p+\mu&\longrightarrow&n+\nu_{\mu}\;.\end{array}$
(1.7)

(Neutrino do not affect the matter thermodynamics so their chemical potential
is set to zero). Eqs. (1.4) supplemented by the constraints (1.3) form a
closed system which determines the equilibrium composition of the $npe\mu$
matter. Once the equilibrium is set, the energy and pressure as a function of
the baryon density can be derived and thus the EoS is obtained.

Given the EoS, the structure of a non rotating NS can be fully determined by
solving the Tolman-Oppenheimer-Volkoff (TOV) equations [95, 96]

$\displaystyle\frac{dP(r)}{dr}$
$\displaystyle=-G\frac{\Bigl{[}\mathcal{E}(r)+P(r)\Bigr{]}\Bigl{[}m(r)+4\pi
r^{3}P(r)\Bigr{]}}{r^{2}\Bigl{[}1-\frac{2Gm(r)}{r}\Bigr{]}}\;,$ (1.8a)
$\displaystyle\frac{dm(r)}{dr}$ $\displaystyle=4\pi r^{2}\mathcal{E}(r)\;,$
(1.8b)

which describe the hydrostatic equilibrium of a static spherically symmetric
star. $\mathcal{E}(r)$ and $P(r)$ are the energy density and the pressure of
the matter, $m(r)$ is the gravitational mass enclosed within a radius $r$, and
$G$ is the Gravitational constant. In the stellar interior $P>0$ and
$dP/dr<0$. The condition $P(R)=0$ fixes the stellar radius $R$. Outside the
star for $r>R$, we have $P=0$ and $\mathcal{E}=0$. Eq. (1.8b) gives thus
$m(r>R)=M=const$, which is total gravitational mass. Starting with a central
energy density $\mathcal{E}_{c}=\mathcal{E}(r=0)$ and using the above
conditions, the TOV equations can be numerically solved and the mass-radius
relation $M=M(R)$ is obtained. It can be shown [1], that the relativistic
corrections to the Newtonian law $dP(r)/dr=-Gm\mathcal{E}(r)/r^{2}$ included
in Eq. (1.8a) give an upperbound to the $M(R)$ relation, i.e. there exists a
maximum mass for a NS in hydrostatic equilibrium. It is important to note
that, given the EoS, the mass-radius relation is univocally determined. Any
modification made on the EoS will lead to a change in the $M(R)$ curve and
thus in the allowed maximum mass.

For $\rho_{b}\gtrsim 2\rho_{0}$, the inner core is thought to have the same
$npe\mu$ composition of the outer core. However, since at high densities the
nucleon gas will be highly degenerate, hyperons with energies lower than a
threshold value will become stable, because the nucleon arising from their
decay cannot find a place in phase space in accordance to the Pauli principle
[2]. Thus, beyond a density threshold we have to take into account the
contribution of hyperons to the $\beta$ equilibrium. Eq. (1.2) becomes a
function of $\rho_{b}$ (baryons: nucleons and hyperons) and $\rho_{l}$
(leptons: electrons and muons). Given the baryon density and imposing
electrical neutrality conditions, the equilibrium equations now read:

$\displaystyle Q_{b}=-1\,:$ $\displaystyle\mu_{b^{-}}$
$\displaystyle=\mu_{n}+\mu_{e}$ $\displaystyle\Rightarrow$
$\displaystyle\mu_{\Omega^{-}}$
$\displaystyle=\mu_{\Xi^{-}}=\mu_{\Sigma^{-}}=\mu_{n}+\mu_{e}\;,$ (1.9a)
$\displaystyle Q_{b}=\phantom{+}0\,:$ $\displaystyle\mu_{b^{0}}$
$\displaystyle=\mu_{n}$ $\displaystyle\Rightarrow$
$\displaystyle\mu_{\Xi^{0}}$
$\displaystyle=\mu_{\Sigma^{0}}=\mu_{\Lambda}=\mu_{n}\;,$ (1.9b)
$\displaystyle Q_{b}=+1\,:$ $\displaystyle\mu_{b^{+}}$
$\displaystyle=\mu_{n}-\mu_{e}$ $\displaystyle\Rightarrow$
$\displaystyle\mu_{\Sigma^{+}}$ $\displaystyle=\mu_{p}=\mu_{n}-\mu_{e}\;,$
(1.9c)

where $Q_{b}$ is the electric charge of a baryon. As soon as the neutron
chemical potential becomes sufficiently large, energetic neutrons can decay
via weak strangeness nonconserving reactions into $\Lambda$ hyperons, leading
to a $\Lambda$ Fermi sea.

We can derive the hyperons threshold densities $\rho_{Y}$ by calculating the
minimum increase of the energy of the matter produced by adding a single
strange particle at a fixed pressure. This can be done by considering the
energy of the matter with an admixture of given hyperons and by calculating
numerically the limit of the derivative

$\displaystyle\lim_{\rho_{Y}\rightarrow
0}\,\frac{\partial\mathcal{E}}{\partial\rho_{Y}}\bigg{|}_{eq}\\!=\mu_{Y}^{0}\;.$
(1.10)

Consider for example the lightest $\Lambda$ hyperon. As long as
$\mu_{\Lambda}^{0}>\mu_{n}$, the strange baryon cannot survive because the
system will lower its energy via an exothermic reaction
$\Lambda+N\longrightarrow n+N$. However, $\mu_{n}$ increases with growing
$\rho_{b}$ and the functions $\mu_{\Lambda}^{0}(\rho_{b})$ and
$\mu_{n}^{0}(\rho_{b})$ intersect at some $\rho_{b}=\rho_{\Lambda}^{th}$ (the
left panel in Fig. 1.4). For $\rho_{b}>\rho_{\Lambda}^{th}$ the $\Lambda$
hyperons become stable in dense matter because their decay is blocked by the
Pauli principle.

Figure 1.4: Threshold chemical potentials of neutral hyperons and neutron
(left panel), and of negatively charged hyperons and the sum $\mu_{n}+\mu_{e}$
(right panel) versus baryon density. Vertical dotted lines mark the thresholds
for the creation of new hyperons. Dashed lines show the minimum chemical
potential $\mu_{Y}^{0}$ of unstable hyperons before the thresholds. Figure
taken from Ref. [1].

Although the $\Lambda$ particle is the lightest among hyperons, one expects
the $\Sigma^{-}$ to appear via

$\displaystyle n+e^{-}\longrightarrow\Sigma^{-}+\nu_{e}$ (1.11)

at densities lower than the $\Lambda$ threshold, even thought the $\Sigma^{-}$
is more massive. This is because the negatively charged hyperons appear in the
ground state of matter when their masses equal $\mu_{n}+\mu_{e}$, while the
neutral hyperon $\Lambda$ appears when its mass equals $\mu_{n}$. Since the
electron chemical potential in matter is typically larger (ultrarelativistic
degenerate electrons $\mu_{e}\sim E_{F_{e}}\sim\hbar
c(3\pi^{2}\rho_{e})^{1/3}>120\leavevmode\nobreak\
\text{MeV\leavevmode\nobreak\ for\leavevmode\nobreak\ }\rho_{e}\sim
5\%\rho_{0}$) than the mass difference
$m_{\Sigma^{-}}-m_{\Lambda}=81.76\leavevmode\nobreak\ \mbox{MeV}$, the
$\Sigma^{-}$ will appear at lower densities. However, in typical neutron
matter calculations with the inclusion of strange degrees of freedom, only
$\Lambda$, $\Sigma^{0}$ and $\Xi^{0}$ hyperons are taken into account due to
charge conservation.

The formation of hyperons softens the EoS because high energy neutrons are
replaced by more massive low energy hyperons which can be accommodated in
lower momentum states. There is thus a decrease in the kinetic energy that
produces lower pressure. The softening of the EoS of the inner core of a NS
induced by the presence of hyperons is generic effect. However, its magnitude
is strongly model dependent.

Calculations based on the extension to the hyperonic sector of the Hartree-
Fock (HF) [97, 98] and Brueckner-Hartree-Fock (BHF) [99, 100] methods, do all
agree that the appearance of hyperons around $2\div 3\rho_{0}$ leads to a
strong softening of the EoS. Consequently, the expected maximum mass is
largely reduced, as shown for instance in Fig. 1.5 and Fig. 1.6. The addition
of the hyperon-nucleon force to the pure nucleonic Hamiltonian, lowers the
maximum mass of a value between $0.4M_{\odot}$ and more than $1M_{\odot}$.
From the pure nucleonic case of $M_{\max}>1.8M_{\odot}$, the limit for
hypernuclear matter is thus reduced to the range
$1.4M_{\odot}<M_{\max}<1.6M_{\odot}$. These results, although compatible with
the canonical limit of $1.4\div 1.5M_{\odot}$, cannot be consistent with the
recent observations of $2M_{\odot}$ millisecond pulsars [6, 7].

Figure 1.5: Mass-radius and mass-central density relations for different NS
EoS obtained in Brueckner-Hartree-Fock calculations of hypernuclear matter.
V18+TBF and V18+UIX’ refer to purely nuclear matter EoS built starting from
two- and three-body nucleon-nucleon potentials (see § 2.1). The other curves
are obtained adding two different hyperon-nucleon forces among the Nijmegen
models to the previous nucleonic EoS. For more details see the original paper
[99]. Figure 1.6: Neutron star mass as a function of the circumferential
radius. QMC700 and MC _i_ -H(F)/N refer to EoS based on quark-meson coupling
model and chiral model in the Hartee(Fock) approximation without hyperons. In
the MC _i_ -H(F)/NY models also hyperons are taken into account. The canonical
maximum mass limit of $\sim 1.45M_{\odot}$ and the mass of the two heavy
millisecond pulsars PSR J1903+0327 ($1.67(2)M_{\odot}$) and PSR J1614-2230
($1.97(4)M_{\odot}$) are shown. Details on the potentials and method adopted
can be found in Ref. [98].

It is interesting to note that the hyperonic $M_{\max}$ weakly depends on the
details of the employed nucleon-nucleon interaction and even less on the
hypernuclear forces. In Ref. [97] the interaction used for the nuclear sector
is an analytic parametrization fitted to energy of symmetric matter obtained
from variational calculations with the Argonne V18 nucleon-nucleon interaction
(see § 2.1) including three-body forces and relativistic boost corrections.
Refs. [99] and [100] adopted the bare $NN$ Argonne V18 supplemented with
explicit three-nucleon forces or phenomenological density-dependent contact
terms that account for the effect of nucleonic and hyperonic three-body
interactions. The hypernuclear forces employed in these work belong to the
class of Nijmegen potentials (see § 2). Finally, in Ref. [98] chiral
Lagrangian and quark-meson coupling models of hyperon matter have been
employed. Despite the differences in the potentials used in the strange and
non strange sectors, the outcomes of these works give the same qualitative and
almost quantitative picture about the reduction of $M_{\max}$ due to the
inclusions of strange baryons. Therefore, the (B)HF results seem to be rather
robust and thus, many doubts arise about the real appearance of hyperons in
the inner core of NSs.

Other approaches, such as relativistic Hartree-Fock [101, 102, 103], standard,
density-dependent and nonlinear Relativistic Mean Field models [104, 105, 106,
107, 108] and Relativistic Density Functional Theory with density-dependent
couplings [109], indicate much weaker effects as a consequence of the presence
of strange baryons in the core of NSs, as shown for example in Fig. 1.7 and
Fig. 1.8. In all these works, it was possible to find a description of
hypernuclear matter, within the models analyzed, that produces stiff EoS,
supporting a $2M_{\odot}$ neutron star. Same conclusion has been reported in
Ref. [110] where the EoS of matter including hyperons and deconfined quark
matter has been constructed on the basis of relativistic mean-field nuclear
functional at low densities and effective Nambu-Jona-Lasinio model of quark
matter. The results of this class of calculations seem to reconcile the onset
of hyperons in the inner core of a NS with the observed masses of order
$2M_{\odot}$.

Figure 1.7: Neutron star mass-radius relations in Hartree (left panel) and
Hartree-Fock (right panel) calculations. CQMC, QMC and QHD+NL denote the
chiral quark-meson coupling, quark-meson coupling and non linear quantum
hadrodynamics employed potentials, with (npY) and without hyperons (np). For
details see Ref. [101]. Figure 1.8: Stellar mass versus circumferential radius
in non linear relativistic mean field model. The purely nucleon case is
denoted with N, the nucleon+hyperon case with NH. In the inset, the effect of
rotation at $f=317$ Hz on the mass-radius relation near $M_{\max}$. The dashed
region refers to the mass of the pulsar PSR J1614-2230. All the details are
reported in Ref. [104].

This inconsistency among different calculations and between the theoretical
results and the observational constraints, at present is still an open
question. For example, given the theoretical evidence about the appearance of
hyperons in the inner core of a NS, the results of all available (B)HF
calculations seem to be in contradiction with the picture drawn by the
relativistic mean field models. On one hand there should be uncontrolled
approximations on the method used to solve the many-body Hamiltonian. On the
other hand the employed hypernuclear interactions might not be accurate enough
in describing the physics of the infinite nuclear medium with strange degrees
of freedom. For instance, as reported in Refs. [111, 106], one of the possible
solutions to improve the hyperon-nucleon interactions might be the inclusion
of explicit three-body forces in the models. These should involve one or more
hyperons (i.e., hyperon-nucleon-nucleon, hyperon-hyperon-nucleon or hyperon-
hyperon-hyperon interactions) and they could eventually provide the additional
repulsion needed to make the EoS stiffer and, therefore the maximum mass
compatible with the current observational limits. On the grounds of this
observation, we decided to revisit the problem focusing on a systematic
construction of a realistic, though phenomenological hyperon-nucleon
interaction with explicit two- and three-body components (§ 2) by means of
Quantum Monte Carlo calculations (§ 3).

Empty page

## Chapter 2 Hamiltonians

The properties of nuclear systems arise from the interactions between the
individual constituents. In order to understand these properties, the starting
point is the determination of the Hamiltonian to be used in the description of
such systems. In principle the nuclear Hamiltonian should be directly derived
from Quantum Chromodynamics (QCD). Many efforts have been done in the last
years [112, 113, 114], but this goal is still far to be achieved.

The problem with such derivation is that QCD is non perturbative in the low-
temperature regime characteristic of nuclear physics, which makes direct
solutions very difficult. Moving from the real theory to effective models, the
structure of a nuclear Hamiltonian can be determined phenomenologically and
then fitted to exactly reproduce the properties of few-nucleon systems. In
this picture, the degrees of freedom are the baryons, which are considered as
non relativistic point-like particles interacting by means of phenomenological
potentials. These potentials describe both short and the long range
interactions, typically via one-boson and two-meson exchanges, and they have
been fitted to exactly reproduce the properties of few-nucleon systems [115].
In more details, different two-body phenomenological forms have been proposed
and fitted on the nucleon-nucleon ($NN$) scattering data of the Nijmegen
database [45, 46] with a $\chi^{2}/N_{data}\simeq 1$. The more diffuse are the
Nijmegen models [116], the Argonne models [117, 118] and the CD-Bonn [119].
Although reproducing the $NN$ scattering data, all these two-nucleon
interactions underestimate the triton binding energy, suggesting that the
contribution of a three-nucleon ($NNN$) interaction (TNI) is essential to
reproduce the physics of nuclei. The TNI is mainly attributed to the
possibility of nucleon excitation in a $\Delta$ resonance and it can be
written as different effective three-nucleon interactions which have been
fitted on light nuclei [120, 121] and on saturation properties of nuclear
matter [122]. The TNIs typically depend on the choice of the two-body $NN$
potential [123], but the final result with the total Hamiltonian should be
independent of the choice.

A different approach to the problem is the realization that low-energy QCD is
equivalent to an Effective Field Theory (EFT) which allows for a perturbative
expansion that is known as chiral perturbation theory. In the last years
modern nucleon-nucleon interaction directly derived from Chiral Effective
Field Theory ($\chi$-EFT) have been proposed, at next-to-next-to-next-to-
leading order (N3LO) in the chiral expansion [124, 125] and recently at
optimized next-to-next-to-leading order (N2LO) [126] (see Ref. [127] for a
complete review). All these potentials are able to reproduce the Nijmegen
phase shifts with $\chi/N_{data}^{2}\simeq 1$. TNIs enter naturally at N2LO in
this scheme, and they play again a pivotal role in nuclear structure
calculations [128]. The contributions of TNIs at N3LO have also been worked
out [129, 130, 131]. The $\chi$-EFT interactions are typically developed in
momentum space, preventing their straightforward application within the
Quantum Monte Carlo (QMC) framework. However, a local version of the
$\chi$-EFT potentials in coordinate space up to N2LO has been very recently
proposed and employed in QMC calculations [132].

Nuclear phenomenological Hamiltonians have been widely used to study finite
and infinite nuclear systems within different approaches. From now on, we will
focus on the Argonne $NN$ potentials and the corresponding TNIs, the Urbana IX
(UIX) [122] and the modern Illinois (ILx) [121] forms. These potentials have
been used to study nuclei, neutron drops, neutron and nuclear matter in
Quantum Monte Carlo (QMC) calculations, such as Variational Monte Carlo (VMC)
[15, 16, 24], Green Function Monte Carlo (GFMC) [133, 118, 134, 17, 135, 136,
19, 20] and Auxiliary Field Diffusion Monte Carlo (AFDMC) [37, 33, 30, 31, 20,
35, 38]. Same bare interactions have been also employed in the Fermi Hyper-
Netted Chain (FHNC) approach [22, 137], both for nuclei and nuclear matter.
With a projection of the interaction onto the model space, these Hamiltonians
are used in Effective Interaction Hyperspherical Harmonics (EIHH) [10, 13] and
Non-Symmetrized Hyperspherical Harmonics (NSHH) [14] calculations. Finally,
same potentials can be also used in Brueckner Hartree Fock (BHF) [138], Shell-
Model (SM) [139], No-Core-Shell-Model (NCSM) [9] and Coupled Cluster (CC) [26]
calculations by means of appropriate techniques to handle the short-range
repulsion of the nucleon-nucleon force, such as Brueckner $G$-matrix approach
[140, 141], $V_{low-k}$ reduction [142, 143, 144], Unitary Correlation
Operator Method (UCOM) [145] or Similarity Renormalization Group (SRG)
evolution [146, 147]. The list of methods that can handle in a successful way
the Argonne+TNIs potentials demonstrates the versatility and reliability of
this class of phenomenological nuclear Hamiltonians.

Moving from the non-strange nuclear sector, where nucleons are the only
baryonic degrees of freedom, to the strange nuclear sector, where also
hyperons enter the game, the picture becomes much less clear. There exists
only a very limited amount of scattering data from which one could construct
high-quality hyperon-nucleon ($YN$) potentials. Data on hypernuclei binding
energies and hyperon separation energies are rather scarce and can only
partially complete the scheme.

After the pioneering work reported in Ref. [148], several models have been
proposed to describe the $YN$ interaction. The more diffuse are the Nijmegen
soft-core models (like NSC89 and NSC97x) [149, 150, 151, 152, 153, 154, 155]
and the Jülic potential (J04) [156, 157, 158]. A recent review of these
interactions, together with Hartree-Fock (HF) calculations have been published
by Ðapo  _et al._ in Ref. [159]. In the same framework, extended soft-core
Nijmegen potentials for strangeness $S=-2$ have been also developed [160,
161]. Very recently, the extended soft-core 08 (ESC08) model has been
completed, which represents the first unified theoretical framework involving
hyperon-nucleon, hyperon-hyperon ($YY$) and also nucleon-nucleon sectors [48].
This class of interaction has been used in different calculations for
hypernuclei [162, 163, 164, 165, 166, 159, 48] and hypermatter [159, 97, 99,
100] within different methods, but the existing data do not constrain the
potentials sufficiently. For example, six different parameterizations of the
Nijmegen $YN$ potentials fit equally well the scattering data but produce very
different scattering lengths, as reported for instance in Ref. [152]. In
addition, these potentials are not found to yield the correct spectrum of
hypernuclear binding energies. For example, the study [166] of
${}^{4}_{\Lambda}$H and ${}^{4}_{\Lambda}$He that uses Nijmegen models, does
not predict all experimental separation energies. Similar conclusions for
single- and double-$\Lambda$ hypernuclei have also been drawn in a study
employing a different many-body technique [165]. Even the most recent ESC08
model produces some overbinding of single-$\Lambda$ hypernuclei and a weakly
repulsive incremental $\Lambda\Lambda$ energy [48], not consistent with the
observed weak $\Lambda\Lambda$ attraction in
${}^{\;\;\,6}_{\Lambda\Lambda}$He.

In analogy with the nucleon-nucleon sector, a $\chi$-EFT approach for the
hyperon-nucleon interaction has been also developed. The first attempt was
proposed by Polinder and collaborators in 2006 [167], resulting in a leading
order (LO) expansion. Only recently the picture has been improved going to
next-to-leading order (NLO) [168, 169, 170]. The $YN$ $\chi$-EFT model is
still far away from the theoretical accuracy obtained in the non-strange
sector, but it is any case good enough to describe the limited available $YN$
scattering data.

As an alternative, a cluster model with phenomenological interactions has been
proposed by Hiyama and collaborators to study light hypernuclei [171, 172,
173, 174, 175, 176]. Interesting results on $\Lambda$ hypernuclei have also
been obtained within a $\Lambda$-nucleus potential model, in which the need of
a functional with a more than linear density dependence was shown, suggesting
the importance of a many-body interaction [177]. While studying $s$-shell
hypernuclei, the $\Lambda N\rightarrow\Sigma N$ coupling as a three-body
$\Lambda NN$ force has been investigated by many authors [166, 178, 179, 180].
Having strong tensor dependence it is found to play an important role,
comparable to the TNI effect in non-strange nuclei.

Finally, starting in the 1980s, a class of Argonne-like interactions has been
developed by Bodmer, Usmani and Carlson on the grounds of quantum Monte Carlo
calculations to describe the $\Lambda$-nucleon force. These phenomenological
interactions are written in coordinates space and they include two- and three-
body hyperon-nucleon components, mainly coming from two-pion exchange
processes and shorter range effects. They have been used in different forms
mostly in variational Monte Carlo calculations for single $\Lambda$
hypernuclei (${}^{3}_{\Lambda}$H [181, 182], ${}^{4}_{\Lambda}$H and
${}^{4}_{\Lambda}$He [183, 181, 182, 184], ${}^{5}_{\Lambda}$He [181, 182,
185, 186, 184, 187, 188, 189], ${}^{9}_{\Lambda}$Be [190, 191],
${}^{13}_{\leavevmode\nobreak\ \Lambda}$C [190],
${}^{17}_{\leavevmode\nobreak\ \Lambda}$O [192, 185]), double $\Lambda$
hypernuclei (${}^{\;\;\,4}_{\Lambda\Lambda}$H,
${}^{\;\;\,5}_{\Lambda\Lambda}$H, ${}^{\;\;\,5}_{\Lambda\Lambda}$He [193] and
${}^{\;\;\,6}_{\Lambda\Lambda}$He [193, 194, 195]) and in the framework of
correlated basis function theory for $\Lambda$ hypernuclei [196], typically in
connection with the Argonne $NN$ potential.

Within the phenomenological interaction scheme, a generic nuclear system
including nucleons and hyperons, can be described by the non relativistic
phenomenological Hamiltonian

$H=H_{N}+H_{Y}+H_{YN}\;,$ (2.1)

where $H_{N}$ and $H_{Y}$ are the pure nucleonic and hyperonic Hamiltonians
and $H_{YN}$ represents the interaction Hamiltonian connecting the two
distinguishable types of baryon:

$\displaystyle H_{N}$
$\displaystyle=\frac{\hbar^{2}}{2m_{N}}\sum_{i}\nabla_{i}^{2}\;+\sum_{i<j}v_{ij}\;\,+\sum_{i<j<k}v_{ijk}\;\;\,+\,\ldots\;,$
(2.2) $\displaystyle H_{Y}$
$\displaystyle=\frac{\hbar^{2}}{2m_{\Lambda}}\sum_{\lambda}\nabla_{\lambda}^{2}\;+\sum_{\lambda<\mu}v_{\lambda\mu}\,+\sum_{\lambda<\mu<\nu}v_{\lambda\mu\nu}\;+\,\ldots\;,$
(2.3) $\displaystyle H_{YN}$ $\displaystyle=\sum_{\lambda i}v_{\lambda
i}\,+\sum_{\lambda,i<j}v_{\lambda ij}\,+\sum_{\lambda<\mu,i}v_{\lambda\mu
i}\,+\,\ldots\;.$ (2.4)

In this context, $A$ is the total number of baryons,
$A=\mathcal{N}_{N}+\mathcal{N}_{Y}$. Latin indices
$i,j,k=1,\ldots,\mathcal{N}_{N}$ label nucleons and Greek symbols
$\lambda,\mu,\nu=1,\ldots,\mathcal{N}_{Y}$ are used for the hyperons. The
Hamiltonians (2.2) and (2.3) contain the kinetic energy operator and two- and
three-body interactions for nucleons and hyperons separately. In principles
they could include higher order many-body forces that however are expected to
be less important. The Hamiltonian (2.4) describes the interaction between
nucleons and hyperons, and it involves two-body ($YN$) and three-body ($YNN$
and $YYN$) forces. At present there is no evidence for higher order terms in
the hyperon-nucleon sector.

As reported in the previous chapter, experimental data are mainly available
for $\Lambda p$ scattering and $\Lambda$ hypernuclei and present experimental
efforts are still mostly concentrated in the study of the $S=-1$ hypernuclear
sector. Information on heavier hyperon-nucleon scattering and on $\Sigma$ or
more exotic hypernuclei are very limited. For these reasons, from now on we
will focus on the phenomenological interactions involving just the $\Lambda$
hyperon. We adopt the class of Argonne-like $\Lambda$-nucleon interaction for
the strange sector and the nucleon-nucleon Argonne force with the
corresponding TNIs (UIX and ILx) for the non-strange sector. An effective
$\Lambda\Lambda$ interaction has been also employed.

###  2.1 Interactions: nucleons

We report the details of the $NN$ Argonne potential [117, 118] and the
corresponding TNIs, the Urbana IX (UIX) [122] and the Illinois (ILx) [121].
These interactions are written in coordinate space and they include different
range components coming from meson (mostly pion) exchange and phenomenological
higher order contributions.

####  2.1.1 Two-body $NN$ potential

The nucleon-nucleon potential Argonne V18 (AV18) [117] contains a complete
electromagnetic (EM) interaction and a strong interaction part which is
written as a sum of a long-range component $v_{ij}^{\pi}$ due to one-pion
exchange (OPE) and a phenomenological intermediate- and short-range part
$v_{ij}^{R}$ :

$v_{ij}=v_{ij}^{\pi}+v_{ij}^{R}\;.$ (2.5)

Ignoring isospin breaking terms, the long-range OPE is given by

$\displaystyle v_{ij}^{\pi}=\frac{f_{\pi
NN}^{2}}{4\pi}\frac{m_{\pi}}{3}\,X_{ij}\,\bm{\tau}_{i}\cdot\bm{\tau}_{j}\;,$
(2.6)

where $\tfrac{f_{\pi NN}^{2}}{4\pi}=0.075$ is the pion-nucleon coupling
constant [197] and

$\displaystyle
X_{ij}=Y_{\pi}(r_{ij})\,\bm{\sigma}_{i}\cdot\bm{\sigma}_{j}+T_{\pi}(r_{ij})\,S_{ij}\;.$
(2.7)

$\bm{\sigma}_{i}$ and $\bm{\tau}_{i}$ are Pauli matrices acting on the spin or
isospin of nucleons and $S_{ij}$ is the tensor operator

$\displaystyle
S_{ij}=3\left(\bm{\sigma}_{i}\cdot\hat{\bm{r}}_{ij}\right)\left(\bm{\sigma}_{j}\cdot\hat{\bm{r}}_{ij}\right)-\bm{\sigma}_{i}\cdot\bm{\sigma}_{j}\;.$
(2.8)

The pion radial functions associated with the spin-spin (Yukawa potential) and
tensor (OPE tensor potential) parts are

$\displaystyle Y_{\pi}(r)$
$\displaystyle=\frac{\operatorname{e}^{-\mu_{\pi}r}}{\mu_{\pi}r}\xi_{Y}(r)\;,$
(2.9) $\displaystyle T_{\pi}(r)$
$\displaystyle=\left[1+\frac{3}{\mu_{\pi}r}+\frac{3}{(\mu_{\pi}r)^{2}}\right]\frac{\operatorname{e}^{-\mu_{\pi}r}}{\mu_{\pi}r}\xi_{T}(r)\;,$
(2.10)

where $\mu_{\pi}$ is the pion reduced mass

$\displaystyle\mu_{\pi}=\frac{m_{\pi}}{\hbar}=\frac{1}{\hbar}\frac{m_{\pi^{0}}+2\,m_{\pi^{\pm}}}{3}\quad\quad\frac{1}{\mu_{\pi}}\simeq
1.4\leavevmode\nobreak\ \text{fm}\;,$ (2.11)

and $\xi_{Y}(r)$ and $\xi_{T}(r)$ are the short-range cutoff functions defined
by

$\displaystyle\xi_{Y}(r)=\xi_{T}^{1/2}(r)=1-\operatorname{e}^{-cr^{2}}\quad\quad
c=2.1\leavevmode\nobreak\ \text{fm}^{-2}\;.$ (2.12)

It is important to note that since $T_{\pi}(r)\gg Y_{\pi}(r)$ in the important
region where $r\lesssim 2$ fm, the OPE is dominated by the tensor part.

The remaining intermediate- and short-range part of the potential is expressed
as a sum of central, $L^{2}$, tensor, spin-orbit and quadratic spin-orbit
terms (respectively labelled as $c$, $l2$, $t$, $ls$, $ls2$) in different $S$,
$T$ and $T_{z}$ states:

$\displaystyle\\!v_{NN}^{R}=v_{NN}^{c}(r)+v_{NN}^{l2}(r)\bm{L}^{2}+v_{NN}^{t}(r)S_{12}+v_{NN}^{ls}(r)\bm{L}\\!\cdot\\!\bm{S}+v_{NN}^{ls2}(r)(\bm{L}\\!\cdot\\!\bm{S})^{2}\;,$
(2.13)

with the radial functions $v_{NN}^{k}(r)$ written in the general form

$\displaystyle
v_{NN}^{k}(r)=I_{NN}^{k}\,T_{\pi}^{2}(r)+\bigg{[}P_{NN}^{k}+(\mu_{\pi}r)\,Q_{NN}^{k}+(\mu_{\pi}r)^{2}\,R_{NN}^{k}\bigg{]}\,W(r)\;,$
(2.14)

where the $T_{\pi}^{2}(r)$ has the range of a two-pion exchange (TPE) force
and $W(r)$ is a Wood-Saxon function which provides the short-range core:

$\displaystyle
W(r)=\Bigl{(}1+\operatorname{e}^{\frac{r-\bar{r}}{a}}\Bigr{)}^{-1}\quad\quad\bar{r}=0.5\leavevmode\nobreak\
\text{fm},\quad a=0.2\leavevmode\nobreak\ \text{fm}\;.$ (2.15)

By imposing a regularization condition at the origin, it is possible to reduce
the number of free parameters by one for each $v_{NN}^{k}(r)$. All the
parameters in the $\xi(r)$ short-range cutoff functions as well as the other
phenomenological constants are fitted on the $NN$ Nijmegen scattering data
[45, 46].

The two-body nucleon potential described above can be projected from $S$, $T$,
$T_{z}$ states into an operator format with 18 terms

$\displaystyle v_{ij}=\sum_{p=1,18}v_{p}(r_{ij})\,\mathcal{O}_{ij}^{\,p}\;.$
(2.16)

The first 14 operators are charge independent and they are the ones included
in the Argonne V14 potential (AV14):

$\displaystyle\mathcal{O}_{ij}^{\,p=1,8}$
$\displaystyle=\Bigl{\\{}1,\bm{\sigma}_{i}\cdot\bm{\sigma}_{j},S_{ij},\bm{L}_{ij}\cdot\bm{S}_{ij}\Bigr{\\}}\otimes\Bigl{\\{}1,\bm{\tau}_{i}\cdot\bm{\tau}_{j}\Bigr{\\}}\;,$
(2.17) $\displaystyle\mathcal{O}_{ij}^{\,p=9,14}$
$\displaystyle=\Bigl{\\{}\bm{L}_{ij}^{2},\bm{L}_{ij}^{2}\;\bm{\sigma}_{i}\cdot\bm{\sigma}_{j},\left(\bm{L}_{ij}\cdot\bm{S}_{ij}\right)^{2}\Bigr{\\}}\otimes\Bigl{\\{}1,\bm{\tau}_{i}\cdot\bm{\tau}_{j}\Bigr{\\}}\;.$
(2.18)

The first eight terms give the higher contribution to the $NN$ interaction and
they are the standard ones required to fit $S$ and $P$ wave data in both
triplet and singlet isospin states. The first six of them come from the long-
range part of OPE and the last two depend on the velocity of nucleons and give
the spin-orbit contribution. In the above expressions, $\bm{L}_{ij}$ is the
relative angular momentum of a couple $ij$

$\displaystyle\bm{L}_{ij}=\frac{1}{2i}({\bf r}_{i}-{\bf
r}_{j})\times(\bm{\nabla}_{i}-\bm{\nabla}_{j})\;,$ (2.19)

and $\bm{S}_{ij}$ the total spin of the pair

$\displaystyle\bm{S}_{ij}=\frac{1}{2}(\bm{\sigma}_{i}+\bm{\sigma}_{j})\;.$
(2.20)

Operators from 9 to 14 are included to better describe the Nijmegen higher
partial waves phase shifts and the splitting of state with different $J$
values. However, the contribution of these operators is small compared to the
total potential energy.

The four last additional operators of the AV18 potential account for the
charge symmetry breaking effect, mainly due to the different masses of charged
and neutral pions, and they are given by

$\displaystyle\mathcal{O}_{ij}^{\,p=15,18}=\Bigl{\\{}T_{ij},(\bm{\sigma}_{i}\cdot\bm{\sigma}_{j})\,T_{ij},S_{ij}\,T_{ij},\tau_{i}^{z}+\tau_{j}^{z}\Bigr{\\}}\;,$
(2.21)

where $T_{ij}$ is the isotensor operator defined in analogy with $S_{ij}$ as

$\displaystyle
T_{ij}=3\,\tau_{i}^{z}\tau_{j}^{z}-\bm{\tau}_{i}\cdot\bm{\tau}_{j}\;.$ (2.22)

The contribution to the total energy given by these four operators is however
rather small.

In QMC calculations reduced versions of the original AV18 potential are often
employed. The most used one is the Argonne V8’ (AV8’) [118] that contains only
the first eight operators and it is not a simple truncation of AV18 but also a
reprojection, which preserves the isoscalar part in all $S$ and $P$ partial
waves as well as in the ${}^{3}D_{1}$ wave and its coupling to ${}^{3}S_{1}$.
AV8’ is about $0.2\div 0.3$ MeV per nucleon more attractive than Argonne V18
in light nuclei [121, 118, 198], but its contribution is very similar to AV18
in neutron drops, where the difference is about 0.06 MeV per neutron [121].
Other common solutions are the Argonne V6’ (AV6’) and V4’ (AV4’) potentials
[118]. AV6’ is obtained by deleting the spin-orbit terms from AV8’ and
adjusting the potential to preserve the deuteron binding. The spin-orbit terms
do not contribute to $S$-wave and ${}^{1}P_{1}$ channel of the $NN$ scattering
and are the smallest contributors to the energy of 4He [21], but they are
important in differentiating between the ${}^{3}P_{0,1,2}$ channels. The AV4’
potential eliminates the tensor terms. As a result, the ${}^{1}S_{0}$ and
${}^{1}P_{1}$ potentials are unaffected, but the coupling between
${}^{3}S_{1}$ and ${}^{3}D_{1}$ channels is gone and the ${}^{3}P_{0,1,2}$
channels deteriorate further. The Fortran code for the AV18 and AVn’
potentials is available at the webpage [199].

####  2.1.2 Three-body $NNN$ potential

The Urbana IX three-body force was originally proposed in combination with the
Argonne AV18 and AV8’ [122]. Although it slightly underbinds the energy of
light nuclei, it has been extensively used to study the equation of state of
nuclear and neutron matter [5, 37, 200, 38, 39, 40]. The Illinois forces
[121], the most recent of which is the Illinois-7 (IL7) [201], have been
introduced to improve the description of both ground- and excited-states of
light nuclei, showing an excellent accuracy [121, 17], but they produce an
unphysical overbinding in pure neutron systems [34].

The three-body Illinois potential consists of two- and three-pion exchange and
a phenomenological short-range component (the UIX force does not include the
three-pion rings):

$\displaystyle V_{ijk}=V_{ijk}^{2\pi}+V_{ijk}^{3\pi}+V_{ijk}^{R}\;.$ (2.23)

The two-pion term, as shown in Fig. 2.1, contains $P$\- and $S$-wave $\pi N$
scattering terms (respectively in Fig. 2.1(a) and Fig. 2.1(b)):

$\displaystyle V_{ijk}^{2\pi}=V_{ijk}^{2\pi,P}+V_{ijk}^{2\pi,S}\;.$ (2.24)

(a) (b)

Figure 2.1: Two-pion exchange processes in the $NNN$ force. 2.1(a) is the
Fujita-Miyazawa $P$-wave term and 2.1(b) the Tucson-Melbourne $S$-wave term.

The $P$-wave component, originally introduced by Fujita-Miyazawa [202],
describes an intermediate excited $\Delta$ resonance produced by the exchange
of two pions between nucleons $i$-$j$ and $j$-$k$, as shown in Fig. 2.1(a),
and it can be written as

$\displaystyle V_{ijk}^{2\pi,P}=A_{2\pi}^{P}\,\mathcal{O}_{ijk}^{2\pi,P}\;,$
(2.25)

where

$\displaystyle A_{2\pi}^{P}$ $\displaystyle=-\frac{2}{81}\frac{f_{\pi
NN}^{2}}{4\pi}\frac{f_{\pi\Delta
N}^{2}}{4\pi}\frac{m_{\pi}^{2}}{m_{\Delta}-m_{N}}\;,$ (2.26a)
$\displaystyle\mathcal{O}_{ijk}^{2\pi,P}$
$\displaystyle=\sum_{cyclic}\left(\phantom{\frac{1}{4}}\\!\\!\\!\\!\Bigl{\\{}X_{ij},X_{jk}\Bigr{\\}}\Bigl{\\{}\bm{\tau}_{i}\cdot\bm{\tau}_{j},\bm{\tau}_{j}\cdot\bm{\tau}_{k}\Bigr{\\}}+\frac{1}{4}\Bigl{[}X_{ij},X_{jk}\Bigr{]}\Bigl{[}\bm{\tau}_{i}\cdot\bm{\tau}_{j},\bm{\tau}_{j}\cdot\bm{\tau}_{k}\Bigr{]}\right)\;,$
(2.26b)

and the $X_{ij}$ operator is the same of Eq. (2.7). The constant
$A_{2\pi}^{P}$ is fitted to reproduce the ground state of light nuclei and
properties of nuclear matter. The $P$-wave TPE term is the longest-ranged
nuclear $NNN$ contribution and it is attractive in all nuclei and nuclear
matter. However it is very small or even slightly repulsive in pure neutron
systems.

The $S$-wave component of TPE three-nucleon force is a simplified form of the
original Tucson-Melbourne model [203], and it involves the $\pi N$ scattering
in the $S$-wave as shown in Fig. 2.1(b). It has the following form:

$\displaystyle V_{ijk}^{2\pi,S}=A_{2\pi}^{S}\,\mathcal{O}_{ijk}^{2\pi,S}\;,$
(2.27)

where

$\displaystyle A_{2\pi}^{S}$ $\displaystyle=\left(\frac{f_{\pi
NN}}{4\pi}\right)^{2}a^{\prime}m_{\pi}^{2}\;,$ (2.28a)
$\displaystyle\mathcal{O}_{ijk}^{2\pi,S}$
$\displaystyle=\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}\,\bm{\tau}_{i}\cdot\bm{\tau}_{k}\;,$
(2.28b)

and the $Z_{\pi}(r)$ function is defined as

$\displaystyle
Z_{\pi}(r)=\frac{\mu_{\pi}r}{3}\Bigl{[}Y_{\pi}(r)-T_{\pi}(r)\Bigr{]}\;.$
(2.29)

The $S$-wave TPE term is required by chiral perturbation theory but in
practice its contribution is only 3%–4% of $V_{ijk}^{2\pi,P}$ in light nuclei.

The three-pion term (Fig. 2.2) was introduced in the Illinois potentials. It
consists of the subset of three-pion rings that contain only one $\Delta$ mass
in the energy denominators.

(a) (b)

Figure 2.2: Three-pion exchange processes in the $NNN$ force.

As discussed in Ref. [121], these diagrams result in a large number of terms,
the most important of which are the ones independent of cyclic permutations of
$ijk$:

$\displaystyle V_{ijk}^{3\pi}=A_{3\pi}\,\mathcal{O}_{ijk}^{3\pi}\;,$ (2.30)

where

$\displaystyle A_{3\pi}$ $\displaystyle=\left(\frac{f^{2}_{\pi
NN}}{4\pi}\frac{m_{\pi}}{3}\right)^{3}\frac{f^{2}_{\pi N\Delta}}{f^{2}_{\pi
NN}}\frac{1}{(m_{\Delta}-m_{N})^{2}}\;,$ (2.31a)
$\displaystyle\mathcal{O}_{ijk}^{3\pi}$
$\displaystyle\simeq\frac{50}{3}S_{ijk}^{\tau}\,S_{ijk}^{\sigma}+\frac{26}{3}A_{ijk}^{\tau}\,A_{ijk}^{\sigma}\;.$
(2.31b)

The letters $S$ and $A$ denote operators that are symmetric and antisymmetric
under the exchange of $j$ with $k$. Superscripts $\tau$ and $\sigma$ label
operators containing isospin and spin-space parts, respectively. The isospin
operators are

$\displaystyle S_{ijk}^{\tau}$
$\displaystyle=2+\frac{2}{3}\left(\bm{\tau}_{i}\cdot\bm{\tau}_{j}+\bm{\tau}_{j}\cdot\bm{\tau}_{k}+\bm{\tau}_{k}\cdot\bm{\tau}_{i}\right)=4\,P_{T=3/2}\;,$
(2.32a) $\displaystyle A_{ijk}^{\tau}$
$\displaystyle=\frac{1}{3}\,i\,\bm{\tau}_{i}\cdot\bm{\tau}_{j}\times\bm{\tau}_{k}=-\frac{1}{6}\Bigl{[}\bm{\tau}_{i}\cdot\bm{\tau}_{j},\bm{\tau}_{j}\cdot\bm{\tau}_{k}\Bigr{]}\;,$
(2.32b)

where $S_{ijk}^{\tau}$ is a projector onto isospin 3/2 triples and
$A_{ijk}^{\tau}$ has the same isospin structure as the commutator part of
$V_{ijk}^{2\pi,P}$. The spin-space operators have many terms and they are
listed in the Appendix of Ref. [121]. An important aspect of this structure is
that there is a significant attractive term which acts only in $T=3/2$
triples, so the net effect of $V_{ijk}^{3\pi}$ is slight repulsion in
$S$-shell nuclei and larger attraction in $P$-shell nuclei. However, in most
light nuclei the contribution of this term is rather small, $\langle
V_{ijk}^{3\pi}\rangle\lesssim 0.1\langle V_{ijk}^{2\pi}\rangle$.

The last term of Eq. (2.23) was introduced to compensate the overbinding in
nuclei and the large equilibrium density of nuclear matter given by the
previous operators. It is strictly phenomenological and purely central and
repulsive, and it describes the modification of the contribution of the TPE
$\Delta$-box diagrams to $v_{ij}$ due to the presence of the third nucleon $k$
(Fig. 2.3). It takes the form:

$\displaystyle
V_{ijk}^{R}=A_{R}\,\mathcal{O}^{R}_{ijk}=A_{R}\sum_{cyclic}T_{\pi}^{2}(r_{ij})\,T_{\pi}^{2}(r_{jk})\;,$
(2.33)

where $T_{\pi}(r)$ is the OPE tensor potential defined in Eq. (2.10).

Figure 2.3: Repulsive short-range contribution included in the $NNN$ force.

Finally, 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}\;.$
(2.34)

###  2.2 Interactions: hyperons and nucleons

We present a detailed description of the $\Lambda N$ and $\Lambda NN$
interaction as developed by Bodmer, Usmani and Carlson following the scheme of
the Argonne potentials [190, 183, 181, 192, 185, 191, 186, 184, 187, 188,
189]. The interaction is written in coordinates space and it includes two- and
three-body hyperon nucleon components with an explicit hard-core repulsion
between baryons and a charge symmetry breaking term. We introduce also an
effective $\Lambda\Lambda$ interaction mainly used in variational [194, 195]
and cluster model [171, 173] calculations for double $\Lambda$ hypernuclei.

####  2.2.1 Two-body $\Lambda N$ potential

#####  $\Lambda N$ charge symmetric potential

The $\Lambda$ particle has isospin $I=0$, so there is no OPE term, being the
strong $\Lambda\Lambda\pi$ vertex forbidden due to isospin conservation. The
$\Lambda$ hyperon can thus exchange a pion only with a $\Lambda\pi\Sigma$
vertex. The lowest order $\Lambda N$ coupling must therefore involve the
exchange of two pions, with the formation of a virtual $\Sigma$ hyperon, as
illustrated in Figs. 2.4(a) and 2.4(b). The TPE interaction is intermediate
range with respect to the long range part of $NN$ force. One meson exchange
processes can only occur through the exchange of a $K,K^{*}$ kaon pair, that
contributes in exchanging the strangeness between the two baryons, as shown in
Fig. 2.4(c). The $K,K^{*}$ potential is short-range and contributes to the
space-exchange and $\Lambda N$ tensor potential. The latter is expected to be
quite weak because the $K$ and $K^{*}$ tensor contributions have opposite sign
[204].

(a) (b) (c)

Figure 2.4: Meson exchange processes in the $\Lambda N$ force. 2.4(a) and
2.4(b) are the TPE diagrams. 2.4(c) represents the kaon exchange channel.

The $\Lambda N$ interaction has been modeled with an Urbana-type potential
[205] with spin-spin and space-exchange components and a TPE tail which is
consistent with the available $\Lambda p$ scattering data below the $\Sigma$
threshold:

$\displaystyle v_{\lambda i}=v_{0}(r_{\lambda
i})(1-\varepsilon+\varepsilon\,\mathcal{P}_{x})+\frac{1}{4}v_{\sigma}T^{2}_{\pi}(r_{\lambda
i})\,{\bm{\sigma}}_{\lambda}\cdot{\bm{\sigma}}_{i}\;,$ (2.35)

where

$\displaystyle v_{0}(r_{\lambda i})=v_{c}(r_{\lambda
i})-\bar{v}\,T^{2}_{\pi}(r_{\lambda i})\;.$ (2.36)

Here,

$\displaystyle
v_{c}(r)=W_{c}\Bigl{(}1+\operatorname{e}^{\frac{r-\bar{r}}{a}}\Bigr{)}^{-1}$
(2.37)

is a Wood-Saxon repulsive potential introduced, similarly to the Argonne $NN$
interaction, in order to include all the short-range contributions and
$T_{\pi}(r)$ is the regularized OPE tensor operator defined in Eq. (2.10). The
term $\bar{v}\,T^{2}_{\pi}(r_{\lambda i})$ corresponds to a TPE mechanism due
to OPE transition potentials $\left(\Lambda N\leftrightarrow\Sigma
N,\Sigma\Delta\right)$ dominated by their tensor components. The $\Lambda p$
scattering at low energies is well fitted with $\bar{v}=6.15(5)$ MeV. The
terms $\bar{v}=(v_{s}+3v_{t})/4$ and $v_{\sigma}=v_{s}-v_{t}$ are the spin-
average and spin-dependent strengths, where $v_{s}$ and $v_{t}$ denote
singlet- and triplet-state strengths, respectively. $\mathcal{P}_{x}$ is the
$\Lambda N$ space-exchange operator and $\varepsilon$ the corresponding
exchange parameter, which is quite poorly determined from the $\Lambda p$
forward-backward asymmetry to be $\varepsilon\simeq 0.1\div 0.38$. All the
parameters defining the $\Lambda N$ potential can be found in Tab. 2.1.

#####  $\Lambda N$ charge symmetry breaking potential

The $\Lambda$-nucleon interaction should distinguish between the nucleon
isospin channels $\Lambda p$ and $\Lambda n$. The mirror pair of hypernuclei
${}^{4}_{\Lambda}$H and ${}^{4}_{\Lambda}$He is the main source of information
about the charge symmetry breaking (CSB) $\Lambda N$ interaction. The
experimental data for $A=4$ $\Lambda$ hypernuclei [77], show indeed a clear
difference in the $\Lambda$ separation energies for the $(0^{+})$ ground state

$\displaystyle B_{\Lambda}\left({}^{4}_{\Lambda}\text{H}\right)$
$\displaystyle=2.04(4)\leavevmode\nobreak\ \text{MeV}\;,$ (2.38a)
$\displaystyle B_{\Lambda}\left({}^{4}_{\Lambda}\text{He}\right)$
$\displaystyle=2.39(3)\leavevmode\nobreak\ \text{MeV}\;,$ (2.38b)

and for the $(1^{+})$ excited state

$\displaystyle B_{\Lambda}^{*}\left({}^{4}_{\Lambda}\text{H}\right)$
$\displaystyle=1.00(6)\leavevmode\nobreak\ \text{MeV}\;,$ (2.39a)
$\displaystyle B_{\Lambda}^{*}\left({}^{4}_{\Lambda}\text{He}\right)$
$\displaystyle=1.24(6)\leavevmode\nobreak\ \text{MeV}\;.$ (2.39b)

The differences in the hyperon separation energies are:

$\displaystyle\Delta B_{\Lambda}$ $\displaystyle=0.35(6)\leavevmode\nobreak\
\text{MeV}\;,$ (2.40a) $\displaystyle\Delta B_{\Lambda}^{*}$
$\displaystyle=0.24(6)\leavevmode\nobreak\ \text{MeV}\;.$ (2.40b)

However, the experimental values $\Delta B_{\Lambda}$ must be corrected to
include the difference $\Delta B_{c}$ due to the Coulomb interaction in order
to obtain the values to be attributed to CSB effects. By means of a
variational calculation, Bodmer and Usmani [183] estimated the Coulomb
contribution to be rather small

$\displaystyle|\Delta B_{c}|$ $\displaystyle=0.05(2)\leavevmode\nobreak\
\text{MeV}\;,$ (2.41a) $\displaystyle|\Delta B_{c}^{*}|$
$\displaystyle=0.025(15)\leavevmode\nobreak\ \text{MeV}\;,$ (2.41b)

and they were able to reproduce the differences in the $\Lambda$ separation
energies by means of a phenomenological spin dependent CSB potential. It was
found that the CSB interaction is effectively spin independent and can be
simply expressed (as subsequently reported in Ref. [186]) by

$\displaystyle v_{\lambda i}^{CSB}=C_{\tau}\,T_{\pi}^{2}\left(r_{\lambda
i}\right)\tau_{i}^{z}\quad\quad C_{\tau}=-0.050(5)\leavevmode\nobreak\
\text{MeV}\;.$ (2.42)

Being $C_{\tau}$ negative, the $\Lambda p$ channel becomes attractive while
the $\Lambda n$ channel is repulsive, consistently with the experimental
results for ${}^{4}_{\Lambda}$H and ${}^{4}_{\Lambda}$He. The contribution of
CSB is expected to be very small in symmetric hypernuclei (if Coulomb is
neglected) but could have a significant effect in hypernuclei with an neutron
(or proton) excess.

####  2.2.2 Three-body $\Lambda NN$ potential

The $\Lambda N$ force as obtained by fitting the $\Lambda p$ scattering does
not provide a good account of the experimental binding energies, as in the
case of nuclei with the bare $NN$ interaction. A three-body $\Lambda NN$ force
is required in this scheme to solve the overbinding. The $\Lambda NN$
potential is at the same TPE order of the $\Lambda N$ force and it includes
diagrams involving two nucleons and one hyperon, as reported in Fig. 2.5.

(a) (b) (c)

Figure 2.5: Two-pion exchange processes in three-body $\Lambda NN$ force.
2.5(a) and 2.5(b) are, respectively, the $P$\- and $S$-wave TPE contributions.
2.5(c) is the phenomenological dispersive term.

The diagrams in Fig. 2.5(a) and Fig. 2.5(b) correspond respectively to the
$P$-wave and $S$-wave TPE

$\displaystyle v^{2\pi}_{\lambda ij}=v^{2\pi,P}_{\lambda
ij}+v^{2\pi,S}_{\lambda ij}\;,$ (2.43)

that can be written in the following form:

$\displaystyle v_{\lambda ij}^{2\pi,P}$
$\displaystyle=\widetilde{C}_{P}\,\mathcal{O}_{\lambda ij}^{2\pi,P}$
$\displaystyle=-\frac{C_{P}}{6}\Bigl{\\{}X_{i\lambda}\,,X_{\lambda
j}\Bigr{\\}}\,{\bm{\tau}}_{i}\cdot{\bm{\tau}}_{j}\;,$ (2.44) $\displaystyle
v_{\lambda ij}^{2\pi,S}$ $\displaystyle=C_{S}\,O_{\lambda ij}^{2\pi,S}$
$\displaystyle=C_{S}\,Z\left(r_{\lambda i}\right)Z\left(r_{\lambda
j}\right)\,{\bm{\sigma}}_{i}\cdot\hat{\bm{r}}_{i\lambda}\,{\bm{\sigma}}_{j}\cdot\hat{\bm{r}}_{j\lambda}\,{\bm{\tau}}_{i}\cdot{\bm{\tau}}_{j}\;.$
(2.45)

The structure of $V_{\lambda ij}^{2\pi}$ is very close to the Fujita-Miyazawa
$P$-wave term and the Tucson-Melbourne $S$-wave term of the nuclear
$V_{ijk}^{2\pi}$ (see Eqs. (2.26) and (2.28)). In the hypernuclear sector,
however, there are simplifications because only two nucleons at a time enter
the picture, so there are no cyclic summations, and the $\Lambda$ particle has
isospin zero, thus there is no $\bm{\tau}_{\lambda}$ operator involved. As
reported in Ref. [121], the strength of $V_{ijk}^{2\pi,S}$ is
$\left|A_{2\pi}^{S}\right|\simeq 0.8$ MeV. However, in other references it is
assumed to have a value of 1.0 MeV. Comparing the Tucson-Melbourne $NNN$ model
with Eq. (2.45) for the $\Lambda NN$ potential, one may write an identical
structure for both $S$-wave $\Lambda NN$ and $NNN$ potentials as follows:

$\displaystyle C_{S}\,\mathcal{O}_{\lambda
ij}^{2\pi,S}=A_{S}^{2\pi}\,\mathcal{O}_{ijk}^{2\pi,S}\;.$ (2.46)

This directly relates $C_{S}$ in the strange sector to $A_{2\pi}^{S}$ in the
non-strange sector. Since the $\Sigma$-$\Lambda$ mass difference is small
compared to the $\Delta$-$N$ mass difference, the $2\pi$ $\Lambda NN$
potential of $S=-1$ sector is stronger than the non-strange $NNN$ potential of
$S=0$ sector. This provides stronger strengths in the case of $\Lambda NN$
potential compared to the $NNN$ potential. It is therefore expected that the
value of $C_{S}$ would be more than 1.0 MeV, and is taken to be 1.5 MeV [189].
However, the $S$-wave component is expected to be quite weak, at least in
spin-zero core hypernuclei, and indeed it has been neglected in variational
calculations for ${}^{17}_{\leavevmode\nobreak\ \Lambda}$O and
${}^{5}_{\Lambda}$He [192, 187, 194].

The last diagram (Fig. 2.5(c)) represents the dispersive contribution
associated with the medium modifications of the intermediate state potentials
for the $\Sigma$, $N$, $\Delta$ due to the presence of the second nucleon.
This term describes all the short-range contributions and it is expected to be
repulsive due to the suppression mechanism associated with the $\Lambda
N$-$\Sigma N$ coupling [206, 207]. The interaction of the intermediate states
$\Sigma$, $N$, $\Delta$ with a nucleon of the medium will be predominantly
through a TPE potential, proportional to $T_{\pi}^{2}(r)$, with an explicit
spin dependence (negligible for spin-zero core hypernuclei):

$\displaystyle v_{\lambda ij}^{D}=W_{D}\,\mathcal{O}_{\lambda
ij}^{D}=W_{D}\,T_{\pi}^{2}\left(r_{\lambda
i}\right)T^{2}_{\pi}\left(r_{\lambda
j}\right)\\!\\!\bigg{[}1+\frac{1}{6}{\bm{\sigma}}_{\lambda}\\!\cdot\\!\left({\bm{\sigma}}_{i}+{\bm{\sigma}}_{j}\right)\bigg{]}\;.$
(2.47)

The radial functions $T_{\pi}(r)$ and $Z_{\pi}(r)$ are the same of the nuclear
potential, see Eq. (2.10) and Eq. (2.29). The operator $X_{\lambda i}$ is the
same of Eq. (2.7), in which the first nucleon is replaced by the $\Lambda$
particle.

It is important to note that the three-body $\Lambda NN$ interaction have been
investigated in variational calculations for ${}_{\Lambda}^{5}$He [185, 187,
189], ${}_{\Lambda\Lambda}^{\;\;\,6}$He [194, 195] and
${}_{\leavevmode\nobreak\ \Lambda}^{17}$O [192, 185], resulting in a range of
values for the $C_{P}$ and $W_{D}$ parameters (see Tab. 2.1) that gives good
description of the properties of the studied hypernuclei. A unique set of
parameters that reproduces all the available experimental energies for single
(and double) $\Lambda$ hypernuclei has not been set yet.

A second crucial observation is that, differently to the nucleon sector, both
two- and three-body lambda-nucleon interactions are at the same TPE order. In
addition, the mass difference between the $\Lambda$ particle and its
excitation $\Sigma$ is much smaller than the mass difference between the
nucleon and the $\Delta$ resonance. Thus, the $\Lambda NN$ interaction can not
be neglected in this framework but it is a key ingredient in addition to the
$\Lambda N$ force for any consistent theoretical calculation involving
$\Lambda$ hyperons.

Constant | Value | Unit
---|---|---
$W_{c}$ | $2137$ | MeV
$\bar{r}$ | $0.5$ | fm
$a$ | $0.2$ | fm
$v_{s}$ | $6.33,6.28$ | MeV
$v_{t}$ | $6.09,6.04$ | MeV
$\bar{v}$ | $6.15(5)$ | MeV
$v_{\sigma}$ | $0.24$ | MeV
$c$ | $2.0$ | fm-2
$\varepsilon$ | $0.1\div 0.38$ | —
$C_{\tau}$ | -0.050(5) | MeV
$C_{P}$ | $0.5\div 2.5$ | MeV
$C_{S}$ | $\simeq 1.5$ | MeV
$W_{D}$ | $0.002\div 0.058$ | MeV
Table 2.1: Parameters of the $\Lambda N$ and $\Lambda NN$ interaction (see
[189] and reference therein). For $C_{P}$ and $W_{D}$ the variational allowed
range is shown. The value of the charge symmetry breaking parameter $C_{\tau}$
is from Ref. [186].

####  2.2.3 Two-body $\Lambda\Lambda$ potential

Due to the impossibility to collect $\Lambda\Lambda$ scattering data,
experimental information about the $\Lambda\Lambda$ interaction can be
obtained only from the $\Lambda\Lambda$ separation energy of the observed
double $\Lambda$ hypernuclei, ${}^{\;\;\,6}_{\Lambda\Lambda}$He [91, 92, 93],
${}^{\;13}_{\Lambda\Lambda}$B [92] and the isotopes of
${}^{\;10}_{\Lambda\Lambda}$Be ($A=10\div 12$) [94, 92]. Evidence for the
production of ${}_{\Lambda\Lambda}^{\;\;\,4}$H has been reported in Ref.
[208], but no information about the $\Lambda\Lambda$ separation energy was
found. On the other hand, there is a theoretical indication for the one-boson
exchange (OBE) part of the $\Lambda\Lambda$ interaction coming from the
$SU(3)$-invariance of coupling constants, but the $\Lambda\Lambda$ force is
still far to be settled.

In the next, we follow the guide line adopted in the three- and four-body
cluster models for double $\Lambda$ hypernuclei [171, 173], which was also
used in Faddeev-Yakubovsky calculations for light double $\Lambda$ hypernuclei
[209] and in variational calculations on ${}^{\;\;\,4}_{\Lambda\Lambda}$H
[193, 210], ${}^{\;\;\,5}_{\Lambda\Lambda}$H and
${}^{\;\;\,5}_{\Lambda\Lambda}$He [193, 211] and
${}^{\;\;\,6}_{\Lambda\Lambda}$He [194, 195, 193, 211], with different
parametrizations. The employed OBE-simulating $\Lambda\Lambda$ effective
interaction is a low-energy phase equivalent Nijmegen interaction represented
by a sum of three Gaussians:

$\displaystyle
v_{\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}}\;.$
(2.48)

The most recent parametrization of the potential (see Tab. 2.2), was fitted in
order to simulate the $\Lambda\Lambda$ sector of the Nijmegen F (NF)
interaction [150, 151, 152]. The NF is the simplest among the Nijmegen models
with a scalar nonet, which seems to be more appropriate than the versions
including only a scalar singlet in order to reproduce the weak binding energy
indicated by the NAGARA event [91]. The components $k=1,2$ of the above
Gaussian potential are determined so as to simulate the $\Lambda\Lambda$
sector of NF and the strength of the part for $k=3$ is adjusted so as to
reproduce the ${}^{\;\;\,6}_{\Lambda\Lambda}$He NAGARA experimental double
$\Lambda$ separation energy of $7.25\pm 0.19^{+0.18}_{-0.11}$ MeV. In 2010,
Nakazawa reported a new, more precise determination of
$B_{\Lambda\Lambda}=6.93\pm 0.16$ MeV for ${}^{\;\;\,6}_{\Lambda\Lambda}$He
[92], obtained via the $\Xi^{-}$ hyperon capture at rest reaction in a hybrid
emulsion. This value has been recently revised to $B_{\Lambda\Lambda}=6.91\pm
0.16$ MeV by the E373 (KEK-PS) Collaboration [93]. No references were found
about the refitting of the $\Lambda\Lambda$ Gaussian potential on the more
recent experimental result, which is in any case compatible with the NAGARA
event. We therefore consider the original parametrization of Ref. [173].

$\mu^{(k)}$ | $0.555$ | $1.656$ |  $8.163$
---|---|---|---
$v_{0}^{(k)}$ | $-10.67$ | $-93.51$ |  $4884$
$v_{\sigma}^{(k)}$ | $0.0966$ | $16.08$ |  $915.8$
Table 2.2: Parameters of the the $\Lambda\Lambda$ interaction. The size
parameters $\mu^{(k)}$ are in fm-2 and the strengths $v_{0}^{(k)}$ and
$v_{\sigma}^{(k)}$ are in MeV [173].

## Chapter 3 Method

In nuclear physics, many-body calculations are used to understand the nuclear
systems in the non-relativistic regime. When interested in low energy
phenomena, a nucleus (or an extensive nucleonic system) can be described as a
collection of particles interacting via a potential that depends on positions,
momenta, spin and isospin. The properties of the system can be determined by
solving a many-body Schrödinger equation. Such calculations can study, for
example, binding energies, excitation spectra, densities, reactions and many
other aspects of nuclei. The equation of state, masses, radii and other
properties are obtained by describing astrophysical objects as a nuclear
infinite medium.

The two main problems related to microscopic few- and many-body calculations
in nuclear physics are the determination of the Hamiltonian and the method
used to accurately solve the Schrödinger equation. In the previous chapter, we
have already seen how to build a realistic nuclear Hamiltonian, including also
strange degrees of freedom. In the next we will focus on the methodological
part presenting a class of Quantum Monte Carlo algorithms, the Diffusion Monte
Carlo (DMC) and, more in detail, the Auxiliary Field Diffusion Monte Carlo
(AFDMC). Such methods are based on evolving a trial wave function in imaginary
time to yield the ground state of the system. The DMC method sums explicitly
over spin and isospin states and can use very sophisticated wave functions.
However, it is limited to small systems. In the AFDMC, in addition to the
coordinates, also the spin and isospin degrees of freedom are sampled. It can
thus treat larger systems but there are some limitations on the trial wave
function and the nuclear potentials that can be handled.

Strangeness can be included in AFDMC calculations by adding hyperons to the
standard nucleons. The interaction between hyperons and nucleons presented in
the previous chapter is written in a suitable form to be treated within this
algorithm. By extending the AFDMC nuclear wave function to the hyperonic
sector, it is possible to study both hypernuclei and hypermatter. A new QMC
approach to strange physics is thus now available.

###  3.1 Diffusion Monte Carlo

The Diffusion Monte Carlo method [212, 136, 213, 214] projects the ground-
state out of a stationary trial wave function $|\psi_{T}\rangle$ not
orthogonal to the true ground state. Consider the many-body time dependent
Schrödinger equation with its formal solution

$\displaystyle i\hbar\frac{\partial}{\partial
t}|\psi(t)\rangle=(H-E_{T})|\psi(t)\rangle\quad\Rightarrow\quad|\psi(t+dt)\rangle=\operatorname{e}^{-\frac{i}{\hbar}(H-E_{T})dt}|\psi(t)\rangle\;,$
(3.1)

and let move to the imaginary time $\tau=it/\hbar$111with this definition
$\tau$ has the dimensions of the inverse of an energy.:

$\displaystyle-\frac{\partial}{\partial\tau}|\psi(\tau)\rangle=(H-E_{T})|\psi(\tau)\rangle\quad\Rightarrow\quad|\psi(\tau+d\tau)\rangle=\operatorname{e}^{-(H-E_{T})d\tau}|\psi(\tau)\rangle\;.$
(3.2)

The stationary states $|\psi(0)\rangle=|\psi_{T}\rangle$ are the same for both
normal and imaginary time Schrödinger equations and we can expand them on a
complete orthonormal set of eigenvectors $|\varphi_{n}\rangle$ of the
Hamiltonian $H$:

$\displaystyle|\psi_{T}\rangle=\sum_{n=0}^{\infty}c_{n}|\varphi_{n}\rangle\;.$
(3.3)

Supposing that the $|\psi_{T}\rangle$ is not orthogonal to the true ground
state, i.e. $c_{0}\neq 0$, and that at least the ground state is non
degenerate, i.e. $E_{n}\geq E_{n-1}>E_{0}$, where $E_{n}$ are the eigenvalues
of $H$ related to $|\varphi_{n}\rangle$, the imaginary time evolution of
$|\psi_{T}\rangle$ is given by

$\displaystyle|\psi(\tau)\rangle$
$\displaystyle=\sum_{n=0}^{\infty}c_{n}\operatorname{e}^{-(E_{n}-E_{T})\tau}|\varphi_{n}\rangle\;,$
$\displaystyle=c_{0}\operatorname{e}^{-(E_{0}-E_{T})\tau}|\varphi_{0}\rangle+\sum_{n=1}^{\infty}c_{n}\operatorname{e}^{-(E_{n}-E_{T})\tau}|\varphi_{n}\rangle\;.$
(3.4)

If the energy offset $E_{T}$ is the exact ground state energy $E_{0}$, in the
limit $\tau\rightarrow\infty$ the components of Eq. (3.4) for $n>0$ vanish and
we are left with

$\displaystyle\lim_{\tau\rightarrow\infty}|\psi(\tau)\rangle=c_{0}|\varphi_{0}\rangle\;.$
(3.5)

Starting from a generic initial trial wave function $|\psi_{T}\rangle$ not
orthogonal to the ground state, and adjusting the energy offset $E_{T}$ to be
as close as possible to $E_{0}$, in the limit of infinite imaginary time, one
can project out the exact ground state $c_{0}|\varphi_{0}\rangle$ giving
access to the lowest energy properties of the system.

Consider the imaginary time propagation of Eq. (3.2) and insert a completeness
on the orthonormal basis $|R^{\prime}\rangle$, where $R$ represents a
configuration $\\{\bm{r}_{1},\ldots,\bm{r}_{\mathcal{N}}\\}$ of the
$\mathcal{N}$ particle system with all its degrees of freedom:

$\displaystyle|\psi(\tau+d\tau)\rangle$
$\displaystyle=\operatorname{e}^{-(H-E_{T})d\tau}|\psi(\tau)\rangle\;,$
$\displaystyle=\int
dR^{\prime}\operatorname{e}^{-(H-E_{T})d\tau}|R^{\prime}\rangle\langle
R^{\prime}|\psi(\tau)\rangle\;.$ (3.6)

Projecting on the coordinates $\langle R|$ leads to

$\displaystyle\langle R|\psi(\tau+d\tau)\rangle=\int dR^{\prime}\,\langle
R|\operatorname{e}^{-(H-E_{T})d\tau}|R^{\prime}\rangle\langle
R^{\prime}|\psi(\tau)\rangle\;,$ (3.7)

where $\langle
R|\operatorname{e}^{-(H-E_{T})d\tau}|R^{\prime}\rangle=G(R,R^{\prime},d\tau)$
is the Green’s function of the operator
$(H-E_{T})+\frac{\partial}{\partial\tau}$. Recalling that $\langle
R|\psi(\tau)\rangle=\psi(R,\tau)$, we can write Eq. (3.2) as

$\displaystyle-\frac{\partial}{\partial\tau}\psi(R,\tau)$
$\displaystyle=(H-E_{T})\psi(R,\tau)\;,$ (3.8)
$\displaystyle\psi(R,\tau+d\tau)$ $\displaystyle=\int
dR^{\prime}\,G(R,R^{\prime},d\tau)\,\psi(R^{\prime},\tau)\;.$ (3.9)

If we consider a non-interacting many-body system, i.e. the Hamiltonian is
given by the pure kinetic term

$\displaystyle
H_{0}=T=-\frac{\hbar^{2}}{2m}\sum_{i=1}^{\mathcal{N}}\nabla_{i}^{2}\;,$ (3.10)

the Schrödinger equation (3.8) becomes a $3\mathcal{N}$-dimensional diffusion
equation. By writing the Green’s function of Eq. (3.9) in momentum space by
means of the Fourier transform, it is possible to show that $G_{0}$ is a
Gaussian with variance proportional to $\tau$

$\displaystyle G_{0}(R,R^{\prime},d\tau)=\left(\frac{1}{4\pi
Dd\tau}\right)^{\frac{3\mathcal{N}}{2}}\\!\operatorname{e}^{-\frac{(R-R^{\prime})^{2}}{4Dd\tau}}\;,$
(3.11)

where $D=\hbar^{2}/2m$ is the diffusion constant of a set of particles in
Brownian motion with a dynamic governed by random collisions. This
interpretation can be implemented by representing the wave function
$\psi(R,\tau)$ by a set of discrete sampling points, called _walkers_

$\displaystyle\psi(R,\tau)=\sum_{k}\delta(R-R_{k})\;,$ (3.12)

and evolving this discrete distribution for an imaginary time $d\tau$ by means
of Eq. (3.9):

$\displaystyle\psi(R,\tau+d\tau)=\sum_{k}G_{0}(R,R_{k},d\tau)\;.$ (3.13)

The result is a set of Gaussians that in the infinite imaginary time limit
represents a distribution of walkers according to the lowest state of the
Hamiltonian, that can be used to calculate the ground state properties of the
system.

Let now consider the full Hamiltonian where the interaction is described by a
central potential in coordinate space:

$\displaystyle
H=T+V=-\frac{\hbar^{2}}{2m}\sum_{i=1}^{\mathcal{N}}\nabla_{i}^{2}+V(R)\;.$
(3.14)

Because $T$ and $V$ in general do not commute, it is not possible to directly
split the propagator in a kinetic and a potential part

$\displaystyle\operatorname{e}^{-(H-E_{T})d\tau}\neq\operatorname{e}^{-Td\tau}\operatorname{e}^{-(V-E_{T})d\tau}\;,$
(3.15)

and thus the analytic solution of the Green’s function $\langle
R|\operatorname{e}^{-(T+V-E_{T})d\tau}|R^{\prime}\rangle$ is not known in most
of the cases. However, by means of the Trotter-Suzuki formula to order
$d\tau^{3}$

$\displaystyle\operatorname{e}^{-(A+B)d\tau}=\operatorname{e}^{-A\frac{d\tau}{2}}\operatorname{e}^{-Bd\tau}\operatorname{e}^{-A\frac{d\tau}{2}}\,+\operatorname{o}\left(d\tau^{3}\right)\;,$
(3.16)

which is an improvement of the standard

$\displaystyle\operatorname{e}^{-(A+B)d\tau}=\operatorname{e}^{-Ad\tau}\operatorname{e}^{-Bd\tau}\,+\operatorname{o}\left(d\tau^{2}\right)\;,$
(3.17)

in the limit of small imaginary time step $d\tau$ it is possible to write an
approximate solution for $\psi(R,\tau+d\tau)$:

$\displaystyle\psi(R,\tau+d\tau)$ $\displaystyle\simeq\int dR^{\prime}\langle
R|\operatorname{e}^{-V\frac{d\tau}{2}}\operatorname{e}^{-Td\tau}\operatorname{e}^{-V\frac{d\tau}{2}}\operatorname{e}^{E_{T}d\tau}|R^{\prime}\rangle\,\psi(R^{\prime},\tau)\;,$
$\displaystyle\simeq\int dR^{\prime}\underbrace{\langle
R|\operatorname{e}^{-Td\tau}|R^{\prime}\rangle}_{G_{0}(R,R^{\prime},d\tau)}\underbrace{\phantom{\langle}\\!\\!\operatorname{e}^{-\left(\frac{V(R)+V(R^{\prime})}{2}-E_{T}\right)d\tau}}_{G_{V}(R,R^{\prime},d\tau)}\psi(R^{\prime},\tau)\;,$
$\displaystyle\simeq\left(\frac{1}{4\pi
Dd\tau}\right)^{\frac{3\mathcal{N}}{2}}\\!\\!\int
dR^{\prime}\operatorname{e}^{-\frac{(R-R^{\prime})^{2}}{4Dd\tau}}\operatorname{e}^{-\left(\frac{V(R)+V(R^{\prime})}{2}-E_{T}\right)d\tau}\psi(R^{\prime},\tau)\;,$
(3.18)

which is the same of Eq. (3.9) with the full Green’s function given by

$\displaystyle G(R,R^{\prime},d\tau)\simeq
G_{0}(R,R^{\prime},d\tau)\,G_{V}(R,R^{\prime},d\tau)\;.$ (3.19)

According to the interacting Hamiltonian, the propagation of $\psi(R,\tau)$
for $d\tau\rightarrow 0$ is thus described by the integral (3.18) and the long
imaginary time evolution necessary to project out the ground state component
of the wave function is realized by iteration until convergence is reached.

The steps of this process, that constitute the Diffusion Monte Carlo
algorithm, can be summarized as follows:

  1. 1.

An initial distribution of walkers $w_{i}$ with $i=1,\ldots,\mathcal{N}_{w}$
is sampled from the trial wave function $\langle
R|\psi_{T}\rangle=\psi_{T}(R)$ and the starting trial energy $E_{T}$ is chosen
(for instance from a variational calculation or close to the expected result).

  2. 2.

The spacial degrees of freedom are propagated for small imaginary time step
$d\tau$ with probability density $G_{0}(R,R^{\prime},d\tau)$, i.e. the
coordinates of the walkers are diffused by means of a Brownian motion

$\displaystyle R=R^{\prime}+\xi\;,$ (3.20)

where $\xi$ is a stochastic variable distributed according to a Gaussian
probability density with $\sigma^{2}=2Dd\tau$ and zero average.

  3. 3.

For each walker, a weight

$\displaystyle\omega_{i}=G_{V}(R,R^{\prime},d\tau)=\operatorname{e}^{-\left(\frac{V(R)+V(R^{\prime})}{2}-E_{T}\right)d\tau}\;,$
(3.21)

is assigned. The estimator contributions (kinetic energy, potential energy,
root mean square radii, densities, …) are evaluated on the imaginary time
propagated configurations, weighting the results according to $\omega_{i}$.

  4. 4.

The _branching_ process is applied to the propagated walkers. $\omega_{i}$
represents the probability of a configuration to multiply at the next step
according to the normalization. This process is realized by generating from
each $w_{i}$ a number of walker copies

$\displaystyle n_{i}=[\omega_{i}+\eta_{i}]\;,$ (3.22)

where $\eta_{i}$ is a random number uniformly distributed in the interval
$[0,1]$ and $[x]$ means integer part of $x$. In such a way, depending on the
potential $V(R)$ and the trial energy $E_{T}$, some configurations will
disappear and some other will replicate, resulting in the evolution of walker
population which is now made of
$\widetilde{\mathcal{N}}_{w}=\sum_{i=1}^{\mathcal{N}_{w}}n_{i}$ walkers. A
simple solution in order to control the fluctuations of walker population is
to multiply the weight $\omega_{i}$ by a factor
$\mathcal{N}_{w}/\widetilde{\mathcal{N}}_{w}$, adjusting thus the branching
process at each time step. This solution is not efficient if the potential
diverges. The corrections applied run-time could generate a lot of copies from
just few good parent walkers and the population will be stabilized but not
correctly represented. A better sampling technique is described in § 3.1.1.

  5. 5.

Iterate from 2 to 4 as long as necessary until convergence is reached, i.e.
for large enough $\tau$ to reach the infinite limit of Eq. (3.5). In this
limit, the configurations $\\{R\\}$ are distributed according to the lowest
energy state $\psi_{0}(R,\tau)$. Therefore, we can compute the ground state
expectation values of observables that commute with the Hamiltonian

$\displaystyle\\!\\!\langle\mathcal{O}\rangle=\frac{\langle\psi_{0}|\mathcal{O}|\psi_{0}\rangle}{\langle\psi_{0}|\psi_{0}\rangle}=\\!\lim_{\tau\rightarrow\infty}\\!\frac{\langle\psi_{T}|\mathcal{O}|\psi(\tau)\rangle}{\langle\psi_{T}|\psi(\tau)\rangle}=\\!\lim_{\tau\rightarrow\infty}\int\\!\\!dR\frac{\langle\psi_{T}|\mathcal{O}|R\rangle\psi(R,\tau)}{\psi_{T}(R)\psi(R,\tau)}\;,$
(3.23)

by means of

$\displaystyle\langle\mathcal{O}\rangle=\frac{\sum_{\\{R\\}}\langle
R|\mathcal{O}|\psi_{T}\rangle}{\sum_{\\{R\\}}\langle
R|\psi_{T}\rangle}=\frac{\sum_{\\{R\\}}\mathcal{O}\psi_{T}(R)}{\sum_{\\{R\\}}\psi_{T}(R)}\;.$
(3.24)

Statistical error bars on expectation values are then estimated by means of
block averages and the analysis of auto-correlations on data blocks. The
direct calculation of the expectation value (3.24) gives an exact result only
when $\mathcal{O}$ is the Hamiltonian $H$ or commutes with $H$, otherwise only
“mixed” matrix elements
$\langle\mathcal{O}\rangle_{m}\neq\langle\mathcal{O}\rangle$ can be obtained.
Among the different methods to calculate expectation values for operators that
do not commute with $H$, the extrapolation method [136] is the most widely
used. Following this method, one has a better approximation to the “pure”
(exact) value by means of a linear extrapolation

$\displaystyle\langle\mathcal{O}\rangle_{p}\simeq
2\,\frac{\langle\psi_{0}|\mathcal{O}|\psi_{T}\rangle}{\langle\psi_{0}|\psi_{T}\rangle}-\frac{\langle\psi_{T}|\mathcal{O}|\psi_{T}\rangle}{\langle\psi_{T}|\psi_{T}\rangle}=2\,\langle\mathcal{O}\rangle_{m}-\langle\mathcal{O}\rangle_{v}\;,$
(3.25)

or, if the operator $\mathcal{O}$ is positive defined, by means of

$\displaystyle\langle\mathcal{O}\rangle_{p}$
$\displaystyle\simeq\frac{\left(\frac{\langle\psi_{0}|\mathcal{O}|\psi_{T}\rangle}{\langle\psi_{0}|\psi_{T}\rangle}\right)^{2}}{\frac{\langle\psi_{T}|\mathcal{O}|\psi_{T}\rangle}{\langle\psi_{T}|\psi_{T}\rangle}}=\frac{\langle\mathcal{O}\rangle_{m}^{2}}{\langle\mathcal{O}\rangle_{v}}\;,$
(3.26)

where $\langle\mathcal{O}\rangle_{v}$ is the variational estimator. The
accuracy of the extrapolation method is closely related to the trial wave
function used in the variational calculation and on the accuracy of the DMC
sampling technique.

For a many-body system, if no constraint is imposed, $H$ has both symmetric
and antisymmetric eigenstates with respect to particle exchange. It can be
proven [215] that the lowest energy solution, and hence the state projected by
imaginary time propagation, is always symmetric. Moreover, in the DMC
algorithm, the walkers distribution is sampled through the wave function, that
must be positive defined in the whole configuration space for the
probabilistic interpretation to be applicable. The projection algorithm
described above is thus referred to Boson systems only. The extension for
Fermion systems is reported in § 3.1.2.

####  3.1.1 Importance Sampling

As discussed in the previous section, the basic version of the DMC algorithm
is rather inefficient because the weight term of Eq. (3.21) could suffer of
very large fluctuations. Indeed, because the Brownian diffusive process
ignores the shape of the potential, there is nothing that prevents two
particles from moving very close to each other, even in presence of an hard-
core repulsive potential.

The _importance function_ techniques [212, 213, 214] mitigates this problem by
using an appropriate importance function $\psi_{I}(R)$ (which is often, but
not necessarily, the same $\psi_{T}(R)$ used for the projection) to guide the
diffusive process. The idea is to multiply Eq. (3.9) by $\psi_{I}(R)$

$\displaystyle\psi_{I}(R)\psi(R,\tau+d\tau)=\int
dR^{\prime}\,G(R,R^{\prime},d\tau)\frac{\psi_{I}(R)}{\psi_{I}(R^{\prime})}\psi_{I}(R^{\prime})\psi(R^{\prime},\tau)\;,$
(3.27)

and define a new propagator

$\displaystyle\widetilde{G}(R,R^{\prime},d\tau)=G(R,R^{\prime},d\tau)\frac{\psi_{I}(R)}{\psi_{I}(R^{\prime})}\;,$
(3.28)

and a new function

$\displaystyle f(R,\tau)=\psi_{I}(R)\psi(R,\tau)\;,$ (3.29)

such that

$\displaystyle f(R,\tau+d\tau)=\int
dR^{\prime}\,\widetilde{G}(R,R^{\prime},d\tau)\,f(R^{\prime},\tau)\;.$ (3.30)

$f(R,\tau)$ represents the new probability density from which sample the
walker distribution. If $\psi_{I}(R)$ is suitably chosen, for example to be
small in the region where the potential presents the hard-core, then
$f(R,\tau)$ contains more information than the original $\psi(R,\tau)$, being
correlated to the potential by construction, and thus there is an improvement
in the quality of the DMC sampling and a reduction of the fluctuations of the
weight (3.21).

By inserting the new propagator $\widetilde{G}(R,R^{\prime},d\tau)$ in Eq.
(3.18) and expanding near $R^{\prime}$, it is possible to show (see for
instance Refs. [214]) that the integration gives an additional drift term in
$G_{0}(R,R^{\prime},d\tau)$

$\displaystyle
G_{0}(R,R^{\prime},d\tau)\rightarrow\widetilde{G}_{0}(R,R^{\prime},d\tau)=\left(\frac{1}{4\pi
Dd\tau}\right)^{\frac{3\mathcal{N}}{2}}\operatorname{e}^{-\frac{(R-R^{\prime}-v_{d}(R^{\prime})Dd\tau)^{2}}{4Dd\tau}}\;,$
(3.31)

where

$\displaystyle\bm{v}_{d}(R)=2\frac{\bm{\nabla}\psi_{I}(R)}{\psi_{I}(R)}\;,$
(3.32)

is a $3\mathcal{N}$ dimensional _drift velocity_ that drives the free
diffusion. The branching factor of Eq. (3.21) modifies in

$\displaystyle\omega_{i}\rightarrow\widetilde{\omega}_{i}=\operatorname{e}^{-\left(\frac{E_{L}(R)+E_{L}(R^{\prime})}{2}-E_{T}\right)d\tau}\;,$
(3.33)

where the potential energy is replaced by the _local energy_

$\displaystyle E_{L}(R)=\frac{H\psi_{I}(R)}{\psi_{I}(R)}\;.$ (3.34)

If the importance function is sufficiently accurate, the local energy remains
close to the ground-state energy throughout the imaginary time evolution and
the population of walkers is not subject to large fluctuations.

Going back to the imaginary time dependent Schrödinger equation, it is
possible to show (details can be found in Refs. [214]) that by multiplying Eq.
(3.8) by $\psi_{I}(R)$ we obtain a non homogenous Fokker-Plank equation for
$f(R,\tau)$

$\displaystyle\\!\\!-\frac{\partial}{\partial\tau}f(R,\tau)=$
$\displaystyle-\frac{\hbar^{2}}{2m}\nabla^{2}f(R,\tau)+\frac{\hbar^{2}}{2m}\bm{\nabla}\\!\cdot\\!\Bigl{[}\bm{v}_{d}(R)f(R,\tau)\Bigr{]}+E_{L}(R)f(R,\tau)\;,$
(3.35)

for which the corresponding Green’s function is given by the two terms of Eqs.
(3.31) and (3.33).

The DMC algorithm including the importance sampling procedure is still the
same described in § 3.1, where now the coordinates of the walkers are diffused
by the Brownian motion and guided by the drift velocity

$\displaystyle R=R^{\prime}+\xi+\bm{v}_{d}Dd\tau\;,$ (3.36)

and the branching process is given by the local energy through the weight
(3.33). The expectation values are still calculated by means of Eq. (3.24) but
now the sampling function $\psi(R,\tau)$ is replaced by $f(R,\tau)$.

####  3.1.2 Sign Problem

As discussed in § 3.1.1, the standard DMC algorithm applies to positive
defined wave function and the result of the imaginary time projection is a
nodeless function. The ground state of a Fermionic system is instead described
by an antisymmetric wave function, to which a probability distribution
interpretation cannot be given. Moreover, the search for an antisymmetric
ground state $|\psi_{0}^{A}\rangle$ corresponds to the search for an excited
state of the many-body Hamiltonian with eigenvalue

$\displaystyle E_{0}^{A}>E_{0}^{S}\;,$ (3.37)

where $E_{0}^{S}$ and $E_{0}^{A}$ are the ground state energies for the
Bosonic and the Fermionic system.

If no constraint is imposed, the Hamiltonian has both eigenstates that are
symmetric and antisymmetric with respect to particle exchange. We can thus
rewrite Eq. (3.4) by separating Bosonic and Fermionic components:

$\displaystyle|\psi(\tau)\rangle=\sum_{n=0}^{\infty}c_{n}^{S}\operatorname{e}^{-(E_{n}^{S}-E_{T})\tau}|\varphi_{n}^{S}\rangle+\sum_{n=0}^{\infty}c_{n}^{A}\operatorname{e}^{-(E_{n}^{A}-E_{T})\tau}|\varphi_{n}^{A}\rangle\;.$
(3.38)

If we want to naively apply the standard DMC algorithm to project out the
Fermionic ground state, we need to propagate the trial wave function for long
imaginary time taking $E_{0}^{A}$ as energy reference. If the Fermionic ground
state is not degenerate, i.e. $E_{n}^{A}\geq E_{n-1}^{A}>E_{0}^{A}$, in the
limit $\tau\rightarrow\infty$ we have

$\displaystyle\lim_{\tau\rightarrow\infty}|\psi(\tau)\rangle=\lim_{\tau\rightarrow\infty}\sum_{n}c_{n}^{S}\operatorname{e}^{-(E_{n}^{S}-E_{0}^{A})\tau}|\varphi_{n}\rangle+c_{0}^{A}|\varphi_{0}^{A}\rangle\;,$
(3.39)

where at least for $E_{0}^{S}$ the Bosonic part diverges due to the condition
(3.37). However, the exponentially growing component along the symmetric
ground state does not affect the expectation of the Hamiltonian. Indeed,
during the evaluation of the integral (3.23) on an antisymmetric trial wave
function $\psi_{T}^{A}(R)$, the symmetric components of $\psi(R,\tau)$ vanish
by orthogonality and in the limit of infinite imaginary time the energy
converges to exact eigenvalue $E_{0}^{A}$. However, the orthogonality
cancellation of the Bosonic terms does not apply to the calculation of the DMC
variance for the antisymmetric energy expectation value $\langle
E_{0}^{A}\rangle$

$\displaystyle\sigma^{2}_{E_{0}^{A}}=\left|\langle
H\rangle_{\psi_{T}^{A}}^{2}-\langle H^{2}\rangle_{\psi_{T}^{A}}\right|\;,$
(3.40)

where the second term diverges. We are left thus with an exact eigenvalue
affected by an exponentially growing statistical error. The signal to noise
ratio exponentially decays. This is the well known _sign problem_ and it
represents the main limit to the straightforward application of the DMC
algorithm to Fermion systems.

In order to extend the DMC method to systems described by antisymmetric wave
functions, it is possible to artificially split the configuration space in
regions where the trial wave function does not change sign. The multi
dimensional surface where the trial wave function vanishes, denoted as _nodal
surface_ , can be used to constrain the diffusion of the walkers: whenever a
walker crosses the nodal surface it is dropped from the calculation. In such a
way only the configurations that diffuse according to region of the wave
function with definite sign are taken into account. The problem reduces thus
to a standard DMC in the subsets of the configuration space delimited by the
nodal surface. This approximate algorithm is called _fixed-node_ [216, 212,
217] and it can be proven that it always provides an upper bound to the true
Fermionic ground state.

The sign problem appears for both real and complex antisymmetric wave
functions. The latter is the case of nuclear Hamiltonians. As proposed by
Zhang _et al._ [218, 219, 220], the _constrained path_ approximation can be
used to deal with the sign problem for complex wave functions. The general
idea is to constraining the path of walkers to regions where the real part of
the overlap with the wave function is positive. If we consider a complex
importance function $\psi_{I}(R)$, in order to keep real the coordinates space
of the system, the drift term in Eq. (3.31) must be real. A suitable choice
for the drift velocity is thus:

$\displaystyle\bm{v}_{d}(R)=2\frac{\bm{\nabla}\operatorname{Re}\left[\psi_{I}(R)\right]}{\operatorname{Re}\left[\psi_{I}(R)\right]}\;.$
(3.41)

Consistently, a way to eliminate the decay of the signal to noise ratio
consists in requiring that the real part of the overlap of each walker with
the importance function must keep the same sign

$\displaystyle\frac{\operatorname{Re}\left[\psi_{I}(R)\right]}{\operatorname{Re}\left[\psi_{I}(R^{\prime})\right]}>0\;,$
(3.42)

where $R$ and $R^{\prime}$ denote the coordinates of the system after and
before the diffusion of a time step. When this condition is violate, i.e. when
the overlap between the importance function and the walker after a diffusive
step changes sign, the walker is dropped. In these scheme, the ground state
expectation value of an observable $\mathcal{O}$ (Eq. (3.24)) is given by

$\displaystyle\langle\mathcal{O}\rangle=\frac{\sum_{\\{R\\}}\mathcal{O}\operatorname{Re}\left[\psi_{T}(R)\right]}{\sum_{\\{R\\}}\operatorname{Re}\left[\psi_{T}(R)\right]}\;.$
(3.43)

Another approach to deal with the complex sign problem is the _fixed phase_
approximation, originally introduced for systems whose Hamiltonian contains a
magnetic field [221]. Let write a complex wave function as

$\displaystyle\psi(R)=\left|\psi(R)\right|\operatorname{e}^{i\phi(R)}\;,$
(3.44)

where $\phi(R)$ is the phase of $\psi(R)$, and rewrite the drift velocity as

$\displaystyle\bm{v}_{d}(R)=2\frac{\bm{\nabla}\left|\psi_{I}(R)\right|}{\left|\psi_{I}(R)\right|}=2\operatorname{Re}\left[\frac{\bm{\nabla}\psi_{I}(R)}{\psi_{I}(R)}\right]\;.$
(3.45)

With this choice, the weight for the branching process becomes

$\displaystyle\widetilde{\omega}_{i}$
$\displaystyle=\exp\Bigg{\\{}-\Bigg{[}\frac{1}{2}\Bigg{(}-\frac{\hbar^{2}}{2m}\frac{\nabla^{2}|\psi_{I}(R)|}{|\psi_{I}(R)|}+\frac{V\psi_{I}(R)}{\psi_{I}(R)}$
$\displaystyle\quad-\frac{\hbar^{2}}{2m}\frac{\nabla^{2}|\psi_{I}(R^{\prime})|}{|\psi_{I}(R^{\prime})|}+\frac{V\psi_{I}(R^{\prime})}{\psi_{I}(R^{\prime})}\Bigg{)}-E_{T}\Bigg{]}d\tau\Bigg{\\}}\times\frac{\left|\psi_{I}(R^{\prime})\right|}{\left|\psi_{I}(R)\right|}\frac{\psi_{I}(R)}{\psi_{I}(R^{\prime})}\;,$
(3.46)

which is the usual importance sampling factor as in Eq. (3.33) multiplied by
an additional factor that corrects for the particular choice of the drift.
Using Eq. (3.44), the last term of the previous equation can be rewritten as

$\displaystyle\frac{\left|\psi_{I}(R^{\prime})\right|}{\left|\psi_{I}(R)\right|}\frac{\psi_{I}(R)}{\psi_{I}(R^{\prime})}=\operatorname{e}^{i[\phi_{I}(R)-\phi_{I}(R^{\prime})]}\;.$
(3.47)

The so called “fixed phase” approximation is then realized by constraining the
walkers to have the same phase as the importance function $\psi_{I}(R)$. It
can be applied by keeping the real part of the last expression. In order to
preserve the normalization, one has to consider an additional term in the
Green’s function due to the phase, that must be added to the weight:

$\displaystyle\exp\Bigg{[}-\frac{\hbar^{2}}{2m}\Bigl{(}\bm{\nabla}\phi(R)\Bigr{)}^{2}d\tau\Bigg{]}\;.$
(3.48)

This factor can be included directly in $\widetilde{\omega}_{i}$ considering
the following relation:

$\displaystyle\operatorname{Re}\left[\frac{\nabla^{2}\psi_{I}(R)}{\psi_{I}(R)}\right]=\frac{\nabla^{2}\left|\psi_{I}(R)\right|}{\left|\psi_{I}(R)\right|}-\Bigl{(}\bm{\nabla}\phi(R)\Bigr{)}^{2}\;.$
(3.49)

Thus, by keeping the real part of the kinetic energy in Eq. (3.46), the
additional weight term given by the fixed phase approximation is automatically
included. The calculation of expectation values is given now by

$\displaystyle\langle\mathcal{O}\rangle=\sum_{\\{R\\}}\operatorname{Re}\left[\frac{\mathcal{O}\psi_{T}(R)}{\psi_{T}(R)}\right]\;,$
(3.50)

i.e. by the evaluation of the real part of a local operator. This is of
particular interest for the technical implementation of the DMC algorithm. As
we will see in § 3.2.4, when dealing with Fermions the wave function can be
written as a Slater determinant of single particle states. It can be shown
(see Appendix A.2) that the evaluation of local operators acting on Slater
determinants can be efficiently implemented by means of the inverse matrix of
the determinant. The fixed phase approximation allows thus to deal with the
Fermion sign problem and also provides a natural scheme to implement the DMC
method. Moreover, the above derivation can be extended to operators other than
the kinetic energy. For example, when dealing with nuclear Hamiltonians like
(2.2), spin and isospin expectation values can be evaluated by taking the real
part of local spin and isospin operators calculated on the Slater determinant.
This is actually the standard way to treat the spin-isospin dependent
components of the nuclear Hamiltonian in the Auxiliary Field Diffusion Monte
Carlo (see § 3.2).

The constrained path and the fixed phase prescriptions are both approximations
introduced to deal with the sign problem for complex wave functions. In
principle they should yield similar results if the importance function is
close enough to the real ground state of the system. Accurate $\psi_{I}(R)$
are thus needed. An additional important observation is that the DMC algorithm
with the constrained path approximation does not necessarily provide an upper
bound in the calculation of energy [222, 223]. Moreover, it has not been
proven that the fixed phase approximation gives an upper bound to the real
energy. Thus, the extension of the DMC algorithm to Fermion systems described
by complex wave functions does not obey to the Rayleigh-Ritz variational
principle. Further details on the fixed node, constrained path and fixed phase
approximations can be found in the original papers and an exhaustive
discussion is reported in the Ph.D. thesis of Armani [224].

####  3.1.3 Spin-isospin degrees of freedom

If we want to study a nuclear many-body system described by the Hamiltonian
(2.2), we need to include also the spin-isospin degrees of freedom in the
picture. In order to simplify the notation, in the next with $A$ we will refer
to the number of nucleons. Starting from § 3.2.4 we will restore the
convention $A=\mathcal{N}_{N}+\mathcal{N}_{\Lambda}$. The typical trial many-
body wave function used in DMC calculation for nuclear systems takes the form
[136, 223]

$\displaystyle|\psi_{T}\rangle=\mathcal{S}\left[\prod_{i<j}\left(1+U_{ij}+\sum_{k}U_{ijk}\right)\right]\prod_{i<j}f_{c}(r_{ij})|\Phi_{A}\rangle\;,$
(3.51)

where $f_{c}(r_{ij})$ is a central (mostly short ranged repulsion)
correlation, $U_{ij}$ are non commuting two-body correlations induced by
$v_{ij}$ (that typically takes the same form of Eq. (2.16) for $p=2,\ldots,6$)
and $U_{ijk}$ is a simplified three-body correlation from $v_{ijk}$.
$|\Phi_{A}\rangle$ is the one-body part of the trial wave function that
determines the quantum numbers of the states and it is fully antisymmetric.
The central correlation is symmetric with respect to particle exchange and the
symmetrization operator $\mathcal{S}$ acts on the operatorial correlation part
of $|\psi_{T}\rangle$ in order to make the complete trial wave function
antisymmetric. The best trial wave function from which (3.51) is derived,
includes also spin-orbit and the full three-body correlations and it is used
in VMC calculations. See Refs. [225, 223].

Given $A$ nucleons ($Z$ protons, $A-Z$ neutrons), the trial wave function is a
complex vector in spin-isospin space with dimension
$\mathcal{N}_{S}\times\mathcal{N}_{T}$, where $\mathcal{N}_{S}$ is the number
of spin states and $\mathcal{N}_{T}$ the number of isospin states:

$\displaystyle\mathcal{N}_{S}=2^{A}\quad\quad\quad\mathcal{N}_{T}=\left(\begin{array}[]{c}A\\\
Z\end{array}\right)=\frac{A!}{Z!(A-Z)!}\;.$ (3.54)

For example, the wave function of an $A=3$ system has 8 spin components and,
considering the physical systems for $Z=1$ $\left({}^{3}\text{H}\right)$ or
$Z=2$ $\left({}^{3}\text{He}\right)$, 3 isospin states, thus a spin-isospin
structure with 24 entries. Using the notation of Ref. [136], we can write the
spin part of an $A=3$ wave function as a complex 8-vector (ignore
antisymmetrization)

$\displaystyle|\Phi_{A=3}\rangle=\left(\begin{array}[]{c}a_{\uparrow\uparrow\uparrow}\\\
a_{\uparrow\uparrow\downarrow}\\\ a_{\uparrow\downarrow\uparrow}\\\
a_{\uparrow\downarrow\downarrow}\\\ a_{\downarrow\uparrow\uparrow}\\\
a_{\downarrow\uparrow\downarrow}\\\ a_{\downarrow\downarrow\uparrow}\\\
a_{\downarrow\downarrow\downarrow}\end{array}\right)\quad\quad\text{with}\quad
a_{\uparrow\downarrow\uparrow}=\langle\uparrow\downarrow\uparrow|\Phi_{A=3}\rangle\;.$
(3.63)

The potentials ($v_{ij}$, $v_{ijk}$) and correlations ($U_{ij}$, $U_{ijk}$)
involve repeated operations on $|\psi_{T}\rangle$ but the many-body spin-
isospin space is closed under the action of the operators contained in the
Hamiltonian. As an example, consider the term
$\bm{\sigma}_{i}\cdot\bm{\sigma}_{j}$:

$\displaystyle\bm{\sigma}_{i}\cdot\bm{\sigma}_{j}$
$\displaystyle=2\left(\sigma_{i}^{+}\sigma_{j}^{-}+\sigma_{i}^{-}\sigma_{j}^{+}\right)+\sigma_{i}^{z}\sigma_{j}^{z}\;,$
$\displaystyle=2\,\mathcal{P}_{ij}^{\sigma}-1\;,$
$\displaystyle=\left(\begin{array}[]{cccc}1&0&0&0\\\ 0&-1&2&0\\\ 0&2&-1&0\\\
0&0&0&1\end{array}\right)\quad\text{acting
on}\quad\left(\begin{array}[]{c}\uparrow\uparrow\\\ \uparrow\downarrow\\\
\downarrow\uparrow\\\ \downarrow\downarrow\end{array}\right)\;.$ (3.72)

The $\mathcal{P}_{ij}^{\sigma}$ exchanges the spin $i$ and $j$, so the
operator $\bm{\sigma}_{i}\cdot\bm{\sigma}_{j}$ does not mix different isospin
components and acts on different, non contiguous, 4-element blocks of
$|\Phi_{A=3}\rangle$. For $i=2$ and $j=3$ we have for example:

$\displaystyle\bm{\sigma}_{2}\cdot\bm{\sigma}_{3}\,|\Phi_{A=3}\rangle=\left(\begin{array}[]{c}a_{\uparrow\uparrow\uparrow}\\\
2a_{\uparrow\downarrow\uparrow}-a_{\uparrow\uparrow\downarrow}\\\
2a_{\uparrow\uparrow\downarrow}-a_{\uparrow\downarrow\uparrow}\\\
a_{\uparrow\downarrow\downarrow}\\\ a_{\downarrow\uparrow\uparrow}\\\
2a_{\downarrow\downarrow\uparrow}-a_{\downarrow\uparrow\downarrow}\\\
2a_{\downarrow\uparrow\downarrow}-a_{\downarrow\downarrow\uparrow}\\\
a_{\downarrow\downarrow\downarrow}\end{array}\right)\;.$ (3.81)

The action of pair operators on $|\psi_{T}\rangle$, that are the most
computationally expensive, results thus in a sparse matrix of (non contiguous)
$4\times 4$ blocks in the $A$-body problem.

In the Green Function Monte Carlo, which slightly differs from the DMC in the
way the propagator is treated, each of the $2^{A}\frac{A!}{Z!(A-Z)!}$ spin-
isospin configurations undergoes to the imaginary time evolution of Eq. (3.9).
The propagation is now acting on the component $a_{\alpha}$, being $\alpha$
the spin-isospin index,

$\displaystyle a_{\alpha}(R,\tau+d\tau)=\sum_{\beta}\int
dR^{\prime}\,G_{\alpha\beta}(R,R^{\prime},d\tau)\,a_{\beta}(R^{\prime},\tau)\;,$
(3.82)

where the Green’s function is a matrix function of $R$ and $R^{\prime}$ in
spin-isospin space, defined as

$\displaystyle G_{\alpha\beta}(R,R^{\prime},d\tau)=\langle
R,\alpha|\operatorname{e}^{-(H-E_{T})d\tau}|R^{\prime},\beta\rangle\;.$ (3.83)

Due to the the factorial growth in the number of components of the wave
function, GFMC cannot deal with systems having a large number of nucleons,
like medium-heavy nuclei or nuclear matter. Standard GFMC calculations are
indeed limited up to 12 nucleons [17, 18, 19] or 16 neutrons [20].

###  3.2 Auxiliary Field Diffusion Monte Carlo

The AFDMC algorithm was originally introduced by Schmidt and Fantoni [29] in
order to deal in an efficient way with spin-dependent Hamiltonians. Many
details on the AFDMC method can be found in Refs. [37, 33, 226, 30, 31, 38,
224, 227, 214]. The main idea is to move from the many particle wave function
of the DMC or GFMC to a single particle wave function. In this representation,
going back to the example of the previous section, the spin part of an $A=3$
wave function becomes a tensor product of 3 single particle spin states
(ignore antisymmetrization):

$\displaystyle|\Phi_{A=3}\rangle$
$\displaystyle=\left(\begin{array}[]{c}a_{1\uparrow}\\\
a_{1\downarrow}\end{array}\right)_{1}\\!\\!\otimes\\!\left(\begin{array}[]{c}a_{2\uparrow}\\\
a_{2\downarrow}\end{array}\right)_{2}\\!\\!\otimes\\!\left(\begin{array}[]{c}a_{3\uparrow}\\\
a_{3\downarrow}\end{array}\right)_{3}\quad\text{with}\quad
a_{k\uparrow}=\,_{{}_{k}\,}\\!\langle\uparrow|\Phi_{A=3}\rangle\;.$ (3.90)

Taking also into account the isospin degrees of freedom, each single particle
state becomes a complex 4-vector and the total number of entries for
$|\Phi_{A=3}\rangle$ is thus 12, half of the number for the full DMC function
of Eq. (3.63). In the general case, the dimension of the multicomponent vector
describing a system with $A$ nucleons scale as $4A$. So, in this picture, the
computational cost for the evaluation of the wave function is drastically
reduced compared to the DMC-GFMC method when the number of particles becomes
large.

The problem of the single particle representation is that it is not closed
with respect to the application of quadratic spin (isospin) operators. As done
in the previous section (Eq. (3.81)), consider the operator
$\bm{\sigma}_{2}\cdot\bm{\sigma}_{3}=2\,\mathcal{P}_{23}^{\sigma}-1$ acting on
$|\Phi_{A=3}\rangle$:

$\displaystyle\bm{\sigma}_{2}\cdot\bm{\sigma}_{3}\,|\Phi_{A=3}\rangle$
$\displaystyle=2\left(\begin{array}[]{c}a_{1\uparrow}\\\
a_{1\downarrow}\end{array}\right)_{1}\\!\\!\otimes\\!\left(\begin{array}[]{c}a_{3\uparrow}\\\
a_{3\downarrow}\end{array}\right)_{2}\\!\\!\otimes\\!\left(\begin{array}[]{c}a_{2\uparrow}\\\
a_{2\downarrow}\end{array}\right)_{3}$ (3.97)
$\displaystyle\phantom{=}-\left(\begin{array}[]{c}a_{1\uparrow}\\\
a_{1\downarrow}\end{array}\right)_{1}\\!\\!\otimes\\!\left(\begin{array}[]{c}a_{2\uparrow}\\\
a_{2\downarrow}\end{array}\right)_{2}\\!\\!\otimes\\!\left(\begin{array}[]{c}a_{3\uparrow}\\\
a_{3\downarrow}\end{array}\right)_{3}\;.$ (3.104)

There is no way to express the result as a single particle wave function of
the form (3.90). At each time step, the straightforward application of the DMC
algorithm generates a sum of single particle wave functions. The number of
these functions will grows very quickly during the imaginary time evolution,
destroying the gain in computational time obtained using a smaller
multicomponent trial wave function.

In order to keep the single particle wave function representation and overcome
this problem, the AFDMC makes use of the Hubbard-Stratonovich transformation

$\displaystyle\operatorname{e}^{-\frac{1}{2}\lambda\mathcal{O}^{2}}=\frac{1}{\sqrt{2\pi}}\int\\!dx\operatorname{e}^{-\frac{x^{2}}{2}+\sqrt{-\lambda}\,x\mathcal{O}}\;,$
(3.105)

to linearize the quadratic dependence on the spin-isospin operators by adding
the integration over a new variable $x$ called _auxiliary field_. It is indeed
possible to show that the single particle wave function is closed with respect
to the application of a propagator containing linear spin-isospin operators at
most:

$\displaystyle\operatorname{e}^{-\mathcal{O}_{j}d\tau}|\Phi_{A}\rangle$
$\displaystyle=\operatorname{e}^{-\mathcal{O}_{j}d\tau}\bigotimes_{i}\left(\begin{array}[]{c}a_{i\uparrow}\\\
a_{i\downarrow}\end{array}\right)_{i}\;,$ (3.108)
$\displaystyle=\left(\begin{array}[]{c}a_{1\uparrow}\\\
a_{1\downarrow}\end{array}\right)_{1}\\!\otimes\cdots\otimes\,\operatorname{e}^{-\mathcal{O}_{j}d\tau}\left(\begin{array}[]{c}a_{j\uparrow}\\\
a_{j\downarrow}\end{array}\right)_{j}\\!\otimes\cdots\otimes\left(\begin{array}[]{c}a_{A\uparrow}\\\
a_{A\downarrow}\end{array}\right)_{A}\;,\quad\quad$ (3.115)
$\displaystyle=\left(\begin{array}[]{c}a_{1\uparrow}\\\
a_{1\downarrow}\end{array}\right)_{1}\\!\otimes\cdots\otimes\left(\begin{array}[]{c}\widetilde{a}_{j\uparrow}\\\
\widetilde{a}_{j\downarrow}\end{array}\right)_{j}\\!\otimes\cdots\otimes\left(\begin{array}[]{c}a_{A\uparrow}\\\
a_{A\downarrow}\end{array}\right)_{A}\;,$ (3.122)

where, working on 2-component spinors, $\mathcal{O}_{j}$ can be a $2\times 2$
spin or isospin matrix. If we are dealing with the full 4-component spinor,
$\mathcal{O}_{j}$ can be an extended $4\times 4$ spin, isospin or
isospin$\,\otimes\,$spin matrix. To get this result we have used the fact that
the operator $\mathcal{O}_{j}$ is the representation in the $A$-body tensor
product space of a one-body operator:

$\displaystyle\mathcal{O}_{j}\equiv\mathbb{I}_{1}\otimes\cdots\otimes\mathcal{O}_{j}\otimes\cdots\otimes\mathbb{I}_{A}\;.$
(3.123)

Limiting the study to quadratic spin-isospin operators and making use of the
Hubbard-Stratonovich transformation, it is thus possible to keep the single
particle wave function representation over all the imaginary time evolution.
This results in a reduced computational time for the propagation of the wave
function compared to GFMC, that allows us to simulate larger systems, from
medium-heavy nuclei to the infinite matter. In the next we will see in detail
how the AFDMC works on the Argonne V6 like potentials (§ 3.2.1), and how it is
possible to include also spin-orbit (§ 3.2.2) and three-body (§ 3.2.3) terms
for neutron systems. Finally the extension of AFDMC for hypernuclear systems
(§ 3.2.5) is presented.

####  3.2.1 Propagator for nucleons: $\bm{\sigma}$,
$\bm{\sigma}\cdot\bm{\tau}$ and $\bm{\tau}$ terms

Consider the first six components of the Argonne $NN$ potential of Eq. (2.16).
They can be conveniently rewritten as a sum of a spin-isospin independent and
a spin-isospin dependent term

$\displaystyle
V_{NN}=\sum_{i<j}\sum_{p=1,6}v_{p}(r_{ij})\,\mathcal{O}_{ij}^{\,p}=V_{NN}^{SI}+V_{NN}^{SD}\;,$
(3.124)

where

$\displaystyle V_{NN}^{SI}$ $\displaystyle=\sum_{i<j}v_{1}(r_{ij})\;,$ (3.125)

and

$\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}+\frac{1}{2}\sum_{i\neq
j}\sum_{\alpha\beta\gamma}\tau_{i\gamma}\,\sigma_{i\alpha}\,A_{i\alpha,j\beta}^{[\sigma\tau]}\,\sigma_{j\beta}\,\tau_{j\gamma}\;$
$\displaystyle\,+\frac{1}{2}\sum_{i\neq
j}\sum_{\gamma}\tau_{i\gamma}\,A_{ij}^{[\tau]}\,\tau_{j\gamma}\;.$ (3.126)

The $3A\times 3A$ matrices $A^{[\sigma]}$, $A^{[\sigma\tau]}$ and the $A\times
A$ matrix $A^{[\tau]}$ are real and symmetric under Cartesian component
interchange $\alpha\leftrightarrow\beta$, under particle exchange
$i\leftrightarrow j$ and fully symmetric with respect to the exchange
$i\alpha\leftrightarrow j\beta$. They have zero diagonal (no self interaction)
and contain proper combinations of the components of AV6 (Latin indices are