Uncertainties in ab initio nuclear structure calculations with chiral interactions

We present theoretical ground state energies and their uncertainties for p-shell nuclei obtained from chiral effective field theory internucleon interactions as a function of chiral order, fitted to two- and three-body data only. We apply a Similary Renormalization Group transformation to improve the numerical convergence of the many-body calculations, and discuss both the numerical uncertainties arising from basis truncations and those from omitted induced many-body forces, as well as chiral truncation uncertainties. With complete Next-to-Next-to-Leading (N2LO) order two- and three-body interactions, we find significant overbinding for the ground states in the upper p-shell, but using higher-order two-body potentials, in combination with N2LO three-body forces, our predictions agree with experiment throughout the p-shell to within our combined estimated uncertainties. The uncertainties due to chiral order truncation are noticeably larger than the numerical uncertainties, but they are expected to become comparable to the numerical uncertainties at complete N3LO.

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 II. In order to assess the numerical uncertainties in our NCSM calculations, we make a detailed comparison with Faddeev-Yakubovsky calculations for 3 H and 4 He using the same interactions; this is described in Section III, together with details about the NCSM. Our results for the binding energies of p-shell nuclei are presented in Section IV. Finally, we give some concluding remarks in Section V.

II. 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 2 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 2 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 4 LO; and the most accurate LENPIC-SMS NN potential, referred to as N 4 LO + , including some contributions from the 6th order in the chiral expansion. The N 4 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 χ 2 datum = 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 2 LO 3NFs exist, implying that the regularization of the 3N interactions is consistent with that of the NN potential, that all relevant symmetries are respected, and that the same LEC values are used in the NN and 3N interactions. The strength of the 2π exchange in the N 2 LO 3NFs (c 1 , c 3 , and c 4 ) 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 3 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 3 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 2 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 3 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 2 LO; and (c) are measured experimentally with sufficiently high accuracy. This can be achieved by e.g. incorporating properties of 4 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 3 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 3 H binding energy fixed with the LENPIC Semilocal Coordinate Space interaction at N 2 LO, and the resulting variation in the 4 He and 12 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 4 He and 16 O due to propagation of the uncertainties in determining the LECs at N 2 LO are much smaller than the chiral truncation errors in those many-body observables at N 2 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 2 LO and thus those uncertainties are included in the uncertainty due to missing higher chiral orders.

