Implementation of Local Chiral Interactions in the Hyperspherical Harmonics Formalism

With the goal of using chiral interactions at various orders to explore the properties of the few-body nuclear systems, we write the recently developed local chiral interactions as spherical irreducible tensors and implement them in the hyperspherical harmonics expansion method. We devote particular attention to three-body forces at next-to-next-to leading order, which play an important role in reproducing experimental data. We check our implementation by benchmarking the ground-state properties of ³H, ³He, and ⁴He against the available Monte Carlo calculations. We then confirm their order-by-order truncation error estimates and further investigate uncertainties in the charge radii obtained by using the precise muonic atom data for single-nucleon radii. Having local chiral Hamiltonians at various orders implemented in our hyperspherical harmonics suites of codes opens up the possibility to test such interactions on other light-nuclei properties, such as electromagnetic reactions.

1. Introduction

In 1935 the seminal idea of Yukawa [1] laid the foundation to the theory of the nuclear forces. His one-pion exchange term is nowadays known as an important contribution to the interaction among nuclei in the long-distance range and is implemented in many nuclear interaction models. In the mid 1990s the first high-precision nucleon-nucleon (NN) potentials able to reproduce at the same time the deuteron properties, the proton-proton, and the proton-neutron scattering data were released. Some notable examples of these interactions are the Argonne v18 (AV18) [2], the Nijmegen (Nijm93) [3], and the charge-dependent Bonn (CD-Bonn) [4]. The subsequent development of three-nucleon (3N) interactions, see for instance references [5, 6], improved the description of the A > 2 nuclear dynamics, initiating a successful theoretical campaign of nuclear structure and reaction predictions, see e.g., references [7–9] and references therein. Despite the great success of the phenomenological interactions, there are still open questions to address, including the difficulty of providing solid uncertainty quantifications in the modeling of the forces, the lack of connection between the NN and 3N interactions and the missing direct link to quantum chromodynamic (QCD), the fundamental theory of the strong force.

An important step forward to address these issues was made when the concept of effective field theory (EFT) was introduced and applied to low-energy QCD. As suggested by Weinberg [10–13], the low-energy nuclear dynamic can be described by a Lagrangian written in terms of pions and nucleons fields and consistent with all the commonly accepted symmetries of QCD, including the (explicitly and spontaneously broken) chiral symmetry which strongly constrains the pion dynamics. The proposed Lagrangian contains an infinite number of terms and a systematic expansion must be introduced to make the theory applicable. Following Weinberg's proposal, in the early 2000s, modern versions of chiral-inspired nuclear interactions were released by many groups—for a compilation of results, see for instance references [14–16] and references therein—each interaction being different by the truncation order of the chiral expansion, by the inclusion or exclusion of the Δ-isobar, by the fitting procedure or by the regularization scheme used. Given that these interactions are derived in field theories written in momentum space, they are highly non-local. One of the consequences is that they are difficult to implement in some of the few- and many-body techniques which are developed in coordinate-space representation.

In recent years a new chiral-inspired set of nuclear interactions at the next-to-next-to-leading order (N2LO) has become available [17–20]. These interactions have a series of interesting properties which make them a promising framework for future nuclear computations. These interactions are completely written in coordinate space and contain only one non-local operator. Furthermore, the NN and 3N terms are regularized consistently, namely the same regulator form and cut-off is used. Here, these interactions are written for the first time as a product of irreducible tensors under space rotations, a required step for the implementation into the hyperspherical harmonics formalism. Using the method of hyperspherical harmonics, we perform benchmark tests in light-nuclear systems, where we compare to available results from the Green's function Monte Carlo (GFMC) and the auxiliary field diffusion Monte Carlo (AFDMC) methods. We note that analogous benchmarks between the hyperspherical harmonic and quantum monte-carlo methods have been successfully performed in recent years for an alternative set of local-chiral interactions that explicitly includes the excitation of the Δ-resonance [21–23].

This paper is summarized as follows. In section 2, we briefly overview the formulation of the hyperspherical harmonics method in coordinate-space representation. In section 3, we present the maximally-local chiral interactions developed in reference [20] and rewrite the 3N force as a product of irreducible tensors under space rotations. In section 4, we show our benchmark results for ³H, ³He, and ⁴He and we discuss uncertainties. Finally, section 5 is reserved for the conclusive remarks and the overview of future prospects.

2. Hyperspherical Harmonics

The hyperspherical harmonic method was firstly introduced in 1935 by Zernike and Brinkman [24], reintroduced later in the 60's by Delves [25], Simonov [26], Zickendraht [27], and Smith [28] and it is extensively applied nowadays to the study of few-body systems. For recent reviews with applications to nuclear physics we refer the reader to the following references [29, 30]. In this work, the hyperspherical harmonic functions are constructed to form irreducible representations of the SO(3) group of space rotations, the O(N) group of dynamical rotations in the space spanned by the N Jacobi vectors, and the S_A permutation group of the A-particle system. The method is briefly reviewed in this section, the formalism introduced follows closely references [31, 32].

We consider a system of A identical nucleons, the Jacobi coordinates {η_i} are commonly introduced in order to separate the internal degrees of freedom from the center of mass. There are several ways to construct the set of N = A − 1 Jacobi coordinates out of the A coordinate vectors {r_i} of the nucleons. One commonly used definition for the relative Jacobi vectors is

$\begin{array}{ll}{\eta }_{j-1}=\sqrt{\frac{j-1}{j}}\left(r_j-\frac{1}{j-1}{\displaystyle \sum _{i=1}^{j-1}}{r}_i\right);j=2,...,A.& (1)\end{array}$

From a given choice of Jacobi coordinates, the hyperspherical coordinates {ρ_N, φ(N), Ω(N)} can be introduced. In this notation, ρ_N is the hyper-radius, Ω(N) ≡ {Ω₁, ..., Ω_N} where Ω_j = (θ_j, ϕ_j) gathers the angular coordinates of the Jacobi vectors, and φ(N) ≡ {φ₂, ..., φ_N} is a set of hyper-angles.

The hyper-radial coordinates ρ₁, ..., ρ_N and the hyper-angular coordinates φ₂, ..., φ_N are constructed recursively. The transformation law for the first two Jacobi coordinates is

$\begin{array}{r}{\eta }_1={\rho }_1={\rho }_{2}cos{\phi }_2,\\ {\eta }_2={\rho }_{2}sin{\phi }_2.& (2)\end{array}$

Assuming that we already know the hyper-radial coordinates ρ₁, ..., ρ_{j − 1} and the hyper-angular coordinates φ₂, ..., φ_j−1 the transformation law for ρ_j and φ_j reads in analogy to Equation (2) as

$\begin{array}{r}{\rho }_{j-1}={\rho }_{j}cos{\phi }_j,\\ {\eta }_j={\rho }_{j}sin{\phi }_j.& (3)\end{array}$

The internal kinetic energy operator for the A-body system is given by the 3N-dimensional Laplace operator Δ(N). In terms of the hyperspherical coordinates it is written as

$\begin{array}{ll}\Delta (N)={\Delta }_{\rho }-\frac{1}{{\rho }^2}{\widehat{K}}_N^2(\phi (N),\Omega (N))& (4)\end{array}$

where the hyper-radial part is

$\begin{array}{ll}{\Delta }_{\rho }=\frac{{\partial }^2}{\partial {\rho }^2}+\frac{3N-1}{\rho }\frac{\partial }{\partial \rho },& (5)\end{array}$

with ρ ≡ ρ_N. ${\widehat{K}}_N^2$ is the grand-angular momentum operator whose eigenfunctions are known as the hyperspherical harmonics.

Denoting ${\widehat{l}}_j$ as the angular momentum operator related to η_j, and ${\widehat{L}}_j^2$ and ${\widehat{M}}_j$ as the total orbital angular momentum operator and z-projection of the system identified by the first j Jacobi coordinates, it is possible to define the grand-angular momentum operator ${\widehat{K}}_N^2$ of the system recursively in terms of ${\widehat{K}}_{N-1}^2$ and ${\widehat{l}}_N$ as [33]

$\begin{array}{l}{\widehat{K}}_N^2=-\frac{{\partial }^2}{\partial {\phi }_N^2}+\frac{3N-6-(3N-2)cos(2{\phi }_N)}{sin(2{\phi }_N)}\frac{\partial }{\partial {\phi }_N}\\ +\frac{1}{{cos}^2{\phi }_N}{\widehat{K}}_{N-1}^2+\frac{1}{{sin}^2{\phi }_N}{\widehat{l}}_N^2& (6)\end{array}$

