# The Hyperspherical Harmonics Method: A Tool for Testing and Improving Nuclear Interaction Models

^{1}Department of Physics “E. Fermi”, University of Pisa, Pisa, Italy^{2}Istituto Nazionale di Fisica Nucleare, Sezione di Pisa, Pisa, Italy^{3}Physique Quantique, and Physique Nucléaire Théorique et Physique Mathématique, Université libre de Bruxelles, Brussels, Belgium^{4}Department of Mathematics and Physics “E. De Giorgi”, University of Salento, Lecce, Italy^{5}Istituto Nazionale di Fisica Nucleare, Sezione di Lecce, Lecce, Italy^{6}Gran Sasso Science Institute, L'Aquila, Italy

The Hyperspherical Harmonics (HH) method is one of the most accurate techniques to solve the quantum mechanical problem for nuclear systems with a number of nucleons *A* ≤ 4. In particular, by applying the Rayleigh-Ritz or Kohn variational principle, both bound and scattering states can be addressed, using either local or non-local interactions. Thanks to this versatility, the method can be used to test the two- and three-nucleon components of the nuclear interaction. In the present review we introduce the formalism of the HH method, both for bound and scattering states. In particular, we describe the implementation of the method to study the *A* = 3 and 4 scattering problems. Second, we present a selected choice of results of the last decade, most representative of the latest achievements. Finally, we conclude with a discussion of what we believe will be the most significant developments within the HH method for the next 5–10 years.

## 1. Introduction

The “standard” picture of a nucleus sees it as a system of *A* nucleons, protons or neutrons, interacting among themselves and eventually with external electroweak probes. The interaction between nucleons, i.e., the nuclear interaction, is the subject of the Research Topic of which this contribution is part. Using the nucleon as the relevant degree of freedom, the system is described by the nuclear Hamiltonian, written as

where the first term is the (non-relativistic) kinetic energy (in the center-of-mass reference frame), *m*_{i} being the *i*th nucleon mass, *V*_{ij} and *V*_{ijk} are, respectively, the two- and three-nucleon interactions, i.e., the interaction between each *ij* pair or *ijk* triple. It has been shown in several studies (for recent ones see references [1, 2]) that even the nuclear systems with medium-large values of *A* can be at least qualitatively described including *V*_{ij} and *V*_{ijk} only: essentially, it seems to be little room [3] for four- or more-nucleon interactions [the dots of Equation (1)]. Therefore, we will neglect from now on the contributions beyond three-nucleon interaction.

There exists a large variety of realistic models for *V*_{ij} and *V*_{ijk}. The most important ones are presented in this Research Topic. These models are very different among themselves, as they can be local, or minimally non-local and expressed in coordinate space, or non-local and expressed in momentum space. Some older models were derived phenomenologically or in a meson-theoretical approach, as the Argonne *v*_{18} [4] or the CDBonn [5] potentials, but, since the seminal work of Weinberg [6], the preferred framework to derive the nuclear interaction is chiral effective field theory. Since the presentation of the different models is assigned to this Research Topic, here we only mention that all the models have a common characteristic: the *V*_{ij} and *V*_{ijk} interactions have an intricate operatorial structure. As a consequence, the solution of the Schrödinger equation for *A* > 2 becomes a challenging problem. The methods which are able to solve the *A*-body quantum mechanical problem without approximations, or with approximations which can be maintained under control, are the so-called *ab initio* methods^{1}. They are clearly essential in order to test a given model for the nuclear interaction against experiment. It is fundamental for these methods to be sufficiently accurate to capture the small details introduced by the complexity of the interaction. As an example, we can quote the case of the triton binding energy. It is nowadays well-known that the triton binding energy calculated just retaining *V*_{ij} in Equation (1) is underestimated by 0.5−1 MeV, depending of the considered model. It is straightforward that the required accuracy of the *ab initio* method to catch this disagreement must be better that 500 keV. Nowadays, the methods for the *A* = 3 bound systems have reached a much higher accuracy, of the order of 1 keV, or even better. And therefore, the presence of *V*_{ijk} is not anymore under discussion. To be noted that (i) all models for the two-nucleon interaction are phase-equivalent, and (ii) each model for *V*_{ijk} is built in conjunction with a given model for *V*_{ij}, and therefore two- and three-nucleon interactions are linked to each other and cannot be uniquely defined.

There are several *ab initio* methods which can solve the *A*-body quantum mechanical problem in different regions of the nuclear chart. A recent review is given in reference [7]. Here we limit ourselves to mention the methods based on Monte Carlo techniques, as the variational Monte Carlo (VMC) or the Green's function Monte Carlo (GFMC) methods (see reference [8], and references therein). There are then the methods linked to the shell model, as the no-core shell model (NCSM) or the realistic shell model (RSM) (see references [9–11]), respectively. All these methods are quite powerful to study medium-mass nuclear bound states, but less accurate, apart from the GFMC and NCSM, for very light nuclei, as those with *A* = 3, 4. Furthermore, their extension to the scattering systems is not so trivial, and, in some cases, still not at reach.

Restricting ourselves to the *A* = 3, 4 nuclear systems, both bound and scattering states, we have at hand very few accurate *ab initio* methods, i.e., the Faddeev (Faddeev-Yakubovsky for *A* = 4) equations (FE) technique, solved in coordinate- or in momentum-space, the method based on the Alt-Grassberger-Sandhas (AGS) equations solved in momentum space, and the Hyperspherical Harmonics (HH) method presented here. We refer the reader to references [12, 13] for the FE method in coordinate space, to references [14, 15] for the FE method in momentum space, to references [16, 17] for recent reviews on the AGS method. Clearly, each method has advantages and drawbacks. For instance, the FE method in momentum space can be applied to *A* = 3, 4 bound and scattering states in a wide energy range. However, the inclusion of the Coulomb interaction for charged particle scattering states is quite problematic. The FE method in coordinate space can handle the Coulomb interaction, but it has not yet been applied to scattering problems at very low-energy, and it has been applied only recently to study systems with larger *A* values [18]. It is though a method with in principle great possibilities of extension [13]. The AGS method, although working in momentum space, can handle the Coulomb interaction and can be applied to a large variety of *A* = 3, 4 scattering states, in a wide energy range. However, the very low energy range, that of interest to nuclear astrophysics, i.e., below about 100 keV, is still not accessible with the AGS method. The method has also not been applied for *A* > 4 yet.

The HH method has a long history, summarized in the introduction of reference [19]. We will concentrate here on the latest developments, essentially those obtained since 2008, year of publication of reference [19]. However, to fully appreciate the major developments of this last decade, it is necessary to briefly outline the state-of-the-art of the HH method at that time. The HH method in 2008 was extensively used by two groups, one formed by some of the present authors, and referred to as the Pisa group, and the other one including, among others, N. Barnea, W. Leidemann, and G. Orlandini. The latter has developed over the years the so-called effective interaction HH method (EIHH), which uses the Suzuki-Lee approach [20–22] to significantly reduce the number of the basis functions needed in the expansion. The method has been applied to study the *A* ≤ 4 nuclear bound systems in references [23, 24] using local realistic interactions, and had been pushed to work up to *A* = 6 with central potential models [23]. In the last decade, the EIHH method has been updated in order to work also with the most recent non-local potential models [25]. Furthermore, the EIHH method has been extensively used in testing the nuclear interaction models, using reactions between light nuclei and electromagnetic probes. For example, a test of the interaction models has been performed studying the 0^{+} resonance of ^{4}He in ^{4}He(*e, e*′), where a large potential model dependence has been found [26, 27]. In this review, however, we have decided to leave out the large body of results for the electromagnetic reactions on light nuclei, which would deserve a review all by itself (see reference [28]).

The HH method as developed by the Pisa group existed in 2008 in two flavors: the correlated HH method, including a pair-correlation function (pair-correlated HH method—PHH) or with a Jastrow type factor (correlated HH method—CHH), and the “pure” HH method. The correlation factor was introduced to describe correlations induced by the strong repulsion of the interaction at short range. The correlation factor describes the particular configuration in which two particles are close to each other and goes to unity for large pair relative distances. Therefore the HH expansion has to take care of reconstructing the wave function outside this range, making the convergence of the expansion much faster. The drawbacks of the PHH and CHH methods are (i) the necessity of performing numerical integrations, which would be instead analytical without correlation factors, reducing the accuracy of the method in the *A* = 4 case; (ii) the not simple extension of the PHH method to work in momentum space. Therefore, it is difficult to apply the PHH method with the non-local potentials mentioned above. This has motivated our group, together with the continuous increasing of computing power, to return to the “pure” HH method. Up to the year 2008, this had been developed and applied to study with great accuracy the *A* = 3, 4 bound states, with both local potentials, expressed in coordinate space, or non-local ones, given in momentum-space. In fact, while the local interactions had been at reach for the HH method from the very beginning [29], the non-local ones were a recent achievement at that time [30]. In 2008, the zero-energy *A* = 3, 4 scattering states were also calculated with both local and non-local interactions [19]. The higher energy scattering states, still below the breakup threshold of the target nucleus, had been studied for both *A* = 3 and 4 systems only with local interactions, in a variety of contributions extensively mentioned in reference [19]. What was still missing in 2008 was the study of the *A* = 3, 4 scattering states, still below the target breakup threshold, with non-local potentials. This has been obtained in references [31–35] for both *A* = 3 and 4, and it is in fact one of the main achievements of the HH method in the last decade. The HH method, in its PHH version, has been applied, including the full electromagnetic interaction, to describe elastic scattering observables in *A* = 3 above the deuteron breakup threshold [36] and in wide energy region [37]. Preliminary studies of the method to treat the breakup channels, as for instance to the process *n* + *d* → *n* + *n* + *p*, can be found in references [38–40]. The application of the method using the Hamiltonian defined in Equation (1) is in progress. In progress is also the further development of the method toward *A* > 4. This has been performed within the so-called non-symmetrized HH method [41] with central potentials, or, as mentioned above, by the EIHH method. The first steps to use the HH method without the Suzuki-Lee procedure have been shown in references [42, 43], and intense research activity is currently underway. The formalism which is presented here is in fact quite general, and can be applied also to the *A* = 5, 6 nuclear systems.

Before concluding this section with the outline of this contribution, we would like to make few remarks: (i) the HH method is extremely powerful, and its application to systems up to *A* ~ 7, 8 is limited essentially by computing power. (ii) The accuracy of the HH method has been tested in a number of benchmark calculations. In particular we quote the benchmark on the *A* = 3 [44] and *A* = 4 [45] bound states, on the *nd* and *pd* scattering phase shifts [46, 47], and, in the last decade, on the *A* = 4 scattering states [34, 48]. (iii) Compared with the other *ab initio* methods, the HH technique seems to be one of the best choices to study low-energy scattering states, in order to obtain accurate predictions for nuclear reactions of astrophysical interest [49, 50].

The present review is organized as follows: in section 2 we discuss the HH formalism, both for bound and scattering states. We will try to keep a somewhat “pedagogical” level, in order to allow the interested reader to perform his/her own algebraic steps and eventually reproduce the already existing results. In section 3 we discuss the most important results obtained within the HH method since the year 2008. In particular, we will show that the HH method has reached such a degree of accuracy for both bound and scattering states, that it has been used in order to construct an accurate model for the three-nucleon interaction, with a procedure similar, in principle, to the one used to derive the nowadays very accurate two-nucleon interaction models. Finally, in section 4 we will give some concluding remarks and an outlook.

## 2. The HH Formalism

We review in this section the HH formalism, focusing in particular on the three- and four-body systems, both bound and scattering states. The approach described below can be used in conjunction with both local and non-local two-nucleon interactions. At present, the method works with only local three-nucleon interactions, but its extension to the non-local case presents no conceptual difficulties.

### 2.1. Hyperspherical Harmonic Functions

