Uncertainties in ab initio nuclear structure calculations with chiral interactions

1 Introduction

An atomic nucleus, consisting of Z protons and N neutrons, is a self-bound quantum many-body system with A = N + Z strongly interacting nucleons. The interactions between these nucleons are in principle governed by QCD–but it is impractical to describe nuclei in terms of quarks and gluons, except for the very lightest systems. Even a microscopic description of nuclei using realistic two-body (NN), three-body (3N) and possibly higher n-body interactions between point-like nucleons remains a formidable task, both in terms of high-performance computing, and in terms of determining realistic nuclear interactions in tractable terms. In order to confront such a description with experimental data, one needs honest assessments of all uncertainties, both those arising from the numerical solution of a many-body problem, and those arising from a necessarily approximate theory of the effective interactions between nucleons.

Any ab initio theory of nuclei in terms of interacting nucleons requires a high-quality NN potential providing an accurate description of NN scattering data. Highly accurate NN potentials have been in existence for several decades now, all incorporating one-pion exchange, and often inspired by one-boson-exchange (OBE) models, adjusted and augmented by phenomenological terms as necessary to fit the available NN data, such as the Argonne [1], (CD)-Bonn [2, 3], and Nijmegen [4] potentials. Although there exist highly accurate NN potentials in terms of describing the NN phase-shift data, that at the same time accurately describe the spectra of light nuclei [5], most realistic NN potentials require more or less phenomenological 3N forces (3NFs) [6–8] for a good description of nuclei in the p-shell and beyond. However, in order to quantify any uncertainties associated with the choice of the NN potential (and 3NFs), we need a more systematic method of arriving at the potential.

Chiral Effective Field Theory (χEFT) allows us to derive nuclear interactions in a systematic way, in terms of an expansion in the pion mass (or the relevant nucleon momentum) over the hadronic or breakdown scale [9–12]; and in principle, it also allows for a quantification of the uncertainties arising from truncating this expansion. However, this chiral expansion is by no means unique, and different choices for e.g., the degrees of freedom to include in the χEFT can lead to very different χEFT interactions, with a different ordering of various types of diagrams, and indeed different orders at which higher n-body forces have to be included. Furthermore, different choices on e.g., how to regulate the various expressions for loop integrals lead to different versions of the NN potentials (and 3NFs) at any given order, even if the ordering of the various diagrams is the same. Each of these different versions of χEFT comes with its own parameters (Low-Energy Constants or LECs) that need to be fitted to data (or, eventually, calculated from e.g., lattice QCD), and with its own uncertainty quantification.

Just like there are different ways to obtain (effective) nuclear interactions, there are different quantum many-body methods being used for ab initio nuclear structure calculations. For up to four nucleons, one can use the Faddeev–Yakubovsky method (the 3- and 4-body reformulation of the Schrödinger equation that permits the incorporation of the appropriate boundary condition for 3- and 4-body systems that are asymptotically clustered), but this has not been extended to A = 6 or beyond. Broadly speaking, the computational methods applicable to nuclei beyond ⁴He, fall into one of three categories: Quantum Monte Carlo simulations (both variational, VMC [13], and Green’s function, GFMC [13]), non-relativistic lattice simulations with nucleons (Nuclear Lattice Effective Field Theory, NLEFT [14, 15]), and Configuration Interaction (CI) methods (No-Core Shell Model (NCSM) [16], Coupled-Cluster (CC) [17], In-Medium Similarity Renormalization Group (IM-SRG) [18]), which are based on an expansion of the many-body wave-functions in terms of basis functions (configurations). Each of these methods has their own uncertainties: Monte Carlo simulations are typically dominated by statistical uncertainties, though there is also a dependence on the variational wave function; lattice simulations have both statistical and systematic (lattice size and lattice spacing) uncertainties; and CI methods are generally dominated by systematic uncertainties due to the truncation of the many-body basis, though one can make use of statistical sampling of the many-body basis [19]. Because each of these methods have different sources of uncertainties, and they are not always easy to identify and quantify, it is very valuable to use two or more of these many-body methods for the same nucleus, using the same interactions.

In this paper we use the NCSM to perform ab initio nuclear structure calculations for the ground state energies of nearly all stable p-shell nuclei (excluding mirror nuclei) from A = 4 to A = 16 using the χEFT interactions from Ref. [20]. We perform a systematic set of order-by-order calculations in the chiral expansion to determine the uncertainties associated with the truncation of the chiral expansion; more details about the χEFT and how we estimate the truncation uncertainty can be found in Section 2. In order to assess the numerical uncertainties in our NCSM calculations, we make a detailed comparison with Faddeev–Yakubovsky calculations for ³H and ⁴He using the same interactions; this is described in Section 3, together with details about the NCSM. Our results for the binding energies of p-shell nuclei are presented in Section 4. Finally, we give some concluding remarks in Section 5.

2 Nuclear interactions from chiral effective field theory

In recent years two different formulations of χEFT have emerged that are being used in ab initio nuclear structure calculations. The most commonly used χEFT is based on only pions and nucleon degrees of freedom [10, 11], for which the Leading Order (LO) and Next-to-Leading Order (NLO) terms consists of just two-body interactions; three-body interactions first appear at Next-to-Next-to-Leading order (N²LO). Alternatively, one can include Δ degrees of freedom into the EFT, in which case three-body interactions appear already at NLO [21, 22]; see Refs. [23, 24] for nuclear structure calculations with these NN plus 3N interactions. The reordering of contributions possibly speeds up the convergence of the chiral expansion.

Here we use the formulation of χEFT based on only pion and nucleon degrees of freedom since high order potentials have already been developed for this approach. This implies that we work with the conventional power-counting scheme, and with only NN potentials at LO and NLO, while 3N interactions arise at N²LO. Specifically, within the Low-Energy Nuclear Physics International Collaboration (LENPIC) we use the semilocal momentum-space (SMS) regulated NN potentials from Ref. [20], which have been developed completely up through N⁴LO; and the most accurate LENPIC-SMS NN potential, referred to as N⁴LO⁺, including some contributions from the 6th order in the chiral expansion. The N⁴LO⁺ potential gives a near-perfect description of the mutually compatible neutron-proton and proton-proton scattering data below E _lab = 300 MeV with a ${\chi }_{\text{\hspace{0.17em}datum}}^2=1.01$. At the moment, the accompanying higher n-body forces have not yet been developed to the same chiral order.

Right now, consistent N²LO 3NFs exist, implying that the regularization of the 3N interactions is consistent with that of the NN potential, all relevant symmetries are respected, and the same LEC values are used in the NN and 3N interactions. The strength of the 2π exchange in the N²LO 3NFs (c₁, c₃, and c₄) has been determined from πN scattering, see Table 1 of Ref. [20]. (Note that, for the 3NF, these values need to be shifted as given in Eq. (2.8) of Ref. [25]). We have not taken uncertainties of these c_i’s into account; this should be part of the N³LO uncertainty estimate given below. These 3NFs have already been used for nucleon-deuteron scattering [26], as well as select light nuclei [27, 28]. Consistent N³LO 3NFs are being developed and tested, and are expected to be available for use in many-body calculations soon; similarly, consistent electroweak operators are also under development.

The 3NFs at N²LO depend on two LECs, generally referred to as c_D and c_E; these two LECs have been determined in Ref. [28] by fitting the ³H binding energy (using the Faddeev approach), as well as the experimental proton-deuteron scattering data [29] for the differential cross-section minimum at the proton beam energy of E = 70 MeV. Note that for the determination of c_D and c_E it is important to identify observables that a) provide sufficiently independent constraints, i.e., are sensitive to the 3NFs and are sufficiently uncorrelated; b) can be predicted accurately at N²LO; and c) are measured experimentally with sufficiently high accuracy. This can be achieved by e.g., incorporating properties of ⁴He (and other nuclei), in addition to A = 3 observables, in the fitting of c_D and c_E [30]. However, here we prefer to only use A = 3 data for the determination of c_D and c_E in order to obtain parameter-free predictions for A > 3, and to avoid interference of 4N (and higher-body) interactions at N³LO and higher. In Ref. [31], it has been observed that the triton binding energy and the proton-deuteron scattering cross section minimum at 70 MeV are fulfilling these requirements.

Note that we keep all LECs in the NN potentials fixed at their values determined from NN scattering; and we do not propagate any uncertainties in these LECs through the many-body calculations. Similarly, we have not explicitly propagated uncertainties in the LECs c_D and c_E for the 3NFs through the many-body calculations. In Ref. [31] we did vary c_D and c_E while keeping the ³H binding energy fixed with the LENPIC Semilocal Coordinate Space interaction at N²LO, and the resulting variation in the ⁴He and ¹²C binding energies, while not negligible and in opposite directions, stayed within the chiral truncation uncertainty estimate for a variation of c_D between 6 and 8, the preferred range based on Nd scattering data for that interaction. Furthermore, in Ref. [32] it was shown that the uncertainties in many-body observables of ⁴He and ¹⁶O due to propagation of the uncertainties in determining the LECs at N²LO are much smaller than the chiral truncation errors in those many-body observables at N²LO. We therefore assume here that any variation of the LECs of the NN and 3N interaction is an effect that is of higher order than N²LO and thus those uncertainties are included in the uncertainty due to missing higher chiral orders.