where ${\widehat{K}}_1^2={\widehat{l}}_1^2$.

The operators ${\widehat{K}}_N^2,...,{\widehat{K}}_2^2$, ${\widehat{L}}_N^2,...,{\widehat{L}}_2^2$, ${\widehat{l}}_N^2,...,{\widehat{l}}_1^2$, and ${\widehat{M}}_N$ commute with each other. As a consequence, it is possible to label hyperspherical states using the set of 3N − 1 quantum numbers {K} ≡ {K_N, ..., K₂, L_N, ..., L₂, l_N, ..., l₁, M_N}. The hyperspherical harmonics functions ${{Y}}_{\left\{K_N\right\}}$ are the eigenfunctions of the grand-angular momentum operator with eigenvalues K_N(K_N + 3N − 2). The explicit expression for the resulting hyperspherical harmonics functions is given by [34]

$\begin{array}{l}{{Y}}_{\left\{K_N\right\}}=[{\displaystyle \sum _{m_1,...,m_N}}{C}_{l_1{m}_1,l_2{m}_2}^{L_2{M}_2}{C}_{L_2{M}_2,l_3{m}_3}^{L_3{M}_3}\times ...\\ \times C_{L_{N-1}{M}_{N-1},l_N{m}_N}^{L_N{M}_N}{\displaystyle \prod _{j=1}^N}{Y}_{l_j{m}_j}({\Omega }_j)]\times \\ \times \left[{\displaystyle \prod _{j=2}^N}{{N}}_j{(sin{\phi }_j)}^{l_j}{(cos{\phi }_j)}^{K_{j-1}}{P}_{n_j}^{\left[l_j+\frac{1}{2},K_{j-1}+\frac{(3j-5)}{2}\right]}(cos2{\phi }_j)\right],& (7)\end{array}$

where $C_{l_i{m}_i,l_j{m}_j}^{LM}$ are the Clebsch-Gordan coefficients, Y_ljmj(Ω_j) are the spherical harmonics associated with η_j, $P_n^{\left[\alpha ,\beta \right]}$ are Jacobi polynomials and

$\begin{array}{ll}{{N}}_j={\left[\frac{(3K_j+3j-2)n_j!\Gamma (n_j+K_{j-1}+l_j+\frac{3j-2}{2})}{\Gamma (n_j+l_j+\frac{3}{2})\Gamma (n_j+K_{j-1}+\frac{3j-3}{2})}\right]}^{\frac{1}{2}}& (8)\end{array}$

is a normalization constant with 2n_j = K_j − K_j−1 − l_j.

In our formulation of the hyperspherical harmonics method we construct hyper-angular functions that form irreducible tensors under the SO(3) group of spatial rotations, the O(N) group of kinematic rotations and the S_A group of permutations of the A nucleons. These symmetry-adapted hyperspherical harmonics, ${{Y}}_{\left[K_N\right]}$, are uniquely identified by the set of quantum numbers [K_N] ≡ {K_N, L_N, M_N, λ_N, α_NY_A, β_A}. For the current purposes, it is enough to specify that λ_N identifies the irreducible representation of O(N), Y_A is the Yamanouchi symbol which specifies the irreducible representations of the group-subgroup chain S₁ ⊂ … ⊂ S_A presented by the appropriate Young diagrams Γ₁, …, Γ_A, while α_N and β_A are additional quantum numbers needed to remove further degeneracies. The O(N) and S_A symmetry-adapted hyperspherical harmonics ${{Y}}_{\left[K_N\right]}$ are constructed recursively. Assuming that ${{Y}}_{\left[K_{N-1}\right]}$ have been already constructed, the Nth Jacobi coordinate is then coupled to this system, so that a state with total angular momentum L_N and grand-angular momentum K_N is formed, let us call this state ${{Y}}_{\left[K_{N-1}\right],K_N{L}_N{M}_N}$. Note that ${{Y}}_{\left[K_{N-1}\right],K_N{L}_N{M}_N}$ is a irreducible tensor under O(N − 1) and S_A−1 but not under O(N) and S_A. The states ${{Y}}_{\left[K_N\right]}$ are obtained as linear combinations of the states ${{Y}}_{\left[K_{N-1}\right],K_N{L}_N{M}_N}$, where the coefficients of the linear combinations are labeled as [(K_N−1, L_N−1, λ_N−1, α_N−1; l_N)K_NL_N|}K_NL_Nλ_Nα_N], [(λ_N−1Γ_Nβ_N)λ_N|}λ_NΓ_Aβ_A] and are known as hyperspherical orthogonal group parentage coefficients (HSOPCs) and orthogonal group coefficients of fractional parentage (OCFPs), respectively.

The full expression of the symmetry-adapted hyperspherical harmonics reads

$\begin{array}{l}{{Y}}_{\left[K_N\right]}={\displaystyle \sum _{{\lambda }_{N-1}{\beta }_N}}\left[\left({\lambda }_{N-1}{\Gamma }_N{\beta }_N\right){\lambda }_N\left|\right\}{\lambda }_N{\Gamma }_A{\beta }_A\right]\times \\ \times {\displaystyle \sum _{K_{N-1},L_{N-1},{\alpha }_{N-1},l_N}}\left[\left(K_{N-1},L_{N-1},{\lambda }_{N-1},{\alpha }_{N-1};l_N\right)K_N{L}_N|\right\}\\ K_N{L}_N{\lambda }_N{\alpha }_N]\times \\ \times {{Y}}_{\left[K_{N-1}\right],K_N{L}_N{M}_N}.& (9)\end{array}$

Nucleons also possess spin and isospin degrees of freedom. Because the nuclear Hamiltonian is rotationally invariant, nuclear states have the total angular momentum J as good quantum number. Furthermore, isospin is an approximate symmetry for the nuclear interaction with the consequence that the total isospin T of a nuclear state is a conserved quantum number. For these reasons we couple the symmetry-adapted hyperspherical harmonics to the S_A symmetry-adapted spin-isospin wavefunction χ of the A-nucleon system

$\begin{array}{ll}{H}_{(K_N)}={\displaystyle \sum _{Y_N}}\frac{{\Lambda }_{{\Gamma }_A,Y_N}}{\sqrt{\left|{\Gamma }_A\right|}}{\displaystyle \sum _{M_N{S}_z}}{C}_{L_N{M}_N,S{S}_z}^{J{J}_z}{{Y}}_{\left[K_N\right]}{\chi }_{\left[S_A\right]}.& (10)\end{array}$

Here $(K_N)\equiv \left\{K_N,L_N,S_N,J_N,{(J_N)}_z,{\lambda }_N,{\alpha }_N^{ST},Y_A,{\beta }_A\right\}$, $\left[S_A\right]\equiv \left\{S,S_z,T,T_z,Y_A,{\alpha }_A^{ST}\right\}$, Λ_ΓA,YN is a phase factor, and |Γ_A| is the dimension of the irreducible representation Γ_A.

Analogously to what has been done with the hyperspherical harmonics, the spin-isospin wavefunctions are constructed recursively. Assuming that the symmetry-adapted wavefunction χ_[Sj−1] have been obtained, the construction of the χ_[Sj] is done by first coupling χ_[Sj−1] to the spin-isospin wavefunction of the jth nucleon, let us call this state χ_{[Sj−1], SjTj}, and then taking linear combinations of χ_{[Sj−1], SjTj} using the coefficients of fractional parentage labeled as $\left[S_{j-1}{S}_j{T}_{j-1}{T}_j{\Gamma }_{j-1}{\alpha }_{j-1}^{ST}|\right\}{S}_j{T}_j{\Gamma }_j{\alpha }_j^{ST}]$. Namely the full expression for χ_[Sj] reads

$\begin{array}{ll}{\chi }_{[S_j]}={\displaystyle \sum _{S_{j-1}{T}_{j-1}{\alpha }_{j-1}^{ST}}}\left[S_{j-1}{S}_j{T}_{j-1}{T}_j{\Gamma }_{j-1}{\alpha }_{j-1}^{ST}\left|\right\}{S}_j{T}_j{\Gamma }_j{\alpha }_j^{ST}\right]{\chi }_{[S_{j-1}],S_j{T}_j}.& (11)\end{array}$

We are finally able to expand the nuclear wavefunction in terms of hyperspherical harmonics. In practice, the expansion is performed up to a maximal value of the grand-angular quantum number K_max as