Let us consider a system of *A* particles with masses *m*_{1}, …, *m*_{A} and spatial coordinates *r*_{1}, …, *r*_{A}, respectively. For separating the internal and center-of-mass (c.m.) motion, it is convenient to introduce another set of coordinates made of *N* = *A* − 1 internal Jacobi coordinates *x*_{1}, …, *x*_{N} and the c.m. coordinate ** X** defined by

where $M=\sum _{i=1}^{A}{m}_{i}$ is the total mass of the system. There are several definitions of the Jacobi coordinates, but a convenient one which will be used through this work is the following

where *m* is a reference mass, ${M}_{j}=\sum _{i=1}^{j}{m}_{i}$, and *j* = 1, …, *N*. In the case where all the particles have the same mass *m*, Equation (3) reduces to

From a given choice of the Jacobi vectors, the hyperspherical coordinates (ρ, Ω_{N}) can be introduced. The hyperradius ρ is defined by

where *x*_{i} is the modulus of the Jacobi vector *x*_{i}. The hyperradius ρ is symmetric with respect to particle exchanges and does not depend on the particular choice of Jacobi coordinates. The set Ω_{N} of hyperangular coordinates,

is made of the angular parts ${\widehat{x}}_{i}=({\theta}_{i},{\varphi}_{i})$ of the spherical components of the Jacobi vectors *x*_{i}, with *i* = 1, …*N*, and of the hyperangles φ_{j}, defined by

where 0 ≤ φ_{j} ≤ π/2 and *j* = 2, …, *N*.

The advantage of using the hyperspherical coordinates can be appreciated noting that the internal kinetic energy operator of the *A*-body system can be decomposed as

where the operator ${\Lambda}_{N}^{2}({\Omega}_{N})$ is the so-called grand-angular momentum operator. Its explicit expression can be found, for instance, in references [19, 51], but it is not essential for our purposes. More important are the eigenfunctions of the grand-angular momentum ${\Lambda}_{N}^{2}({\Omega}_{N})$, the so-called hyperspherical harmonics (HH). They can be defined as

Here ${Y}_{{l}_{i}}({\widehat{x}}_{i})$ is a spherical harmonic function for *i* = 1, …, *N*, *L* is the total orbital angular momentum, *M*_{L} its projection on the *z* axis, and

with *n*_{1} = 0, *j* = 1, …, *N*, and *K*_{N} ≡ *K* is the so-called grand-angular momentum. The notation [*K*] stands for the collection of all the quantum numbers [*l*_{1}, …, *l*_{N}, *L*_{2}, …, *L*_{N − 1}, *n*_{2}, …, *n*_{N}]. The functions ${{}^{(j)}{P}}_{{n}_{j}}^{{K}_{j-1},{l}_{j}}({\phi}_{j})$ in Equation (9) are defined by

where ${P}_{{n}_{j}}^{{\nu}_{j-1},{l}_{j}+1/2}(cos2{\phi}_{j})$ are Jacobi polynomials [52], with

and the normalization factors ${{N}}_{n}^{l,\nu}$ are given by

with Γ indicating the standard Gamma function [52]. To be noticed that for *A* = 3, *j* = 1, 2, and since *n*_{1} = 0, there is only one index *n*_{2} ≡ *n*. In this case, *K* = *l*_{1} + *l*_{2} + 2*n*. For the convergence on *n* or alternatively *K*, see the discussion in section 3. With the definition of Equation (9), the HH functions are eigenvectors of the grand-angular momentum operator ${\Lambda}_{N}^{2}({\Omega}_{N})$, the square of the total orbital angular momentum ** L**, its

*z*component

*L*_{z}, and the parity operator Π. Therefore we have

We remark here two useful properties of the HH functions. First of all, the HH functions are orthonormal with respect to the volume element *dΩ*_{N}, i.e.,

with

and

Therefore, the number of HH functions for a given *K* increases fast with *K*, but is always finite. In fact, according with Equation (10), $K=\sum _{i=1,N}({l}_{i}+2{n}_{i})$. Furthermore, independently of the specific choice of Jacobi coordinates used to define the hyperspherical ones or of the order of the coupling of the spherical harmonics in Equation (9), the HH functions constitute a complete basis.

Secondly, in order to evaluate matrix elements of a given many-body operator between HH functions, it is often useful to determine the effect of a particles permutation on an HH function. Since the grand-angular and the total orbital angular momenta are fully symmetric, and since the HH functions constitute a complete basis, the permuted HH functions ${{Y}}_{\left[K\right]}^{KL{M}_{L}}({\Omega}_{N}^{p})$ can be written as linear combinations of unpermuted HH functions ${{Y}}_{\left[{K}^{\prime}\right]}^{KL{M}_{L}}({\Omega}_{N})$ with same *K*, *L*, and *M*_{L} values. Therefore, we can write

The transformation coefficients ${a}_{\left[K\right];\left[{K}^{\prime}\right]}^{KL,p}$ do not depend on the quantum number *M*_{L}. For *A* = 3, they are called the Raynal-Revai coefficients [53]. To be remarked that $\left[{K}^{\prime}\right]\equiv \left[{l}_{1}^{\prime},\dots ,{l}_{N}^{\prime},{L}_{2}^{\prime},\dots ,{L}_{N-1}^{\prime},{n}_{2}^{\prime},\dots ,{n}_{N}^{\prime}\right]$, but such that *K*′ = *K*. Note that *L* is conserved. For *A* > 3, see references [42, 54].

Let us consider more specifically a system of *A* nucleons described within the isospin formalism. The *A*-nucleon wave function contains spatial, spin, and isospin parts. We can define the spin functions ${\chi}_{\left[S\right]}^{S{M}_{S}}$ with total spin *S* and total spin projection *M*_{S} and the isospin functions ${\xi}_{\left[T\right]}^{T{M}_{T}}$ with total isospin *T* and total isospin projection *M*_{T} by coupling the individual spin functions χ_{1/2, ±1/2} or isospin functions ξ_{1/2, ±1/2}, respectively, of each nucleon, as

So now [*S*] stands for [*S*_{2}, …, *S*_{N − 1}] and [*T*] for [*T*_{2}, …, *T*_{N − 1}].

Including the spin and isospin functions, the HH basis functions read

where *J* is the total angular momentum, *J*_{z} its projection, and [*KST*] stands for [*K*][*S*][*T*]. To be noticed that also the spin-isospin part of ${Y}_{\left[KST\right]}^{KLSJ{J}_{z}T{M}_{T}}({\Omega}_{N})$ constructed with a given ordering of the particles, can be rewritten in terms of a different permutation, using the Wigner 6j coefficients [55].

We conclude by noting that the HH functions can also be built in momentum space instead of configuration space. They can be obtained by replacing the hyperspherical coordinates (ρ, Ω_{N}) associated with the Jacobi coordinates {*x*_{i}}_{i = 1, …, N} by the hyperspherical coordinates $(Q,{\Omega}_{N}^{(q)})$ associated with the *N* Jacobi conjugate momenta {*q*_{i}}_{i = 1, …, N}. The rest of the formalism remains unchanged. For more details, see references [19, 30, 56].

### 2.2. The HH Method for *A* = 3 and 4

We discuss in some detail the method for systems with *A* = 3, 4 nucleons within the isospin formalism for both bound and scattering states in sections 2.2.1 and 2.2.2, respectively. The extension to *A* > 4 is straightforward, but leads to more lengthy expressions.

#### 2.2.1. The *A* = 3 and 4 Bound States

The wave function of an *A*-body bound state, with *A* = 3, 4, having total angular momentum *J, J*_{z} and parity π, and third component of the total isospin *M*_{T}, can be decomposed as a sum of Faddeev-like amplitudes as:

Here the sum on *p* runs up to *N*_{p} = 3 or 12 even permutations of the *A* particles, with *A* = 3 or 4, respectively, and the coordinates ${x}_{1}^{(p)},\cdots \phantom{\rule{0.3em}{0ex}},{x}_{N}^{(p)}$ are the Jacobi coordinates as defined in Equation (3), but here we show explicitly the dependence on the permutation *p*. To be noticed that, increasing the number of particles, different arrangements of them in sub-clusters allow for different definitions of the Jacobi coordinates. For example, in *A* = 4 two different sets exist corresponding to have a 3+1 or a 2+2 asymptotic configuration. However in the sub-space defined by the grand-angular momentum *K*, HH functions defined in different sets of Jacobi coordinates result to be linearly dependent. In the following we always refer to the set defined in Equation (3).

The coordinate-space hyperspherical coordinates are given in Equations (5)–(7), and the hyperangular variables are φ_{2} for *A* = 3 and φ_{2}, φ_{3} for *A* = 4.

We rewrite here the HH basis of Equation (24) for the *A* = 3 and 4 cases. Historically, the angular, spin and isospin quantum numbers have been collected in the so-called channels α, defined explicitly by

so that we can write

for *A* = 3, and

for *A* = 4. To be noticed that, in order to ensure the antisymmetry of the wave function, the Faddeev-like amplitudes have to change sign under exchange of particle 1 and 2. Therefore, the sum *l*_{2α} + *S*_{aα} + *T*_{aα} for *A* = 3 and *l*_{3α} + *S*_{aα} + *T*_{aα} for *A* = 4 must be odd. Furthermore, *l*_{1α} + *l*_{2α} for *A* = 3 and *l*_{1α} + *l*_{2α} + *l*_{3α} for *A* = 4 must be an even or odd number in correspondence to a positive or negative parity state. Even with these restrictions, there is an infinite number of channels. However, the contributions of the channels with higher and higher values for *l*_{1α} + *l*_{2α} for *A* = 3 and *l*_{1α} + *l*_{2α} + *l*_{3α} for *A* = 4 should become less and less important, due to the centrifugal barrier. Therefore, it is found that the number of channels with a significant contribution is relatively small for bound and low-energy scattering states. The most important ones for *A* = 3 and for *A* = 4 are listed, respectively, in Tables 1, 2 of reference [19].

By using Equations (28) and (29), the *A*-body wave function Ψ_{A} of Equation (25) can be written in coordinate-space as

for *A* = 3, and

for *A* = 4. The sum over *n*_{2} in Equation (30) and *n*_{2}, *n*_{3} in Equation (31) is restricted to independent states, see below. The hyperradial functions *u*_{αn2}(ρ) (*u*_{αn2n3}(ρ) for *A* = 4) are themselves expanded in terms of known functions. It is common to use Laguerre polynomials multiplied by an exponential function, as they have been found to give a nice convergence of this expansion. Therefore,

where the sum is truncated at *N*_{L}, and the functions *f*_{m}(ρ) are written as

Here *D* ≡ 3*N* − 1, ${L}_{m}^{(D-1)}(\gamma \rho )$ is a Laguerre polynomial [52], and γ is a non-linear parameter, to be variationally optimized. The exponential factor e^{−γρ/2} ensures that *f*_{m}(ρ) → 0 for ρ → ∞. The optimal value of γ depends on the potential model, and it is typically in the interval 2.5–4.5 fm^{−1} for local and 4–8 fm^{−1} for non-local potentials. Also *N*_{L} depends on the potential models, but typically with *N*_{L} ~ 20 − 30 a convergence at the 1 keV (10 keV) level for the *A* = 3 (*A* = 4) binding energies is achieved.

When working in momentum space, the *A*-body wave function Ψ_{A} is written as in Equations (30) and (31), with *u*_{αn2}(ρ) and *u*_{αn2n3}(ρ) replaced with *w*_{αn2}(*Q*) and *w*_{αn2n3}(*Q*), i.e., functions of the hypermomentum *Q*, while the HH functions are expressed in terms of conjugate Jacobi momenta. The *w*-functions are related to the *u*-functions as

where ${J}_{K+\frac{D}{2}-1}(Q\rho )$ are Bessel functions of the first kind [19], and *K* is again the grand-angular momentum.

At the end, the *A*-body wave function of Equations (30)–(34) can be cast in the form