2.1 Chiral truncation uncertainty estimates

Assuming that the chiral expansion of the nuclear interactions translates into a similar expansion for the physical observables, one expects that an observable X follows a similar expansion pattern. Consider therefore an observable X, and write it as

$\begin{array}{lll}\hfill X& =\hfill & \hfill X^{\left(0\right)}+\mathrm{\Delta }{X}^{\left(2\right)}+\mathrm{\Delta }{X}^{\left(3\right)}+\dots ,\end{array}(1)$

where X⁽⁰⁾ is the LO term, ΔX⁽²⁾ = X⁽²⁾ − X⁽⁰⁾ and ΔX⁽³⁾ = X⁽³⁾ − X⁽²⁾ are the NLO and N²LO correction terms, respectively, and the dots represent higher-order corrections. If this observable indeed follows the same expansion pattern as the nuclear interaction itself, then the correction terms ΔX⁽ⁱ⁾ behave like Qⁱ for increasing i, where Q = max(p, M_π)/Λ_B is the chiral expansion parameter (typically the maximum of the relevant momentum p and the pion mass M_π over the breakdown scale Λ_B). Note that there is no term linear in Q in this expansion: the first correction, at NLO, is quadratic in the expansion parameter Q, at least for observables governed solely by the strong interaction. For electroweak observables, the power-counting is different.

For the purpose of a Bayesian analysis, it is more convenient to rewrite this in terms of dimensionless expansion coefficient c_i, with the scale set by an overall reference value X_ref. Thus we can rewrite the expansion for X as

$\begin{array}{lll}\hfill X& =\hfill & \hfill X_{\text{ref}}\left(c_0+c_2{Q}^2+c_3{Q}^3+\dots \right).\end{array}(2)$

Now we can use Bayesian analysis on the coefficients c_i to estimate the chiral truncation uncertainties. Here we follow the Bayesian model of Ref. [33] for pointwise truncation errors with hyperparameters ν₀ = 1.5 and τ₀ = 1.5 [28]. We apply this to the ground state energy of the p-shell nuclei, with the experimental value as our reference value X_ref. Furthermore, we use an effective pion mass of $M_{\pi }^{\text{\hspace{0.17em}eff}}\approx 200$ MeV and a breakdown scale of Λ_B ≈ 650 MeV [27, 28], and therefore a dimensionless expansion parameter Q ≈ 0.31. Note that in Ref. [34] it was observed that the average momentum of the nucleons inside a nucleus increases with A, and one might therefore have to increase Q with A as well; but up to ¹⁶O this average momentum remains below 200 MeV so we use the same value for Q throughout the p-shell. Nevertheless, a Bayesian analysis of correlated uncertainties for ground states and excited states of a subset of p-shell nuclei does suggest a slightly larger value of Q for the upper p-shell [25].

Finally, although we use the LENPIC-SMS NN potentials from LO up to N⁴LO⁺, we only have the corresponding 3NFs at N²LO. We therefore perform our chiral truncation uncertainty analysis for the N²LO through N⁴LO⁺ NN potentials, all in combination with the N²LO 3NFs, as if they were all N²LO interactions; that is, we include only the coefficients c₀, c₂, and c₃ in Eq. 2 (and again, there is no term linear in Q).

3 No-core shell model

3.1 Numerical method

In the No-Core Shell Model (NCSM) [16], the wavefunction Ψ of a nucleus consisting of Z protons and N neutrons is expanded in a finite A = Z + N-body basis of Slater determinants Φ_k of single-particle wavefunctions ${\varphi }_{\mathit{nljm}}(\vec{r})$

$\begin{array}{lll}\hfill \mathrm{\Psi }\left({\vec{r}}_1,\dots ,{\vec{r}}_A\right)& =\hfill & \hfill \sum a_k{\mathrm{\Phi }}_k\left({\vec{r}}_1,\dots ,{\vec{r}}_A\right).\end{array}(3)$

With such an expansion, the many-body Schrödinger equation

$\begin{array}{lll}\hfill \widehat{\mathbf{H}}\mathrm{\Psi }\left({\vec{r}}_1,\dots ,{\vec{r}}_A\right)& =\hfill & \hfill E\mathrm{\Psi }\left({\vec{r}}_1,\dots ,{\vec{r}}_A\right)\end{array}(4)$

becomes an eigenvalue problem

$\begin{array}{lll}\hfill H_{ik}{a}_k& =\hfill & \hfill E{a}_i,\end{array}(5)$

for the coefficients a_k of the expansion in Eq. 3. The matrix H_ik consists of matrix elements ${\mathrm{\Phi }}_i\widehat{\mathbf{H}}{\mathrm{\Phi }}_k$ (where integration over all spatial degrees of freedom is understood) of the many-body hamiltonian

$\begin{array}{lll}\hfill \widehat{\mathbf{H}}& =\hfill & \hfill {\widehat{\mathbf{T}}}_{\text{\hspace{0.17em}rel}}+{\widehat{\mathbf{V}}}_{NN}+{\widehat{\mathbf{V}}}_{3N}+\cdots \end{array}(6)$

consisting of the relative kinetic-energy operator, a two-body potential, and, in general three-body and higher n-body interaction terms. If the interaction is limited to n-body terms, the matrix H_ik for a nucleus with A > n becomes a sparse matrix; in practice, the NCSM is generally applied with up to three-body interactions, and the corresponding Hamiltonian matrices are extremely sparse for A ≥ 6. For any finite basis expansion, the obtained eigenvalue E gives a strict upper bound for the energy in the complete (though infinitely large) basis, at least for the lowest states of a given Z, N, and spin-parity quantum numbers J^P; and the corresponding eigenvector $\vec{a}$ gives an approximation to the A-body wavefunction $\mathrm{\Psi }({\vec{r}}_1,\dots ,{\vec{r}}_A)$. As one increases the basis size, the obtained eigenvalues E of the matrix H_ik approach the exact eigenvalues for a given Hamiltonian $\widehat{\mathbf{H}}$.

In the conventional NCSM one uses a harmonic oscillator (HO) basis for the single-particle wavefunctions ${\varphi }_{\mathit{nljm}}(\vec{r})$, characterised by its scale parameter ℏω. One particular advantage of a HO basis is that one can treat the center-of-mass motion exactly: the Talmi–Moshinksy brackets [35, 36] can be used to convert between HO matrix elements in single-particle coordinates and relative plus center-of-mass coordinates; furthermore, with a many-body truncation on the total number of oscillator quanta in the many-body basis, the obtained wavefunctions factorize exactly into a center-of-mass wavefunction and a relative wavefunction [37, 38]. The single-particle wavefunctions ${\varphi }_{\mathit{nljm}}(\vec{r})$ are labelled by their radial quantum number n, orbital motion quantum number l, total single-particle spin $j=l\pm \frac{1}{2}$, and magnetic projection m which satisfies −j ≤ m ≤ j. In a HO basis, the combination (2n + l) gives the number of HO quanta for each state; thus, in a HO basis with a truncation on ∑_i(2n_i + l_i) over all A nucleons, the factorization of the center-of-mass wavefunction is guaranteed. We add a Lagrange multiplier acting on the center-of-mass coordinates of the many-body system to the Hamiltonian $\widehat{\mathbf{H}}$ to remove center-of-mass excited states from the low-lying spectrum [37, 38]; thus all low-lying states will have a 0s HO center-of-mass wavefunction. Note that this does not alter the eigenvalues nor the eigenvectors for these states, it merely separates the center-of-mass excited states from the states with the lowest center-of-mass motion.

All NCSM calculations presented here were performed using the code Many-Fermion Dynamics–nuclear physics [39–41]. It solves the eigenvalue problem Eq. (5) for the lowest eigenvalues, starting from two- and three-body matrix elements in a HO basis. MFDn is a platform-independent Fortran 90 code using a hybrid MPI+OpenMP programming model. The actual calculations have been performed on Theta at the Argonne Leadership Computing Facility (ALCF) and Cori at the National Energy Research Scientific Computing center (NERSC). For each nucleus and interaction, we performed a series of calculations, using a range of different values of ℏω and the truncation parameter N_max, which is defined as the number of HO quanta above the minimal number of HO quanta in the many-body basis for that nucleus. That is, an N_max = 0 calculation corresponds to a calculation in the lowest oscillator configuration. Here we are only interested in the normal or natural parity states (the parity of the N_max = 0 space), and we increase N_max in steps of 2 starting from N_max = 0 up to at least N_max = 8. Some of the largest calculations for this study were for ¹⁴N and ¹⁵N at N_max = 8, both with dimensions of over one billion, and about 76 × 10¹² nonzero matrix elements, i.e., less than 1 in 10,000 matrix elements is nonzero with three-body interactions for these largest computations.

