Implementation of local chiral interactions in the hyperspherical harmonics formalism

With the goal of using chiral interactions at various orders to explore 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 $^3$H, $^3$He and $^4$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.


INTRODUCTION
In 1935 the seminal idea of Yukawa [1] laid the foundation to the theory of the nuclear forces. His onepion 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 Refs. [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., Refs. [7,8,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.
The hyper-radial coordinates ρ 1 , ..., ρ N and the hyper-angular coordinates ϕ 2 , ..., ϕ N are constructed recursively. The transformation law for the first two Jacobi coordinates is Assuming that we already know the hyper-radial coordinates ρ 1 , ..., ρ j−1 and the hyper-angular coordinates ϕ 2 , ..., ϕ j−1 the transformation law for ρ j and ϕ j reads in analogy to Eq. (2) as 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 where the hyper-radial part is whileK 2 N is the grand-angular momentum operator whose eigenfunctions are known as the hyperspherical harmonics.
Denotingl j as the angular momentum operator related to η j , andL 2 j andM 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 operatorK 2 N of the system recursively in terms ofK 2 N −1 andl N as [29] whereK 2 1 =l 2 1 . The operatorsK 2 N , ...,K 2 2 ,L 2 N , ...,L 2 2 ,l 2 N , ...,l 2 1 andM N commute with each others. As a consequence, it is possible to label hyperspherical states using the set of 3N − 1 quantum numbers {K} ≡ {K N , ..., K 2 , L N , ..., L 2 , l N , ..., l 1 , M N }. The hyperspherical harmonics functions Y {K N } 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 [30] where C LM l i m i ,l j m j are the Clebsch-Gordan coefficients, Y l j m j (Ω j ) are the spherical harmonics associated with η j and 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, 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 1 ⊂ ... ⊂ S A presented by the appropriate Young diagrams Γ 1 , ..., Γ 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 [K N ] are constructed recursively. Assuming that Y [K N −1 ] have been already constructed, the N th 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 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 This is a provisional file, not the final typeset article Nucleons posses also 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 Here Analogously to what has been done with the hyperspherical harmonics, the spin-isospin wavefunctions are constructed recursively. Assuming that the symmetry-adapted wavefunction χ [S j−1 ] have been obtained, the construction of the χ [S j ] is done by first coupling χ [S j−1 ] to the spin-isospin wavefunction of the jth nucleon, let us call this state χ [S j−1 ],S j T j , and then taking linear combinations of χ [S j−1 ],S j T j using the coefficients of fractional parentage labeled as 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 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 v n (ρ N ). Again, the model space is truncated to a given maximum number of Laguerre polynomials n max 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 [31], 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 Ref. [32], and also in the more recent review [33].

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 whereT is the sum of the non-relativistic kinetic energy operators of the individual nucleons,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 and when working with antisymmetrized wavefunctions, the expectation values of the NN and 3N terms become where only the first two (or three) particles are involved 1 .
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,11,12,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/Λ χ ) ν , 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 [34,35,36,37,38,39], 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,18,19].
The local chiral NN forces are composed by contact (ct) terms and pion-exchange (π) terms so that the interaction between particle 1 and 2 can be written as 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 where r 12 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 2 and δ is the delta function.
At NLO, the following new terms enter 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 3 , 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 5 is the only non-local operator appearing in this maximally local chiral interaction.
Following Refs. [17,40], all the pion-exchange interactions up to N2LO can be written in a complete local form as where, S 12 is the well known tensor operator, defined as wherer 12 is the unitary vector related to the relative distance r 12 . The local functions V C (r 12 ), W C (r 12 ), V S (r 12 ), W S (r 12 ), V T (r 12 ) and W T (r 12 ) 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 Ref. [40]. In Ref. [19] pion loops are regularized using the spectral-function regularization (SFR) with an ultraviolet cut-offΛ = 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' [41]. 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 Fig. 1. On the one hand, the 2π-term comes with the LECs c 1 , c 3 and c 4 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.
With respect to Ref. [19], 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 Eq. ( 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 where the coupling constants are A = c 1 The 2π-terms include the following functions with analogous expressions for Y 13 and U 13 . The operator χ 12 (and analogously χ 13 ) is defined as with In the last expression, r 0 is the cut-off and following Refs. [19,18] n is taken to be equal to 4.
For the 1π-interaction terms there are two options and W 1π,c D 2 While the difference between the two is due to regulator artifacts, in this work only the second choice is implemented, namely W 1π,c D 2 1,23 .

Frontiers
For the contact term there are different options on the operator structure, which come from different choices in the Fierz rearrangement. In this work the following two are considered and  Table 1. Fit values for the couplings c D and c E for different choices of 3N cut-offs as reported in [18,19]. The constants c 1,3,4 are tuned in the pion-nucleon sector, see Ref. [15].
The value of all LECs entering the 3N forces at N2LO are shown in Table 1. In Refs. [18,19] c D and c E have been fitted in order to reproduce the 4 He binding energy and the n-α P -wave phase shift.

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 can not 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 Σ λ ij , Σ λ,Λ ij,k and configuration space X λ ij , X (λ,λ )Λ ij,ij irreducible tensor operators as where i, j, k are generic particle indexes andr 1i is the rank 1 normalized spherical tensor associated to the relative distance between particle 1 and particle i. With the notation [r 1i ×r 1j ] λ we intend the two rank-one coordinate space tensors coupled into a rank-λ tensor, and analogously for [σ i × σ j ] λ and [τ i × τ j ] λ in spin and isospin space, respectively. Furthermore, we define wherer 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 Eqs. (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 1 becomes the 2π-exchange term that depends on c 3 becomes while the term that depends on c 4 can be expressed as To write the above expression in a compact form, we have introduced the following definitions The 1π-exchange contribution takes the following form while the contact terms become and W E1
We have implemented these expressions in our hypershperical 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 Ref. [42]. Below we present the benchmark results we obtained with these local-chiral forces on few-body systems such as 3 H, 3 He and 4 He.

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 were r 2 pt is the calculated point-proton radius, r p = 0.8751(61) fm [43] is the root-mean-square (rms) charge radius of the proton, r 2 n = −0.1161(22) fm 2 [43] 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 [44] which depends on the proton mass m p . We neglect the spin-orbit relativistic contribution, since it is negligible in s-shell nuclei [45], 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 value for r p and r n as in Ref. [19], which follows the CODATA-2014 recommendations [43]. Hence, in a first stage we will not be using the more modern results for r p/n from Refs. [46,47].
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 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 to have a clean test on the NN interactions. Then we move to the N2LO, where the three-body forces are included.

Benchmarks at LO and NLO
We study the maximally-local chiral interactions for two different regulator cut-offs, indicated by r 0 , namely exploring the two possibilities of r 0 = 1.0 fm and r 0 = 1.2 fm. The latter gives rise to a softer interaction compared to the first one. For the benchmarks at LO and NLO, the 4 He nucleus is used as a testing ground. We compute point-proton charge radii, r 2 pt , and ground-state energies, E 0 , for the two different cut-off choices at increasing values of the grand-angular momentum quantum number and compare to the GFMC calculations.  Table 2. Ground-state energies and point-proton radii for the 4 He nuclear system at LO and NLO computed with the EIHH method. For comparison we report the GFMC results and the experimental values taken from Ref. [48,49].
The final results are shown in Table 2, where the uncertainty is computed as explained above using Eq. (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 agree with the GFMC calculations within uncertainties. 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 0 = 1.2 fm is much softer than the other, since the relative observables converge with a smaller model space. Finally, it is to note 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.

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 can not be fitted in the NN sector. In Ref. [18] these couplings have been fitted to reproduce the 4 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 saturation properties in neutron matter [18].
As a testing ground for our N2LO Hamiltonian expressed in terms of spherical tensors outlined in the previous section, we study the three-body 3 He and 3 H and the four-body 4 He nuclear systems. We compute ground-state energies, E 0 , and charge radii, r 2 c , for the two different cut-off choices r 0 = 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 Fig. 2, Fig. 3, and Fig. 4, where a comparison to the GFMC method is made. As it can be seen from Fig. 2 and 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-q interactions [50].    For the 4 He nucleus shown in Figure 4, we obtain a very nice agreement with the GFMC method for the cut-off value r 0 = 1.2 fm, while for the cut-off r 0 = 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 Eq. (38) with K max = 22 are shown in Table 3 in comparisons with the GFMC, AFDMC and the experimental data. We note that the small difference found for 4 He ground-state energy is just at the level of 0.03 MeV in the central values and is non-significant when the full uncertainty of the EIHH method is considered. Similar kind of sub-percentage differences between EIHH and GFMC were also observed in other benchmarks [54]   and can be found at this level of precision. It is to note that the cut-off r 0 = 1.0 fm leads to a harder force, where in fact, quite a large discrepancy is seen also between the GFMC and the AFDMC computations. Therefore, 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.
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 4 He, given that 3N forces are fit to reproduce the 4 He binding energy, however a better agreement is also found for 3 He and 3 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 Ref. [55] the proton-charge radius r p = 0.8751(61) fm and the neutron-charge radius r 2 n = −0.1161(22) fm 2 recommended by CODATA-2014 were used in the evaluation of nuclear charge radii using Eq. (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 Eq. (37). To address this point in Table 4 we compare our results for the charge radii of 3 He, 3 H and 4 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) [46] and the new value of the rms charge radius of the neutron r 2 n = −0.106(7) fm 2 [47]. We denote the first choice with e − r c and the second with µ − r c . The general effect of using this choices 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 Ref. [56] and is included in Table 4.
For a graphical representation of our findings, in Fig 5 we show the 4 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 oder-by-order convergence patter, already discussed in Refs. [18,19], but there shown only for the binding energy and the point-proton radius, which does not include the single nucleon input.
Interestingly, when comparing the 4 He theoretical charge radius with the newest muonic atom measurement from Ref. [57], we see that the µ − r c results are still consistent with experiment, leaving however space for meson exchange currents to help improving the theoretical precision, which is by far lower than the experimental one.  Table 4. Nuclear rms charge radii for 3 He, 3 H and 4 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 ). The theoretical results are compared to data from the electron-nucleus system [51,53] and, when available, to data obtained from the muon-nucleus system [57].

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 uncertainties due to the chiral order expansion is higher. 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% uncertainties 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.  [51] and and from muonic atoms (dashed line) [57].
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 [33] and of the 4 He monopole transition form factor [58] in an order-by-order chiral expansions. We reserve these applications to future studies.