where

in coordinate-space (a similar expression holds in momentum-space). The decomposition proposed in Equation (25) ensures the complete antisymmetrization of the state through the sum on the permutations as indicated in Equation (36). Indeed, the hyperangular-spin-isospin basis state |*K, m*〉 is completely antisymmetric. However, the sum over the permutations for fixed values of *K* produces linear dependent states that have to be individuated and eliminated from the basis set [42, 54, 57]. This procedure could be delicate from a numerical point of view as the number of *K* increases. In such a case, one needs a robust orthonormal procedure capable to deal with the presence of large numerical cancellations. However, if one is successful in this step, at the end one can work with a basis of independent antisymmetrical states, whose number is noticeably less than the degeneracy of the full HH basis. Attempts to use the HH basis without symmetrization has been recently proposed [41]. The idea is then to use the complete HH basis in which all symmetries are represented to describe a particular state. The diagonalization of the Hamiltonian produces eigenvectors with well-defined permutation properties reflecting the symmetries in it. Different applications followed this procedure for bosons as well as for fermions (see references [41, 58–61]). The advantage of eliminating the orthonormalization of the states has to be balanced by the fact that in this case one has to work with the full basis of HH functions, whose degeneracy rapidly increases with *K* and the number of particles *A*.

Once the antisymmetric state |*K, m*〉 is constructed, what is left is to obtain the unknown coefficients *c*_{K; m} of the expansion. In order to do so, we apply the Rayleigh-Ritz variational principle, which states that the quantity 〈Ψ_{A}|*H* − *E*|Ψ_{A}〉 is stationary with respect to the variation of any unknown coefficient. Here *H* is the nuclear Hamiltonian and *E* = −*B* the energy of the state, which, in the case of a bound state, is negative and opposite to the binding energy *B*.

When differentiating respect to *c*_{K; m} we obtain the following equation

where the matrix elements of the Hamiltonian *H* and of the identity operator can be calculated with standard numerical techniques (see reference [19] for more details). Equation (37) represents a generalized eigenvalue-eigenvector problem, which can be solved with a variety of numerical algorithms. Widely used within the HH method is the Lanczos algorithm [62], since the HH basis can become quite large (up to about 10,000 terms for *A* = 3 and about one order of magnitude larger for *A* = 4 are used in practice).

The results obtained solving Equation (37) for a variety of nuclear interaction models will be presented in section 3.

#### 2.2.2. The *A* = 3 and 4 Scattering States

The HH method has been also applied to the scattering problem. In particular, the method can study the elastic *N* + *Y* → *N* + *Y* process, where *N* is a nucleon and *Y* a bound system (*A*_{Y} + 1 ≡ *A* = 3, 4), both below and above the *Y* nucleus breakup threshold. The extension of the HH method to the full breakup problem, i.e., for *A* = 3 the process *n* + *d* → *n* + *n* + *p*, is currently underway and will not be discussed here.

The wave function ${\text{\Psi}}_{NY}^{LSJ{J}_{z}}$ describing the *N* + *Y* scattering state with incoming orbital angular momentum *L*, channel spin $\overrightarrow{S}\equiv \overrightarrow{\frac{1}{2}}+{\overrightarrow{S}}_{Y}$, parity π = (−1)^{L}, and total angular momentum *J, J*_{z}, is written as

Here we have introduced ${\text{\Psi}}_{C}^{LSJ{J}_{z}}$, which is the so-called “core” wave function, describing the system in the region where all the particles are close to each other and their mutual interaction is strong, and ${\text{\Psi}}_{A}^{LSJ{J}_{z}}$, which is the so-called “asymptotic” wave function, describing the relative motion between nucleon *N* and nucleus *Y* in the asymptotic region, where the *N* − *Y* interaction is negligible or reduces to the Coulomb interaction in the case of *N* ≡ *p*. The core function ${\text{\Psi}}_{C}^{LSJ{J}_{z}}$ has to vanish at large *N* − *Y* distances, and can be expanded in terms of the HH basis as for the bound state. Therefore, using Equation (35), we can write

The asymptotic wave function ${\text{\Psi}}_{A}^{LSJ{J}_{z}}$ is the solution of the Schrödinger equation of the relative *N* + *Y* motion. It is written as a linear combination of the following functions

Here we have indicated with ${C}$ a normalization factor [to be explained below, see Equation (49)]. The sum runs over the *N*_{p} even permutations of the *A* nucleons necessary to antisymmetrize the function ${\Omega}_{LSJ{J}_{z}}^{\lambda}$, χ_{1/2}(*N*) and ϕ_{SY}(*Y*) are the nucleon *N* and nucleus *Y* wave functions, respectively, and *y*_{p} is the relative distance between *N* and the c.m. of nucleus *Y* and is proportional to *x*_{N − j+1} of Equation (3). Furthermore, ${Y}_{L}({\widehat{y}}_{p})$ is the standard spherical harmonic function, and the functions ${R}_{L}^{\lambda}({y}_{p})$ for λ = *R, I* are respectively the regular and irregular solutions of the two-body *N* + *Y* Schrödinger equation without the nuclear interaction. They are explicitly written as [19, 31]

where *q* is the modulus of the *N* − *Y* relative momentum, such that the total kinetic energy in the c.m. frame is ${T}_{c.m.}={q}^{2}/2\mu $, μ being the *N* − *Y* reduced mass, $\eta ={Z}_{N}{Z}_{Y}\mu {e}^{2}/q$ is the Coulomb parameter, where *Z*_{N} and *Z*_{Y} are the charge numbers of *N* and *Y*, and *F*_{L}(η, *qy*_{p}) and *G*_{L}(η, *qy*_{p}) are the regular and irregular Coulomb functions defined in the standard way [52]. The factor *C*_{L}(η) is defined in reference [52] as

The factor $(2L+1)!!{q}^{L}{C}_{L}(\eta )$ has been introduced so that the functions ${R}_{L}^{R}({y}_{p})$ and ${R}_{L}^{I}({y}_{p})$ have a finite limit for *q* → 0. Finally, the function *f*(*b, y*_{p}) in Equation (42) is given by

so that the divergent behavior of *G*_{L}(η, *qy*_{p}) for small values of *y*_{p} is cured, and ${R}_{L}^{I}({y}_{p})$ is well-defined also in this limit. The trial parameter *b* is determined by requiring *f*(*b, y*_{p}) → 1 for large values of *y*_{p}, leaving therefore unchanged the asymptotic behavior of the scattering wave function. A value of *b* ~ 0.25 fm^{−1} has been found appropriate in all the considered cases. The non-Coulomb case of Equations (41) and (42) is obtained if either *Z*_{N} or *Z*_{Y} = 0, so that the functions *F*_{L}(η, *qy*_{p})/(*qy*_{p}) and *G*_{L}(η, *qy*_{p})/(*qy*_{p}) are replaced by the regular and irregular Riccati-Bessel functions *j*_{L}(*qy*_{p}) and *n*_{L}(*qy*_{p}) as defined in reference [52], and the factor (2*L*+1)*!!C*_{L}(η) reduces to 1 for η → 0 [52].

With these definitions, ${\text{\Psi}}_{A}^{LSJ{J}_{z}}$ can be cast in the form

where the parameters ${{R}}_{LS,{L}^{\prime}{S}^{\prime}}^{J}(q)$ give the relative weight between the regular and irregular components of the wave function. These parameters can be written in terms of the reactance matrix (${K}$-matrix) elements as [19, 31]

The ${K}$-matrix, by definition, is such that its eigenvalues are tanδ_{LSJ}, δ_{LSJ} being the phase shifts. The sum over *L*′ and *S*′ in Equation (45) is over all values compatible with a given *J* and parity π, and therefore the sum over *L*′ is limited to include either even or odd values since (−1)^{L}′ = π.

Using Equations (39) and (45), the full scattering wave functions is written as

where the unknown quantities are the coefficients *c*_{K; m} and ${{R}}_{LS,{L}^{\prime}{S}^{\prime}}^{J}(q)$. In order to determine their values, we use the Kohn variational principle [63], which states that the functional

has to be stationary with respect to variations of the trial parameters *c*_{K; m} and ${{R}}_{LS,{L}^{\prime}{S}^{\prime}}^{J}(q)$ in ${\text{\Psi}}_{NY}^{LSJ{J}_{z}}$. Here *E* is the total energy of the system, and the normalization coefficients ${C}$ of the asymptotic functions ${\Omega}_{LSJ{J}_{z}}^{\lambda}$ in Equation (40) are chosen so that

The variation of the diagonal functionals of Equation (48) with respect to the linear parameters *c*_{K; m} leads to a system of linear inhomogeneous equations,

where the two terms *D*^{λ} corresponding to λ ≡ *R, I* are defined as

Therefore, two sets of the coefficients ${c}_{K;m}^{\lambda}$ are obtained, depending on λ ≡ *R, I*, and consequently, we can introduce two core functions, defined as

The matrix elements ${{R}}_{LS,{L}^{\prime}{S}^{\prime}}^{J}(q)$ are obtained varying the diagonal functionals of Equation (48) with respect to them. This leads to the following set of algebraic equations

with the coefficients *X* and *Y* defined as

The solution of Equation (53) provides a first-order estimate of the matrix elements ${{R}}_{LS,{L}^{\prime}{S}^{\prime}}^{J}(q)$. A second-order estimate of ${{R}}_{LS,{L}^{\prime}{S}^{\prime}}^{J}(q)$, and consequently of ${{K}}_{LS,{L}^{\prime}{S}^{\prime}}^{J}(q)$, is given by the quantities $[{{R}}_{LS,{L}^{\prime}{S}^{\prime}}^{J}(q)]$, obtained by substituting in Equation (48) the first order results of Equations (50) and (53). Such second-order calculation provides then a symmetric ${K}$-matrix. This condition is not imposed *a priori*, and therefore it is a useful test of the numerical accuracy reached by the method.

The Kohn variational principle as explained so far is particularly useful in the case of *q* = 0 (zero-energy scattering). For *q* = 0 the scattering can occur only in the *L* = 0 channel, and the observables of interest are the scattering lengths. Within the present approach, they can be easily obtained from the relation

from which

An alternative version of the Kohn variational principle is the so-called complex Kohn variational principle for the ${S}$-matrix, quite convenient when *q* ≠ 0 and especially above the *Y* nucleus breakup threshold, as explained in reference [64]. In this case, the Kohn variational principle of Equation (48) becomes

where

${\text{\Psi}}_{C}^{LSJ{J}_{z}}$ being expanded as in Equation (39) and

The functions ${\Omega}_{LSJ{J}_{z}}^{\lambda}$ have been given in Equation (40). Note that, with the above definition, the reactance ${K}$-matrix elements can be related to the ${S}$-matrix elements as

The differentiation of the complex Kohn variational principle of Equation (57) leads to a set of equations for *c*_{K; m} and ${{S}}_{LS,{L}^{\prime}{S}^{\prime}}^{J}(q)$ similar to those given in Equations (50) and (53), where now λ stands for λ = +, −.