Of course, for two- and three-body systems it is more efficient and straightforward to work with wavefunctions in relative coordinates, rather than in single-particle wavefunctions. However, beyond four nucleons, the necessary anti-symmetrization becomes increasingly cumbersome in relative coordinates, whereas the NCSM in single-particle coordinates is straightforward to implement for an arbitrarily large number of nucleons; however, the size of the matrix does grow dramatically with the number of nucleons. Nevertheless, in recent years the NCSM has been implemented in Jacobi coordinates (J-NCSM) [42] and applied to (hyper)nuclei with up to eight (hyper)nucleons [43]. The codes MFDn and J-NCSM have been benchmarked against each other, and generally agree to within 10–20 keV for A = 3 and 4, and to within about 30 keV for A = 6, i.e., to within 0.1% of the obtained eigenvalues. The differences of up to about 0.1% have been attributed to differences in the implementations of transforming the three-body forces from their momentum-space expressions to HO matrix elements, including differences in the implementations of the Similarity Renormalization Group (SRG) transformations discussed next.

3.2 Convergence and similarity renormalization group evolution

In the left panel of Figure 1 we show the obtained ground state energy of ⁴He at NLO and N²LO for two different values of the regulator Λ as a function of N_max at ℏω = 24 MeV; as illustration of the effect of the 3NFs, we also include results using only the NN potential at N²LO, without the 3NFs (while at NLO, there are no 3NFs, so there is only the NN potential). Even at N_max = 16, the NCSM results are still several MeV above the corresponding Yakubovsky results, and far from being converged with N_max; and for the upper half of the p-shell nuclei, for A ≥ 10 we are restricted to N_max = 8 in the presence of 3NFs due to computational limitations. Clearly, we have to improve the numerical convergence while keeping the computational needs under control in order to obtain meaningful results for the ground state energies and other observables. There are several methods to do so, which generally fall into four categories (and of course one can also use a combination of these techniques!)

• modify the underlying single-particle wavefunctions to improve the numerical convergence, e.g., start with a Hartree–Fock basis, and/or use natural orbitals [44–46];

• modify the truncation scheme, e.g., select only the most important basis states at each step in N_max (importance-truncated NCSM) [47], or use symmetries to reduce the number of basis states as N_max increases (symmetry-adapted NCSM) [48–50];

• reduce the 3N interaction to an effective NN interaction by normal-ordering the 3N interaction, which typically gains one step in N_max in terms of computational needs [51, 52];

• apply a unitary transformation on the Hamiltonian to improve the convergence at (relatively) small values of N_max [53, 54].

FIGURE 1

FIGURE 1. Convergence of NCSM calculations for ⁴He without (A) and with (B) SRG evolution of the interactions. Induced 3NFs are included in the calculations in the right-hand panel; for comparison results obtained with Yakubovsky calculations in momentum space are shown for the same interactions as well.

Each of these methods has its advantages and drawbacks; and each of them is likely to obfuscate any uncertainty quantification of the numerical results; furthermore, with the first two methods listed above one might lose the exact factorization of the center-of-mass motion. Here we choose to improve the numerical convergence in finite bases by applying a suitable SRG transformation on the Hamiltonian.

The SRG approach [54–57] provides a robust framework for consistently evolving (softening) the Hamiltonian, including three-body terms [58–61], as well as operators for other observables, by applying a unitary transformation on the operator(s) of interest. This unitary transformation is formulated in terms of a flow equation

$\begin{array}{lll}\hfill \frac{d}{d\alpha }{\widehat{\mathbf{H}}}_{\alpha }& =\hfill & \hfill \left[{\eta }_{\alpha },{\widehat{\mathbf{H}}}_{\alpha }\right],\end{array}(7)$

with a continuous flow parameter α. The physics of the SRG evolution is governed by the anti-hermitian generator η_α. A specific form widely used in nuclear physics [54] is given by

$\begin{array}{lll}\hfill {\eta }_{\alpha }& =\hfill & \hfill m_N^2\left[{\widehat{\mathbf{T}}}_{\text{rel}},{\widehat{\mathbf{H}}}_{\alpha }\right],\end{array}(8)$

where m_N is the (average) nucleon mass and ${\widehat{\mathbf{T}}}_{\text{rel}}$ is the relative kinetic-energy operator. This generator drives the Hamiltonian towards a diagonal form in a basis of eigenstates of the intrinsic kinetic energy, i.e., towards a diagonal in momentum space. The initial (or ‘bare’) Hamiltonian provides the initial condition at α = 0 for this flow equation; at NLO, this is just an NN-potential, but at N²LO (and higher orders) it also includes the explicit 3NFs. The width of the diagonal of the potential matrix elements in momentum space is proportional to λ _SRG = 1/α⁴ [62]. For a typical value of α = 0.04 fm⁴, λ _SRG ≈ 2.24 fm⁻¹, which can be considered as an effective cutoff in momentum space; lowering this cutoff improves the convergence of NCSM calculations.

Along with a decoupling of low-momentum and high-momentum components, this SRG induces many-body operators beyond the rank of the initial Hamiltonian. In principle, all induced terms up to the A-body level should be retained to ensure that the transformation is a unitary transformation, such that the spectrum of the Hamiltonian is independent of the flow parameter α. In practice however, one has to truncate these many-body forces induced by the SRG evolution; here we follow the common practice of truncating the SRG evolution at the 3N level, omitting induced four-nucleon (and higher) induced interactions. Of course, this violates unitarity, and therefore introduces a fictitious dependence on the SRG parameter α for A ≥ 4 which we have to monitor, and include in our uncertainty budget. Unfortunately, it is as of yet unclear how to identify an expansion parameter that allows for an estimate of uncertainties due to missing higher-body induced interactions in A ≥ 4 nuclei.

The flow equation for the three-nucleon system is solved numerically using a HO basis in Jacobi-coordinates [60] at a fixed HO basis parameter of ℏω = 36 MeV. The intermediate sums in this three-body Jacobi basis are truncated at N_max = 40 for channels with J < 9/2, N_max = 38 for J = 9/2, and N_max = 36 for all J > 9/2. (Note that the flow equation at the two-nucleon level is solved numerically to a much higher numerical accuracy.) The SRG evolution and transformations first from ℏω = 36 MeV to the desired ℏω value in the range from 14 to 32 MeV, and subsequently from Jacobi coordinates to single-particle coordinates, were all performed on a single multicore CPU node using an efficient OpenMP parallelized code.

In the right panel of Figure 1 we show the ground state energy of ⁴He for the same initial interactions, and the same ℏω = 24 MeV, as in the left panel, but after first performing an SRG evolution of the Hamiltonian to α = 0.04 fm⁴, and including induced 3N interactions, but omitting any induced 4N interactions. Comparing these two panels it is immediately evident that the SRG evolution has indeed dramatically improved the convergence: after the SRG evolution to α = 0.04 fm⁴, the obtained ground state energies at N_max = 4 are already closer to the Yakubovsky results than the ground states energies at N_max = 16 without any SRG evolution, and at N_max = 14 they appear to be almost converged and in agreement with the Yakubovsky results to within a fraction of an MeV. Empirically, α in the range of 0.04 ≤ α ≤ 0.08 fm⁴ appears to be a good compromise between the convergence of the NCSM calculations and minimizing the contributions of the SRG-induced four- and higher-body forces.

3.3 Extrapolating to the complete basis

Although the right-hand panel of Figure 1 looks converged to well within 100 keV, it is not completely converged; furthermore, in the upper half of the p-shell we are limited to at most N_max = 8, at which point the results are clearly not yet converged. However, the approach to convergence appears to be smooth, and if we plot the difference between our results at successive N_max values, at fixed ℏω values near the variational minimum in the largest basis, see Figure 2, it is evident that these differences decrease almost exponentially with increasing N_max. Inspired by this behavior, we therefore use exponential extrapolation in N_max at fixed ℏω,

$\begin{array}{lll}\hfill E^{\hslash \omega }\left(N_{\text{max}}\right)& \approx \hfill & \hfill E_{\infty }^{\hslash \omega }+a{\mathrm{e}}^{\left(-b{N}_{\text{max}}\right)}\end{array}(9)$

based on three consecutive values of N_max at or slightly above the variational minimum in ℏω to extract binding energies in the complete (but infinitely large) basis. Indeed, such an empirical exponential has been widely used for a range of different interactions [63–65], and appears to be reasonably reliable and accurate, at least for true bound states. Furthermore, an exponential approach to convergence for the binding energy is also suggested by various analytic investigations into the asymptotic behavior [66–69]. Insights into the approach to convergence allows one to improve the extrapolation [70], but these analytic expressions generally depend on the underlying structure of the state. Here, we restrict ourselves to the simple ansatz Eq. (9) since it works well for all ground state energies considered, without the need to adapt the extrapolation to the specific structure of each nucleus.

