The Hyperspherical Harmonics method: a tool for testing and improving nuclear interaction models

The Hyperspherical Harmonics (HH) method is one of the most accurate techniques to solve the quantum mechanical problem for nuclear systems with $A\le 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 problem. 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 five-to-ten years.


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 ith 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 mesontheoretical 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 Abody 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][10][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][21][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][32][33][34][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][39][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.

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 threenucleon interactions, but its extension to the non-local case presents no conceptual difficulties.

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 = A i=1 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 = j i=1 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 partsx i = (θ i , φ 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 2 N ( 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 2 N ( N ), the so-called hyperspherical harmonics (HH). They can be defined as Here Y l i (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 grandangular 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 K j−1 ,l j n j (ϕ j ) in Equation (9) are defined by (j) P K j−1 ,l j n j = N l j ,ν j n j (cos ϕ j ) l j (sin ϕ j ) K j−1 P ν j−1 ,l j +1/2 n j (cos 2ϕ j ), (11) where P ν j−1 ,l j +1/2 n j (cos 2ϕ j ) are Jacobi polynomials [52], with and the normalization factors N l,ν n 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 + 2n. 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 2 N ( 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 dx 1 · · · dx N = ρ 3N−1 dρ d N (19) and sin θ j dθ j dφ j (cos ϕ j ) 2 (sin ϕ j ) 3j−4 dϕ j . (20) 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 = i=1,N (l i + 2n 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 The transformation coefficients a KL,p [K];[K ′ ] do not depend on the quantum number M L . For A = 3, they are called the Raynal-Revai coefficients [53]. To be remarked that . . , n ′ N ], 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 χ SM S [S] with total spin S and total spin projection M S and the isospin functions ξ TM T [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 Including the spin and isospin functions, the HH basis functions read where J is the total angular momentum, J z its projection, and ( 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, (q) N ) 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].

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.

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 N 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).
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 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 lowenergy 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 αn 2 (ρ) (u αn 2 n 3 (ρ) 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 [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 αn 2 (ρ) and u αn 2 n 3 (ρ) replaced with w αn 2 (Q) and w αn 2 n 3 (Q), i.e., functions of the hypermomentum Q, while the HH functions are expressed in terms of conjugate Jacobi momenta. The wfunctions are related to the u-functions as where J K+ D 2 −1 (Qρ) 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 momentumspace). 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][59][60][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 (37) where the matrix elements of the Hamiltonian H and of the identity operator 1 1 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.

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 LSJJ z NY describing the N + Y scattering state with incoming orbital angular momentum L, channel spin S ≡ 1 2 + S Y , parity π = (−1) L , and total angular momentum J, J z , is written as Here we have introduced LSJJ z C , 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 LSJJ z A , 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 LSJJ z C 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 LSJJ z A 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 λ LSJJ z , χ 1/2 (N) and φ S Y (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 (ŷ p ) is the standard spherical harmonic function, and the functions R λ L (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µ, µ being the N − Y reduced mass, η = Z N Z Y µ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 (η) has been introduced so that the functions R R L (y p ) and R I L (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 I L (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 (2L + 1)!!C L (η) reduces to 1 for η → 0 [52].
With these definitions, LSJJ z A can be cast in the form where the parameters R J LS,L ′ S ′ (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 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 are obtained, depending on λ ≡ R, I, and consequently, we can introduce two core functions, defined as The matrix elements R J LS,L ′ S ′ (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) (50) and (53). Such secondorder 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 +,LSJJ z NY LSJJ z C being expanded as in Equation (39) and The functions λ LSJJ z 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 J LS,L ′ S ′ (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 LSJJ z C 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 λ LSJJ z , i.e., K, m|V| λ LSJJ z present in Equation (51), and λ ′ L ′ S ′ JJ z + L ′ S ′ JJ Z ,λ ′ C |V| λ LSJJ z of Equation (54), with λ, λ ′ = R, I. In particular, we note that, being λ LSJJ z 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 momentumspace 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].

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-toleading 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 lowenergy 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 twonucleon 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 twonucleon 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  3 He, and 4  The underlined values are used in the LECs fitting procedure. In the last row, we show the 3 H ( 3 He) experimental binding energy of 8.482 MeV (7.718 MeV), lowered (increased) by 7 keV in order to take into account the n − p mass difference. See text for more details. All the results presented here are in very good agreement with the values reported in the literature. The experimental binding energies are taken from reference [78]. The experimental uncertainty is well below the 1 keV level, and therefore it is not quoted in the table.
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.

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 threenucleon 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. 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 twoand 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 The experimental value for 2 a nd is from reference [82], that for 4 a nd is from reference [83], while those for 2 a pd and 4 a pd are from reference [84]. 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 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][80][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]. 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.

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 nonlocal 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  (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.  The different number N of experimental data is also indicated. The data are from references [89,90] at Ec.m. = 0.666 MeV, and from reference [86] at Ec.m. = 1.33, 1.66, and 2.0 MeV.
Frontiers in Physics | www.frontiersin.org 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 energydependent 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 phaseshift 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 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). 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]).
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 twoand 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 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 twoand three-nucleon interaction (dashed-dotted blue lines), and to the experimental data of reference [90]. 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.
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 (61) 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 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].
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 centerof-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 4 π ), E i = e i (F 4 π 3 ), i = 5, 7, F π = 92.4 MeV being the pion decay constant, are of order 1 as expected.
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.

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

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 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][96][97]. See text for more details. 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.
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 twonucleon interaction. Therefore, it could be used for a more stringent test of the two-nucleon force.

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 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][107][108][109][110].
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 threenucleon 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.