We conclude this section with the following remarks: (i) the calculation of the matrix elements between the core functions ${\text{\Psi}}_{C}^{LSJ{J}_{z}}$ can be performed with the HH expansion either in coordinate- or in momentum-space, depending on what is more convenient. Therefore, regarding this part, we can apply the method with any potential model, both local or non-local. (ii) Some difficulties arise with the calculations of the potential matrix elements which involve ${\Omega}_{LSJ{J}_{z}}^{\lambda}$, i.e., $\langle K,m|V|{\Omega}_{LSJ{J}_{z}}^{\lambda}\rangle $ present in Equation (51), and $\langle {\Omega}_{{L}^{\prime}{S}^{\prime}J{J}_{z}}^{{\lambda}^{\prime}}+{\text{\Psi}}_{C}^{{L}^{\prime}{S}^{\prime}J{J}_{Z},{\lambda}^{\prime}}|V|{\Omega}_{LSJ{J}_{z}}^{\lambda}\rangle $ of Equation (54), with λ, λ′ = *R*, *I*. In particular, we note that, being ${\Omega}_{LSJ{J}_{z}}^{\lambda}$ given in coordinate-space, which is particularly suitable when the Coulomb interaction is considered, as for *p* − *d* scattering, the non-local potential expressed in momentum-space is Fourier transformed to work in coordinate-space. The consequent integration on the momentum transfer are easily performed for the recent chiral and *V*_{low − k} potential models, but not for the non-local meson-theoretic CDBonn potential model, which has a high-momentum tail. Therefore, the CDBonn potential model has not been used in the study of the scattering processes presented here. We further refer to reference [31] for all technical details. (iii) The three-nucleon interaction models which at the moment have been implemented with the HH method are only the local ones, like the Urbana IX potential (UIX) of reference [65] and the N2LO model of reference [66]. The models used so far, besides being local, have a well-defined operatorial structure. In this case, the projection procedure as used for the two-nucleon interaction is not needed and the approach follows well-established footsteps, as explained in references [67, 68].

## 3. Selected Results

We present in this section selected results obtained with the HH method described above. The method has been applied widely since many years, and therefore a selection is mandatory. We have followed these criteria: (i) we focus on the results obtained after 2008, year of the publication of the review of reference [19] on the same method. (ii) We restrict ourselves to the potential models, mostly discussed in the Research Topic of which this contribution is part. They are the most widely used models. (iii) We concentrate on the results obtained for the *A* = 3, 4 elastic scattering observables, but we present briefly also the corresponding bound state results.

The aim of this section is two-fold: first of all we wish to show the effectiveness of the HH method for few-nucleon systems; secondly, we want to emphasize that the HH method, as well as any *ab initio* method, is an essential tool for testing and eventually improving nuclear interaction models.

All the results presented here are obtained at convergence, i.e., the HH expansion and the expansion on the Laguerre polynomials [see Equations (32) and (33)] has been pushed so that an accuracy of 1 keV (10 keV) is reached for the *A* = 3 (*A* = 4) binding energies, and the numerical accuracy on the scattering lengths is of the order of 0.001 fm. For a discussion on the convergence of the expansion see, for instance, references [30, 31].

The potentials which will appear in the following subsections include both two- and three-nucleon interactions. They are the phenomenological two-nucleon interaction Argonne *v*_{18} (AV18) [4], augmented by the three-nucleon Urbana IX (UIX) model [65], the meson-theoretic CDBonn potential [5] (CDB), together with the three-nucleon Tucson-Melbourne [69, 70] (TM) model, and the *V*_{low − k} potential [71], obtained from the AV18 with Λ = 2.2 fm^{−1}, so that the triton binding energy is reproduced. We consider in addition also chiral potentials, in particular the two-nucleon interaction models of the Idaho group of reference [72], obtained at next-to-next-to-next-to-leading order (N3LO), and here labeled with N3LO-I, and the more recent models derived by the same group in reference [73], here labeled according to the chiral order, i.e., from leading order (LO) up to next-to-next-to-next-to-next-to leading order (N4LO). All these two-nucleon models have been augmented with a (local) three-nucleon interaction derived up to N2LO as in reference [66]. The momentum-cutoff value is set equal to Λ = 500 MeV, unless differently specified. Note that the low-energy constants (LECs) *c*_{D} and *c*_{E} are those of reference [66] when the N2LO three-nucleon interaction is used in conjunction with the N3LO-I two-nucleon potential, while the LECs are those of reference [74] when the N2LO three-nucleon interaction is used in conjunction with the N2LO, N3LO, and N4LO two-nucleon interactions of reference [73] (no three-nucleon interaction is present at lower chiral order). To be remarked that the LECs *c*_{D} and *c*_{E}, and more generally the parameters entering the three-nucleon interaction model, depend on which two-nucleon interaction is used, as well as on which set of observables is used for their determination. This is why, for instance, the N2LO three-nucleon interaction in conjunction with the N2LO two-nucleon interaction has different values for the LECs compared to those present in the same three-nucleon interaction considered together with the N3LO or N4LO two-nucleon interaction. Finally, we will present results obtained also with the minimally non-local chiral potentials of the Norfolk group, as derived in reference [75] for the two-nucleon, and in references [1, 76] for the three-nucleon interaction. The two-nucleon models are labeled NVIa, NVIIa, NVIb, and NVIIb depending on the cutoff value and the maximum laboratory energy of the considered *NN* database. When the three-nucleon interaction are included, we will refer to NV2+3/Ia, NV2+3/IIa, and so on, corresponding to the fitting procedure of reference [1], and NV2+3/Ia*, NV2+3/IIa*, and so on, corresponding to the fitting procedure of reference [76]. We discuss in more details these fitting procedures below, and we refer the reader to the original references, or to the contributions present in this Research Topic. To be noticed that when the HH method is used to study the bound states, the local AV18, AV18/UIX, NV, and NV2+3 potentials have been all augmented by the full electromagnetic interaction, which includes corrections up to α^{2} (α is the fine-structure constant) [77]. On the other hand, the non-local CDB, CDB/TM, and all the non-local chiral potentials retain only the point-Coulomb interaction. The point-Coulomb interaction, and not the full electromagnetic one, is also used when studying the scattering states presented below.

### 3.1. *A* = 3, 4 Bound States

The results for the trinucleon and ^{4}He binding energies, obtained using all the above mentioned potentials, are given in Table 1. To be noticed that in many cases, the experimental trinucleon binding energy is used for the LECs fitting procedure performed applying the HH method. When this occurs, the corresponding HH results is underlined in the table. The results not underlined are obtained using the “original” two- and three-nucleon interactions, whose parameters are usually fitted to the triton binding energy and other observables, applying *ab initio* methods different than the HH method. The HH results are therefore not necessarily in perfect agreement with the experimental data.

**Table 1**. The binding energies in MeV for ^{3}H, ^{3}He, and ^{4}He, calculated with the HH technique using different Hamiltonian models.

We briefly outline the fitting procedure for the LECs *c*_{D} and *c*_{E} in order to better understand the results, and we refer to references [1, 76, 79] for more details. The ^{3}H and ^{3}He ground state wave functions are calculated using a given two- and three-nucleon potential, and the corresponding LECs *c*_{D} and *c*_{E} are determined by fitting the *A* = 3 experimental binding energies, corrected for a small contribution (+7 keV in ^{3}H and −7 keV in ^{3}He), due to the *n* − *p* mass difference [44], since in the present HH method this effect is neglected. This procedure generates two trajectories, one for ^{3}H and one for ^{3}He, in the {*c*_{D}, *c*_{E}} plane, so that each point of the trajectory corresponds to the correct binding energy. The two trajectories are typically extremely close to each other and the average can be safely considered, since the points of the average trajectory typically lead to *A* = 3 binding energies within 10 keV of the experimental values. A second observable is needed in the fitting procedure. In reference [1] the *n* − *d* doublet scattering length ^{2}*a*_{nd} has been used, which leads in the {*c*_{D}, *c*_{E}} plane to another trajectory, which is very close to the one corresponding to the ^{3}H binding energy, but not exactly overlapping. This is a well-known fact, that the ^{3}H binding energy and ^{2}*a*_{nd} are correlated observables. However, it is possible to find an intersection point of the two trajectories, which allows to determine the LECs. This procedure has been used for the NV2+3/Ia, NV2+3/Ib, NV2+3/IIa, and NV2+3/IIb potential models. The corresponding {*c*_{D}, *c*_{E}} values, as given in Table 1 of reference [1], are {3.666, −1.638}, {−2.061, −0.982}, {1.278, −1.029}, {−4.480, −0.412}, respectively. Alternatively we can choose as the second observable the Gamow-Teller matrix element of tritium β-decay, to take advantage of the fact that the LEC *c*_{D} enters also in the two-nucleon axial current operator at N2LO [76, 79–81]. This second procedure has been used for the N2LO/N2LO, N3LO/N2LO, and N4LO/N2LO potentials of reference [74], and the NV2+3/Ia*, NV2+3/Ib*, NV2+3/IIa*, and NV2+3/IIb* potential models of reference [76]. In this last case, we report the corresponding {*c*_{D}, *c*_{E}} values for completeness, which are {−0.635, −0.090}, {−4.710, 0.550}, {−0.610, −0.350}, {−5.250, 0.050}, respectively.

We can now proceed with some comments regarding the binding energies results of Table 1. (i) The large variety of models for the nuclear interaction which the HH method can handle is an indication of how strong and reliable this method has become. Furthermore, we should mention that the theoretical uncertainty is of 1 keV (10 keV) for the *A* = 3 (^{4}He) binding energies. The HH method is therefore extremely accurate. Furthermore, all the HH results are in very good agreement with the values reported in the literature, when available. (ii) In order to reproduce the experimental binding energies the inclusion of three-nucleon force is essential. In all cases, the triton binding energy is well-reproduced, within few keV. On the other hand, the ^{4}He binding energies can differ from the experimental value of even up to 700 keV (in the CDB/TM case). (iii) In the case of the NV2+3 potential models, when the observables used to fit the LECs are the triton binding energy and ^{2}*a*_{nd}, we notice a systematic overestimation of the ^{3}He binding energy. (iv) All the results for the *A* = 3 (*A* = 4) binding energies obtained with any model for the two- and three-nucleon interaction are within 10 (400) keV from the experimental values. Therefore we can conclude that any of the constructed model is essentially able to reproduce these very light nuclei.

### 3.2. *N* − *d* Scattering

One of the remarkable features of the HH method resides in its capability of dealing with local as well as with non-local potentials, formulated in either coordinate or momentum space, not only for the bound states, as we have seen above, but also for *N* − *d* scattering observables. This has been demonstrated in reference [31], in which the local AV18 and the non-local chiral N3LO-I potential models were used to calculate the *N* − *d* elastic scattering observables below the deuteron breakup threshold. Here we present results with a subset of all the potential models mentioned above, and in particular with the AV18, AV18/UIX, the N3LO-I, N3LO-I/N2LO, and some of the NV and NV2+3 models. A further class of nuclear interactions that has been tested using the HH method is represented by the so-called *V*_{low − k} potential obtained from the AV18 with Λ = 2.2 fm^{−1}, so chosen to reproduce the triton binding energy when the complete electromagnetic interaction is used [71]. We do not report here detailed investigations on the convergence of the HH expansion, but we can mentioned that this convergence is faster for the non-local potentials as compared to the local ones, due to the much softer behavior at small distances. For instance, for *N* − *d* elastic scattering in the channel *J*^{π} = 1/2^{+}, the HH basis can be of the order of 12000 (7000) elements with the NV (N3LO-I) potential to get convergence.

We first consider the converged results for the *n* − *d* and *p* − *d* doublet and quartet scattering lengths, which are given in Table 2, together with the very precise experimental result from reference [82] for ^{2}*a*_{nd}, and the older experimental results from reference [83] for ^{4}*a*_{nd}. Though no experimental data are available for ^{2}*a*_{pd} and ^{4}*a*_{pd}, the results of the energy-dependent phase-shift analysis of reference [84], using very low *p* − *d* data, is reported. All the results are obtained using the pure Coulomb electromagnetic interaction. When the full electromagnetic interaction is used, ^{4}*a*_{nd} remains practically unchanged, while ^{2}*a*_{nd} becomes smaller. For the NVIa and NVIb potentials, for instance, ^{2}*a*_{nd} = 1.103 fm and 1.293 fm, respectively, with the full electromagnetic interaction. As it is clear from inspection of Table 2, while ^{4}*a*_{nd} is very little model-dependent and in good agreement with experiment, the same is not true for ^{2}*a*_{nd}. In particular, the inclusion of a three-nucleon force appears necessary to bring the results closer to the experimental datum. However, not every model agrees with the experiment. The disagreement is more pronounced for the *V*_{low − k} interaction, showing that this observable cannot be simply reproduced by increasing the attraction of the two-nucleon interaction, as is done in this case by choosing the right value for Λ to describe the triton; instead, a subtle balance between attraction and repulsion in the three-nucleon system has to be reached. Indeed, being the zero-energy *n* − *d* scattering state orthogonal to the triton, the associated wave function presents a node in the relative distance, whose precise position, which is related to the scattering length, depends on the interplay between attraction and repulsion. The results of the *p* − *d* phase-shift analysis given in reference [84] are a first tentative to determine the *p* − *d* scattering lengths from *p* − *d* data. Very few experimental data exist for center-of-mass energies below 500 keV, introducing large uncertainties in the quartet case. In the case of the doublet scattering length, difficulties arise from the particular pole structure of the doublet *p* − *d* effective range expansion close to threshold (see, for example, Figure 1 of reference [84]).