FIGURE 2

FIGURE 2. Convergence of NCSM calculations for ⁴He (A) and ⁶Li (B): Difference in obtained ground energies for successive N_max values for different SRG evolution parameters α.

Following Refs. [25, 28, 34], we take as our best estimate for E_∞ in the complete basis the value of $E_{\infty }^{\hslash \omega }$ near the variational minimum in ℏω for which $|E_{\infty }^{\hslash \omega }-E^{\hslash \omega }(N_{\text{max}})|$ is minimal. Of course, this extrapolation is not exact, and will depend (slightly) on the ℏω value; furthermore, we have to include an extrapolation uncertainty in our uncertainty budget. Again, we resort to an empirical estimate of this uncertainty based on the variation with ℏω and N_max, and our estimate of the corresponding extrapolation uncertainty is the maximum of

• the difference in $E_{\infty }^{\hslash \omega }$ for two successive extrapolations using data for (N_max − 6, N_max − 4, N_max − 2) and (N_max − 4, N_max − 2, N_max) respectively;

• half the variation in $E_{\infty }^{\hslash \omega }$ over a 8 MeV interval in ℏω around the variational minimum;

• 20% of $|E_{\infty }^{\hslash \omega }-E^{\hslash \omega }(N_{\text{max}})|$.

Note that this empirical uncertainty estimate is a conservative estimate, based on calculations with several different interactions [65, 71], and has been shown to give decreasing uncertainties with increasing N_max, with the higher-N_max results generally within the uncertainty estimates of the lower-N_max results. However, these uncertainty estimates cannot be interpreted statistically; for that one should use e.g., the Bayesian analysis of [72].

3.4 Combined numerical uncertainties

In Table 1 we give our extrapolated NCSM ground state energies for ³H and ⁴He, with our extrapolation uncertainty estimates, together with results obtained in momentum space with the Faddeev–Yakubovsky equations [25]. For ³H the NCSM and Faddeev results agree very well, comfortably within the estimated extrapolation uncertainties of the NCSM calculations; and the results obtained with SRG α = 0.08 fm⁴ are more precise than those obtained with α = 0.04 fm⁴, judging by their smaller uncertainties. For the ground state energy of ⁴He we do see differences between the NCSM results and the Yakubovsky calculations, beyond the 10 keV uncertainty in the Yakubovsky calculations and the estimated NCSM uncertainties. These differences can be attributed to the omitted induced 4NFs. They are generally larger with the Λ = 500 MeV regulator than with the Λ = 450 MeV regulator, as one might expect, given that the Λ = 500 MeV interactions are converging slower than the Λ = 450 MeV interactions (see Figure 1); and it may be counter-intuitive that the effects of omitted 4NFs are larger for α = 0.04 fm⁴ than for α = 0.08 fm⁴, but this accidental, as can be seen from Figure 3. This figure also shows good agreement between the Yakubovsky and NCSM calculations.

TABLE 1

TABLE 1. SRG dependence of ground state energies in MeV for A = 3 and 4, compared to Faddeev–Yakubovsky calculations [25]. Explicit N²LO 3NFs are included in the N²LO through N⁴LO⁺ calculations. Quoted uncertainties are the estimated NCSM extrapolation uncertainties only.

FIGURE 3

FIGURE 3. SRG depenence of the ground state energy of ⁴He, with the N⁴LO⁺ NN potentials plus the N²LO 3NFs and Λ = 450 MeV.

For A ≥ 6 we do not have any calculations without SRG evolution for comparisons–or rather, NCSM calculations for these interactions without SRG evolution are so far from convergence for A ≥ 6 that they are not very useful for comparison. However, we can gain insight in effects of omitted induced many-body forces by comparing results obtained with different values for the SRG parameter α, see Tables 2, 3. Table 2 shows a differences of about 0.2 MeV in the binding energies due to the two different SRG parameters α, both for ⁶He and ⁶Li, and almost independent of the chiral order of the NN potential; though at NLO the difference is somewhat smaller, probably due to the lack of an explicit 3N interaction at NLO. However, this difference of about 0.2 MeV is the same order of magnitude as our estimated extrapolation uncertainties at α = 0.04 fm⁴, which prevents one from making firm conclusions.

TABLE 2

TABLE 2. SRG dependence for A = 6 ground state energies in MeV for SRG parameter α = 0.04 fm⁴ and α = 0.08 fm⁴, together with their difference. Quoted uncertainties are the estimated NCSM extrapolation uncertainties only.

TABLE 3

TABLE 3. SRG dependence for select 6 < A < 16 ground state energies in MeV at N²LO for SRG parameter α = 0.04 fm⁴ and α = 0.08 fm⁴, together with their difference. Quoted uncertainties are the estimated NCSM extrapolation uncertainties only.

Note that the extrapolation uncertainties for ⁶He and ⁶Li at α = 0.08 fm⁴ are a factor of two to three smaller than those obtained with α = 0.04 fm ⁴, clearly indicating the improved convergence as the interaction is further evolved with SRG. The exception are the results for A = 6 at NLO; this is most likely caused by the fact that the obtained binding energies at NLO are actually above threshold for ⁶He, and right around threshold for ⁶Li, as was already observed in Ref. [28]. Indeed, for states above threshold, the simple exponential extrapolation may not be very reliable since neglected continuum effects could be significant.

In Table 3 we show the ground state energies for selected stable nuclei with 6 < A < 16 using the N²LO interaction (including 3NFs), SRG evolved to α = 0.04 fm⁴ and α = 0.08 fm⁴. As is the case for A = 6, the convergence improves with the SRG evolution: the extrapolation uncertainty estimates are a factor of 3–8 smaller at α = 0.08 fm⁴ than at α = 0.04 fm⁴. This effect of improved convergence becomes more pronounced as A increases, and is consistent with what we saw for A = 6 in Table 2. Furthermore, at α = 0.04 fm⁴, starting from A = 10, the Λ = 450 MeV interaction converges noticeably better than the Λ = 500 MeV interaction; in qualitative agreement with the picture for ⁴He (see Figure 2); however, at α = 0.08 fm⁴ this difference has washed away, and both regulators give a similar level of convergence for the ground state energies.

Somewhat surprisingly, the difference in ground state energies between the two SRG values, α = 0.04 fm⁴ and α = 0.08 fm⁴, remains almost constant, at around 0.5 MeV, from A = 10 to A = 15, at least for the nuclei considered in Table 3, with only a slight tendency to increase with A, and with similar tendencies for both regulator values. Furthermore, this difference is similar to the estimated NCSM uncertainty at α = 0.08 fm⁴ (which also increases slowly with A), but a factor of 3–8 smaller than the NCSM uncertainty at α = 0.04 fm⁴, which makes it hard to draw a firm conclusion. Nevertheless, based on these observations, we conclude that, for the calculations described here, it is realistic to include an SRG uncertainty that is equal to the NCSM extrapolation uncertainty estimate for stable A = 10 to A = 16 p-shell nuclei.

On the other hand, for A = 6 and ⁸Li the difference in ground state energies between the two SRG values is noticeably larger than the NCSM extrapolation uncertainty for α = 0.08 fm⁴: for ⁸Li it is a factor of two larger; and for A = 6 it is approximately a factor of three larger. Note that this coincides with a larger N_max value used for the lower half of the p-shell: for A = 6 and 7 we can perform our calculations up to N_max = 12 (and are in fact limited by the size of the input files with three-body HO matrix elements) and the calculations for A = 8 and 9 extend up to N_max = 10, whereas for A ≥ 10 we are limited to N_max = 8. Of course this upper limit in N_max also determines the level of numerical convergence that can be achieved, and hence the order of magnitude of the NCSM extrapolation uncertainty. Again, based on these observations we estimate the SRG uncertainty in the binding energy to be about 0.2 MeV for A = 6 and 7, and about 0.3 MeV for A = 8 and 9.

4 Ground state energies of p-shell nuclei

In Tables 4–8 we present our results for the ground state energies of most stable p-shell nuclei, excluding mirror nuclei. We also include ⁸Be, despite it being above the 2α threshold. All calculations were done in the NCSM approach, extrapolated to the (infinitely-large) complete basis, using NN (and 3N) potentials, SRG evolved to α = 0.08 fm⁴, including induced 3N interactions, but omitting higher-body induced interactions. The first set of uncertainties in these tables is our estimate of the combined numerical uncertainties; the second is our estimate of the chiral truncation uncertainty; both as described in the previous section. The numerical uncertainty is estimated strictly based on the numerical convergence pattern and the SRG dependence, and cannot be interpreted statistically. The chiral truncation uncertainty is based on a Bayesian model. We give here the 68% degree of belief (DoB) values.

TABLE 4