A. 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 where (2) are the NLO and N 2 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 (i) behave like Q i 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 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 ν 0 = 1.5 and τ 0 = 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 eff π ≈ 200 MeV and a breakdown scale of Λ B ≈ 650 MeV [26,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 16 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 [27].
Finally, although we use the LENPIC-SMS NN potentials from LO up to N 4 LO + , we only have the corresponding 3NFs at N 2 LO. We therefore perform our chiral truncation uncertainty analysis for the N 2 LO through N 4 LO + NN potentials, all in combination with the N 2 LO 3NFs, as if they were all N 2 LO interactions; that is, we include only the coefficients c 0 , c 2 , and c 3 in Eq. (2) (and again, there is no term linear in Q).

A. 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 ϕ nljm (⃗ r) With such an expansion, the many-body Schrödinger equation becomes an eigenvalue problem for the coefficients a k of the expansion in Eq. (3). The matrix H ik consists of matrix elements Φ i HΦ k (where integration over all spatial degrees of freedom is understood) of the many-body hamiltonian 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 ⃗ a gives an approximation to the A-body wavefunction Ψ(⃗ r 1 , . . . , ⃗ 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 H.
In the conventional NCSM one uses a harmonic oscillator (HO) basis for the single-particle wavefunctions ϕ nljm (⃗ 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 ϕ nljm (⃗ r) are labelled by their radial quantum number n, orbital motion quantum number l, total single-particle spin j = l ± 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 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][40][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 14 N and 15 N at N max = 8, both with dimensions of over one billion, and about 76 × 10 12 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 to 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.

B. Convergence and Similarity Renormalization Group evolution
In the left panel of Fig. 1 we show the obtained ground state energy of 4 He at NLO and N 2 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 2 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][45][46]; • modify the truncation scheme, e.g. select only the most important basis states at each step in N max (importancetruncated NCSM) [47], or use symmetries to reduce the number of basis states as N max increases (symmetryadapted NCSM) [48][49][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].
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][55][56][57] provides a robust framework for consistently evolving (softening) the Hamiltonian, including three-body terms [58][59][60][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 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 where m N is the (average) nucleon mass and T 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 2 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/α 4 [62]. For a typical value of α = 0.04 fm 4 , λ SRG ≈ 2.24 fm −1 , 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 Fig. 1 we show the ground state energy of 4 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 4 , 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 4 , 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 4 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.

C. Extrapolating to the complete basis
Although the right-hand panel of Fig. 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 Fig. 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 ℏω, 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][64][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][67][68][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. Following Refs. [27,28,34], we take as our best estimate for E ∞ in the complete basis the value of E ℏω ∞ near the variational minimum in ℏω for which |E ℏω ∞ − E ℏω (N 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 ℏω ∞ 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 ℏω ∞ over a 8 MeV interval in ℏω around the variational minimum; 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].

D. Combined numerical uncertainties
In Table I we give our extrapolated NCSM ground state energies for 3 H and 4 He, with our extrapolation uncertainty estimates, together with results obtained in momentum space with the Faddeev-Yakubovsky equations [27]. For 3 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 4 are more precise than those obtained with α = 0.04 fm 4 , judging by their smaller uncertainties. For the ground state energy of 4 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 Fig. 1); and it may be counter-intuitive that the effects of omitted 4NFs are larger for α = 0.04 fm 4 than for α = 0.08 fm 4 , but this accidental, as can be seen from Fig. 3. This figure also shows good agreement between the Yakubovsky and NCSM calculations.
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 II and III. Table II shows a differences of about 0.2 MeV in the binding energies due to the two different SRG parameters α, both for 6 He and 6 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 4 , which prevents one from making firm conclusions.   Note that the extrapolation uncertainties for 6 He and 6 Li at α = 0.08 fm 4 are a factor of two to three smaller than those obtained with α = 0.04 fm 4 , 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 6 He, and right around threshold for 6 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 III    (including 3NFs), SRG evolved to α = 0.04 fm 4 and α = 0.08 fm 4 . As is the case for A = 6, the convergence improves with the SRG evolution: the extrapolation uncertainty estimates are a factor of 3 to 8 smaller at α = 0.08 fm 4 than at α = 0.04 fm 4 . This effect of improved convergence becomes more pronounced as A increases, and is consistent with what we saw for A = 6 in Table II. Furthermore, at α = 0.04 fm 4 , starting from A = 10, the Λ = 450 MeV interaction converges noticeably better than the Λ = 500 MeV interaction; in qualitative agreement with the picture for 4 He (see Fig. 2); however, at α = 0.08 fm 4 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 4 and α = 0.08 fm 4 , remains almost constant, at around 0.5 MeV, from A = 10 to A = 15, at least for the nuclei considered in Table III, 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 4 (which also increases slowly with A), but a factor of 3 to 8 smaller than the NCSM uncertainty at α = 0.04 fm 4 , 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 8 Li the difference in ground state energies between the two SRG values is noticeably larger than the NCSM extrapolation uncertainty for α = 0.08 fm 4 : for 8 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.

IV. GROUND STATE ENERGIES OF P-SHELL NUCLEI
In Tables IV-VIII we present our results for the ground state energies of most stable p-shell nuclei, excluding mirror nuclei. We also include 8 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 4 , 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.
The NCSM calculations for the nuclei presented in Table IV were performed up to N max = 14 for 4 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 4 . Therefore, the numerical uncertainties (the first set of quoted uncertainties in Table IV) 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 the chiral truncation    uncertainties, an approximate bound state, or rather, resonance energy is sufficient for our purpose. Similarly, 6 He appears to be unbound at NLO; but again, for estimating the chiral uncertainty at N 2 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 2 LO they are, as expected, more than a factor of three smaller. Remember that we use N 2 LO 3NFs in combination with the higher-order NN potentials, and we therefore apply the N 2 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 2 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 2 LO in combination with the N 2 LO 3NFs. This brings the ground state energy of 4 He in closer agreement with experiment, whereas the ground state energies of 6 He, 6 Li, and 7 Li are reasonably close to experiment with the N 2 LO NN potential, and move away from their experimental values when including higher-order for the NN potential, to the point that 6 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 IV. This can be easily understood in terms of their structure: 6 He, 6 Li, and 7 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 4 He ground state energy (remember, the deuteron binding energy is fitted exactly, and the triton binding energies is fitted at N 2 LO and up). However, there are subtle but important details that can make a difference: whereas the ground state energy of 4 He changes only by about 150 to 200 keV going from the N 2 LO to the N 3 LO NN potential, the difference between these two potentials for the A = 6 ground state energies is about 600 to 700 keV, and that for the 7 Li ground state is about 1 MeV, for both regulators.
In Table V we give our results the ground state energies of 8 He, 8 Li, 8 Be, and 9 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 IV, 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 V. 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, 8 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 2 LO they are about a factor of three smaller. The agreement with  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 V. Still, given the combined numerical and chiral uncertainty estimates, our results for the 8 Be ground state energy are in good agreement with experiment. Furthermore it is interesting to note that the chiral uncertainty estimates for 8 Be are approximately twice that of 4 He, whereas the other three nuclei have chiral uncertainty estimates that are quite similar to those in Table IV. This can be easily understood by realizing that 8 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 8 Li and 9 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 8 He is smaller than that of any of the other p-shell nuclei.
Moving to the middle of the p-shell, in Table VI we have our results for the ground state energies of 9 He, 10 Be, 10 B, and 11 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 9 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 2 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 2 LO NN plus 3N interaction, vs. using an NN potential at N 3 LO or higher in combination with the N 2 LO 3NFs (and remember, we are using the N 2 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. 9 Li, see Table V.) Around A = 10, the agreement with experiment is no longer uniformly better with the N 2 LO NN plus 3N interaction than with the higher-order NN potentials. The N 2 LO NN plus 3N interaction is the only combination for which 9 Be is truly bound with respect to two α particles plus a neutron, for both regulator values, whereas with the higher-order NN potentials 9 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 2 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 3 LO NN plus N 2 LO 3N interaction that gives the best agreement with experiment. In fact, the 10 B ground  Tables VII and VIII. At LO, the ground states are still unbound, except for 12 B at Λ = 450 MeV; at NLO they are bound and in agreement with experiment, given the (granted, rather large) uncertainty estimates; and at N 2 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 3 LO and higher NN potentials, in combination with the N 2 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 15   experimental values just at the edge of our uncertainty intervals. We have visually summarized our findings in Fig. 4, which clearly shows that with the N 2 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 2 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.
Clearly, for A > 12 our uncertainty estimates for the N 2 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 2 LO fit to the NN scattering data not being sufficiently accurate, and that discrepancies between the N 2 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 2 LO NN potential was fitted only up to E lab = 125 MeV, whereas the N 4 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 208 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 3 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 2 LO, 8 He and 9 Li clearly stand out among the N 4 LO + results in Fig. 4; and also our predictions for 8 Li do not agree, to within their estimated uncertainties, with experiment. Interestingly, 8 He and 9 Li are two of the most neutron-rich nuclei, out of the 20 nuclei shown in Fig. 4, with N − Z = 4 and 3, respectively; and also 8 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.

V. 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 2 LO, and with NN potentials up to N 4 LO + in combination with N 2 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 Fig. 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 2 LO NN plus 3N interaction do not agree with the experimental binding energies for the upper p-shell, our results with the N 4 LO + NN potential plus N 2 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 3 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][75][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.