**Table 2**. *n* − *d* and *p* − *d* doublet and quartet scattering lengths in fm calculated with the HH technique using different Hamiltonian models.

With the purpose of investigating the capability of some widely used models of three-nucleon interaction to reproduce ^{2}*a*_{nd}, a sensitivity study was conducted in reference [85] taking the AV18 as the reference two-nucleon interaction. Three different models of the three-nucleon interactions were considered: the UIX, the TM and the chiral N2LO of reference [66]. Their parameters were adjusted, constraining them to reproduce simultaneously ^{2}*a*_{nd} and the triton binding energy, and the resulting value for the ^{4}He binding energy was calculated. For the UIX model, a reasonable description of these three observables was possible, at the cost of a sizable increase of the repulsive term, as compared to the original parameterization. A similar conclusion held for the TM model, where a repulsive short-range term was found to be necessary. Finally, for the N2LO three-nucleon interaction, the relative importance of the parameters involving the *P*-wave pion rescattering had to be changed. This is not surprising, due to the mismatch between the physics underlying the adopted models for two- and three-nucleon interactions. Also in this case, a repulsive short-range interaction was preferred. Then, a set of polarization observables on elastic *p* − *d* scattering were computed using the AV18 augmented by the modified versions of the three-nucleon interactions models as described above. These led to three classes of interaction models. As an interesting result, all models within a given class led to very similar predictions, but for some observables, namely the proton *A*_{y} and the deuteron *iT*_{11}. These predictions were different from class to class, and all in disagreement with the data. This is shown in Figure 1. Since the three classes of models mostly differ in their short-distance behavior, it follows that an improvement in this component of the three-nucleon interaction is needed to explain the data. Indeed, no sensible improvement was obtained as compared to the original AV18/UIX model.

**Figure 1**. The vector analyzing powers *A*_{y} and *iT*_{11} for *p* − *d* elastic scattering at center-of-mass energy *E*_{c.m.} = 2 MeV, using models in the AV18/TM class (cyan bands), AV18/UIX (violet bands), and AV18/N2LO (red bands). The predictions of the original AV18/UIX model (solid lines) and the experimental points from reference [86] are also shown.

In order to be more quantitative, as to the accuracy of the existing models of two- and three-nucleon interaction, we show in Table 2 the χ^{2}/datum for all *p* − *d* elastic scattering observables at different center-of-mass energies, as obtained with the AV18 and N3LO-I two-body interactions, without or with the inclusion of the UIX and N2LO three-nucleon interaction models [31]. It is clear that all considered models fail to give a satisfactory description of all polarization observables, especially for *A*_{y} and *iT*_{11}. From the previous discussion, there are strong hints that the improvement should come from a more accurate modeling of the short distance structure of the three-nucleon interaction. Therefore, in reference [87] all the possible short-distance (contact) structures for the three-nucleon interaction have been classified up to the subleading order of a systematic low-energy expansion. It has been found that the short-distance component of the three-nucleon interaction can be parameterized by ten LECs, denoted by *E*_{i} with *i* = 1, …, 10. The corresponding three-nucleon potential in configuration space can be written as

where **σ**_{i} (**τ**_{i}) are the Pauli spin (isospin) matrices of particle *i*, **r**_{ij} is the relative distance between particles *i* and *j*, and *S*_{ij} and (**L**·** S**)

_{ij}are, respectively, the tensor and spin-orbit operators. The profile functions

*Z*

_{0}(

*r*; Λ) are written as

with *F*(**k**^{2}; Λ) a suitable cutoff function which suppresses the momentum transfers **k** above a given short-distance cutoff Λ. In Equation (61), the basis of operators has been chosen so that most terms in the potential can be viewed as an ordinary interaction of particles *ij* with a further dependence on the coordinate of the third particle *k*. In reference [88], elastic *p* − *d* scattering data at *E*_{c.m.} = 2 MeV center-of-mass energy have been used to fit the *E*_{i} LECs, when the subleading three-nucleon interaction given in Equation (61) is considered in addition to the AV18/UIX interaction. Also ^{2}*a*_{nd} and ^{4}*a*_{nd} and the triton binding energy have been included in the fit.

The results of reference [88] can be summarized as follows. First of all, we noticed that the operators which play a leading role in reducing the large χ^{2}/datum of Table 3 are the spin-orbit and tensor interactions, which depend on the LECs *E*_{5} and *E*_{7}. We present in Table 4 the results of a fit where only the terms proportional to *E*_{0}, *E*_{5}, and *E*_{7} are kept. The LEC *E*_{0} is used to fix the triton binding energy. Then the experimental data for the doublet and quartet *n* − *d* scattering lengths of references [82, 83], and those of several *p* − *d* scattering observables at 2 MeV center-of-mass energy of reference [86] are used for the determinations of the LECs. As it is shown in Table 4, the χ^{2}/datum is drastically reduced to ~2 for the short distance cutoff Λ of Equation (62) between 200 and 500 MeV. More sophisticated fits, including all the involved LECs, lead to only slightly better χ^{2}/datum ~1.6. In Figure 2 we show the corresponding fitted curves compared to the AV18 and AV18/UIX predictions. It is clear that a very accurate description can be obtained with only the spin-orbit and tensor subleading operators. We also note that the values of the LECs *e*_{0}, *e*_{5}, *e*_{7}, defined in terms of *E*_{0}, *E*_{5}, *E*_{7} as ${E}_{0}={e}_{0}/({F}_{\pi}^{4}\Lambda )$, ${E}_{i}={e}_{i}({F}_{\pi}^{4}{\Lambda}^{3})$, *i* = 5, 7, *F*_{π} = 92.4 MeV being the pion decay constant, are of order 1 as expected.

**Table 3**. χ^{2}/datum of the *p* − *d* elastic scattering observables at center-of-mass energies *E*_{c.m.} = 0.666, 1.33, 1.66 and 2.0 MeV, calculated with the N3LO-I or AV18 two-nucleon only, and the N3LO-I/N2LO or AV18/UIX two- and three-nucleon Hamiltonian models.

**Table 4**. χ^{2}/datum of the two-parameter fit to *p* − *d* elastic scattering data at *E*_{c.m.} = 2 MeV, obtained neglecting in Equation (61) all the subleading operators except the leading contact term proportional to the LEC *E*_{0}, and the tensor *S*_{ij} and spin-orbit (**L**·** S**)

_{ij}operators, proportional to the LECs

*E*

_{5}and

*E*

_{7}, considered on top of the AV18/UIX potential model.

**Figure 2**. Curves obtained including only the tensor and spin-orbit subleading contact operator on the top of the AV18/UIX interaction, fitted to a set of cross section and polarization observables in *p* − *d* elastic scattering at 2 MeV center-of-mass energy [86], for Λ = 200 − 500 MeV (red bands), are compared to the purely two-body AV18 interaction (dashed black lines) and to the AV18/UIX two- and three-nucleon interaction (dashed-dotted blue lines).

With the interaction fitted using the *E*_{c.m.} = 2 MeV data of reference [86], we can perform a study at lower energies, where experimental data exist. As a representative example we show in Figure 3 the results corresponding to *E*_{c.m.} = 0.666 MeV, from which we can observe that the adopted interaction captures quite nicely the energy dependence of the data. In reference [88], a fit including all the subleading operators of Equation (61) leads to predictions in even better agreement with the data. However, in order to obtain further improvements, a global fit at multiple energies should be performed.

**Figure 3**. Predictions obtained with the three-nucleon interaction models discussed in the text with Λ = 200 − 500 MeV (red bands) for a set of cross section and polarization *p* − *d* observables at 0.666 MeV center-of-mass energy, as compared to the purely two-body AV18 interaction (dashed black lines), to the AV18/UIX two- and three-nucleon interaction (dashed-dotted blue lines), and to the experimental data of reference [90].

### 3.3. *p*−^{3}He and *n* − ^{3}H Scattering

The study of *N* − *d* scattering to constrain the three-nucleon force has the limitation of being mostly restricted to the isospin *T* = 1/2 channel. From this perspective, *A* = 4 systems open new possibilities, besides being of direct relevance for the role they play in many reactions of astrophysical and cosmological interest. The HH method has been used in this context to address first of all the *n* − ^{3}H [32] and *p*−^{3}He [33, 35] elastic scattering at low energies. The HH method applied to these systems has been benchmarked in reference [34] with the only two other *ab initio* methods which can study low-energy scattering states, with full inclusion of the Coulomb interaction. They are the AGS equations solved in momentum space (see for a review references [16, 17] and references therein), and the FE method in configuration space (see reference [13]. This topic is also covered in the present Research Topic). All these methods differ by <1%, which is smaller than the experimental uncertainties of the available data. The agreement found using softer potentials of the *V*_{low − k}-type is even better.

The *n* − ^{3}H elastic scattering total cross section is shown in Figure 4. From inspection of the figure, we can see a sizable dependence on the three-nucleon interaction, both in the very low-energy region and in the peak region (for neutron laboratory energy *E*_{n} ~ 3.5 MeV). Indeed, at very low energies, it is crucial to have a correct description of the triton binding energy in order to reproduce the data, whereas in the peak region there is more model dependence. The HH calculations of Figure 4 have been performed using the non-local chiral N3LO-I two-nucleon potential, also supplemented by the chiral N2LO three-nucleon interaction of reference [66] with the LECs fixed to reproduce the *A* = 3, 4 binding energies. This leads to a remarkable agreement with the available experimental data in the low-energy region. The chiral N3LO-I model seems to perform better than the AV18 one also in the peak region.

**Figure 4**. *n* − ^{3}H total cross sections calculated with the AV18 (dashed black line), AV18/UIX (solid black line), N3LO-I (dashed blue line), and the N3LO-I/N2LO (solid red line) potential models as function of the incident neutron laboratory energy *E*_{n}. The experimental data are from reference [91].

In Figure 5 we show the *n* − ^{3}H differential cross section compared to the experimental data at three different neutron laboratory energies. As it is clear from inspection of the figure, the N3LO-I/N2LO results are in nice agreement with the data. A further study of convergence with respect to chiral orders and of cutoff dependence would be highly desirable, and it is currently underway.

**Figure 5**. *n*−^{3}H differential cross sections calculated with the N3LO-I (dashed blue lines) and the N3LO-I/N2LO (solid red lines) interaction models for three different incident neutron energies. The experimental data are from reference [92].