TABLE 4. Ground state energies in MeV of ⁴He, ⁶He, ⁶Li, and ⁷Li, for LO through N⁴LO⁺ NN potentials, with 3NFs at N²LO, for N²LO through N⁴LO⁺, for Λ = 450 MeV (top) and Λ = 500 MeV (bottom). Both our estimated numerical uncertainties (first set of uncertainties) and chiral truncation uncertainty estimates (second set of uncertainties, not available for LO) are given. Entries with an * indicate energies above threshold, indicating a resonance, rather than a bound state.

TABLE 5

TABLE 5. Ground state energies in MeV of A = 8 and ⁹Li, for LO through N⁴LO⁺ NN potentials, with 3NFs at N²LO, for N²LO through N⁴LO⁺, for Λ = 450 MeV (top) and Λ = 500 MeV (bottom). Both our estimated numerical (first) and chiral truncation (second) uncertainties are given; and an * indicates ground states with energies above threshold.

TABLE 6

TABLE 6. Ground state energies in MeV of ⁹Be, ¹⁰Be, ¹⁰B, and ¹¹B, for LO through N⁴LO⁺ NN potentials, with 3NFs at N²LO, for N²LO through N⁴LO⁺, for Λ = 450 MeV (top) and Λ = 500 MeV (bottom). Both our estimated numerical (first) and chiral truncation (second) uncertainties are given; and an * indicates ground states with energies above threshold.

TABLE 7

TABLE 7. Ground state energies in MeV of A = 12 and A = 13 nuclei, for LO through N⁴LO⁺ NN potentials, with 3NFs at N²LO, for N²LO through N⁴LO⁺, for Λ = 450 MeV (top) and Λ = 500 MeV (bottom). Both our estimated numerical (first) and chiral truncation (second) uncertainties are given; and an * indicates ground states with energies above threshold.

TABLE 8

TABLE 8. Ground state energies in MeV of ¹⁴C, ¹⁴N, ¹⁵N, and ¹⁶O, for LO through N⁴LO⁺ NN potentials, with 3NFs at N²LO, for N²LO through N⁴LO⁺, for Λ = 450 MeV (top) and Λ = 500 MeV (bottom). Both our estimated numerical (first) and chiral truncation (second) uncertainties are given.

The NCSM calculations for the nuclei presented in Table 4 were performed up to N_max = 14 for ⁴He, and up to N_max = 12 for A = 6 and 7, which is generally sufficient to reach convergence for the ground state energies to within 0.1% (or even better) for a given set of input HO two- and three-body matrix elements, thanks to the interaction being SRG evolved to α = 0.08 fm⁴. Therefore, the numerical uncertainties (the first set of quoted uncertainties in Table 4) are dominated by the uncertainties in the SRG evolution, which is mostly coming from the omission of induced 4-body forces (as well as higher-body forces for A > 4), as well as from the numerical implementation of the SRG evolution and transformations from momentum space expressions to HO matrix elements.

The exceptions are the A = 6 and 7 ground state energies at LO, because it turns out that at LO, these states are not bound, as indicated by the * in the tables: they are all above threshold for decays into α plus two neutrons, or plus a deuteron or a triton, respectively. Hence the numerical convergence of the NCSM calculations is poor (at best it would converge to a quasi-bound state), and neither the extrapolation to the complete basis (nor its uncertainty estimate) is likely to be accurate, which is why the extrapolation uncertainties at LO are noticeably higher than at higher chiral orders. Since we only include the LO results to improve our estimate of the chiral truncation uncertainties, an approximate bound state, or rather, resonance energy is sufficient for our purpose. Similarly, ⁶He appears to be unbound at NLO; but again, for estimating the chiral uncertainty at N²LO that is not a real problem.

The estimated chiral truncation uncertainties are significantly larger than any of the numerical uncertainty estimates. At NLO, these uncertainties are too large to draw any meaningful conclusions, but at N²LO they are, as expected, more than a factor of three smaller. Remember that we use N²LO 3NFs in combination with the higher-order NN potentials, and we therefore apply the N²LO power-counting rules for estimating the chiral uncertainties for these higher-order NN potentials. It should therefore not be surprising that the obtained chiral uncertainty estimates are the same at these higher order as those with the N²LO NN potential. The central values however do change: all of the A = 4, 6, and 7 nuclei become less bound when using NN potentials beyond N²LO in combination with the N²LO 3NFs. This brings the ground state energy of ⁴He in closer agreement with experiment, whereas the ground state energies of ⁶He, ⁶Li, and ⁷Li are reasonably close to experiment with the N²LO NN potential, and move away from their experimental values when including higher-order for the NN potential, to the point that ⁶He appears to be barely bound, or maybe even slightly unbound, with these higher-order NN potentials. However, they all still agree with their corresponding experimental values, to within our combined numerical and chiral uncertainty estimates, and it is therefore too early to draw firm conclusions.

Finally, it is interesting to note that the estimated chiral truncation uncertainties are very similar for each of the four nuclei in Table 4. This can be easily understood in terms of their structure: ⁶He, ⁶Li, and ⁷Li can be described as bound states of an α plus two neutrons, an α plus a deuteron, and an α plus a triton, respectively. It is therefore not surprising that the chiral uncertainties of these states follow that of the ⁴He ground state energy (remember, the deuteron binding energy is fitted exactly, and the triton binding energies is fitted at N²LO and up). However, there are subtle but important details that can make a difference: whereas the ground state energy of ⁴He changes only by about 150–200 keV going from the N²LO to the N³LO NN potential, the difference between these two potentials for the A = 6 ground state energies is about 600–700 keV, and that for the ⁷Li ground state is about 1 MeV, for both regulators.

In Table 5 we give our results the ground state energies of ⁸He, ⁸Li, ⁸Be, and ⁹Li. The NCSM calculations for these nuclei were performed up to N_max = 10, making the extrapolation uncertainties for these nuclei somewhat larger than those in Table 4, and they become of the same order as the estimated numerical uncertainties coming from the SRG evolution. This increased numerical uncertainty is reflected in the first set of error estimates in Table 5. And just as for the A = 6 and 7 nuclei, at LO none of these nuclei are actually bound–leading to larger extrapolation uncertainties. Again, the main purpose for the LO calculations is to set the scale for the estimate of the chiral truncation uncertainties; and for that purpose, an approximate ground state energy is sufficient. Furthermore, ⁸Be is unbound at all chiral orders considered here, in agreement with experiment.

As for the A = 4, 6, and 7 nuclei, the estimated chiral truncation uncertainties for these A = 8 and 9 nuclei is significantly larger than the estimated numerical uncertainties. Again, at NLO, these uncertainties are too large to draw any meaningful conclusions, but at N²LO they are about a factor of three smaller. The agreement with experiment is best with the N²LO NN plus 3N interaction; moving to higher orders for the NN potential while retaining the N²LO 3NFs leads to significant underbinding for ⁸He, ⁸Li, and ⁹Li, with the experimental values outside combined numerical and chiral uncertainty estimates for the N³LO and higher NN potentials. It will be very interesting to see whether or not consistent 3NFs at N³LO can restore or improve on the level of agreement for these ground state energies obtained with N²LO NN plus 3N interactions.

Beryllium-8 remains unbound according to our calculations, for all of these interactions, in qualitative agreement with experiment; and the extracted ground state energies may therefore be not as precise as for the other three nuclei in Table 5. Still, given the combined numerical and chiral uncertainty estimates, our results for the ⁸Be ground state energy are in good agreement with experiment. Furthermore it is interesting to note that the chiral uncertainty estimates for ⁸Be are approximately twice that of ⁴He, whereas the other three nuclei have chiral uncertainty estimates that are quite similar to those in Table 4. This can be easily understood by realizing that ⁸Be is a loosely bound state, or rather, slightly unbound state, of two α particles, so the uncertainty is simply twice that of one α particle. It may be more surprising that the chiral uncertainties of ⁸Li and ⁹Li, neither of which are α-cluster states, are qualitatively similar to that of the A = 4, 6, and 7; and it is also surprising that the estimated chiral uncertainties of ⁸He is smaller than that of any of the other p-shell nuclei.

Moving to the middle of the p-shell, in Table 6 we have our results for the ground state energies of ⁹He, ¹⁰Be, ¹⁰B, and ¹¹B. Starting from A = 10, the NCSM calculations are limited to N_max = 8, and therefore the extrapolation uncertainties become a significant factor in the uncertainty budget. Nevertheless, qualitatively, the overall picture remains the same: at LO all nuclei are unbound, but at NLO and beyond, they are generally bound, with the exception of ⁹Be. Also, the estimated chiral truncation uncertainties for these nuclei remains significantly larger than the estimated numerical uncertainties; at NLO, these uncertainties are too large to draw any meaningful conclusions, but at N²LO and beyond they are about a factor of three smaller. Here, we also start to see significant differences between the chiral uncertainties with the N²LO NN plus 3N interaction, vs. using an NN potential at N³LO or higher in combination with the N²LO 3NFs (and remember, we are using the N²LO counting rules for all these calculations with higher-order NN potential)—the effect of the higher-order NN potentials is becoming more pronounced with increasing A, and more so with Λ = 500 MeV than with Λ = 450 MeV. (This trend is already noticeable for e.g., ⁹Li, see Table 5.).