$\begin{array}{ll}\text{Ψ}={\displaystyle \sum _{(K_N)}^{K_{\text{max}}}}{{R}}_{(K_N)}({\rho }_N)H_{(K_N)}({\Omega }_N).& (12)\end{array}$

When we insert this wavefunction into the Schrödinger equation, an eigenvalues equation is obtained for the hyper-radial wavefunction ${{R}}_{(K_N)}$, the eigenvalue equation is then solved by expanding the hyper-radial wavefunction in terms of an orthogonal set of functions. In this work the set is taken as the generalized Laguerre polynomials $L_n^v({\rho }_N)$. Again, the model space is truncated to a given maximum number of Laguerre polynomials n_max

$\begin{array}{ll}{{R}}_{(K_N)}={\displaystyle \sum _{n=0}^{n_{max}}}{C}_{(K_N)}^n{L}_n^v({\rho }_N).& (13)\end{array}$

With the introduction of this further model space, the resulting eigenvalue equation is solved with direct diagonalization routines, or with the Lanczos method when the model space is too big for a direct diagonalization. In essence, the hyperspherical harmonics method is a powerful technique that allows for an exact solution of the Schrödinger equation for few-body systems. In the limit where n_max → ∞ and K_max → ∞ the solution correspond to the exact solution to the Schrödinger equation. While we observe that good convergence can be reached with n_max ≤ 50, the convergence in terms of K_max will be carefully investigated. The uncertainty coming from the truncation of the model space, in particular of K_max, can be estimated by looking at the convergence pattern of the observables of interest, for instance the binding energy and the radius. As a consequence, the method is an excellent candidate for uncertainty quantifications in nuclear physics, with the possibility of performing tests over commonly accepted nuclear Hamiltonians or making precise predictions for few-nucleon systems. Because the formulation we present here is developed in coordinate space, the method benefits from having local forces, such as the AV18 potential. While one can formulate hyperspherical harmonics also in momentum space [35], the goal of this paper is to work in coordinate space and implement local-chiral interactions. To further improve the convergence with respect to the model space, we make use of the effective interaction hyperspherical harmonics (EIHH) method. The interested reader can find more details on this approach in reference [36], and also in the more recent review [37].

3. Nuclear Hamiltonians

Nuclear physics is mainly formulated in the framework of non-relativistic quantum mechanics. The relevant degrees of freedom are represented by the nucleons, whose interactions are remnants of the color forces among the quarks. In this picture, the nucleus is a compound object of A non-relativistic nucleons and the dynamic of the system is specified by the nuclear Hamiltonian operator

$\begin{array}{ll}Ĥ=\widehat{T}+\widehat{V}+Ŵ+...={\displaystyle \sum _{i=1}^A}{\widehat{T}}_i+{\displaystyle \sum _{i>j=1}^A}{\widehat{V}}_{ij}+{\displaystyle \sum _{i>j>k=1}^A}{Ŵ}_{ijk}+...,& (14)\end{array}$

where $\widehat{T}$ is the sum of the non-relativistic kinetic energy operators of the individual nucleons, $\widehat{V}$ is a sum of NN interactions, and Ŵ is a sum of 3N interactions. The dots stand for higher order forces not explicitly included in this work.

Our goal is to solve the Schödinger equation

$\begin{array}{ll}Ĥ|\text{Ψ}\rangle =E|\text{Ψ}\rangle & (15)\end{array}$

and when working with antisymmetrized wavefunctions, the expectation values of the NN and 3N terms become

$\begin{array}{l}\langle \text{Ψ}|\widehat{V}|\text{Ψ}\rangle =\frac{A(A-1)}{2}\langle \text{Ψ}|{\widehat{V}}_{12}|\text{Ψ}\rangle ,\\ \langle \text{Ψ}|Ŵ|\text{Ψ}\rangle =\frac{A(A-1)(A-2)}{6}\langle \text{Ψ}|{Ŵ}_{123}|\text{Ψ}\rangle ,& (16)\end{array}$

where only the first two (or three) particles are involved.¹

In the modern theory of nuclear forces, interactions are derived from the chiral effective field theory (ChEFT). In this theory, proposed first by Weinberg [10–13], the chiral Lagrangian is constructed in terms of pion and nucleon fields and is consistent with the commonly accepted symmetries of QCD, including the explicitly and spontaneously broken chiral symmetry. This effective Lagrangian has infinitely many terms, therefore one needs to introduce an ordering scheme to render the theory predictive.

In ChEFT, the terms in the chiral Lagrangian are analyzed counting powers of a small external momentum over the large scale: ${(Q/{\Lambda }_{\chi })}^{\nu }$, where Q stands for an external momentum or a pion mass and Λ_χ is the chiral symmetry breaking scale, whose value is approximately given by the mass of the ρ-meson Λ_χ ~ m_ρ = 770 MeV. Determining systematically the power of ν has become known as power counting. The lowest possible value of ν is conventionally referred to as the leading order (LO), the second lowest is the next-to-leading order (NLO), the third lowest is the next-to-next-to leading order (N2LO), and so on. While there are many proposed power counting schemes [38–43], in this work we adopt the Weinberg power counting, which makes use of naive dimensional analysis [11, 12].

Given that ChEFT is naturally formulated in momentum space, the derived nuclear interactions are strongly non-local, which is a disadvantage for methods that are formulated in coordinate space. However, it has been recently found that it is possible to construct maximally local chiral interactions by regularizing in coordinate space and exploiting Fierz ambiguities to remove non-localities in the short-distance interactions [17–20].

The local chiral NN forces are composed of contact (ct) terms and pion-exchange (π) terms so that the interaction between particle 1 and 2 can be written as

$\begin{array}{ll}{V}_{12}=V_{12}^{ct}+V_{12}^{\pi }.& (17)\end{array}$

When working with totally anti-symmetric systems, it is possible to exploit Fierz ambiguities for removing the non-local operators contributing to the contact NN interactions. This means that the interactions can be chosen to have the following operator structure [17] at LO

$\begin{array}{ll}{V}_{12}^{ct,\text{LO}}=\left(C_{\text{S}}+C_{\text{T}}{\sigma }_1·{\sigma }_2\right)\delta (r_{12}),& (18)\end{array}$

where r₁₂ is the relative distance between nucleon 1 and nucleon 2, σ_1/2 are the vector-spin Pauli matrices operating in the space of the first/second nucleon² and δ is the delta function.

At NLO, the following new terms enter

$\begin{array}{l}{V}_{12}^{ct,\text{NLO}}=-\left(C_1+C_2{\tau }_1·{\tau }_2\right)\Delta \delta (r_{12})\\ -\left(C_3+C_4{\tau }_1·{\tau }_2\right){\sigma }_1·{\sigma }_2\Delta \delta (r_{12})\\ +\frac{C_5}{2}\frac{{\partial }_{r_{12}}\delta (r_{12})}{r_{12}}\text{L}·\text{S}+\left(C_6+C_7{\tau }_1·{\tau }_2\right)\\ \times [\left({\sigma }_1·{\widehat{r}}_{12}\right)\left({\sigma }_2·{\widehat{r}}_{12}\right)\left[\frac{{\partial }_{r_{12}}\delta (r_{12})}{r_{12}}-{\partial }_{r_{12}}^2\delta (r_{12})\right]\\ -{\sigma }_1·{\sigma }_2\frac{{\partial }_{r_{12}}\delta (r_{12})}{r_{12}}],& (19)\end{array}$

where τ_1/2 are the vector-isospin Pauli matrices, Δ is the Laplace operator, L and S are the total orbital angular momentum and spin operator in the two-body system represented by the two interacting nucleons³, ${\widehat{r}}_{12}$ is the unit vector related to the relative distance r₁₂, and the δ-function will have to be regularized. The {C_i} are a set of low energy constants (LECs). The term proportional to the LEC C₅ is the only non-local operator appearing in this maximally local chiral interaction.

Following references [17, 44], all the pion-exchange interactions up to N2LO can be written in a complete local form as

$\begin{array}{l}{V}_{12}^{\pi }=V_{\text{C}}(r_{12})+W_{\text{C}}(r_{12}){\tau }_1·{\tau }_2\\ +\left[V_{\text{S}}(r_{12})+W_{\text{S}}(r_{12}){\tau }_1·{\tau }_2\right]{\sigma }_1·{\sigma }_2\\ +\left[V_{\text{T}}(r_{12})+W_{\text{T}}(r_{12}){\tau }_1·{\tau }_2\right]S_{12},& (20)\end{array}$