Much more accurate data are available for *p*−^{3}He elastic scattering, whose polarization observables have also been accurately measured [93]. Similarly to the *p* − *d* case, there is a strong discrepancy between theory and experiment for the proton analyzing power *A*_{y}. In reference [35] the HH method has been applied with the N3LO-I/N2LO chiral potential model, in this case obtained with two different values of the momentum cutoff Λ = 500, 600 MeV [94], and two different procedures to fix the LECs entering the three-nucleon interaction, i.e., either reproducing the *A* = 3, 4 binding energies [66], or reproducing the triton binding energy and Gamow-Teller matrix element in tritium β-decay [79]. We show in Figure 6 the corresponding results for proton laboratory energy of 5.54 MeV, compared to experimental data. The two bands reflect the cutoff dependence and the model dependence introduced by the LECs determinations. As it is clear, the *A*_{y} discrepancy is largely reduced down to the 8–10% level. Note that these asymmetries are 10 times larger in the *A* = 4 systems than for *p* − *d* and *n* − *d*. The remaining discrepancy, although it appears small, is of the order of 0.05, the size of *A*_{y} for *p* − *d*. Therefore, we expect that the subleading components of the three-nucleon interactions discussed in section 3.2 could give a correction of the necessary order of magnitude to solve the remaining discrepancy. Work is in progress in this direction.

**Figure 6**. *p*−^{3}He differential cross section, analyzing powers and various spin correlation coefficients at proton laboratory energy *E*_{p} = 5.54 MeV, calculated with only the two-nucleon N3LO-I (light cyan band) or with two- and three-nucleon interaction N3LO-I/N2LO (darker blue band). The experimental data are from references [95–97]. See text for more details.

### 3.4. *p*−^{3}H and *n* − ^{3}He Scattering

The treatment of *p*−^{3}H and *n* − ^{3}He scattering, even below the *d* + *d* threshold, is more challenging due to the coupling between these two channels and to the presence of both isospin 0 and 1 states. Also in this case, recently, in reference [48], a benchmark calculation has been performed with the HH, AGS and FE methods, using the N3LO-I interaction. Good agreement among the three methods has been found, with discrepancies smaller than the uncertainties in the experimental data. In references [98, 99], we have studied with the HH method the effect of the inclusion of the N2LO three-nucleon interaction, with the LECs fixed from the triton binding energy and the Gamow-Teller matrix element in the tritium β-decay [79]. We show in Figure 7 the *p*−^{3}H differential cross section, for which, only at very low energies, below the opening of the *n* − ^{3}He channel, some sizable effects are visible. Otherwise, the three-nucleon interaction contributions are found very small. The *p*−^{3}H analyzing power at three values of the laboratory beam energy are shown in Figure 8. Also for this observable, the three-nucleon interaction effect is found too small to improve the agreement with the available experimental data.

**Figure 7**. *p*−^{3}H differential cross section at several values of the proton laboratory beam energy *E*_{p}, calculated with the N3LO-I (dashed blue lines) and with the N3LO-I/N2LO (solid red lines) interactions. The experimental data are from references [100–105].

**Figure 8**. *p*−^{3}H proton analyzing power at three values of the proton laboratory beam energy *E*_{p} calculated with the N3LO-I (dashed blue lines) and with the N3LO-I/N2LO (solid red lines) interactions. The experimental data are from reference [105].

We conclude showing in Figure 9 the HH results for the differential cross section and proton analyzing power of the charge-exchange reaction *p*+^{3}H → *n*+^{3}He at three different proton laboratory energies, compared with the experimental data. By inspection of the figure, we can see that also in this case the effects of the three-nucleon interaction are quite small, and sometimes go in the wrong direction as compared to the experimental data, as for the analyzing power *A*_{y0}. It is important to notice that this observable is mostly sensitive to the two-nucleon interaction. Therefore, it could be used for a more stringent test of the two-nucleon force.

**Figure 9**. *p*+^{3}H → *n*+^{3}He differential cross section and proton analyzing power at three values of the proton laboratory beam energy *E*_{p} calculated with the N3LO-I (dashed blue lines) and with the N3LO-I/N2LO (solid red lines) interactions. The experimental data are from references [106–110].

## 4. Conclusions and Outlook

In this work we have presented a review of the HH method, focusing on the most significant achievements after the year 2008, when the previous review on the HH method [19] was published. We have also included a presentation of the HH formalism with some detail, in order to make the reader appreciate the main concepts of the method and to provide him/her the instruments needed to implement the method by him/herself. We have then focused on the latest results obtained within the HH method. We can summarize the situation as follows: the HH method can solve the three- and four-body bound-state problem with great accuracy and with essentially any (local and non-local) model for the two-nucleon interaction available in the literature. The three-nucleon interaction models used so far are however only local. The *A* = 3, 4 scattering states have been studied with several local and non-local potentials below the target nucleus breakup threshold. Using local potentials, also the elastic channel above the breakup threshold have been investigated. The HH method has then a wide range of applications: it has been used not only to test the models for the two- and three-nucleon interactions, but also to determine the parameters entering in the subleading three-nucleon contact interaction, derived in reference [87]. This has allowed one to construct a model for the three-nucleon interaction able to solve, at least within the (preliminary) hybrid framework of reference [88], some long-standing puzzles, as the *A*_{y}-puzzle. Furthermore, the HH method has been widely used in the study of nuclear reactions of astrophysical interest, as well as the electroweak structure of light nuclei [50, 111, 112].

The HH method has still a lot of potentialities, which will be explored in the near future. First of all, we will implement the method to the case of a non-local three-nucleon interaction. This is widely requested, in order to have consistency in the two- and three-nucleon cutoff functions which appear in the models derived within chiral effective field theory, for instance in references [72, 73]. Once the LECs *c*_{D} and *c*_{E} will be determined using the non-local three-nucleon interaction with the same procedure outlined in section 3, they will be used in *fully consistent* studies of other systems, as nuclear and neutron matter.

Secondly, we can mention only preliminary applications of the HH method to describe breakup reactions in *A* = 3 [40]. Work on the implementation of the HH method to the breakup channels in *A* = 3, 4 is currently underway. It does not require significant modifications of the method, but still it has not been performed yet. Once done, the three- and four-body nuclear systems will be completely covered by the method.

As mentioned above, the extension of the method to the *A* = 5, 6 nuclear systems has been investigated and the first results obtained using a *V*_{low − k} interaction will appear soon and are indeed very promising. This is a major step for the HH method, as it will allow us to tackle a large number of interesting subjects, and especially a large number of nuclear reactions of astrophysical interest. From a first investigation, the further extension of the method to even larger values of *A*, i.e., *A* = 7, 8, seems feasible.

Finally, in order to have access to higher mass nuclear systems, both bound and scattering states, we could take advantage of the strong clusterization present in some of them, as, for instance, in ^{9}Be, which can be studied as a α−α−*n* system. In order to do so, the HH method must then be extended to the case of non-equal mass systems. And this, in turn, will allow to study also more exotic systems, as hypernuclei, where one nucleon is replaced with an hyperon. Works along this line have started in reference [61], and are conducted also by other groups [113, 114].

In conclusion, the HH method has quite a “glorious” history, and has fulfilled its service in the continuous test of the nuclear interaction models. However, this service is not yet at an end, and we expect to see the HH method playing a protagonist role also in the next years.

## Author Contributions

All the authors have contributed essentially equally to this work, read, and corrected the whole manuscript. In particular, JD-E, AG, AK, and LM have worked in writing the section on the formalism. LG, AK, LM, and MV have taken care of collecting the presented results. AK and LM have taken care of writing the introductory and concluding sections.

## Funding

This work of JD-E was supported by the Fonds de la Recherche Scientifique (FNRS) under Grant Number 4.45.10.08.

## Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

## Acknowledgments

The computational resources of the Istituto Nazionale di Fisica Nucleare (INFN), Sezione di Pisa, were gratefully acknowledged.

## Footnotes

1. ^The expression “*ab initio* method” has been quite widely used in the literature, with a somewhat less strict definition, than the one used here. Here we follow the definition of reference [7].

## References

1. Piarulli M, Baroni A, Girlanda L, Kievsky A, Lovato A, Lusk E, et al. Light-nuclei spectra from chiral dynamics. *Phys Rev Lett*. (2018) **120**:052503. doi: 10.1103/PhysRevLett.120.052503

2. Epelbaum E, Golak J, Hebeler K, Hüther T, Kamada H, Krebs H, et al. Few- and many-nucleon systems with semilocal coordinate-space regularized chiral two- and three-body forces. *Phys Rev C*. (2019) **99**:024313. doi: 10.1103/PhysRevC.99.024313

3. Rozpedzik D, Golak J, Skibinski R, Witala H, Glöckle W, Epelbaum E, et al. A first estimation of chiral four-nucleon force effects in ^{4}He. *Acta Phys Polon B*. (2006) **37**:2889.

4. Wiringa RB, Stoks VGJ, Schiavilla R. Accurate nucleon-nucleon potential with charge-independence breaking. *Phys Rev C*. (1995) **51**:38. doi: 10.1103/PhysRevC.51.38

5. Machleidt R. High-precision, charge-dependent Bonn nucleon-nucleon potential. *Phys Rev C*. (2001) **63**:024001. doi: 10.1103/PhysRevC.63.024001

6. Weinberg S. Nuclear forces from chiral lagrangians. *Phys Lett B*. (1990) **251**:288. doi: 10.1016/0370-2693(90)90938-3

7. Leidemann W, Orlandini G. Modern *ab initio* approaches and applications in few-nucleon physics with *A*≥4. *Prog Part Nucl Phys*. (2013) **68**:158. doi: 10.1016/j.ppnp.2012.09.001

8. Carlson J, Gandolfi S, Pederiva F, Pieper SC, Schiavilla R, Schmidt KE, et al. Quantum Monte Carlo methods for nuclear physics. *Rev Mod Phys*. (2015) **87**:1067. doi: 10.1103/RevModPhys.87.1067

9. Barrett BR, Navrátil P, Vary JP. *Ab initio* no core shell model. *Prog Part Nucl Phys*. (2013) **69**:131. doi: 10.1016/j.ppnp.2012.10.003

10. Fukui T, De Angelis L, Ma YZ, Coraggio L, Gargano A, Itaco N, et al. Realistic shell-model calculations for *p*-shell nuclei including contributions of a chiral three-body force. *Phys Rev C*. (2018) **98**:044305. doi: 10.1103/PhysRevC.98.044305

11. Coraggio L, Covello A, Gargano A, Itaco N, Kuo TTS. Effective shell-model hamiltonians from realistic nucleon-nucleon potentials within a perturbative approach. *Ann Phys*. (2012) **327**:2125. doi: 10.1016/j.aop.2012.04.013

12. Lazauskas R, Carbonell J. Testing nonlocal nucleon-nucleon interactions in four-nucleon systems. *Phys Rev C*. (2004) **70**:044002. doi: 10.1103/PhysRevC.70.044002

13. Lazauskas R, Carbonell J. The Faddeev-Yakubovsky Symphony. *Few Body Syst*. (2019) **60**:62. doi: 10.1007/s00601-019-1529-5

14. Glöckle W, Witala H, Hüber D, Kamada H, Golak J. The three-nucleon continuum: achievements, challenges and applications. *Phys Rep*. (1996) **274**:107. doi: 10.1016/0370-1573(95)00085-2

15. Nogga A, Hüber D, Kamada H, Glöckle W. Triton binding energies for modern NN forces and the π−π exchange three-nucleon force. *Phys Lett B*. (1997) **409**:19. doi: 10.1016/S0370-2693(97)00841-1

16. Fonseca AC, Deltuva A. Numerical exact *ab initio* four-nucleon scattering calculations: from dream to reality. *Few Body Syst*. (2017) **58**:46. doi: 10.1007/s00601-017-1211-8

17. Deltuva A, Fonseca AC, Sauer PU. Coulomb force effects in few-nucleon systems. *Few Body Syst*. (2019) **60**:29. doi: 10.1007/s00601-019-1496-x

18. Lazauskas R. Solution of the *n* − ^{4} He elastic scattering problem using the Faddeev-Yakubovsky equations. *Phys Rev C*. (2018) **97**:044002. doi: 10.1103/PhysRevC.97.044002

19. Kievsky A, Rosati S, Viviani M, Marcucci LE, Girlanda L. A High-precision variational approach to three- and four-nucleon bound and zero-energy scattering states. *J Phys G Nucl Part Phys*. (2008) **35**:063101. doi: 10.1088/0954-3899/35/6/063101