Around A = 10, the agreement with experiment is no longer uniformly better with the N²LO NN plus 3N interaction than with the higher-order NN potentials. The N²LO NN plus 3N interaction is the only combination for which ⁹Be is truly bound with respect to two α particles plus a neutron, for both regulator values, whereas with the higher-order NN potentials ⁹Be becomes unbound or right at threshold in contrast with experiment where it is bound by about 1.6 MeV. Within the combined uncertainty estimates however, it is still in agreement with experiment for all of these interactions. For the A = 10 ground state energies the situation is different: at Λ = 450 MeV, the N²LO NN plus 3N interaction gives slightly better agreement with experiment than the higher-order NN potentials, but at Λ = 500 MeV, it is the N³LO NN plus N²LO 3N interaction that gives the best agreement with experiment. In fact, the ¹⁰B ground state energy with the N²LO NN plus 3N potential is just outside the combined uncertainty estimates. Furthermore, the N²LO NN potential plus 3NFs give a ground state energy for ¹¹B that is just outside the combined uncertainty estimates with both the Λ = 450 MeV and Λ = 500 MeV regulators.

This trend becomes more pronounced for A ≥ 12, see Tables 7, 8. At LO, the ground states are still unbound, except for ¹²B at Λ = 450 MeV; at NLO they are bound and in agreement with experiment, given the (granted, rather large) uncertainty estimates; and at N²LO they are all significantly overbound, with the experimental values outside the combined numerical and chiral uncertainty estimates. Increasing the chiral order of the NN potential improves the agreement with experiment again: for A = 12 and 13 our results with the N³LO and higher NN potentials, in combination with the N²LO 3NFs, agree with experiment, well within our uncertainty estimates, with both the Λ = 450 MeV and Λ = 500 MeV regulators. For A = 14 this is also the case with Λ = 450 MeV, but Λ = 500 MeV leads to modest overbinding, though still within our uncertainty estimates. Also for ¹⁵N and ¹⁶O the ground state energies agree with experiment with Λ = 450 MeV, but with Λ = 500 MeV there is significant overbinding, with the experimental values just at the edge of our uncertainty intervals.

We have visually summarized our findings in Figure 4, which clearly shows that with the N²LO NN plus 3N interaction one finds good agreement with experiment for the ground state energies of nuclei up to about A = 9, but significant overbinding starting from about A = 11, more than the estimated uncertainties for A = 13 and beyond. On the other hand, using higher-order NN potentials, in combination with N²LO 3NFs, reduces this overbinding in the upper half of the p-shell, while maintaining reasonable agreement, taking into account both numerical and chiral truncation uncertainties, in the lower half of the p-shell, with only a few exceptions, out of the 20 nuclei considered here.

FIGURE 4

FIGURE 4. Comparison of ground state energies of p-shell nuclei between chiral EFT calculations at N²LO and N⁴LO⁺, each for two values of the regulator Λ, and experiment. Both numerical uncertainty estimates (dark colored) and chiral truncation uncertainties (light colored, corresponding to 68% DoB) are shown.

Clearly, for A > 12 our uncertainty estimates for the N²LO NN plus 3N interaction are noticeably smaller than the deviation from both the experimental data and from the calculations with higher-order NN potentials. We speculate that this may be caused by the N²LO fit to the NN scattering data not being sufficiently accurate, and that discrepancies between the N²LO fit and NN data should be taken into account as uncertainties in the LECs, whereas at higher orders in the NN potential, the NN scattering data are described much more accurately, and this is therefore not necessary. (Note that the N²LO NN potential was fitted only up to E_lab = 125 MeV, whereas the N⁴LO⁺ potential was fitted to 260 MeV in Ref. [20].) Of course, one should then also incorporate the uncertainties in the 3NFs, c_D and c_E [30], and propagate all these uncertainties through the many-body bound state calculations [32, 72].

Another possible explanation could be that NN (and 3N) systems cannot sufficiently constrain the LECs–in which case one necessarily has to include properties of A ≥ 4 nuclei for fitting some (or even all) of the LECs. Indeed, impressive progress has been made in recent years along this way, extending ab initio calculations all the way to ²⁰⁸Pb [73], but one loses some of the predictive power of χEFT by incorporating select many-body observables in the fitting procedures, and the results will depend on exactly which observables are included in the fitting. Yet another cause could be that the actual expansion parameter increases with A, as suggested in Ref. [34]. Calculations with consistent 3NFs at N³LO, propagation of the uncertainties in the LECs through the many-body calculations, and Bayesian inference for both the chiral truncation uncertainties and the numerical uncertainties should help to resolve this issue.

Besides this general trend of increasing deviations with increasing A at N²LO, ⁸He and ⁹Li clearly stand out among the N⁴LO⁺ results in Figure 4; and also our predictions for ⁸Li do not agree, to within their estimated uncertainties, with experiment. Interestingly, ⁸He and ⁹Li are two of the most neutron-rich nuclei, out of the 20 nuclei shown in Figure 4, with N − Z = 4 and 3, respectively; and also ⁸Li is a neutron-rich nucleus. This could be an indication of some deficiencies in the neutron-neutron (or three-neutron) part of the interactions. Unfortunately, there are no accurate neutron-neutron data, let alone three-neutron data, to constrain the LECs; the LECs of the interactions were all fitted to 2- and 3-body data involving at least one proton.

5 Concluding remarks and outlook

We have performed systematic calculations for the binding energies of p-shell nuclei using LENPIC-SMS χEFT NN and 3N interactions complete up through N²LO, and with NN potentials up to N⁴LO⁺ in combination with N²LO 3NFs. We have made a careful analysis of all sources of uncertainties, and incorporated our best estimates of these uncertainties in our comprehensive tables with order-by-order results and in Figure 4. Note that all LECs in the χEFT had been fitted to A = 2 and A = 3 data prior to these many-body calculations, and the obtained binding energies are therefore parameter-free predictions. Although our results with the N²LO NN plus 3N interaction do not agree with the experimental binding energies for the upper p-shell, our results with the N⁴LO⁺ NN potential plus N²LO 3NFs do agree with experiment throughout the p-shell within the combined numerical uncertainty estimates and the chiral truncation uncertainty estimates at the 68% DoB.

In future work we plan to extend these calculations to include consistent N³LO 3NFs, which should bring the chiral truncation uncertainties down, and they may become comparable to the estimated numerical uncertainties. We therefore also intend to further reduce our numerical uncertainties; promising new developments include, among others, the use of Artificial Neural Networks [74–76] and Bayesian inference [72] for extrapolating NCSM binding energies to the complete basis. The latter is particularly interesting, since with Bayesian methods for both the numerical and the chiral truncation uncertainties one can consider correlated uncertainties of different states. This naturally leads to reduced uncertainties for excitation energies (compared to the uncertainties on the binding energies themselves), as well as e.g., neutron separation energies and various cluster thresholds.

Last but not least, we plan to use the obtained wavefunctions, in combination with consistent χEFT operators, to evaluate other observables, in particular radii, charge densities, magnetic and quadrupole moments, and electroweak transitions.

Data availability statement

The original contributions presented in the study are included in the article, further inquiries can be directed to the corresponding author.

Author contributions

All authors contributed to the calculations presented in this article, utilizing several different codes, developed at Iowa State University, the Forschungszentrum Jülich, and the Technische Universität Darmstadt. The initial draft for this article was written by PM, and all authors contributed to the subsequent discussions and final version of this article.

Funding

This work was supported by the US Department of Energy under Grants DE-FG02-87ER40371, DE-SC0018223 and DE-SC0023495, by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) Project-ID 279384907, SFB 1245, by DFG and NSFC through funds provided to the Sino-German CRC 110 “Symmetries and 329 the Emergence of Structure in QCD” (NSFC Grant No. 12070131001, Project-ID 196253076–TRR 110), by the BMBF through Verbundprojekt 05P2021 (ErUM-FSP T07, Contract No. 05P21RDFNB), and by the MKW NRW under the funding code NW21-024-A. This research used resources of the National Energy Research Scientific Computing Center (NERSC) and the Argonne Leadership Computing Facility (ALCF), which are US Department of Energy Office of Science user facilities, supported under Contracts No. DE-AC02-05CH11231 and No. DE-AC02-06CH11357, and computing resources provided under the INCITE award ‘Nuclear Structure and Nuclear Reactions’ from the US Department of Energy, Office of Advanced Scientific Computing Research. Further computing resources were provided on LICHTENBERG II at the TU Darmstadt and on JURECA and the JURECA Booster of the Jülich Supercomputing Center, Jülich, Germany.

Acknowledgments