where, S₁₂ is the well-known tensor operator, defined as

$\begin{array}{ll}{S}_{12}=3({\sigma }_1·{\widehat{r}}_{12})({\sigma }_2·{\widehat{r}}_{12})-({\sigma }_1·{\sigma }_2).& (21)\end{array}$

The local functions V_C(r₁₂), W_C(r₁₂), V_S(r₁₂), W_S(r₁₂), V_T(r₁₂), and W_T(r₁₂) have dependencies on the axial-vector coupling constant of the nucleon g_A, on the pion decay constant F_π and on the pion mass m_π. These functions are evaluated at each order in ChEFT (LO, NLO, and N2LO) and details can be found in reference [44]. In reference [20] pion loops are regularized using the spectral-function regularization (SFR) with an ultraviolet cut-off $\tilde{\Lambda }=1$ GeV and we follow this prescription.

The local chiral NN interactions up to N2LO are already written or can be written with minimal modifications as irreducible tensors under space rotations. Thus, they can be easily implemented in the hyperspherical harmonics formalism in coordinate space. In fact, they have pretty much the same structure as the Argonne potential AV8' [45] with the only difference being that, compared to this interaction, there is no isospin-dependent spin-orbit operator. The same does not apply to 3N interactions.

Three-body interactions arise at NLO in Weinberg power counting. However, at this order their contribution is canceled out. The first non-zero contributions start at N2LO. The 3N force at this order is composed of a two-pion (2π) exchange, a one-pion (1π) exchange and a 3N contact (ct) interaction (see Figure 1). On the one hand, the 2π-term comes with the LECs c₁, c₃, and c₄ that already appear at the subleading two-pion-exchange interaction at the NN level at the same chiral order which highlights the consistency of the NN and 3N interactions in ChEFT. On the other hand, the one-pion exchange and the 3N contact diagrams introduce two new LECs, c_D and c_E, which must be fitted on A ≥ 3 observables.

FIGURE 1

Figure 1. Feynman diagrams of the chiral 3N force at N2LO, from the left to the right: 2π-term, 1π-term, and ct-term.

With respect to reference [20], here the 3N interaction is written for a given triplet of nucleons, since at the end we use the fact that the wavefunction is anti-symmetric to compute the expectations values as in Equation (16). The 3N interaction reads

$\begin{array}{l}{W}_{123}={\displaystyle \sum _{\text{cyc}}}{W}_{1,23}\\ ={\displaystyle \sum _{\text{cyc}}}\left[W_{1,23}^{2\pi ,c_1}+W_{1,23}^{2\pi ,c_3}+W_{1,23}^{2\pi ,c_4}+W_{1,23}^{1\pi ,c_D}+W_{1,23}^{ct,c_E}\right],& (22)\end{array}$

where the sum runs over the cyclic permutations of the particle triplet and the notation has the intention to highlight the symmetry of the interaction over the exchange of particles 2 and 3. Each term is denoted with a label that includes the associated LEC.

The 2π exchange terms are given by

$\begin{array}{l}{W}_{1,23}^{2\pi ,c_1}=A{U}_{12}{Y}_{12}{U}_{13}{Y}_{13}({\tau }_2·{\tau }_3)({\sigma }_2·{\widehat{r}}_{12})({\sigma }_3·{\widehat{r}}_{13}),\\ W_{1,23}^{2\pi ,c_3}=B\left\{{\tau }_1·{\tau }_2,{\tau }_1·{\tau }_3\right\}\left\{{\chi }_{12},{\chi }_{13}\right\},\\ W_{1,23}^{2\pi ,c_4}=-C\left[{\tau }_1·{\tau }_2,{\tau }_1·{\tau }_3\right]\left[{\chi }_{12},{\chi }_{13}\right],& (23)\end{array}$

where the coupling constants are $A=c_1\frac{g_A^2{m}_{\pi }^4{(\hslash c)}^2}{16{\pi }^2{F}_{\pi }^4},B=c_3\frac{g_A^2{m}_{\pi }^4{(\hslash c)}^2}{1152{\pi }^2{F}_{\pi }^4},$ and $C=c_4\frac{g_A^2{m}_{\pi }^4{(\hslash c)}^2}{2304{\pi }^2{F}_{\pi }^4}$. The 2π-terms include the following functions

$\begin{array}{l}{Y}_{12}=Y(r_{12})=\frac{e^{-m_{\pi }{r}_{12}}}{r_{12}},\\ U_{12}=U(r_{12})=1+\frac{1}{m_{\pi }{r}_{12}},& (24)\end{array}$

with analogous expressions for Y₁₃ and U₁₃. The operator χ₁₂ (and analogously χ₁₃) is defined as

$\begin{array}{ll}{\chi }_{12}=X_{12}-\frac{4\pi }{m_{\pi }^2}{\delta }_{12}{\sigma }_1·{\sigma }_2=T_{12}{S}_{12}+{Ỹ}_{12}{\sigma }_1·{\sigma }_2& (25)\end{array}$

with

$\begin{array}{l}{X}_{12}=T_{12}{S}_{12}+Y_{12}{\sigma }_1·{\sigma }_2,\\ {Ỹ}_{12}=Y_{12}-\frac{4\pi }{m_{\pi }^2}{\delta }_{12},\\ T_{12}=T(r_{12})=(1+\frac{3}{m_{\pi }{r}_{12}}+\frac{3}{{(m_{\pi }{r}_{12})}^2})Y_{12},\\ {\delta }_{12}={\delta }_{r_0}(r_{12})=\frac{1}{\frac{4\pi }{n}\Gamma (\frac{3}{n})r_0^3}{e}^{-{(r_{12}/r_0)}^n}.& (26)\end{array}$

In the last expression, r₀ is the cut-off and following references [18–20] n is taken to be equal to 4.

For the 1π-interaction terms there are two options

$\begin{array}{l}{W}_{1,23}^{1\pi ,c_{D}1}=D({\tau }_2·{\tau }_3)[X_{23}(r_{12}){\delta }_{13}+X_{23}(r_{13}){\delta }_{12}\\ -\frac{8\pi }{m_{\pi }^2}{\delta }_{12}{\delta }_{13}{\sigma }_2·{\sigma }_3]\\ \text{\hspace{1em}\hspace{1em} \hspace{1em}and}\\ W_{1,23}^{1\pi ,c_{D}2}=D({\tau }_2·{\tau }_3){\chi }_{23}({\delta }_{12}+{\delta }_{13}),& (27)\end{array}$

with $D=c_D\frac{g_A{m}_{\pi }^2{(\hslash c)}^4}{96\pi {\Lambda }_{\chi }{F}_{\pi }^4}$. While the difference between the two is due to regulator artifacts, in this work only the second choice is implemented, namely $W_{1,23}^{1\pi ,c_{D}2}$.

For the contact term there are different options on the operator structure, which come from different choices in the Fierz rearrangement. In this work only the following one is considered

$\begin{array}{ll}{W}_{1,23}^{ct,c_E\tau }=E({\tau }_2·{\tau }_3)\delta (r_{12})\delta (r_{13}),& (28)\end{array}$

with $E=c_E\frac{{(\hslash c)}^6}{{\Lambda }_{\chi }{F}_{\pi }^4}$.

The value of all LECs entering the 3N forces at N2LO are shown in Table 1. In references [19, 20] c_D and c_E have been fitted in order to reproduce the ⁴He binding energy and the n-α P-wave phase shift.

TABLE 1

Table 1. Fit values for the couplings c_D and c_E for different choices of 3N cut-offs as reported in references [19, 20].

3.1. Three-Nucleon Forces as Spherical Tensors

The above expressions for the 3N force are not written in terms of irreducible spherical tensors, so that they cannot be implemented directly into the hyperspherical formalism. In this section we address this point and write the interaction in terms of irreducible spherical tensors, both in coordinate-spin space and in isospin space.

For convenience, we denote the general spin space ${\Sigma }_{ij}^{\lambda }$, ${\Sigma }_{ij,k}^{\lambda ,\Lambda }$ and configuration space $X_{ij}^{\lambda }$, $X_{ij,ij}^{(\lambda ,{\lambda }^{\prime })\Lambda }$ irreducible tensor operators as