20. Okubo S. Diagonalization of Hamiltonian and Tamm-Dancoff equation. *Progr Theor Phys*. (1954) **12**:603. doi: 10.1143/PTP.12.603

21. Suzuki K. Construction of Hermitian effective interaction in nuclei: general relation between Hermitian and non-Hermitian forms. *Progr Theor Phys*. (1982) **68**:246. doi: 10.1143/PTP.68.246

22. Suzuki K, Okamoto R. Degenerate perturbation theory in quantum mechanics. *Progr Theor Phys*. (1983) **70**:439. doi: 10.1143/PTP.70.439

23. Barnea N, Leidemann W, Orlandini G. State dependent effective interaction for the hyperspherical formalism. *Phys Rev C*. (2000) **61**:054001. doi: 10.1103/PhysRevC.61.054001

24. Barnea N, Leidemann W, Orlandini G. State-dependent effective interaction for the hyperspherical formalism with noncentral forces. *Nucl Phys A*. (2001) **693**:565. doi: 10.1016/S0375-9474(01)00794-1

25. Barnea N, Leidemann W, Orlandini G. Hyperspherical effective interaction for nonlocal potentials. *Phys Rev C*. (2010) **81**:064001. doi: 10.1103/PhysRevC.81.064001

26. Bacca S, Barnea N, Leidemann W, Orlandini G. Isoscalar monopole resonance of the alpha particle: a prism to nuclear Hamiltonians. *Phys Rev Lett*. (2013) **110**:042503. doi: 10.1103/PhysRevLett.110.042503

27. Bacca S, Barnea N, Leidemann W, Orlandini G. Examination of the first excited state of ^{4} He as a potential breathing mode. *Phys Rev C*. (2015) **91**:024303. doi: 10.1103/PhysRevC.91.024303

28. Bacca S, Pastore S. Electromagnetic reactions on light nuclei. *J Phys G Nucl Part Phys*. (2014) **41**:123002. doi: 10.1088/0954-3899/41/12/123002

29. Kievsky A, Marcucci LE, Rosati S, Viviani M. High-precision calculation of the triton ground state within the hyperspherical-harmonics method. *Few Body Syst*. (1997) **22**:1. doi: 10.1007/s006010050049

30. Viviani M, Marcucci LE, Rosati S, Kievsky A, Girlanda L. Variational calculation on *A* = 3 and 4 nuclei with non-local potentials. *Few Body Syst*. (2006) **39**:159. doi: 10.1007/s00601-006-0158-y

31. Marcucci LE, Kievsky A, Girlanda L, Rosati S, Viviani M. *N*-*d* elastic scattering using the hyperspherical harmonics approach with realistic local and nonlocal interactions. *Phys Rev C*. (2009) **80**:034003. doi: 10.1103/PhysRevC.80.034003

32. Viviani M, Kievsky A, Girlanda L, Marcucci LE. Neutron-triton elastic scattering. *Few Body Syst*. (2009) **45**:119. doi: 10.1007/s00601-009-0036-5

33. Viviani M, Girlanda L, Kievsky A, Marcucci LE, Rosati S. Proton-^{3}He elastic scattering at low energies and the ‘*A*_{y} Puzzle'. *EPJ Web Conf*. (2010) **3**:05011. doi: 10.1051/epjconf/20100305011

34. Viviani M, Deltuva A, Lazauskas R, Carbonell J, Fonseca AC, Kievsky A, et al. Benchmark calculation of n-^{3}H and p-^{3}He scattering. *Phys Rev C*. (2011) **84**:054010. doi: 10.1103/PhysRevC.84.054010

35. Viviani M, Girlanda L, Kievsky A, Marcucci LE. Effect of three-nucleon interaction in *p*−^{3}He elastic scattering. *Phys Rev Lett*. (2013) **111**:172302. doi: 10.1103/PhysRevLett.111.172302

36. Kievsky A, Viviani M, Rosati S. Polarization observables in *p* − *d* scattering below 30 MeV. *Phys Rev C*. (2001) **64**:024002. doi: 10.1103/PhysRevC.64.024002

37. Kievsky A, Viviani M, Marcucci LE. *N*-*d* scattering including electromagnetic forces. *Phys Rev C*. (2004) **69**:014002. doi: 10.1103/PhysRevC.69.014002

38. Kievsky A, Rosati S, Viviani M. Proton-deuteron elastic scattering above the deuteron breakup. *Phys Rev Lett*. (1999) **82**:3759. doi: 10.1103/PhysRevLett.82.3759

39. Viviani M, Kievsky A, Rosati S. The Kohn variational principle for elastic proton deuteron scattering above deuteron breakup threshold. *Few Body Syst*. (2001) **30**:39. doi: 10.1007/s006010170017

40. Garrido E, Kievsky A, Viviani M. Breakup of three particles within the adiabatic expansion method. *Phys Rev C*. (2014) **90**:014607. doi: 10.1103/PhysRevC.90.014607

41. Gattobigio M, Kievsky A, Viviani M. Nonsymmetrized hyperspherical harmonic basis for an *A*-body system. *Phys Rev C*. (2011) **83**:024001. doi: 10.1103/PhysRevC.83.024001

42. Dohet-Eraly J, Viviani M. Computing an orthonormal basis of symmetric or antisymmetric hyperspherical harmonics. *Comp Phys Comm*. (2020) 107183. doi: 10.1016/j.cpc.2020.107183

43. Dohet-Eraly J, Viviani M, Kievsky A, Marcucci LE, Gnech A. Five-nucleon systems with an hyperspherical harmonic method. In: Orr NA, Ploszajczak M, Marqués FM, Carbonell J, editors. *Recent Progress in Few-Body Physics*. Cham: Springer International Publishing (2020). p. 163–7. doi: 10.1007/978-3-030-32357-8_29

44. Nogga A, Kievsky A, Kamada H, Glöckle W, Marcucci LE, Rosati S, et al. Three-nucleon bound states using realistic potential models. *Phys Rev C*. (2003) **67**:034004. doi: 10.1103/PhysRevC.67.034004

45. Kamada H, Nogga A, Glöckle W, Hiyama E, Kamimura M, Varga K, et al. Benchmark test calculation of a four nucleon bound state. *Phys Rev C*. (2001) **64**:044001. doi: 10.1103/PhysRevC.64.044001

46. Kievsky A, Viviani M, Rosati S, Hüber D, Glöckle W, Kamada H, et al. Benchmark calculations for polarization observables in three nucleon scattering. *Phys Rev C*. (1998) **58**:3085. doi: 10.1103/PhysRevC.58.3085

47. Deltuva A, Fonseca AC, Kievsky A, Rosati S, Sauer PU, Viviani M. Benchmark calculation for proton-deuteron elastic scattering observables including Coulomb. *Phys Rev C*. (2005) **71**:064003. doi: 10.1103/PhysRevC.71.064003

48. Viviani M, Deltuva A, Lazauskas R, Fonseca AC, Kievsky A, Marcucci LE. Benchmark calculation of *p*−^{3}H and *n* − ^{3}He scattering. *Phys Rev C*. (2017) **95**:034003. doi: 10.1103/PhysRevC.95.034003

49. Marcucci LE, Nollett KM, Schiavilla R, Wiringa RB. Modern theories of low-energy astrophysical reactions. *Nucl Phys A*. (2006) **777**:111. doi: 10.1016/j.nuclphysa.2004.09.008

50. Marcucci LE, Mangano G, Kievsky A, Viviani M. Implication of the proton-deuteron radiative capture for Big Bang nucleosynthesis. *Phys Rev Lett*. (2016) **116**:102501. [Erratum: (2016) Phys. Rev. Lett. 117, 049901]. doi: 10.1103/PhysRevLett.116.102501

51. Fabre de la Ripelle M. The potential harmonic expansion method. *Ann Phys*. (1983) **147**:281. doi: 10.1016/0003-4916(83)90212-9

52. Abramowitz M, Stegun IA. *Handbook of Mathematical Functions*. New York, NY: Dover Publications (1964).

53. Raynal J, Revai J. Transformation coefficients in the hyperspherical approach to the three-body problem. *Nuovo Cim A*. (1970) **68**:612. doi: 10.1007/BF02756127

54. Viviani M. Transformation coefficients of hyperspherical harmonic functions of an *A*-body system. *Few Body Syst*. (1998) **25**:177. doi: 10.1007/s006010050101

55. Edmonds AR. *Angular Momentum in Quantum Mechanics*. Princeton, NJ: Princeton University Press (1954).

56. Rosati S. The hyperspherical harmonic method: a review and some recent developments. *Ser Adv Quant Many Body Theor*. (2002) **7**:339. doi: 10.1142/9789812777072_0009

57. Viviani M, Kievsky A, Rosati S. Calculation of the α-particle ground state within the hyperspherical harmonic basis. *Phys Rev C*. (2005) **71**:024006. doi: 10.1103/PhysRevC.71.024006

58. Gattobigio M, Kievsky A, Viviani M, Barletta P. Harmonic hyperspherical basis for identical particles without permutational symmetry. *Phys Rev A*. (2009) **79**:032513. doi: 10.1103/PhysRevA.79.032513

59. Gattobigio M, Kievsky A, Viviani M. Spectra of helium clusters with up to six atoms using soft-core potentials. *Phys Rev A*. (2011) **84**:052503. doi: 10.1103/PhysRevA.84.052503

60. Deflorian S, Barnea N, Leidemann W, Orlandini G. Nonsymmetrized hyperspherical harmonics with realistic NN potentials. *Few Body Syst*. (2013) **54**:1879. doi: 10.1007/s00601-013-0717-y

61. Nannini A, Marcucci LE. Non-symmetrized hyperspherical harmonics method for non-equal mass three-body systems. *Front Phys*. (2018) **6**:122. doi: 10.3389/fphy.2018.00122

62. Chen CR, Payne GL, Friar JL, Gibson BF. Faddeev calculations of the 2π−3*N* force contribution to the ^{3}*H* binding energy. *Phys Rev C*. (1986) **33**:1740. doi: 10.1103/PhysRevC.33.1740

63. Kohn W. Variational methods in nuclear collision problems. *Phys Rev*. (1948) **74**:1763. doi: 10.1103/PhysRev.74.1763

64. Kievsky A. The complex Kohn variational method applied to *N* − *d* scattering. *Nucl Phys A*. (1997) **624**:125. doi: 10.1016/S0375-9474(97)81832-5

65. Pudliner BS, Pandharipande VR, Carlson J, Wiringa RB. Quantum Monte Carlo calculations of *A* ≤ 6 nuclei. *Phys Rev Lett*. (1995) **74**:4396. doi: 10.1103/PhysRevLett.74.4396

66. Navrátil P. Local three-nucleon interaction from chiral effective field theory. *Few Body Syst*. (2007) **41**:117. doi: 10.1007/s00601-007-0193-3

67. Kievsky A, Viviani M, Rosati S. Study of bound and scattering states in three-nucleon systems. *Nucl Phys A*. (1994) **577**:511. doi: 10.1016/0375-9474(94)90931-8

68. Kievsky A, Rosati S, Tornow W, Viviani M. Critical comparison of experimental data and theoretical predictions for N-d scattering below the breakup threshold. *Nucl Phys A*. (1996) **607**:402. doi: 10.1016/0375-9474(96)00240-0

69. Coon SA, Glöckle W. The two pion exchange three nucleon potential: partial wave analysis in momentum space. *Phys Rev C*. (1981) **23**:1790. doi: 10.1103/PhysRevC.23.1790

70. Coon SA, Han HK. Reworking the Tucson-Melbourne three-nucleon potential. *Few Body Syst*. (2001) **30**:131. doi: 10.1007/s006010170022

71. Bogner SK, Furnstahl RJ, Ramanan S, Schwenk A. Low-momentum interactions with smooth cutoffs. *Nucl Phys A*. (2007) **784**:79. doi: 10.1016/j.nuclphysa.2006.11.123