We thank all members of the LENPIC collaboration for useful and inspiring discussions. In particular, we would like to thank Evgeny Epelbaum, Dick Furnstahl, Kai Hebeler, Hermann Krebs and Jordan Melendez.

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.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

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

CrossRef Full Text | Google Scholar

2. Machleidt R, Holinde K, Elster C. The Bonn meson exchange model for the nucleon nucleon interaction. Phys Rept (1987) 149:1–89. doi:10.1016/S0370-1573(87)80002-9

CrossRef Full Text | Google Scholar

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

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

CrossRef Full Text | Google Scholar

5. Shirokov AM, Shin IJ, Kim Y, Sosonkina M, Maris P, Vary JP. N3LO NN interaction adjusted to light nuclei in ab exitu approach. Phys Lett B (2016) 761:87–91. doi:10.1016/j.physletb.2016.08.006

CrossRef Full Text | Google Scholar

6. Coon SA, Scadron MD, McNamee PC, Barrett BR, Blatt DWE, McKellar BHJ. The two pion exchange, three nucleon potential and nuclear matter. Nucl Phys A (1979) 317:242–78. doi:10.1016/0375-9474(79)90462-7

CrossRef Full Text | Google Scholar

7. Pudliner BS, Pandharipande VR, Carlson J, Pieper SC, Wiringa RB. Quantum Monte Carlo calculations of nuclei with A<∼7. Phys Rev C (1997) 56:1720–50. doi:10.1103/PhysRevC.56.1720

CrossRef Full Text | Google Scholar

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

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

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

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

CrossRef Full Text | Google Scholar

12. Hammer HW, König S, van Kolck U. Nuclear effective field theory: Status and perspectives. Rev Mod Phys (2020) 92:025004. doi:10.1103/RevModPhys.92.025004

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

14. Lähde TA, Meißner UG. Nuclear lattice effective field theory: An introduction. Cham: Springer (2019). doi:10.1007/978-3-030-14189-9

CrossRef Full Text | Google Scholar

15. Lee D. Recent progress in nuclear lattice simulations. Front Phys (2020) 8:174. doi:10.3389/fphy.2020.00174

CrossRef Full Text | Google Scholar

16. Barrett BR, Navratil P, Vary JP. Ab initio no core shell model. Prog Part Nucl Phys (2013) 69:131–81. doi:10.1016/j.ppnp.2012.10.003

CrossRef Full Text | Google Scholar

17. Hagen G, Papenbrock T, Hjorth-Jensen M, Dean DJ. Coupled-cluster computations of atomic nuclei. Rept Prog Phys (2014) 77:096302. doi:10.1088/0034-4885/77/9/096302

CrossRef Full Text | Google Scholar

18. Hergert H, Bogner SK, Morris TD, Schwenk A, Tsukiyama K. The in-medium similarity renormalization group: A novel ab initio method for nuclei. Phys Rept (2016) 621:165–222. doi:10.1016/j.physrep.2015.12.007

CrossRef Full Text | Google Scholar

19. Otsuka T, Honma M, Mizusaki T, Shimizu N, Utsuno Y. Monte Carlo shell model for atomic nuclei. Prog Part Nucl Phys (2001) 47:319–400. doi:10.1016/S0146-6410(01)00157-0

CrossRef Full Text | Google Scholar

20. Reinert P, Krebs H, Epelbaum E. Semilocal momentum-space regularized chiral two-nucleon potentials up to fifth order. Eur Phys J A (2018) 54:86. doi:10.1140/epja/i2018-12516-4

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

22. Piarulli M, Girlanda L, Schiavilla R, Navarro Pérez R, Amaro JE, Ruiz Arriola E. 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

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

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

25. Bernard V, Epelbaum E, Krebs H, Meißner UG. Subleading contributions to the chiral three-nucleon force: Long-range terms. Phys Rev C (2008) 77:064004. doi:10.1103/PhysRevC.77.064004

CrossRef Full Text | Google Scholar

26. Epelbaum E, Golak J, Hebeler K, Kamada H, Krebs H, Meißner UG, et al. Towards high-order calculations of three-nucleon scattering in chiral effective field theory. Eur Phys J A (2020) 56:92. doi:10.1140/epja/s10050-020-00102-2

CrossRef Full Text | Google Scholar

27. Maris P, Roth R, Epelbaum E, Furnstahl RJ, Golak J, Hebeler K, et al. Nuclear properties with semilocal momentum-space regularized chiral interactions beyond N2LO. Phys Rev C (2022) 106:064002. doi:10.1103/PhysRevC.106.064002

CrossRef Full Text | Google Scholar

28. Maris P, Epelbaum E, Furnstahl RJ, Golak J, Hebeler K, Hüther T, et al. Light nuclei with semilocal momentum-space regularized chiral interactions up to third order. Phys Rev C (2021) 103:054001. doi:10.1103/PhysRevC.103.054001

CrossRef Full Text | Google Scholar

29. Sekiguchi K, Sakai H, Witala H, Glöckle W, Golak J, Hatano M, et al. Complete set of precise deuteron analyzing powers at intermediate energies: Comparison with modern nuclear force predictions. Phys Rev C (2002) 65:034003. doi:10.1103/PhysRevC.65.034003

CrossRef Full Text | Google Scholar

30. Wesolowski S, Svensson I, Ekström A, Forssén C, Furnstahl RJ, Melendez JA, et al. Rigorous constraints on three-nucleon forces in chiral effective field theory from fast and accurate calculations of few-body observables. Phys Rev C (2021) 104:064001. doi:10.1103/PhysRevC.104.064001

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

32. Carlsson BD, Ekström A, Forssén C, Strömberg DF, Jansen GR, Lilja O, et al. Uncertainty analysis and order-by-order optimization of chiral nuclear interactions. Phys Rev X (2016) 6:011019. doi:10.1103/PhysRevX.6.011019

CrossRef Full Text | Google Scholar

33. Melendez JA, Furnstahl RJ, Phillips DR, Pratola MT, Wesolowski S. Quantifying correlated truncation errors in effective field theory. Phys Rev C (2019) 100:044001. doi:10.1103/PhysRevC.100.044001

CrossRef Full Text | Google Scholar

34. Binder S, Calci A, Epelbaum E, Furnstahl RJ, Golak J, Hebeler K, et al. Few-nucleon and many-nucleon systems with semilocal coordinate-space regularized chiral nucleon-nucleon forces. Phys Rev C (2018) 98:014002. doi:10.1103/PhysRevC.98.014002

CrossRef Full Text | Google Scholar

35. Talmi I. Nuclear spectroscopy with harmonic oscillator wave-functions. Helv Phys Acta (1952) 25:185–234. doi:10.3929/ethz-a-000089257

CrossRef Full Text | Google Scholar

36. Moshinsky M. Transformation brackets for harmonic oscillator functions. Nucl Phys (1959) 13:104–16. doi:10.1016/0029-5582(59)90143-9

CrossRef Full Text | Google Scholar

37. Lipkin HJ. Center-of-Mass motion in the nuclear shell model. Phys Rev (1958) 110:1395–7. doi:10.1103/PhysRev.110.1395

CrossRef Full Text | Google Scholar

38. Gloeckner DH, Lawson RD. Spurious center-of-mass motion. Phys Lett B (1974) 53:313–8. doi:10.1016/0370-2693(74)90390-6

CrossRef Full Text | Google Scholar

39. Aktulga HM, Yang C, Ng EG, Maris P, Vary JP. Improving the scalability of a symmetric iterative eigensolver for multi-core platforms. Concurrency Comput Pract Experience (2014) 26:2631–51. doi:10.1002/cpe.3129

CrossRef Full Text | Google Scholar

40. Shao M, Aktulga HM, Yang C, Ng EG, Maris P, Vary JP. Accelerating nuclear configuration interaction calculations through a preconditioned block iterative eigensolver. Comput Phys Commun (2017) 222:1–13. doi:10.1016/j.cpc.2017.09.004

CrossRef Full Text | Google Scholar

41. Cook B, Maris P, Shao M, Wichmann N, Wagner M, O’Neill J, et al. High performance optimizations for nuclear physics MFDn on KNL. In: M Taufer, B Mohr, and JM Kunkel, editors. High performance computing. Cham: Springer International Publishing (2016). p. 366–77.

CrossRef Full Text | Google Scholar

42. Liebig S, Meißner UG, Nogga A. Jacobi no-core shell model for p-shell nuclei. Eur Phys J A (2016) 52:103. doi:10.1140/epja/i2016-16103-5

CrossRef Full Text | Google Scholar

43. Le H, Haidenbauer J, Meißner UG, Nogga A. Jacobi no-core shell model for p-shell hypernuclei. Eur Phys J A (2020) 56:301. doi:10.1140/epja/s10050-020-00314-6

CrossRef Full Text | Google Scholar

44. Constantinou C, Caprio MA, Vary JP, Maris P. Natural orbital description of the halo nucleus ⁶He. Nucl Sci Tech (2017) 28:179. doi:10.1007/s41365-017-0332-6