$\begin{array}{l}{\Sigma }_{ij}^{\lambda }={\left[{\sigma }_i\times {\sigma }_j\right]}^{\lambda },\\ {\Sigma }_{ij,k}^{\lambda ,\Lambda }=[{\left[{\sigma }_k\times {\left[{\sigma }_i\times {\sigma }_j\right]}^{\lambda }\right]}^{\Lambda },\\ X_{ij}^{\lambda }={\left[{\widehat{r}}_{1i}\times {\widehat{r}}_{1j}\right]}^{\lambda },\\ X_{ij,ij}^{(\lambda ,{\lambda }^{\prime })\Lambda }={\left[{\left[{\widehat{r}}_{1i}\times {\widehat{r}}_{1j}\right]}^{\lambda }\times {\left[{\widehat{r}}_{1i}\times {\widehat{r}}_{1j}\right]}^{{\lambda }^{\prime }}\right]}^{\Lambda },& (29)\end{array}$

where i, j, k are generic particle indexes and ${\widehat{r}}_{1i}$ is the rank 1 normalized spherical tensor associated to the relative distance between particle 1 and particle i. With the notation ${\left[{\widehat{r}}_{1i}\times {\widehat{r}}_{1j}\right]}^{\lambda }$ we intend the two rank-one coordinate space tensors coupled into a rank-λ tensor, and analogously for ${\left[{\sigma }_i\times {\sigma }_j\right]}^{\lambda }$ and ${\left[{\tau }_i\times {\tau }_j\right]}^{\lambda }$ in spin and isospin space, respectively. Furthermore, we define

$\begin{array}{ll}{X}^{\lambda }({\widehat{r}}_{ij},{\widehat{r}}_{ij})={\left[{\widehat{r}}_{ij}\times {\widehat{r}}_{ij}\right]}^{\lambda },& (30)\end{array}$

where ${\widehat{r}}_{ij}$ is the rank-1 normalized spherical tensor associated to the relative distance between particles i and j.

At this point, after rearranging the couplings with a few Racah algebra steps and by using the previously introduced notation, one can rewrite the 3N interactions of Equations (23), (27), (28) in terms of irreducible tensors in isospin space and in the coupled spin-configuration space.

The 2π-exchange term depending on c₁ becomes

$\begin{array}{l}{W}_{1,23}^{2\pi ,c_1}=A{U}_{12}{Y}_{12}{U}_{13}{Y}_{13}({\tau }_2·{\tau }_3)({\sigma }_2·{\widehat{r}}_{12})({\sigma }_3·{\widehat{r}}_{13})\\ =-\sqrt{3}A{\left[{\tau }_2\times {\tau }_3\right]}^0{F}_{UU}({\Sigma }_{23}^0·X_{23}^0-{\Sigma }_{23}^1·X_{23}^1\\ +{\Sigma }_{23}^2·X_{23}^2),& (31)\end{array}$

the 2π-exchange term that depends on c₃ becomes

$\begin{array}{l}{W}_{1,23}^{2\pi ,c_3}=B\left\{{\tau }_1·{\tau }_2,{\tau }_1·{\tau }_3\right\}\left\{{\chi }_{12},{\chi }_{13}\right\}\\ =-2\sqrt{3}B{\left[{\tau }_2\times {\tau }_3\right]}^0(\\ {\Sigma }_{23}^0·(+F_{TT}{X}_{23}^0-\frac{1}{\sqrt{3}}(3F_{YY}+F_{TY}+F_{YT}))\\ +{\Sigma }_{23}^1·(-F_{TT}{X}_{23}^1)\\ +{\Sigma }_{23}^2·(+F_{TT}{X}_{23}^2+F_{TY}{X}_{22}^2+F_{YT}{X}_{33}^2)),& (32)\end{array}$

while the term that depends on c₄ can be expressed as

$\begin{array}{l}{W}_{1,23}^{2\pi ,c_4}=-C\left[{\tau }_1·{\tau }_2,{\tau }_1·{\tau }_3\right]\left[{\chi }_{12},{\chi }_{13}\right]\\ =4\sqrt{3}C{\left[{\tau }_1\times {\left[{\tau }_2\times {\tau }_3\right]}^1\right]}^0[\\ {\Sigma }_{23,1}^{1,0}·\left(-F_{TT}{X}_{23,23}^{(1,1)0}-\frac{1}{\sqrt{3}}(3F_{YY}+F_{TY}+F_{YT})\right)\\ +{\Sigma }_{23,1}^{0,1}·\left(+F_{TT}{X}_{23,23}^{(1,0)1}\right)\\ +{\Sigma }_{23,1}^{2,1}·\left(+F_{TT}{X}_{23,23}^{(1,2)1}\right)\\ +{\Sigma }_{23,1}^{1,2}·\left(-F_{TT}{X}_{23,23}^{(1,1)2}-\frac{1}{2}(F_{TY}{X}_{22}^2+F_{YT}{X}_{33}^2)\right)\\ +{\Sigma }_{23,1}^{2,2}·\left(-F_{TT}{X}_{23,23}^{(1,2)2}-\frac{\sqrt{3}}{2}(F_{TY}{X}_{22}^2-F_{YT}{X}_{33}^2)\right)\\ +{\Sigma }_{23,1}^{2,3}·\left(+F_{TT}{X}_{23,23}^{(1,2)3}\right)].& (33)\end{array}$

To write the above expression in a compact form, we have introduced the following definitions

$\begin{array}{l}{F}_{UU}=U_{12}{Y}_{12}{U}_{13}{Y}_{13},\\ F_{TT}=18T_{12}{T}_{13},\\ F_{YY}=2({Ỹ}_{12}-T_{12})({Ỹ}_{13}-T_{13}),\\ F_{TY}=6T_{12}({Ỹ}_{13}-T_{13}),\\ F_{YT}=6({Ỹ}_{12}-T_{12})T_{13}.& (34)\end{array}$

The 1π-exchange contribution takes the following form

$\begin{array}{l}{W}_{1,23}^{1\pi ,c_{D}2}=D({\tau }_2·{\tau }_3){\chi }_{23}({\delta }_{12}+{\delta }_{13})\\ =-\sqrt{3}D{\left[{\tau }_2\times {\tau }_3\right]}^0[\\ {\Sigma }_{23}^0·\left(-\sqrt{3}({\delta }_{12}+{\delta }_{13}){Ỹ}_{23}\right)\\ +{\Sigma }_{23}^2·\left(3({\delta }_{12}+{\delta }_{13})T_{23}{X}^2({\widehat{r}}_{23},{\widehat{r}}_{23})\right)],& (35)\end{array}$

while the contact term becomes

$\begin{array}{ll}{W}_{1,23}^{ct,c_E\tau }=E({\tau }_2·{\tau }_3){\delta }_{12}{\delta }_{13}=-\sqrt{3}E{\left[{\tau }_2\times {\tau }_3\right]}^0{\delta }_{12}{\delta }_{13}.& (36)\end{array}$

We have implemented these expressions in our hyperspherical harmonics codes. Since the interaction is now written in terms of irreducible tensors, the spin and isospin matrix elements can be computed analytically. For the calculation of the spatial matrix elements one can reduce the six-dimensional integration in the two Jacobi coordinates to a two-dimensional numerical quadrature, as explained in details in reference [46]. Below we present the benchmark results we obtained with these local-chiral forces on few-body systems, such as ³H, ³He, and ⁴He.

4. Results

In this section we show the benchmark tests of the maximally-local-chiral interactions using the EIHH method. We compute ground-state energies and charge radii in three- and four-nucleon systems and compare to two Monte Carlo methods, namely the GFMC and AFDMC methods.

In the computations of nuclear charge radii, we use

$\begin{array}{ll}\langle r_{\text{c}}^2\rangle =\langle r_{\text{pt}}^2\rangle +r_{\text{p}}^2+\frac{A-Z}{Z}{r}_{\text{n}}^2+\frac{3{\hslash }^2}{4m_p^2{c}^2},& (37)\end{array}$

where $\sqrt{\langle r_{\mathrm{\text{pt}}}^2\rangle }$ is the calculated point-proton radius, r_p = 0.8751(61) fm [47] is the root-mean-square (rms) charge radius of the proton, $r_{\mathrm{\text{n}}}^2=-0.1161(22)$ fm² [47] is the squared charge radius of the neutron, and Z is the number of protons in the nucleus. The last term is the Darwin-Foldy correction to the proton-charge radius [48] which depends on the proton mass m_p. We neglect the spin-orbit relativistic contribution, since it is negligible in s-shell nuclei [49], as well as meson exchange currents.

Keeping in mind that the goal of this work is to benchmark our expressions for the 3N forces at N2LO by comparing to the Monte Carlo results, we have used the same numerical values for r_p and r_n as in reference [20], which follows the CODATA-2014 recommendations [47]. Hence, in a first stage we will not be using the more modern results for r_p/n from references [50, 51].

A few words addressing the estimation of the numerical uncertainties are in line. As already said, the EIHH method allows for an exact solution of the Schrödinger equation, the computed wavefunction converges to the true eigenfunction of the Hamiltonian operator in the limit of infinite model space. The model space is mostly given by the maximal number, n_max, of Laguerre polynomials and the choice of the maximal value of the grand-angular momentum quantum number, K_max, in the construction of the hyperspherical harmonics functions. It has been practically found that beyond a value n_max = 50, the expectation values are negligibly modified. The convergence in terms of K_max is more delicate, so that in order to estimate the uncertainty coming from the truncation of the model space, we analyze the converging pattern at increasing values of K_max.

To quantify our numerical uncertainty we proceed as follows: denoting with O(K_max) the expectation value of an observable Ô computed by setting a given maximal value of the grand-angular momentum quantum number, K_max, in the wavefunction, our uncertainty in this observable is estimated by

$\begin{array}{l}\delta (O)=|O(K_{\text{max}})-O(K_{\text{max}}-2)|\\ +|O(K_{\text{max}}-2)-O(K_{\text{max}}-4)|+{\delta }_{\text{res}},& (38)\end{array}$

where δ_res is the residual uncertainty (not due to the K_max behavior) obtained by varying: the number of radial grid points (from 70 to 90), the maximal values of the angular momentum in the construction of the two-body effective interaction (from 60 to 120) and the maximal number of three-body angular momentum (from 5/2 to 7/2) in the partial wave expansion of the 3N force.

First, we address and discuss the benchmark of the interactions at LO and NLO, so as to have a clean test on the NN interactions. Then we move to the N2LO, where the three-body forces are included.

4.1. Benchmarks at LO and NLO

We study the maximally-local chiral interactions for two different regulator cut-offs, indicated by r₀, namely exploring the two possibilities of r₀ = 1.0 fm and r₀ = 1.2 fm. The latter gives rise to a softer interaction compared to the first one. For the benchmarks at LO and NLO, the ⁴He nucleus is used as a testing ground. We compute point-proton charge radii, $\sqrt{\langle {\mathrm{\text{r}}}_{\mathrm{\text{pt}}}^2\rangle }$, and ground-state energies, E₀, for the two different cut-off choices at increasing values of the grand-angular momentum quantum number and compare to the GFMC calculations.

The final results are shown in Table 2, where the uncertainty is computed as explained above using Equation (38) with K_max = 22. An extended table with all the various K_max can be found in the Supplementary Material. We observe that as we enlarge the model space a nice converging pattern is obtained and our final EIHH results basically agree with the GFMC calculations when the uncertainties are taken into account. By looking at the converging pattern of the studied observables as the model space is increased (see Supplementary Material), we clearly observe that the interaction with r₀ = 1.2 fm is much softer than the other, since the relative observables converge with a smaller model space. Finally, it is notable that, as shown in Table 2, the LO and NLO results do not reproduce the measured values, but the discrepancy decreases in going from LO to NLO.

TABLE 2

Table 2. Ground-state energies and point-proton radii for the ⁴He nuclear system at LO and NLO computed with the EIHH method.

4.2. Benchmarks at N2LO

We now turn to the benchmark at the next order. At N2LO we have the first appearance of 3N forces, so this will serve as a check of our irreducible tensor representation. The 3N interaction involves two new LECs, c_D and c_E, coming from the 1π-term and from the ct-term of the 3N forces, respectively, that cannot be fitted in the NN sector. In reference [19] these couplings have been fitted to reproduce the ⁴He binding energy and the n-α scattering P-wave phase shift, for which the values reported in Table 1 were obtained. We use the same values in this work, as our goal is to perform a benchmark. In particular, here we implement only the (D2, Eτ) 3N interactions, which we chose since the Eτ term has a more general isospin structure. Different choices of the 3N contact term have been shown to lead to different properties in neutron matter [18, 19].

As a testing ground for our N2LO Hamiltonian expressed in terms of spherical tensors outlined in the previous section, we study the three-body ³He and ³H and the four-body ⁴He nuclear systems. We compute ground-state energies, E₀, and charge radii, $\sqrt{\langle {\mathrm{\text{r}}}_{\mathrm{\text{c}}}^2\rangle }$, for the two different cut-off choices r₀ = 1 and 1.2 fm and carefully study the convergence at increasing K_max values. A complete table of our data is shown in the Supplementary Material. The K_max convergence is also explicitly shown in a graphical manner in Figures 2–4, where a comparison to the GFMC method is made.

FIGURE 2

Figure 2. The ground-state energy and the charge radius of the nuclear ³He system as a function of the grand-angular momentum quantum number K_max. The green and blue error-bands are the GFMC results with the relative statistical uncertainty.

FIGURE 3

Figure 3. The ground-state energy and the charge radius of the nuclear ³H system as a function of the grand-angular momentum quantum number K_max. The green and blue error-bands are the GFMC results with the relative statistical uncertainty.

FIGURE 4

Figure 4. The ground state energy and the charge radius of the nuclear ⁴He system as a function of the grand-angular momentum quantum number K_max. The green and blue error-bands are the GFMC results with the relative statistical uncertainty.

As it can be seen from Figures 2, 3, the EIHH method is in excellent agreement with the GFMC computations for the three-body nuclei, for both the ground-state energies and the charge radii. The typical non-monotonic convergence patter of the EIHH method is observed, and a very good convergence is reached already at K_max = 12. This shows that these forces are softer than the AV18 potential, but harder than the low-k interactions [54].

For the ⁴He nucleus shown in Figure 4, we obtain a very nice agreement with the GFMC method for the cut-off value r₀ = 1.2 fm, while for the cut-off r₀ = 1.0 fm, we perfectly reproduce the charge radius, but we observe a small deviation for the ground-state energy with respect to the GFMC.

Our final EIHH results with uncertainties quantified as explained above using Equation (38) with K_max = 22 are shown in Table 3 in comparison with the GFMC, AFDMC, and the experimental data. The theoretical computations agree with each other at a similar level as observed in other benchmarks [59]. In particular, for the cut-off r₀ = 1.0 fm, which leads to a harder force, EIHH and GFMC are in agreement within error bars, while some discrepancy is seen with respect to the AFDMC results. Given that the latter method is known to be less accurate, we do not think that this difference is significant and we consider all these results to constitute a successful benchmark of our implementation of 3N forces.

TABLE 3

Table 3. Ground-state energies and charge radii for the nuclear ³He, ³H, and ⁴He systems at N2LO in the chiral expansion computed with the EIHH, GFMC, and AFDMC methods.

As can be seen in Table 3, at N2LO a much improved agreement with experiment is obtained. In fact, if one compares the experimental binding energies to the LO and NLO calculations in Table 2 one observes that these low orders overbind (LO) or underbind (NLO) the few-body nuclei, while at N2LO nice agreement is observed. This is expected for ⁴He, given that 3N forces are fit to reproduce the ⁴He binding energy, however a better agreement is also found for ³He and ³H due to the strong correlation between the three- and four-body binding energy. Interestingly, a nice converging pattern is also found for the nuclear charge radii.

From a careful look at Table 3, one can appreciate that our EIHH calculations are more precise than the GFMC and AFDMC results in the three-nucleon sector and that our numerical uncertainty is comparable to the experimental uncertainties for the radii. While this may be an advantage of our method, it is important to note that the error bars quoted in this table do not include the uncertainties coming from the ChEFT expansion, so they do not constitute the full uncertainty of the theory.

We conclude this section with a further investigation on the charge radii of light-nuclear systems. In reference [55] the proton-charge radius r_p = 0.8751(61) fm and the neutron-charge radius $r_{\mathrm{\text{n}}}^2=-0.1161(22)$ fm² recommended by CODATA-2014 were used in the evaluation of nuclear charge radii using Equation (37). Such single-nucleon data come from experiments that study the electron-nucleon system. Recently, these quantities were measured more precisely by investigating muonic atoms, and one could ask what is the effect of this increased precision in the nuclear charge radius when applying Equation (37). To address this point in Table 4 we compare our results for the charge radii of ³He, ³H, and ⁴He at N2LO using the CODATA-2014 single-nucleon input with the results obtained using the rms proton-charge radius coming from the muonic-hydrogen r_p = 0.84087(39) [50] and the new value of the rms charge radius of the neutron $r_{\mathrm{\text{n}}}^2=-0.106(7)$ fm² [51]. We denote the first choice with e − r_c and the second with μ − r_c. The general effect of using this choice of the proton and neutron charge radii amounts to a systematic reduction of roughly 1% of the charge radii of these light nuclei. This has to be contrasted with the full uncertainty of the theory that includes not only the EIHH numerical error, but also considers the uncertainty coming from the order-by-order chiral expansion. The latter is estimated using the algorithm proposed first in reference [60] and is included in Table 4.

TABLE 4

Table 4. Nuclear rms charge radii for ³He, ³H, and ⁴He systems at N2LO computed using either the single-nucleon CODATA-2014 values (columns e-r_c) or the more precise muonic atoms data (columns μ-r_c).

For a graphical representation of our findings, in Figure 5 we show the ⁴He nuclear charge radius at increasing chiral orders computed for different choices for the proton and neutron charge radii. We observe that the chiral order uncertainty is of the order of 2%, hence larger than the effect of the more precise single-nucleon input. Overall, we confirm the chiral order-by-order convergence patter, already discussed in references [19, 20], but there shown only for the binding energy and the point-proton radius, which does not include the single nucleon input.

FIGURE 5

Figure 5. The ⁴He charge radius computed at increasing orders of the chiral expansion. The uncertainty bars include the numerical uncertainty of the EIHH method as well as the uncertainties coming from the truncation of the chiral expansion. The horizontal lines are the experimental values from electron scattering (solid line) [56] and from muonic atoms (dashed line) [61].

Interestingly, when comparing the ⁴He theoretical charge radius with the newest muonic atom measurement from reference [61], we see that the μ − r_c results are still consistent with the experiment, leaving, however, space for meson exchange currents to help improve the theoretical precision, which is by far lower than the experimental one.

5. Conclusion and Outlook

In this work, the maximally local chiral interactions are implemented for the first time in the hyperspherical harmonic formalism. The benchmark tests performed in light nuclei show general agreement between hyperspherical harmonic results and the previously available Monte Carlo calculations. As expected, at N2LO with the inclusion of the 3N forces the experimental results are much better reproduced with respect to the LO and NLO calculations. With this study we thus confirm the nice order-by-order convergence in the ground-state energies and in the radii that was already observed in the Monte Carlo studies.

While our numerical precision of the EIHH calculations lies in the sub-percent range, we find that the uncertainty due to the chiral order expansion is larger. In case of the charge radius, we observed that using the most updated values of the proton and nucleon radii instead of the CODATA-2014 values leads to a variation of 1%, which is smaller than the 2% uncertainty found in the chiral order-by-order truncation at N2LO. Addressing first the latter by going to N3LO should be the priority if the goal is to reduce theoretical uncertainties.

Having these new interactions implemented in our formalisms opens up the possibility of investigating other few-body observables in the future. Our most immediate goals include the investigation of muonic atoms [37], of the ⁴He monopole transition form factor [62], and the longitudinal response function [63] in order-by-order chiral expansions. We reserve these applications for future studies.

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author Contributions

NB derived the first expressions for the spherical tensors, which were then checked by SSL and SB. The two-body force was implemented by SB while SSL implemented the expressions of the 3N spherical tensors in the hyperspherical harmonics code and run the calculations, and results were discussed in the group at every step. All authors contributed in equal parts to this paper and contributed to the writing of the text.

Funding

This work was supported by the Deutsche Forschungsgemeinschaft (DFG) through the Collaborative Research Center [The Low-Energy Frontier of the Standard Model (SFB 1044)], and through the Cluster of Excellence Precision Physics, Fundamental Interactions, and Structure of Matter (PRISMA⁺ EXC 2118/1) funded by the DFG within the German Excellence Strategy (Project ID 39083149). Calculations were performed on the mogon2 cluster in Mainz.

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

SSL and SB would like to acknowledge Joel Lynn, Ingo Tews, and Diego Lonardoni for useful discussions.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphy.2021.671869/full#supplementary-material

Footnotes

1. ^This property will be used later when we will write explicitly the form of the nuclear forces between (among) two (three) particles.

2. ^Even though they are operators in spin space, we do not use the hat in our notation, as they are vectors, whose components are operators.

3. ^We drop the hat from vectors whose components are operators.

References

1. Yukawa H. On the interaction of elementary particles. I. Prog Theor Phys Suppl. (1955) 1:1–10. doi: 10.1143/PTPS.1.1

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Stoks VGJ, Klomp RAM, Terheggen CPF, de Swart JJ. Construction of high-quality NN potential models. Phys Rev C. (1994) 49:2950–62. doi: 10.1103/PhysRevC.49.2950

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

5. Pudliner BS, Pandharipande VR, Carlson J, Wiringa RB. Quantum Monte Carlo Calculations of A ≤ 6 Nuclei. Phys Rev Lett. (1995) 74:4396–9. doi: 10.1103/PhysRevLett.74.4396

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Pieper SC, Pandharipande VR, Wiringa RB, Carlson J. Realistic models of pion-exchange three-nucleon interactions. Phys Rev C. (2001) 64:014001. doi: 10.1103/PhysRevC.64.014001

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

9. Rocco N. Ab initio calculations of lepton-nucleus scattering. Front Phys. (2020) 8:116. doi: 10.3389/fphy.2020.00116

CrossRef Full Text | Google Scholar

10. Weinberg S. Phenomenological Lagrangians. Physica. (1979) 96A:327–40. doi: 10.1016/0378-4371(79)90223-1

CrossRef Full Text | Google Scholar

11. Weinberg S. Nuclear forces from chiral Lagrangians. Phys Lett B. (1990) 251:288–92. doi: 10.1016/0370-2693(90)90938-3

CrossRef Full Text | Google Scholar

12. Weinberg S. Effective chiral Lagrangians for nucleon-pion interactions and nuclear forces. Nucl Phys B. (1991) 363:3–18. doi: 10.1016/0550-3213(91)90231-L

CrossRef Full Text | Google Scholar

13. Weinberg S. Three-body interactions among nucleons and pions. Phys Lett B. (1992) 295:114–21. doi: 10.1016/0370-2693(92)90099-P

CrossRef Full Text | Google Scholar

14. Epelbaum E, Hammer HW, Meißner UG. Modern theory of nuclear forces. Rev Mod Phys. (2009) 81:1773–825. doi: 10.1103/RevModPhys.81.1773

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

16. Epelbaum E, Krebs H, Reinert P. High-precision nuclear forces from chiral EFT: state-of-the-art, challenges, and outlook. Front Phys. (2020) 8:98. doi: 10.3389/fphy.2020.00098

CrossRef Full Text | Google Scholar

17. Gezerlis A, Tews I, Epelbaum E, Freunek M, Gandolfi S, Hebeler K, et al. Local chiral effective field theory interactions and quantum Monte Carlo applications. Phys Rev C. (2014) 90:054323. doi: 10.1103/PhysRevC.90.054323

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Tews I, Gandolfi S, Gezerlis A, Schwenk A. Quantum Monte Carlo calculations of neutron matter with chiral three-body forces. Phys Rev C. (2016) 93:024305. doi: 10.1103/PhysRevC.93.024305

CrossRef Full Text | Google Scholar

19. Lynn J, Tews I, Carlson J, Gandolfi S, Gezerlis A, Schmidt KE, et al. Chiral three-nucleon interactions in light nuclei, neutron-α scattering and neutron matter. Phys Rev Lett. (2016) 116:062501. doi: 10.1103/PhysRevLett.116.062501

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Lynn J, Tews I, Carlson J, Gandolfi S, Gezerlis A, Schmidt KE, et al. Quantum Monte Carlo calculations of light nuclei with local chiral two- and three-nucleon interactions and neutron matter. Phys Rev C. (2017) 96:054007. doi: 10.1103/PhysRevC.96.054007

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

22. Piarulli M, Girlanda L, Schiavilla R, Kievsky A, Lovato A, Marcucci LE, et al. Local chiral potentials with Δ-intermediate states and the structure of light nuclei. Phys Rev C. (2016) 94:054007. doi: 10.1103/PhysRevC.94.054007

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Zernike F, Brinkman HC. Hyperspharische Funktionen und die in sphärischen Bereichen orthogonalen Polynome. Proc K Ned Akad Wett. (1935) 38:161.

25. Delves LM. Tertiary and general-order collisions. Nucl Phys. (1959) 9:391–399. doi: 10.1016/0029-5582(58)90372-9

CrossRef Full Text | Google Scholar

26. Simonov YA. The three-body problem. Complete system of angular functions. Sov J Nucl Phys. (1966) 3:461.

Google Scholar

27. Zickendraht W. Construction of a complete orthogonal system for the quantum-mechanical three-body problem. Ann Phys. (1965) 35:18. doi: 10.1016/0003-4916(65)90067-9

CrossRef Full Text | Google Scholar

28. Smith FT. Generalized angular momentum in many-body collisions. Phys Rev. (1960) 120:1058. doi: 10.1103/PhysRev.120.1058

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

30. Marcucci LE, Dohet-Eraly J, Girlanda L, Gnech A, Kievsky A, Viviani M. The hyperspherical harmonics method: a tool for testing and improving nuclear interaction models. Front Phys. (2020) 8:69. doi: 10.3389/fphy.2020.00069

CrossRef Full Text | Google Scholar

31. Barnea N. Exact Solution of the Schrödinger and Faddeev Equations for Few-Body Systems. Jerusalem: Hebrew University (1997).

Google Scholar

32. Bacca S. Study of Electromagnetic Reactions on Light Nuclei With the Lorentz Integral Transform Method. Trento; Mainz: Universita' degli Studi di Trento; Johannes Gutenberg Universität Mainz (2005).

Google Scholar

33. Efros VD. To K harmonics method in problem of several nucleons. Sov J Nucl Phys. (1972) 15:128.

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

37. Ji C, Bacca S, Barnea N, Hernandez OJ, Nevo-Dinur N. Ab initio calculation of nuclear structure corrections in muonic atoms. J Phys G. (2018) 45:093002. doi: 10.1088/1361-6471/aad3eb

CrossRef Full Text | Google Scholar

38. Kaplan DB, Savage MJ, Wise MB. A new expansion for nucleon-nucleon interactions. Phys Lett B. (1998) 424:390–6. doi: 10.1016/S0370-2693(98)00210-X

CrossRef Full Text | Google Scholar

39. Kaplan DB, Savage MJ, Wise MB. Two-nucleon systems from effective field theory. Nucl Phys B. (1998) 534:329–55. doi: 10.1016/S0550-3213(98)00440-4

CrossRef Full Text | Google Scholar

40. Nogga A, Timmermans RGE, Kolck U. Renormalization of one-pion exchange and power counting. Phys Rev C. (2005) 72:054006. doi: 10.1103/PhysRevC.72.054006

CrossRef Full Text | Google Scholar

41. Pavón Valderrama M, Arriola ER. Renormalization of the NN interaction with a chiral two-pion-exchange potential: central phases and the deuteron. Phys Rev C. (2006) 74:054001. doi: 10.1103/PhysRevC.74.054001

CrossRef Full Text | Google Scholar

42. Long B, Yang CJ. Renormalizing chiral nuclear forces: triplet channels. Phys Rev C. (2012) 85:034002. doi: 10.1103/PhysRevC.85.034002

CrossRef Full Text | Google Scholar

43. van Kolck U. Few-nucleon forces from chiral Lagrangians. Phys Rev C. (1994) 49:2932–41. doi: 10.1103/PhysRevC.49.2932

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Piarulli M, Tews I. Local nucleon-nucleon and three-nucleon interactions within chiral effective field theory, and neutron matter. Front Phys. (2020) 7:245. doi: 10.3389/fphy.2019.00245

CrossRef Full Text | Google Scholar

45. Wiringa RB, Pieper SC. Evolution of nuclear spectra with nuclear forces. Phys Rev Lett. (2002) 89:182501. doi: 10.1103/PhysRevLett.89.182501

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Barnea N, Efros VD, Leidemann W, Orlandini G. Incorporation of three-nucleon force in the effective interaction hyperspherical harmonics approach. Few Body Syst. (2004) 35:155–67. doi: 10.1007/s00601-004-0066-y

CrossRef Full Text | Google Scholar

47. Mohr PJ, Newell DB, Taylor BN. CODATA recommended values of the fundamental physical constants: 2014. Rev Mod Phys. (2016) 88:035009. doi: 10.1103/RevModPhys.88.035009

CrossRef Full Text | Google Scholar

48. Friar JL, Martorell J, Sprung DWL. Nuclear sizes and the isotope shift. Phys Rev A. (1997) 56:4579–86. doi: 10.1103/PhysRevA.56.4579

CrossRef Full Text | Google Scholar

49. Ong A, Berengut JC, Flambaum VV. Effect of spin-orbit nuclear charge density corrections due to the anomalous magnetic moment on halonuclei. Phys Rev C. (2010) 82:014320. doi: 10.1103/PhysRevC.82.014320

CrossRef Full Text | Google Scholar

50. Antognini A, Nez F, Schuhmann K, Amaro FD, Biraben F, Cardoso JMR, et al. Proton structure from the measurement of 2S-2P transition frequencies of muonic hydrogen. Science. (2013) 339:417–20. doi: 10.1126/science.1230016

PubMed Abstract | CrossRef Full Text | Google Scholar

51. Filin AA, Baru V, Epelbaum E, Krebs H, Möller D, Reinert P. Extraction of the neutron charge radius from a precision calculation of the deuteron structure radius. Phys Rev Lett. (2020) 124:082501. doi: 10.1103/PhysRevLett.124.082501

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Lynn JE, Carlson J, Epelbaum E, Gandolfi S, Gezerlis A, Schwenk A. Quantum Monte Carlo calculations of light nuclei using chiral potentials. Phys Rev Lett. (2014) 113:192501. doi: 10.1103/PhysRevLett.113.192501

PubMed Abstract | CrossRef Full Text | Google Scholar

53. NuDat. Available online at: https://www.nndc.bnl.gov/nudat2 (accessed January 22, 2021).

54. Bacca S, Schwenk A, Hagen G, Papenbrock T. Helium halo nuclei from low-momentum interactions. Eur Phys J A. (2009) 42:553. doi: 10.1140/epja/i2009-10815-5

CrossRef Full Text | Google Scholar

55. Lonardoni D, Gandolfi S, Lynn JE, Petrie C, Carlson J, Schmidt KE, et al. Auxiliary field diffusion Monte Carlo calculations of light and medium-mass nuclei with local chiral interactions. Phys Rev C. (2018) 97:044318. doi: 10.1103/PhysRevC.97.044318

CrossRef Full Text | Google Scholar

56. Sick I. Zemach moments of ³He and ⁴He. Phys Rev C. (2014) 90:064002. doi: 10.1103/PhysRevC.90.064002

CrossRef Full Text | Google Scholar

57. Purcell JE, Kelley JH, Kwan E, Sheu CG, Weller HR. Energy levels of light nuclei A=3. Nucl Phys A. (2010) 848:1–74. doi: 10.1016/j.nuclphysa.2010.08.012

CrossRef Full Text | Google Scholar

58. Angeli I, Marinova KP. Table of experimental nuclear ground state charge radii: an update. Atomic Data Nucl Data Tables. (2013) 99:69–95. doi: 10.1016/j.adt.2011.12.006

CrossRef Full Text | Google Scholar

59. Nevo Dinur N, Hernandez OJ, Bacca S, Barnea N, Ji C, Pastore S, et al. Zemach moments and radii of ², 3H and ³, 4He. Phys Rev C. (2019) 99:034004. doi: 10.1103/PhysRevC.99.034004

CrossRef Full Text | Google Scholar

60. Epelbaum E, Krebs H, Meißner UG. Improved chiral nucleon-nucleon potential up to next-to-next-to-next-to-leading order. Eur Phys J A. (2015) 51:53. doi: 10.1140/epja/i2015-15053-8

CrossRef Full Text | Google Scholar

61. Krauth J, Schuhmann K, Ahmed MA, Amaro FD, Amaro P, Biraben F, et al. Measuring the α-particle charge radius with muonic helium-4 ions. Nature. (2022) 589:527–31. doi: 10.1038/s41586-021-03183-1

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

63. Bacca S, Barnea N, Leidemann W, Orlandini G. Role of the final-state interaction and three-body force on the longitudinal response function of ⁴He. Phys Rev Lett. (2009) 102:162501. doi: 10.1103/PhysRevLett.102.162501

PubMed Abstract | CrossRef Full Text | Google Scholar