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

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 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 meson-theoretical approach, as the Argonne v₁₈ [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¹. 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 ⁴He in ⁴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₁, …, m_A and spatial coordinates r₁, …, 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₁, …, 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₁ = 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₁, …, l_N, L₂, …, L_{N − 1}, n₂, …, 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₁ = 0, there is only one index n₂ ≡ n. In this case, K = l₁ + l₂ + 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 ${\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+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 ${{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₂, …, S_{N − 1}] and [T] for [T₂, …, 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 ,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 φ₂ for A = 3 and φ₂, φ₃ 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₂ in Equation (30) and n₂, n₃ 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 ≡ 3N − 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⁻¹ for local and 4–8 fm⁻¹ 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{Ψ}}_{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{Ψ}}_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{Ψ}}_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{Ψ}}_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{Ψ}}_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⁻¹ 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, ${\text{Ψ}}_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{Ψ}}_{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{Ψ}}_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{Ψ}}_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{Ψ}}_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₁₈ (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⁻¹, 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 α² (α 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 ⁴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.

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 ³H and ³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 ³H and −7 keV in ³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 ³H and one for ³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 ²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 ³H binding energy, but not exactly overlapping. This is a well-known fact, that the ³H binding energy and ²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 (⁴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 ⁴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 ²a_nd, we notice a systematic overestimation of the ³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⁻¹, 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 ²a_nd, and the older experimental results from reference [83] for ⁴a_nd. Though no experimental data are available for ²a_pd and ⁴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, ⁴a_nd remains practically unchanged, while ²a_nd becomes smaller. For the NVIa and NVIb potentials, for instance, ²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 ⁴a_nd is very little model-dependent and in good agreement with experiment, the same is not true for ²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]).

With the purpose of investigating the capability of some widely used models of three-nucleon interaction to reproduce ²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 ²a_nd and the triton binding energy, and the resulting value for the ⁴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₁₁. 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.

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 χ²/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₁₁. 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₀(r; Λ) are written as

with F(k²; Λ) 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 ²a_nd and ⁴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 χ²/datum of Table 3 are the spin-orbit and tensor interactions, which depend on the LECs E₅ and E₇. We present in Table 4 the results of a fit where only the terms proportional to E₀, E₅, and E₇ are kept. The LEC E₀ 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 χ²/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 χ²/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₀, e₅, e₇, defined in terms of E₀, E₅, E₇ 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.

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.

3.3. p−³He and n − ³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 − ³H [32] and p−³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 − ³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 − ³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−³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.

3.4. p−³H and n − ³He Scattering

The treatment of p−³H and n − ³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−³H differential cross section, for which, only at very low energies, below the opening of the n − ³He channel, some sizable effects are visible. Otherwise, the three-nucleon interaction contributions are found very small. The p−³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+³H → n+³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.

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

References