CrossRef Full Text | Google Scholar

45. Tichai A, Müller J, Vobig K, Roth R. Natural orbitals for ab initio no-core shell model calculations. Phys Rev C (2019) 99:034321. doi:10.1103/PhysRevC.99.034321

CrossRef Full Text | Google Scholar

46. Fasano PJ, Constantinou C, Caprio MA, Maris P, Vary JP. Natural orbitals for the ab initio no-core configuration interaction approach. Phys Rev C (2022) 105:054301. doi:10.1103/PhysRevC.105.054301

CrossRef Full Text | Google Scholar

47. Roth R. Importance truncation for large-scale configuration interaction approaches. Phys Rev C (2009) 79:064324. doi:10.1103/PhysRevC.79.064324

CrossRef Full Text | Google Scholar

48. Dytrych T, Sviratcheva KD, Bahri C, Draayer JP, Vary JP. Evidence for symplectic symmetry in ab initio no-core shell model results for light nuclei. Phys Rev Lett (2007) 98:162503. doi:10.1103/PhysRevLett.98.162503

PubMed Abstract | CrossRef Full Text | Google Scholar

49. McCoy AE, Caprio MA, Dytrych T. Symplectic no-core configuration interaction framework for ab initio nuclear structure. Ann Acad Rom Sci Ser Phys Chem Sci (2018) 3:17. arXiv.1802.01771

Google Scholar

50. Launey KD, Mercenne A, Dytrych T. Nuclear dynamics and Reactions in the ab initio symmetry-adapted framework. Ann Rev Nucl Part Sci (2021) 71:253–77. doi:10.1146/annurev-nucl-102419-033316

CrossRef Full Text | Google Scholar

51. Hagen G, Papenbrock T, Dean DJ, Schwenk A, Nogga A, Wloch M, et al. Coupled-cluster theory for three-body Hamiltonians. Phys Rev C (2007) 76:034302. doi:10.1103/PhysRevC.76.034302

CrossRef Full Text | Google Scholar

52. Roth R, Binder S, Vobig K, Calci A, Langhammer J, Navratil P. Medium-mass nuclei with normal-ordered Chiral NN+3N Interactions. Phys Rev Lett (2012) 109:052501. doi:10.1103/PhysRevLett.109.052501

PubMed Abstract | CrossRef Full Text | Google Scholar

53. Roth R, Reinhardt S, Hergert H. Unitary correlation operator method and similarity renormalization group: Connections and differences. Phys Rev C (2008) 77:064003. doi:10.1103/PhysRevC.77.064003

CrossRef Full Text | Google Scholar

54. Bogner SK, Furnstahl RJ, Schwenk A. From low-momentum interactions to nuclear structure. Prog Part Nucl Phys (2010) 65:94–147. doi:10.1016/j.ppnp.2010.03.001

CrossRef Full Text | Google Scholar

55. Glazek SD, Wilson KG. Renormalization of Hamiltonians. Phys Rev D (1993) 48:5863–72. doi:10.1103/PhysRevD.48.5863

CrossRef Full Text | Google Scholar

56. Wegner F. Flow-equations for Hamiltonians. Annalen der Physik (1994) 506:77–91. doi:10.1002/andp.19945060203

CrossRef Full Text | Google Scholar

57. Bogner SK, Furnstahl RJ, Maris P, Perry RJ, Schwenk A, Vary JP. Convergence in the no-core shell model with low-momentum two-nucleon interactions. Nucl Phys A (2008) 801:21–42. doi:10.1016/j.nuclphysa.2007.12.008

CrossRef Full Text | Google Scholar

58. Jurgenson ED, Navratil P, Furnstahl RJ. Evolution of nuclear many-body forces with the similarity renormalization group. Phys Rev Lett (2009) 103:082501. doi:10.1103/PhysRevLett.103.082501

PubMed Abstract | CrossRef Full Text | Google Scholar

59. Roth R, Langhammer J, Calci A, Binder S, Navratil P. Similarity-transformed chiral NN+3N interactions for the ab initio description of 12-C and 16-O. Phys Rev Lett (2011) 107:072501. doi:10.1103/PhysRevLett.107.072501

PubMed Abstract | CrossRef Full Text | Google Scholar

60. Roth R, Calci A, Langhammer J, Binder S. Evolved chiral NN+3N Hamiltonians for ab initio nuclear structure calculations. Phys Rev C (2014) 90:024325. doi:10.1103/PhysRevC.90.024325

CrossRef Full Text | Google Scholar

61. Binder S, Langhammer J, Calci A, Roth R. Ab initio path to heavy nuclei. Phys Lett B (2014) 736:119–23. doi:10.1016/j.physletb.2014.07.010

CrossRef Full Text | Google Scholar

62. Bogner SK, Furnstahl RJ, Perry RJ, Schwenk A. Are low-energy nuclear observables sensitive to high-energy phase shifts? Phys Lett B (2007) 649:488–93. doi:10.1016/j.physletb.2007.04.048

CrossRef Full Text | Google Scholar

63. Maris P, Vary JP, Shirokov AM. Ab initio no-core full configuration calculations of light nuclei. Phys Rev C (2009) 79:014308. doi:10.1103/PhysRevC.79.014308

CrossRef Full Text | Google Scholar

64. Jurgenson ED, Maris P, Furnstahl RJ, Navratil P, Ormand WE, Vary JP. Structure of p-shell nuclei using three-nucleon interactions evolved with the similarity renormalization group. Phys Rev C (2013) 87:054312. doi:10.1103/PhysRevC.87.054312

CrossRef Full Text | Google Scholar

65. Maris P, Shin IJ, Vary JP. Ab initio structure of p-shell nuclei with chiral effective field theory and Daejeon16 interactions. In: 6th International Conference Nuclear Theory in the Supercomputing Era (2019). arXiv.1908.00155

Google Scholar

66. Coon SA, Avetian MI, Kruse MKG, van Kolck U, Maris P, Vary JP. Convergence properties of ab initio calculations of light nuclei in a harmonic oscillator basis. Phys Rev C (2012) 86:054002. doi:10.1103/PhysRevC.86.054002

CrossRef Full Text | Google Scholar

67. Furnstahl RJ, Hagen G, Papenbrock T. Corrections to nuclear energies and radii in finite oscillator spaces. Phys Rev C (2012) 86:031301. doi:10.1103/PhysRevC.86.031301

CrossRef Full Text | Google Scholar

68. More SN, Ekström A, Furnstahl RJ, Hagen G, Papenbrock T. Universal properties of infrared oscillator basis extrapolations. Phys Rev C (2013) 87:044326. doi:10.1103/PhysRevC.87.044326

CrossRef Full Text | Google Scholar

69. Wendt KA, Forssén C, Papenbrock T, Sääf D. Infrared length scale and extrapolations for the no-core shell model. Phys Rev C (2015) 91:061301. doi:10.1103/PhysRevC.91.061301

CrossRef Full Text | Google Scholar

70. Forssén C, Carlsson BD, Johansson HT, Sääf D, Bansal A, Hagen G, et al. Large-scale exact diagonalizations reveal low-momentum scales of nuclei. Phys Rev C (2018) 97:034328. doi:10.1103/PhysRevC.97.034328

CrossRef Full Text | Google Scholar

71. Maris P, Vary JP. Ab initio nuclear structure calculations of p-shell nuclei with JISP16. Int J Mod Phys E (2013) 22:1330016. doi:10.1142/S0218301313300166

CrossRef Full Text | Google Scholar

72. Gazda D, Htun TY, Forssén C. Nuclear physics uncertainties in light hypernuclei. Phys Rev C (2022) 106:054001. doi:10.1103/PhysRevC.106.054001

CrossRef Full Text | Google Scholar

73. Hu B, Jiang W, Miyagi T, Sun Z, Ekstrom A, Forssen C, et al. Ab initio predictions link the neutron skin of ²⁰⁸Pb to nuclear forces. Nat Phys (2022) 18:1196–200. doi:10.1038/s41567-022-01715-8

PubMed Abstract | CrossRef Full Text | Google Scholar

74. Negoita GA, Vary JP, Luecke GR, Maris P, Shirokov AM, Shin IJ, et al. Deep learning: Extrapolation tool for ab initio nuclear theory. Phys Rev C (2019) 99:054308. doi:10.1103/PhysRevC.99.054308

CrossRef Full Text | Google Scholar

75. Jiang WG, Hagen G, Papenbrock T. Extrapolation of nuclear structure observables with artificial neural networks. Phys Rev C (2019) 100:054326. doi:10.1103/PhysRevC.100.054326

CrossRef Full Text | Google Scholar

76. Knöll M, Wolfgruber T, Agel ML, Wenz C, Roth R. Machine learning for the prediction of converged energies from ab initio nuclear structure calculations. Phys Lett B (2023) 839:137781. doi:10.1016/j.physletb.2023.137781

CrossRef Full Text | Google Scholar