72. Entem DR, Machleidt R. Accurate charge-dependent nucleon-nucleon potential at fourth order of chiral perturbation theory. *Phys Rev C*. (2003) **68**:041001. doi: 10.1103/PhysRevC.68.041001

73. Entem DR, Machleidt R, Nosyk Y. High-quality two-nucleon potentials up to fifth order of the chiral expansion. *Phys Rev C*. (2017) **96**:024004. doi: 10.1103/PhysRevC.96.024004

74. Marcucci LE, Sammarruca F, Viviani M, Machleidt R. Momentum distributions and short-range correlations in the deuteron and ^{3}He with modern chiral potentials. *Phys Rev C*. (2019) **99**:034003. doi: 10.1103/PhysRevC.99.034003

75. Piarulli M, Girlanda L, Schiavilla R, Pérez RN, Amaro JE, Arriola ER. Minimally nonlocal nucleon-nucleon potentials with chiral two-pion exchange including Δ-resonances. *Phys Rev C*. (2015) **91**:024003. doi: 10.1103/PhysRevC.91.024003

76. Baroni A, Schiavilla R, Marcucci LE, Girlanda L, Kievsky A, Lovato A, et al. Local chiral interactions, the tritium Gamow-Teller matrix element, and the three-nucleon contact term. *Phys Rev C*. (2018) **98**:044003. doi: 10.1103/PhysRevC.98.044003

77. Stoks VGJ, de Swart JJ. Magnetic moment interaction in nucleon-nucleon phase-shift analyses. *Phys Rev C*. (1990) **42**:1235. doi: 10.1103/PhysRevC.42.1235

78. Wang M, Audi G, Wapstra AH, Kondev FG, MacCormick M, Xu X, et al. The Ame2012 atomic mass evaluation. *Chin Phys C*. (2012) **36**:1603. doi: 10.1088/1674-1137/36/12/003

79. Marcucci LE, Kievsky A, Rosati S, Schiavilla R, Viviani M. Chiral effective field theory predictions for muon capture on deuteron and ^{3}He. *Phys Rev Lett*. (2012) **108**:052502. [Erratum: (2018) Phys. Rev. Lett. 121, 049901]. doi: 10.1103/PhysRevLett.108.052502

80. Gårdestig A, Phillips DR. How low-energy weak reactions can constrain three-nucleon forces and the neutron-neutron scattering length. *Phys Rev Lett*. (2006) **96**:232301. doi: 10.1103/PhysRevLett.96.232301

81. Gazit D, Quaglioni S, Navrátil P. Three-nucleon low-energy constants from the consistency of interactions and currents in chiral effective field theory. *Phys Rev Lett*. (2009) **103**:102502. [Erratum: (2019) Phys. Rev. Lett. 122, 029901]. doi: 10.1103/PhysRevLett.103.102502

82. Schoen K, Jacobson DL, Arif M, Huffman PR, Black TC, Snow WM, et al. Precision neutron interferometric measurements and updated evaluations of the *n* − *p* and *n* − *d* coherent neutron scattering lengths. *Phys Rev C*. (2003) **67**:044005. doi: 10.1103/PhysRevC.67.044005

83. Dilg W, Koester L, Nistler W. The neutron-deuteron scattering lengths. *Phys Lett B*. (1971) **36**:208. doi: 10.1016/0370-2693(71)90070-0

84. Black TC, Karwowski HJ, Ludwig EJ, Kievsky A, Rosati S, Viviani M. Determination of proton-deuteron scattering lengths. *Phys Lett B*. (1999) **471**:103. doi: 10.1016/S0370-2693(99)01366-0

85. Kievsky A, Viviani M, Girlanda L, Marcucci LE. Comparative study of three-nucleon force models in *A* = 3, 4 systems. *Phys Rev C*. (2010) **81**:044003. doi: 10.1103/PhysRevC.81.044003

86. Shimizu S, Sagara K, Nakamura H, Maeda K, Miwa T, Nishimori N, et al. Analyzing powers of *p* + *d* scattering below the deuteron breakup threshold. *Phys Rev C*. (1995) **52**:1193. doi: 10.1103/PhysRevC.52.1193

87. Girlanda L, Kievsky A, Viviani M. Subleading contributions to the three-nucleon contact interaction. *Phys Rev C*. (2011) **84**:014001. doi: 10.1103/PhysRevC.84.014001

88. Girlanda L, Kievsky A, Viviani M, Marcucci LE. Short-range three-nucleon interaction from A=3 data and its hierarchical structure. *Phys Rev C*. (2019) **99**:054003. doi: 10.1103/PhysRevC.99.054003

89. Kievsky A, Wood MH, Brune CR, Fisher BM, Karwowski HJ, Leonard DS, et al. Evidence for three nucleon force effects in proton–deuteron elastic scattering. *Phys Rev C*. (2001) **63**:024005. doi: 10.1103/PhysRevC.63.024005

90. Wood MH, Brune CR, Fisher BM, Karwowski HJ, Leonard DS, Ludwig EJ, et al. Low-energy *p* − *d* scattering: High precision data, comparisons with theory, and phase shift analyses. *Phys Rev C*. (2002) **65**:034002. doi: 10.1103/PhysRevC.65.034002

91. Phillips TW, Berman BL, Seagrave JD. Neutron total cross section for tritium. *Phys Rev C.* (1980) **22**:384. doi: 10.1103/PhysRevC.22.384

92. Seagrave JD, Cranberg L, Simmons JE. Elastic scattering of fast neutrons by tritium and he^{3}. *Phys Rev.* (1960) **119**:1981. doi: 10.1103/PhysRev.119.1981

93. Daniels TV, Arnold CW, Cesaratto JM, Clegg TB, Couture AH, Karwowski HJ, et al. Spin-correlation coefficients and phase-shift analysis for *p*+^{3}He elastic scattering. *Phys Rev C*. (2010) **82**:034002. doi: 10.1103/PhysRevC.82.034002

94. Machleidt R, Entem DR. Chiral effective field theory and nuclear forces. *Phys Rep*. (2011) **503**:1. doi: 10.1016/j.physrep.2011.02.001

95. Fisher BM, Brune CR, Karwowski HJ, Leonard DS, Ludwig EJ, Black TC, et al. Proton-He^{3} elastic scattering at low energies. *Phys Rev C.* (2006) **74**:034001. doi: 10.1103/PhysRevC.74.034001

96. Alley MT, Knutson LD. Spin correlation measurements for p-^{3}He elastic scattering between 4.0 and 10.0 MeV. *Phys Rev C.* (1993) **48**:1890. doi: 10.1103/PhysRevC.48.1890

97. Viviani M, Kievsky A, Rosati S, George EA, Knutson LD. The Ay problem for p-^{3}He elastic scattering. *Phys Rev Lett.* (2001) **86**:3739. doi: 10.1103/PhysRevLett.86.3739

98. Viviani M, Girlanda L, Kievsky A, Marcucci LE. Three-nucleon force effects in $\overrightarrow{p}{-}^{3}H$ and $\overrightarrow{n}{-}^{3}$ He scattering. *Few Body Syst*. (2017) **58**:110. doi: 10.1007/s00601-017-1264-8

99. Viviani M, Girlanda L, Kievsky A, Marcucci LE, Dohet-Eraly J. Three-nucleon force effects in *A* = 4 scattering. *Few Body Syst*. (2018) **59**:73. doi: 10.1007/s00601-018-1379-6

100. Hemmendinger A, Jarvis GA, Taschek RF. Scattering of protons by tritons. *Phys Rev*. (1949) **76**:1137. doi: 10.1103/PhysRev.76.1137

101. Claassen RS, Brown RJS, Freier GD, Stratton WR. The scattering of protons by tritons. *Phys Rev*. (1951) **82**:589. doi: 10.1103/PhysRev.82.589

102. Balashko YG, Barit IY, Dulkova LS, Kurepin AB. Elastic scattering protons on the tritium by energy below (p,n) reactions. *Bull Russ Acad Sci Phys*. (1965) **28**:1028.

103. Manduchi C, Moschini G, Tornielli G, Zannoni G. Experimental study on low-energy levels of ^{4}He. *Nuovo Cim B*. (1968) **57**:340. doi: 10.1007/BF02710206

104. Ivanovich M, Young PG, Ohlsen GG. Elastic scattering of several hydrogen and helium isotopes from tritium. *Nucl Phys A*. (1968) **110**:441. doi: 10.1016/0375-9474(68)90552-6

105. Kankowsky R, Fritz JC, Kilian K, Neufert A, Fick D. Elastic scattering of polarized protons on tritons between 4 and 12 MeV. *Nucl Phys A*. (1976) **263**:29. doi: 10.1016/0375-9474(76)90181-0

106. Willard HB, Bair JK, Kington JD. The reactions T(*p*, *n*)He^{3} and T(*p*, γ)He^{4}. *Phys Rev*. (1953) **90**:865. doi: 10.1103/PhysRev.90.865

107. Jarvis GA. Charged particle cross sections, the reaction T(*p*, *n*)He^{3}. *Los Alamos Sci Rep*. (1956) **2014**:35.

108. Drosg M. The ^{3}H(*p*, *n*)^{3}He differential cross sections below 5 MeV and the *n* − ^{3}He cross section. *Los Alamos Sci Rep*. (1980) **1980**:8215. doi: 10.2172/5207437

109. Doyle MA, Clark HW, Dries LJ, Regner JL, Donoghue TR, Hale GM. Comparison of analyzing power and polarization in the reaction ^{3}H(*p*, *n*)^{3} He. *Nucl Phys A*. (1981) **371**:225. doi: 10.1016/0375-9474(81)90065-8

110. Tornow W, Byrd RC, Lisowski PW, Walter RL, Donoghue TR. Comparison of analyzing power and polarization in the reaction ^{3}H(*p*, *n*)^{3} He. *Nucl Phys A*. (1981) **371**:235. doi: 10.1016/0375-9474(81)90066-X

111. Marcucci LE, Viviani M, Schiavilla R, Kievsky A, Rosati S. Electromagnetic structure of *A* = 2 and 3 nuclei and the nuclear current operator. *Phys Rev C*. (2005) **72**:014001. doi: 10.1103/PhysRevC.72.014001

112. Marcucci LE, Gross F, Peña MT, Piarulli M, Schiavilla R, Sick I, et al. Electromagnetic structure of few-nucleon ground states. *J Phys G Nucl Part Phys*. (2016) **43**:023002. doi: 10.1088/0954-3899/43/2/023002

113. Theeten M, Matsumura H, Orabi M, Baye D, Descouvemont P, Fujiwara Y, et al. Three-body model of light nuclei with microscopic nonlocal interactions. *Phys Rev C*. (2007) **76**:054003. doi: 10.1103/PhysRevC.76.054003

Keywords: hyperspherical harmonics method, ab initio methods, nuclear interactions, few-nucleon systems, light nuclei, *A* = 3, 4 scattering

Citation: Marcucci LE, Dohet-Eraly J, Girlanda L, Gnech A, Kievsky A and Viviani M (2020) The Hyperspherical Harmonics Method: A Tool for Testing and Improving Nuclear Interaction Models. *Front. Phys.* 8:69. doi: 10.3389/fphy.2020.00069

Received: 29 November 2019; Accepted: 28 February 2020;

Published: 07 April 2020.

Edited by:

Nunzio Itaco, University of Campania Luigi Vanvitelli, ItalyReviewed by:

Sonia Bacca, Johannes Gutenberg University Mainz, GermanyAndreas Nogga, Julich Research Centre, Germany

Winfried Leidemann, University of Trento, Italy

Copyright © 2020 Marcucci, Dohet-Eraly, Girlanda, Gnech, Kievsky and Viviani. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Laura E. Marcucci, laura.elisa.marcucci@unipi.it