Wavefunction matching for solving quantum many-body problems

Ab initio calculations have an essential role in our fundamental understanding of quantum many-body systems across many subfields, from strongly correlated fermions1–3 to quantum chemistry4–6 and from atomic and molecular systems7–9 to nuclear physics10–14. One of the primary challenges is to perform accurate calculations for systems where the interactions may be complicated and difficult for the chosen computational method to handle. Here we address the problem by introducing an approach called wavefunction matching. Wavefunction matching transforms the interaction between particles so that the wavefunctions up to some finite range match that of an easily computable interaction. This allows for calculations of systems that would otherwise be impossible owing to problems such as Monte Carlo sign cancellations. We apply the method to lattice Monte Carlo simulations15,16 of light nuclei, medium-mass nuclei, neutron matter and nuclear matter. We use high-fidelity chiral effective field theory interactions17,18 and find good agreement with empirical data. These results are accompanied by insights on the nuclear interactions that may help to resolve long-standing challenges in accurately reproducing nuclear binding energies, charge radii and nuclear-matter saturation in ab initio calculations19,20.

Ab initio calculations play an essential role in our fundamental understanding of quantum many-body systems across many subfields, from strongly correlated fermions [1][2][3] to quantum chemistry [4][5][6] and from atomic and molecular systems [7][8][9] to nuclear physics [10][11][12][13][14].One of the primary challenges is to perform accurate calculations for systems where the interactions may be complicated and difficult for the chosen computational method to handle.Here we address the problem by introducing a new approach called wavefunction matching.Wavefunction matching transforms the interaction between particles so that the wavefunctions up to some finite range match that of an easily computable interaction.This allows for calculations of systems that would otherwise be impossible due to problems such as Monte Carlo sign cancellations.We apply the method to lattice Monte Carlo simulations [15,16] of light nuclei, medium-mass nuclei, neutron matter, and nuclear matter.We use high-fidelity chiral effective field theory interactions [17,18] and find good agreement with empirical data.These results are accompanied by new insights on the nuclear interactions that may help to resolve long-standing challenges in accurately reproducing nuclear binding energies, charge radii, and nuclear matter saturation in ab initio calculations [19,20].
Quantum Monte Carlo simulations are a powerful and efficient method that can describe strong correlations in quantum manybody systems.No assumptions about the nature of the system are necessary, and the computational effort grows only as a low power of the number of particles.For many problems of interest, a simple Hamiltonian H S can be found that describes the energies and other observables of the many-body system in fair agreement with empirical data and is easily computable using Monte Carlo methods.On the other hand, realistic high-fidelity Hamiltonians usually suffer from severe sign problems with positive and negative contributions to the averages cancelling each other, so that Monte Carlo calculations become impractical.In this work, we solve this problem using a new approach called wave function matching.While keeping the observable physics unchanged, wave function matching creates a new high-fidelity Hamiltonian H ′ such that wave functions at short distances match that of a simple Hamiltonian H S which is easily computed.This allows for a rapidly converging expansion in powers of the difference H ′ − H S .Wave function matching can be used with any computational scheme.In the following analysis, we focus on the case of quantum Monte Carlo simulations, where the method presents a promising and practical strategy for evading the sign problem in realistic calculations of nuclear quantum many-body systems.
A unitary transformation U is a linear transformation that maps normalized orthogonal states to other normalized orthogonal states.Starting from a realistic high-fidelity Hamiltonian H, wave function matching defines a new Hamiltonian H ′ = U † HU , where U † is the Hermitian conjugate of U .In each two-body angular momentum channel, the unitary transformation U is active only when the separation distance between two particles is less than some chosen distance R. Let us write ψ 0 (r), ψ ′ 0 (r), and ψ S 0 (r) for the ground state wave functions of H, H ′ , and the simple Hamiltonian H S , respectively.The transformation U is defined such that ψ ′ 0 (r) is proportional to ψ S 0 (r) for r < R. The simple Hamiltonian is chosen so that the constant of proportionality is close to 1.For r > R, however, U is not active and so ψ ′ 0 (r) remains equal ψ 0 (r).This is illustrated in the left panel of Fig. 1.The simple Hamiltonian H S is an easily computable Hamiltonian while the high-fidelity Hamiltonian H is not.A unitary transformation on the two-nucleon interaction with finite range R is used to produce a new Hamiltonian H ′ that is close to H S .In each two-body channel, the ground state wave function of H ′ matches the ground state wave function of H for r > R and is proportional to the ground state wave function of H S for r < R. Right Panel: The Tjon band correlation between the binding energies of 3 H and 4 He.The gray band is the predicted result from Ref. [21].The black open box shows the empirical point.The green-diamond, blue-circle, and read-square points show the results at LO, NLO, and N3LO in chiral effective field theory, respectively.The open points show the results from the first-order perturbative calculations using the Hamiltonian H, while the filled points are the results of the first-order perturbative calculations using the Hamiltonian H ′ .
Wave function matching will now be applied to ab initio Monte Carlo nuclear lattice simulations [15,[22][23][24][25] using the framework of chiral effective field theory (χEFT) [17,26].For our realistic Hamiltonian H, we use χEFT at next-to-next-tonext-to-leading order (N3LO) with lattice spacing a = 1.32 fm.For our simple Hamiltonian H S , we employ a χEFT interaction at leading order.Details of the interactions can be found in Methods.In the following, we use the term "local" for interactions that do not change the positions of particles while "nonlocal" refers to interactions that do change the relative positions of particles.The "range" of the interaction refers to the separation distance at which the interaction between particles becomes negligible.
We calculate all quantities up to first order in perturbation theory, which corresponds to one power in the difference H ′ − H S .As a first test, we consider the energy of the deuteron, 2 H.The wave function matching calculation gives a binding energy of 2.02 MeV, in comparison with 2.21 MeV for the true binding energy of H and 2.22 MeV for the experimentally observed value.The residual error of 0.1 MeV per nucleon is due to corrections beyond first order in powers of H ′ − H S .If one does not use wave function matching and instead performs the analogous calculation to first order in H − H S , the result is a much less accurate binding energy of 0.68 MeV.
As a second test of wave function matching, we calculate the binding energies of 3 H and 4 He.The Tjon band describes the universal correlations between the 3 H and 4 He binding energies [21,27].Provided that there are no long-range nonlocal interactions, any realistic two-nucleon interaction produces binding energies that lie on the Tjon band.The inclusion of any shortrange three-nucleon interaction also preserves this universal relation.In Fig. 1 we show wave function matching calculations using two-nucleon interactions only.At leading order (LO) the calculated point falls outside the Tjon band since the Coulomb interaction is not included, while the next-to-leading order (NLO) and N3LO results lie squarely in the middle of the band.We are using a low-energy scheme where the two-nucleon interaction is the same at NLO and next-to-next-to-leading order (NNLO).[28].The empirical point is also shown in Fig. 1.The good agreement with the Tjon band suggests a residual error of 0.1 MeV per nucleon or less.This can be compared with the substantial deviation from the Tjon line if one does not use wave function matching and performs the analogous calculation to first order in H − H S .
Before proceeding to larger nuclei and many-body systems, we first comment on the current status of ab initio calculations of nuclear structure using χEFT.There has been tremendous progress in the last few years towards producing accurate results for nuclear structure across much of the nuclear chart using a variety of different computational approaches [29][30][31][32][33][34][35][36][37][38].But there is also ample evidence that the calculations are sensitive to the manner in which the short-distance features of the interactions are regulated [39][40][41][42], a warning sign that systematic errors are not fully under control.Furthermore, no ab initio calculations have been able to provide an accurate description of nuclear binding energies, charge radii, and the saturation properties of symmetric nuclear matter with equal numbers of protons and neutrons.Our goal here is to identify the problem and point to a viable solution.
The results in Ref. [43,44] showed that the range and locality of the nuclear interactions have a strong influence on nuclear binding.The 4 He nucleus, also called an α particle, is a spatially compact object whose radius is comparable to the range of the interactions between nucleons.The central density of the α particle is about 50% greater than the saturation density of nuclear matter.As a result, the αα interaction is highly sensitive to the range and locality of the nucleonic interactions.These same arguments apply to other interactions involving α particles and nucleons, which we denote as N .Using the formalism of cluster effective field theory (also called halo effective field theory) [45][46][47][48] for α particles and nucleons, the two-cluster interactions are αα and αN , and the three-cluster interactions are ααα, ααN , and αN N , with αN N having two possible isospin channels.In addition to the αα interaction, there will be some dependence on short-distance physics arising from these five cluster interactions in their most attractive channels.While cluster effective field theory is not well suited for describing the structure of heavy nuclei and nuclear matter, it does provides a useful framework for categorizing systematic errors in nuclear many-body systems due to short-distance features of the interactions.
We will remove the unwanted dependence on short-distance features of the nucleonic interactions by fitting the binding energies of several light and medium-mass nuclei using three-nucleon interactions that go beyond N3LO order.In χEFT, three nucleon forces first appear at order NNLO.These include terms associated with the exchange of two pions and whose coefficients are determined from pion-nucleon scattering.There are also two interactions with singular short-distance properties that must be regulated and the corresponding couplings fitted to empirical data.As shown in Fig. 2, c D corresponds to the short-range interaction of two nucleons linked to a third nucleon through the exchange of a pion, and c E corresponds to the short-range interaction of all three nucleons.At order N3LO, there are additional terms associated with the exchange of two pions as well as readjustments of the c D and c E coefficients at N3LO order [49][50][51].The higher-order three-nucleon interactions we consider are short-range modifications to the c D and c E terms that are designed to be strongly correlated with the cluster interactions described above.In Methods we present the details of these new interactions along with a detailed uncertainty analysis.We find that with just one parameter, the root-mean-square-deviation (RMSD) for the energy per nucleon drops from about 2 MeV down to 0.4 MeV.With the addition of five additional parameters, the RMSD per nucleon drops further to 0.1 MeV.These results are consistent with the hypothesis that the αα interaction plays a key role in nuclear binding and that there are five additional cluster interactions that are sensitive to short distance physics.
In the right panel of Fig. 2, we present the results for the nuclear binding energies using wave function matching.We show ground state and excited state energies of selected nuclei with up to A = 40 nucleons and comparison with experimental data.The symbols with a black border indicate nuclei with unequal numbers of protons and neutrons.The nuclei used in the fit of the higher-order three-nucleon interactions are labelled with open squares, while the other nuclei are predictions denoted with filled diamonds.The one-standard-deviation error bars shown in Fig. 2 represent uncertainties due to computational uncertainties from Monte Carlo errors, infinite volume extrapolation, and infinite projection time extrapolation.We estimate the additional systematic errors due to truncation of the expansion of the H ′ − H S to be approximately 0.1 MeV per nucleon.The theoretical uncertainty due to the fitting of N3LO interactions have been empirically reduced by fitting the higher-order three-nucleon interactions to the binding energies of several nuclei.We estimate the residual uncertainty due to the interactions to be less than 0.1 MeV per nucleon.
In the left panel of Fig. 3, we present the results for the charge radii of nuclei with up to A = 40 nucleons.No charge radii data were used to fit any interaction parameters.The one-standard-deviation error bars shown in Fig. 3 [52], Auxiliary Field Diffusion MC calculations (GCR) [53], calculated with N3LO/NNLO (two-nucleon/three-nucleon) chiral interactions (EM 500 MeV, EGM 450/500 MeV and EGM 450/700 MeV) [54] and NNLO chiral interactions with explicit delta degrees of freedom (∆NNLO) [55].

represent computational
uncertainties due to Monte Carlo errors, infinite volume extrapolation, and infinite projection time extrapolation.The theoretical systematic uncertainties for the charge radii calculations are more difficult to predict without a comprehensive study of the dependence of charge radii upon the individual interactions.However, the agreement with empirical results is quite good, with an RMSD of about 0.03 fm.Note that the larger errors for the heaviest nuclei are statistical and can be decreased by utilizing greater computational resources.
In the right panel of Fig. 3, we present lattice results for the energy per nucleon versus density for pure neutron matter and symmetric nuclear matter.None of the neutron matter and symmetric nuclear matter data were used to fit any interaction parameters.The density is expressed as a fraction of the saturation density for nuclear matter, ρ 0 = 0.16 fm −3 .For the neutron matter calculations, we consider 14 to 80 neutrons in periodic box lengths ranging from 6.58 fm to 13.2 fm.For the symmetric nuclear matter calculations, we use system sizes from 12 to 160 nucleons in a periodic box of length 9.21 fm.The comparisons with several other published work are shown and detailed in the figure caption.We see that the neutron matter calculations are in good agreement with previous calculations, and the symmetric nuclear matter calculations pass through the empirical saturation point.The one-standard-deviation error bars represent computational uncertainties due to Monte Carlo errors and infinite projection time extrapolation.Uncertainties associated with extrapolation to large system size are not included here and will be explored in future work.These lattice simulations of symmetric nuclear matter are qualitatively different from other theoretical calculations that assume a homogeneous phase.The lattice simulations exhibit phase separation and cluster formation, just as in the real physical system.Due to the finite number of nucleons in these calculations, some oscillations due to nuclear shell effects can be seen in the energy per nucleon.
Another interesting feature of the lattice results is that symmetric nuclear matter without three-nucleon forces is underbound rather than overbound.This is opposite to what is found in other calculations using renormalization group (RG) methods to "soften" the nucleonic interactions [56][57][58].The term "soft" refers to the reduction of the strength of the interaction at large momentum values.The difference arises from the fact that wave function matching produces a two-nucleon interaction that is more nonlocal and also less "soft" when compared with two-nucleon interactions produced using RG flow equations.While there are some parallels in design, wave function matching is qualitatively different from RG methods.Wave function matching implements a unitary transformation that has finite range, and so the transformation can be viewed as redefining the χEFT twonucleon interaction.One does not need to explicitly compute the many-body forces induced by the unitary transformation and needs only to determine the new coefficients of the short-range χEFT three-nucleon and higher-nucleon forces.This has been demonstrated in a toy model [59].The approach used in wave function matching is similar to that of the Unitary Correlation Operator Method (UCOM) [60][61][62], though the operator-based unitary transformations in UCOM are significantly different in nature and purpose.
In summary, we have presented a new approach for solving quantum many-body systems called wave function matching.Wave function matching uses a transformation of the particle interactions to allow for calculations of systems that would otherwise be difficult or impossible.We have applied the method to lattice Monte Carlo simulations of light nuclei, medium-mass nuclei, neutron matter, and nuclear matter at N3LO in χEFT and found good agreement with empirical data.We have discussed the connection between systematic errors in nuclear many-body systems at short distances and interactions of α particles with themselves and nucleons.These new developments may help resolve long-standing challenges in accurately reproducing nuclear binding energies, charge radii, and nuclear matter saturation in ab initio calculations.While we have focused on Monte Carlo simulations for nuclear physics here, wave function matching can be used with any computational method applied to any quantum many-body system.This also includes quantum computing algorithms where wave function matching can be used to reduce the number of quantum gates required.All that is needed is a simple Hamiltonian H S that produces fair agreement with empirical data for the many-body system of interest and is easily computable using the method of choice.Further details on the implementation of wave function matching are given in Methods.
We are grateful for discussions with members and partners of the Nuclear Lattice Effective Field Theory Collaboration (Joaquín Drut, Gustav Jansen, Stefan Krieg, Zhengxue Ren, Avik Sarkar, and Qian Wang) as well as Scott Bogner, Andreas Ekström, Heiko Hergert, Morten Hjorth-Jensen, Daniel Phillips, Achim Schwenk, and Witek Nazarewicz.We

Hamiltonian Translators
Before describing how wave function matching is implemented in practice, we first discuss a class of transformations called Hamiltonian translators.Let H A and H B be two Hamiltonians acting on the same linear space.Suppose that U AB is a unitary transformation mapping eigenvectors of H B to the eigenvectors of H A .We then call U AB a Hamiltonian translator from AB is then a Hamiltonian translator from H A to H B .We note the curious fact that H ′ A = U BA H A U AB is a Hamiltonian with energy eigenvalues the same as H A but with eigenvectors corresponding to H B .Likewise, H ′ B = U AB H B U BA is a Hamiltonian with energy eigenvalues the same as H B but with the eigenvectors corresponding to H A .
Since H ′ A and H B share the same eigenvectors, H ′ A and H B commute with each other.In order to compute any energy eigenvalue of H ′ A , it suffices to prepare the corresponding eigenvector of H B and compute the energy expectation value of H ′ A .We can express these facts using the language of perturbation theory.If we write , then the zeroth-order expansion of the eigenvectors is exact and the first-order expansion of the energies is exact.
We can construct Hamiltonian translators using quantum adiabatic evolution [63].Let f (t) be a smooth function such that f (0) = 0 and f (1) = 1.Then for any T > 0, we can define the time-dependent Hamiltonian We also define the unitary transformation where ← − T is the time ordering symbol placing operators at later times on the left.In the limit of large T , U T is a Hamiltonian translator from H B to H A .In the limit of large T , U T maps every eigenvector of H B to an eigenvector of H A .Within each symmetry subspace that is invariant under both Hamiltonians H A and H B , the mapping U T preserves the ordering of energy eigenvalues.

Wave Function Matching
As discussed in the main text, we let H S be a simple Hamiltonian that is easily computable and H be a Hamiltonian with realistic interactions.Wave function matching can be viewed as an approximate Hamiltonian translator from H S to H.It is only an approximate translator because the unitary transformation U will be restricted to the space of two nucleons up to some maximum separation distance R and will only map the lowest eigenvector of H S to the lowest eigenvector of H in each scattering channel.
In the following, we restrict our focus to the space of two nucleons in some angular-momentum scattering channel.We impose hard wall boundary conditions at some very large separation distance R wall .The energy eigenstates of H for our chosen two-nucleon scattering channel will be denoted |ψ n ⟩ for n = 0, 1, • • • , and E n will be the corresponding energy eigenvalues.We let the corresponding energy eigenstates for H S be denoted |ψ S n ⟩.We now define a short-distance projection operator P R that projects out the portion of the two-nucleon state with separation distance less than or equal to R. We let |m⟩ for m = 1, • • • , m R be an orthogonal basis spanning the set of two-nucleon channel states at short distances so that P R = m R m=1 |m⟩ ⟨m|.In the lattice calculations presented here, we take |m⟩ to be radial position eigenstates, sorted according to increasing radial distance.See, for example, Ref. [64] for a discussion of radial position eigenstates on the lattice.For continuum space calculations, it is more convenient to work with smooth basis functions.For example, one can take |m⟩ to be energy eigenstates sorted in order of increasing energy for a steep harmonic oscillator potential.The projection operator P R = m R m=1 |m⟩ ⟨m| does not have strictly finite range, but nevertheless diminishes very rapidly as a function of distance.We take R to be the distance at which the interaction is negligible.
Let |ψ 0 ⟩ R and |ψ S 0 ⟩ R be the normalized short-distance portions of the ground states of H and H S respectively, Let us define a unitary transformation U such that |ψ S 0 ⟩ R is mapped to |ψ 0 ⟩ R .There are several ways to define the remaining properties of U .In this work, we use Gram-Schmidt orthogonalization to define an ordered sequence of orthonormal basis states, We require that U maps each basis vector |j⟩ S ⊥ to the corresponding basis vector |j⟩ ⊥ for each j We proceed by defining the transformed Hamiltonian H ′ = U † HU .Let |ψ ′ 0 ⟩ be the ground state of H ′ and let We note that |ψ ′ 0 ⟩ = U † |ψ 0 ⟩, ignoring any irrelevant overall phases.Since U is negligible at distances greater than R, |ψ ′ 0 ⟩ must equal |ψ 0 ⟩ at distances greater than R.At distances less than R, |ψ ′ 0 ⟩ R must equal |ψ S 0 ⟩ R .Hence the short distance behavior of |ψ ′ 0 ⟩ must be proportional to |ψ S 0 ⟩, with constant of proportionality κ.If |ψ 0 ⟩ and |ψ S 0 ⟩ have approximately equal normalizations at distances greater than R, then κ will be numerically close to 1, ignoring any irrelevant overall phases.
We note that wave function matching introduces some nonlocality to the two-nucleon interactions in H ′ .In principle, the process of wave function matching could increase the range of the interactions up to a distance scale equal to R. In practice, however, the range of interactions in H ′ is determined by the radial distances for which the difference between the wave functions |ψ 0 ⟩ and |ψ S 0 ⟩ is significant.For the lattice calculations presented here we take R to be 3.7 fm, while the difference between |ψ 0 ⟩ and |ψ S 0 ⟩ is appreciable at distances less than 1.9 fm.

Simple Hamiltonian
We now present the details of our simple Hamiltonian H S .We construct the Hamiltonian using a χEFT interaction at leading order, where K is the kinetic energy term with nucleon mass m = 938.92MeV, the :: symbol indicates normal ordering, and ρ(d) and ρ(d) I are density operators that are smeared both locally and non-locally, The smeared annihilation and creation operators, ã and ã † , have with spin i = 0, 1 (up, down) and isospin j = 0, 1 (proton, neutron) indices, Through our calculations we use local smearing parameter s L = 0.07 and nonlocal smearing parameter s NL = 0.5.In addition to the short-range SU(4) symmetric interaction, we also have a long-range one-pion-exchange (OPE) potential at leading order χEFT interaction.We define our one-pion-exchange potential following a recently developed regularization method [65], Here f π is a local regulator in momentum space defined as f S ′ ,S is the locally-regulated pion correlation function, and with g A = 1.287 the axial-vector coupling constant, f π = 92.2MeV the pion decay constant and M π = 134.98 MeV the pion mass.The term given in Eq. ( 9) is a counterterm introduced to remove the short-distance admixture in the one-pionexchange potential [65].In our simple Hamiltonian, we set Λ π = 180 MeV and C π = 0, and we compute the difference and the OPEP counterterm V Λπ Cπ perturbatively.Here we use the notation and for the density operators.
Hamiltonian at N3LO (next-to-next-to-next-to-leading order) We now give the details of our realistic Hamiltonian H. Let us define the functions and We define the Hamiltonian H using χEFT interactions at N3LO, where K is the kinetic energy term.V Λπ OPE and V Λπ Cπ are defined in Eqs. ( 8) and ( 9) with 2N,WFM is the wave function matching interaction defined as H ′ − H, and 2N,WFM is the GIR correction of the wave function matching interaction.
For the details of the Coulomb interaction and the 2N short-range interactions we refer the reader to Ref. [28].The threenucleon interactions at Q 3 consists of a contact potential, one-pion exchange potential, and two-pion exchange potential [66][67][68], and in this work we defined two additional SU(4) symmetric potentials denoted by V (l) c E and V (t) c E .Therefore, the three-nucleon interactions at Q 3 has the form where The V (TPE) 3N potential can be separated into the following three parts, We perform our calculations using lattice spacing a = 1.32 fm, and we determine the low-energy constants (LECs) of the 2N short-range interaction up to N3LO of χEFT by reproducing the neutron-proton scattering phase shifts and mixing angles of the Nijmegen partial wave analysis (PWA) [69].The lattice spacing of a = 1.32 fm corresponds to the momentum space cut-off of 470 MeV, which corresponds to the resolution scale, at which the hidden spin-isospin symmetry of the NN interactions is best fulfilled [70].In Fig. 4 we plot the calculated neutron-proton scattering phase shifts up to N3LO of χEFT as functions of relative momenta with comparison to the Nijmegen PWA.Only the statistical errors and not systematic errors are included in the Nijmegen PWA.In the next subsection we discuss the approach used to estimate the uncertainties in our calculations and the determination of the LECs of the three-nucleon interactions.

Uncertainty Analysis
We discuss the analysis and quantify the errors which include relevant sources of uncertainty in the calculations presented in this paper.To perform such a complete analysis to estimate the uncertainties for the calculations using chiral interactions, a global parameter search for the LECs of the chiral interactions is required.This task could be extremely difficult as one has to perform a search over a high-dimensional parameter space.Nevertheless, by performing some prior analyses we can reduce both the dimension and the volume of the parameter space to be searched.For instance, the LECs of the two-pion exchange three-body potentials are already fixed from pion-nucleon scattering data, c 1 = −1.10(3),c 3 = −5.54(6) and c 4 = 4.17 (4) all in GeV −1 [71], and the LECs of two-nucleon potentials can be constrained using the empirical partial wave phase shifts and mixing angles, which help to shrink the parameter space.However, the main difficulty is the determination of the unknown  FIG. 4. Plots of the neutron-proton scattering phase shifts as functions of relative momenta with lattice spacing a = 1.32 fm.For comparison we also plot the phase shifts extracted from the Nijmegen partial wave analysis [69].
LECs of the three-nucleon forces.As we discuss in the following, we treat these LECs as unknown regression coefficients of an emulator employed with history matching, which has been shown to be an effective approach for parameter searches [72][73][74].
In our calculations we determine the LECs of the 2N chiral interactions by fitting the calculated neutron-proton scattering phase shifts and mixing angles on the lattice to the Nijmegen PWA [69].To estimate the systematic theoretical uncertainties, we use the method introduced in Refs.[75,76] and calculate the theoretical errors for the neutron-proton scattering phase shifts and mixing angles as a function of the relative momenta, and the results are shown in Fig. 5.One can find detailed discussions about the convergence of the chiral effective field theory expansion on the lattice and theoretical errors in Ref. [28].As seen in Fig. 5 the theoretical error bands at N3LO give a considerable-sized parameter space of LECs.Therefore, we use history matching together with an emulator to filter the 2N LECs as well as to identify the 3N LECs which provide an acceptable matching between ab initio calculations and the experimental data for some selected nuclei.
The first step in our analysis is to use the theoretical errors shown in Fig. 5 and get posterior distributions of the 2N LECs sampling with Markov Chain Monte Carlo (MCMC) and the results are shown in Fig. 6.For each interaction, we approximately collect 10 3 samples.Now we can use 2N LECs from these distributions as inputs into our ab-initio model, f (x), for which we define the emulator equation given by The first term on the right hand side is the model output of only 2N interactions using 2N inputs, which is denoted by the vector x, from the distributions shown in Fig. 6.The second term on the right hand side is the output of only 3N interactions and it is a calibration term where β k are unknown regression coefficients which will correspond to the 3N LECs at the end of the analysis.
For this term we consider five different modified forms of the SU(4) symmetric three-nucleon forces and three different modified forms of the one-pion-exchange three-nucleon forces at short distances, V c D and V c D .In the emulator Eq. ( 26), the derivatives of the energies with respect to x and β stand for the calculated operators on the lattice using first-order perturbation theory.Also, in the emulator equation all inputs are active variables, therefore, we do not make any such distinction.The emulator works for any number of outputs, f o (x) with o = 1, 2, . . ., q.When we emulate our model's behavior for x i inputs, with i = 1, • • • , d, we use the experimental data corresponding to the outputs and the calibration coefficients β k so that there is acceptable agreement between our model and the experimental data.Therefore, the calibration coefficients β k are obtained from a nonlinear least-squares fitting by minimizing the expression for each input set in the vector x, where z n is the experimental data.When we evaluate Eq. ( 27), we impose a constraint on the least square problem so that only natural sized parameters are allowed for the coefficient β k .
In the second step, we use Eq. ( 27) to link our model to reality by evaluating inputs from the distributions shown in Fig. 6 and determining the 3N LECs.To achieve this, we use history matching to perform an efficient parameter search.History matching is an iterative process that identifies and eliminates implausible parts of the input space by measuring implausibility of inputs, shrinks the input space in every iteration, and repeats the non-implausible input search in the smaller input space.Here we use Q w to denote the volume of non-implausible input space.The implausibility measure of a given d inputs x = (x 1 , x 2 , • • • , x d ) for an n th output is written as, where the numerator is the distance between the experimental data z n and the calculated output f n (x), and the denominator is the square root of the sum of the variance of the distance r n (x) and the variance of the errors of prediction ϵ n .Any particular value of x, say x i , that give a large value of I n (x) can be considered as implausible which means that for this input it is unlikely to get an acceptable match between the calculated outputs and the experimental data.Therefore, we can use a maximum implausibility measure I M (x) and discard the input from the input space provided that I M (x) > c where c is a cutoff for non-implausibility.
The plot of the posterior distribution of the 2N LECs at N3LO sampled using MCMC method.
In our analysis we define the maximum implausibility measure by maximizing over all outputs, where Z w is the set of outputs that we consider in the w th iteration, known as wave, and we use c = 3 by appealing to Pukelsheim's 3-sigma rule [77] which indicates that 95% of the all population lies within ±3σ.Even though the maximum implausibility measure is defined by Eq. ( 29), in our analysis we look at the first three maximum implausibilities in every iteration, and it provides us a tool to select three-nucleon forces which work as a calibrator in Eq. ( 28).This is a crucial step both as it gives the information about three-nucleon forces which are inevitable to obtain an acceptable agreement between the calculated outputs and the experimental data with some reasonable credibility and as it guides us to not discard any useful information in the input space in early iterations.
For the lattice Monte Carlo simulations, the general strategy is to use the expansion H ′ = H S + (H ′ − H S ) and apply corrections up to first order in perturbation theory.In order to accelerate the convergence of perturbation theory further, we find it helpful to consider the more general partition H ′ = H ′S + (H ′ − H ′S ).The modified simple Hamiltonian H ′S has the same form as H S in Eq. ( 4), but we allow for different coupling strengths c SU(4) and c I .We then minimize the energy and use the variational principle to optimize the parameters.
We refer the readers to Ref. [78] for a general description of and a recipe for history matching and to Ref. [74] for a nuclear physics application.In the following we discuss history matching method and outline the detailed procedure of its application to our model.
Our choice for the collection of inputs in the volume Q 0 is based on the fact that only interactions with zero angular momentum are giving significant contributions to the nuclei used as informative observables for this first wave.Therefore, we define the current sets of non-implausible inputs x from the distributions shown in Fig. 6, and we calculate the model outputs from Eq. ( 26) by performing a calibration using least square regression over the second term on the right hand side with the simplest 3N interactions V (0) Then, the implausibility measures I M (x) are calculated using the first maximum implausibility over outputs, and implausibility cutoffs are imposed to define the new non-implausible volume denoted by Q 1 .The results of the fitted (predicted) observables are shown by black open (filled) down-triangle points in Fig. 7 where we plots the deviations between the model outputs and experimental data given in percentage.From Fig. 7 it can be seen clearly that when we run the emulator for outputs Z 1 , deviations for all predicted outputs are negative, which is important information for choosing outputs to be included in the emulator in the next iteration.WAVE 2: At this wave, we run the implausibility measure over outputs Z 2 = Z 1 { 8 Be, 12 C, 16 O} by considering the fact that there exist strong evidences that the structures of these additional nuclei are in the form of distinct geometrical configuration of alpha clusters, then we evaluate the inputs for all 2N interactions.The current non-implausible parameter space for the 1 S 0 and 3 S 1 − 3 D 1 channels are set by WAVE 1, while for the channels that have not been evaluated yet the current sets of non-implausible parameter space is defined from the distributions shown in Fig. 6.At this wave we perform the least square regression given in Eq. ( 27) by using two additional 3N interactions.To make selection for these additional 3N interactions we run the emulator for all possible sets of 3N interactions and compute the first three non-zero maximum implausibility measures.Sets of 3N interaction used as calibrators and the corresponding results for the ratio of the non-implausible space volumes, Q 2 /Q 1 , are shown in Table I.
As seen from the results, in the most cases one of the calculated outputs has a large deviation with small variance, and as a result the first maximum implausibility measures are less than 1%, and in the cases where two additional 3N interactions are c E we obtain the largest non-implausible space volume.Therefore, from Table I we find the most promising set of 3N interaction as c E shown in the first row of Table I.Therefore, we define the new nonimplausible space volume Q 2 at this iteration.

WAVE 3:
The results of the fitted (predicted) observables at WAVE 2 are shown by cyan open (filled) up-triangle points in Fig. 7.We find that using four 3N interactions in WAVE 3 already gives a decent description for nuclei with N = Z = 8, while deviations for neutron-rich nuclei with Z ≤ 8 are large.Therefore, at this wave we run the the implausibility measure over outputs Z 3 = Z 2 { 13 C, 14 C, 17 O, 18 O} and we evaluate the inputs for all 2N interactions.The current non-implausible input volume is set by WAVE 2, Q 2 .Also, we run Eq.( 27) by using two more 3N interactions in addition to c E .To select two additional 3N interactions out of four we run the emulator for all possible sets of 3N interactions and compute the first three non-zero maximum implausibility measures.Sets of 3N interaction used as calibrators and the corresponding results for the ratio of the non-implausible space volumes, Q 3 /Q 2 , are shown in Table II.

Reading the results from Table II we determine the most promising set of 3N interaction as
c E , and we define new non-implausible volume for inputs denoted as Q 3 .TABLE I. Results for the first three non-zero maximum implausibility measures at WAVE 2. We run the implausibility measure over outputs 3 H, 4 He, 8 Be, 12 C, and 16 O.

3N interactions
ratio of non-implausible input volumes Q2/Q1 (%) in addition to 3N interactions ratio of non-implausible input volumes Q3/Q2 (%) in addition to M (x) I (3) WAVE 4: We show the results of the fitted (predicted) observables at WAVE 3 in Fig. 7 by green open (filled) right-triangle points.At this wave we use the same 3N interactions used in WAVE 3 we run the the implausibility measure over outputs Z 4 = Z 3 { 7 Li, 9 Be, 10 Be, 10 B, 11 B, 14 N, 15 N} where we included odd-odd nuclei with N = Z as well as neutron-rich nuclei with Z ≤ 8. Therefore, we define new non-implausible volume for inputs denoted as Q 4 .The results of the fitted (predicted) observables are shown by blue open (filled) circle points in Fig. 7.

WAVE 5:
The results in Fig. 7 show that the agreements between the calculated and experimental data for light nuclei are already in an acceptable level.Therefore, by taking into consideration the level of agreement for medium-mass nuclei and the fact that we are left with two more additional 3N interactions to be used as a calibrator, at this iteration we perform a numerical experiment to determine whether these additional 3N interactions have any considerable influence on improving the agreement.To stabilize the calculated results for medium-mass nuclei at this wave we run the implausibility measure over outputs Z 5 = Z 4 { 40 Ca} and run the emulator for all possible sets of 3N interactions to compute the first three non-zero maximum implausibility measures.The results are shown in Table III.
Based on the numerical results given in Table III, at this iteration we run the implausibility measure over outputs Z 5 using the 3N interactions After performing five waves described above we discard implausible regions from the parameter spaces of the 2N LECs and we determine the minimum number of 3N interaction needed to get a descent agreement between the calculated and experimental data for the nuclear binding energies of nuclei with A ≤ 40.Finally we use the final parameter spaces and perform bootstrap re-sampling of the calculated results to quantify the systematic theoretical uncertainties due to the 2N LECs.TABLE III.Results for the first three non-zero maximum implausibility measures at WAVE 5. We run the implausibility measure over outputs 3 H, 4 He, 7 Li, 8 Be, 9 Be, 10 Be, 10 B, 11 B, 12 C, 13 C, 14 C, 14 N, 15 N, 16 O, 17 O, 18 O, and 40 Ca.

3N interactions
ratio of non-implausible input volumes Q5/Q4 (%) in addition to In the next step, using the calculated results with quantified uncertainties we perform a second analysis for the 3N interactions to show the importance of regulators and the connection with cluster EFT as well as nuclear binding energies for nuclei as discussed in the main text.To end this, we use the central values of the 2N LECs, which give the scattering phase shifts shown by blue dashed-lines in Fig. 5, with the quantified uncertainties, then we fine-tune the 3N interactions using some selected nuclear binding energies (to be explained below).In this analysis we calculate the RMSD (root mean square deviation) over all calculated nuclear binding energies shown in Fig. 7, and we use these results to assess the extent to which set of 3N interactions is the best for accurately describing some certain nuclear binding energies.The results are shown in Table IV-X where we also present the energy differences E 4  8 − 2E .From the results shown in Tables I, II and III, it can be seen that nuclear binding energies are highly sensitive to the range of the three-nucleon interactions.The first few rows of Table I clearly indicates that the nearest-neighbor smearing and the next-to-nearest-neighbor smearing are playing significantly different roles in describing the light and α-like nuclei.Therefore, we start the RMSD analysis using the simpliest two 3N interactions, V   As seen in Table IV, the RMSD of the binding energies per nucleon over all calculated nuclear binding energies is 1.236, and it corresponds to 17.3 % of the average nuclear binding energy per nucleon over all calculated nuclear binding energies 01, where N is the number of nuclei, z i is the binding energy of the i th nucleus and A i is the number of nucleon of the i th nucleus.Now we use one more 3N interaction in addition to c D and fit the nuclear binding energies for 4 He, 8 Be, 12 C, and 16 O.The results are given in Table V.We compute the RMSD over all calculated nuclear binding energies, and we find the lowest value as 0.403 MeV which corresponds to 5.6 % of the average nuclear binding energy per nucleon and is obtained using the 3N interactions In the next analysis, we use two additional 3N interactions and fit the nuclear binding energies for 4 He, 8 Be, 12 C, and 16 O.We compute the RMSD over all calculated nuclear binding energies and the results are given in Table VI.By using c E we find that the lowest value for the RMSD is 0.293 MeV, and this means that roughly we obtain the average nuclear binding energy per nucleon of all calculated nuclei with 4.1 % errors.The RMSD analysis and the results shown in Table VI are consistent with results of history matching given in Table I and cyan-colored up-triangle points in FIG. 7.
After finding four 3N interactions which give a good description for light α-like nuclei, we include neutron-rich nuclei in RMSD analyses.We use three 3N interactions in addition to V (0) c E and V (0) c D and fit the nuclear binding energies for 3 H, 4 He, 7 Li, 8 Be, 9 Be, 10 Be, 10 B, 11 B, 12 C, 13 C, 14 C, 14 N, 15 N, 16 O, 17 O, and 18 O.We compute the RMSD over all calculated nuclear binding energies and the results are given in Table VII.We find that the lowest value for the RMSD is 0.131 MeV and is obtained by using c E , which means that the average nuclear binding energy per nucleon of  IX, respectively.The lowest value for the RMSD is found as 0.109(0.102)MeV when six(seven) 3N interactions are used in the fit, and this value corresponds to 1.54 %(1.43%) of the average nuclear binding energy per nucleon.
Finally, we fit the nuclear binding energies for 3 H, 4 He, 7 Li, 8 Be, 9 Be, 10 Be, 10 B, 11 B, 12 C, 13 C, 14 C, 14 N, 15 N, 16 O, 17 O, 18 O, and 40 Ca using six 3N interactions in addition to V  X.The lowest value for the RMSD is found as 0.079 MeV, which corresponds to 1.11 % of the average nuclear binding energy per nucleon.

Results
In this subsection, we give numerical details of the results presented in the main text.We have calculated results for the ground state and excited state energies as well as the charge radii of some selected nuclei up to A = 40, pure neutron matter, and symmetric nuclear matter at N3LO in chiral effective field theory using wave function matching.In Table XI we present the numerical details of the results for the calculated ground state and excited state energies and charge radii as well as experimental data for comparison.The quoted errors includes all the statistic and systematic uncertainties such as computational uncertainties from Monte Carlo errors, infinite volume extrapolation, infinite projection time extrapolation, truncation error of the expansion of the H ′ − H S and the systematic uncertainty due to the fitting the 2N short-range interactions.The errors on the charge radii grow for the heaviest nuclei, and the main sources of these errors are Monte Carlo statistics.Therefore, they can be reduced by using larger computer resources.In Table XII we show the calculated pure neutron matter energies for various numbers of neutrons and box sizes as well as corresponding densities.Similarly, large scale calculations corresponding to high densities have larger errors which can be reduced by larger computer resources.In Table XIII we show the results for symmetric nuclear matter energies for various numbers of nucleons as well as corresponding densities.The large errors on the results with large numbers of nucleon are the statistical errors of Monte Carlo calculations and can be reduced by more computer power.

Nuclear Lattice Effective Field Theory
Nuclear lattice effective field theory (NLEFT) [15,22] combines the frameworks of chiral effective field theory (χEFT) for the forces between nucleons, lattice field theory, and stochastic Monte Carlo algorithms.Each of these forms a cornerstone of the modern approach to strongly interacting many-fermion systems in many fields of study, notably nuclear, particle, condensed matter and atomic physics.By a unique merger of these cornerstones, NLEFT has matured into a leading framework for the investigation of nuclear structure properties.Building upon the early developments [79], NLEFT has extended its reach to light and medium-mass nuclei, mostly of eveneven type.This has enabled detailed predictions [80][81][82] of nuclear structure properties, in particular the binding energy for 3 H and 4 He at NNLO [83], as well as binding energies for isotopic chains with Z = 1, 2, 4, 6 and 8 [84], together with root-meansquare radii and proton and neutron densities.Moreover, NLEFT has given further impulse to the investigation of α-clustering in nuclear matter.The analyses of the Hoyle state [85,86] and further low-lying excited states [25,86] of 12 C, as well as the study of the structure and EM properties of the 0 + 1 , the 0 + 2 and 2 + 1 states of 16 O [87], and the prediction of ground-state properties of even-even self-conjugate nuclei with Z ≤ 14 [88], are noteworthy in this respect.
An immediate strength of NLEFT is the favorable computational scaling with nucleon number A, as the A-body Hamiltonian in NLEFT is diagonalized stochastically by means of an auxiliary field Monte Carlo (AFMC) algorithm.The two-body, three-body, and pion exchange interactions are described by auxiliary fields, which are sampled at each step in the AFMC algorithm.Energy eigenvalues are obtained using the adiabatic projection method, whereby transition amplitudes mediated by the Schrödinger time evolution operator between trial (Slater-determinant) states are constructed.In practical AFMC calculations of NLEFT, the temporal evolution of the trial states is discretized into L t equal steps, separated by a temporal lattice spacing a t .Hence, the time-evolution operator is subdivided into a sequence of L t transfer matrices (Trotter decomposition).
On the one hand, NLEFT is firmly grounded in low-energy QCD, as the interactions between nucleons are described by the NNLO or N3LO Lagrangians of χEFT.The resulting Hamiltonian is formulated with a finite cubic lattice regulator, and avoids the need for pre-diagonalization techniques such as the similarity renormalization group (SRG) [89].On the other hand, it is known that the AFMC implementation of the full NNLO or N3LO χEFT Hamiltonian produces a severe sign problem, which eventually compromises the effectiveness of the numerical technique.
A firmly-rooted remedy in the NLEFT literature [15] entails the replacement of the χEFT transfer matrices in the initial and the final L to time steps by pionless SU (4)-symmetric transfer matrices.It follows that the full χEFT Hamiltonian acts only in the middle L ti ≡ L t − L to time steps through the transfer matrices, whereas Wigner's SU(4) one acts as a low-energy filter [15] at the boundaries of the time interval.This procedure extends to the expectation value of any operator representing a physical observable, inserted in the midpoint of the chain of transfer-matrices.The positiveness of the transfer matrices, where Wigner's SU(4) action is used, is the major advantage of the method, and is guaranteed for systems with even number of nucleons and either spin-singlet or isospin-singlet quantum numbers [15].However, distortions in the energy eigenvalues as a result of the usage of an isospin-preserving SU(4) Lagrangian in the action do appear in the form of lower bounds [15].Although remedies such as eigenvector continuation [90] or symmetry-sign extrapolation [91] would allow for a broader use of the χEFT transfer matrix, in this work we tackle the problem by using the method of wave function matching.
c E and V (0) c D to fit the nuclear binding energies for 3 H, 4 He, 7 Li, 8 Be, 9 Be, 10 Be, 10 B, 11 B, 12 C, 13 C, 14 C, 14 N, 15   Other ab initio methods Over the past two decades, nuclear ab initio methods have had great success.The widespread adoption of nuclear forces from Effective Field Theory (EFT) [17,18,[92][93][94][95][96][97] laid the foundation for these developments.More specifically, the inclusion of chiral two-nucleon to higher orders [17,18,95,98] and the leading three-nucleon [36,94,99,100] forces have shown to be key elements for high quality calculations.The strong short-range repulsion in nuclear potentials is often soften by a renormalization procedure [58].Besides Nuclear Lattice Effective Field Theory (NLEFT), there are many other promising ab initio approaches being used to calculate the properties of few-and many-nucleon systems.Full configuration methods like no-core shell model (NCSM) [101][102][103][104], symmetry-adapted NCSM [105,106] and quantum Monte Carlo (QMC) in several different varieties [107][108][109][110][111] were able to extend their reach into the lower sd-shell.

FIG. 1 .
FIG.1.Left Panel: Pictoral representation of wave function matching.The simple Hamiltonian H S is an easily computable Hamiltonian while the high-fidelity Hamiltonian H is not.A unitary transformation on the two-nucleon interaction with finite range R is used to produce a new Hamiltonian H ′ that is close to H S .In each two-body channel, the ground state wave function of H ′ matches the ground state wave function of H for r > R and is proportional to the ground state wave function of H S for r < R. Right Panel: The Tjon band correlation between the binding energies of 3 H and 4 He.The gray band is the predicted result from Ref.[21].The black open box shows the empirical point.The green-diamond, blue-circle, and read-square points show the results at LO, NLO, and N3LO in chiral effective field theory, respectively.The open points show the results from the first-order perturbative calculations using the Hamiltonian H, while the filled points are the results of the first-order perturbative calculations using the Hamiltonian H ′ .

FIG. 2 .
FIG. 2. Left Panel: Short-range three nucleon forces at NNLO.The first is the one-pion exchange term cD shown on the left.The other is the purely short-range term cE shown on the right.At order N3LO there are additional three-nucleon interactions associated with the exchange of two pions, as well as the corrections from the renormalization of the cD and cE terms at N3LO order.Right Panel: Results for nuclear binding energies using wave function matching.Calculated ground state and excited state energies of some selected nuclei with up to A = 40 at N3LO in chiral effective field theory and comparison with experimental data.The symbols with a black border indicate nuclei with unequal numbers of protons and neutrons.The nuclei used in the fit of the higher-order three-nucleon interactions are labelled with open squares, while the other nuclei are predictions denoted with filled diamonds.
acknowledge funding by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) and the NSFC through the funds provided to the Sino-German Collaborative Research Center TRR110 "Symmetries and the Emergence of Structure in QCD" (DFG Project ID 196253076 -TRR 110, NSFC Grant No. 12070131001), the Chinese Academy of Sciences (CAS) President's International Fellowship Initiative (PIFI) (Grant No. 2018DM0034), Volkswagen Stiftung (Grant No. 93562), the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (ERC AdG EXOTIC, grant agreement No. 101018170, and ERC AdG NuclearTheory, grant agreement No. 426661267), the Scientific and Technological Research Council of Turkey (TUBITAK project no.120F341), the National Natural Science Foundation of China (Grants No. 12105106 and No. 12275259), U.S. National Science Foundation (PHY-1913620 and PHY-2209184), U.S. Department of Energy (DE-SC0013365 and DE-SC0021152), the Nuclear Computational Low-Energy Initiative (NUCLEI) SciDAC-4 project (DE-SC0018083), the Rare Isotope Science Project of the Institute for Basic Science funded by the Ministry of Science and ICT (MSICT) and by the National Research Foundation of Korea (2013M7A1A1075764), and the Institute for Basic Science (IBS-R031-D1-2022-a00).Computational resources provided by the Gauss Centre for Supercomputing e.V. (www.gauss-centre.eu) for computing time on the GCS Supercomputer JUWELS at Jülich Supercomputing Centre (JSC) and special GPU time allocated on JURECA-DC as well as the Oak Ridge Leadership Computing Facility through the INCITE award "Ab-initio nuclear structure and nuclear reactions", and partially provided by TUBITAK ULAKBIM High Performance and Grid Computing Center (TRUBA resources).vuComputational resources were also partly provided by the National Supercomputing Center of Korea with supercomputing resources including technical support ( KSC-2021-CRE-0429 , KSC-2022-CHA-0003).The work of Y.H.S. and Y.K. was supported by the Rare Isotope Science Project of Institute for Basic Science, funded by Ministry of Science and ICT (MSICT) and by National Research Foundation of Korea (2013M7A1A1075764).M.K. is supported by the Institute for Basic Science (IBS-R031-D1-2022-a00).

FIG. 5 .
FIG.5.The plot of the theoretical error bands for the neutron-proton scattering phase shifts and mixing angles versus the relative momenta.

c E and V ( 2 )
c D as calibrators.The results of the fitted (predicted) observables are shown by red open (filled) square points in Fig. 7.

c
E and V (0) c D , and the nuclear binding energies for light nuclei whose structure are in the form of distinct geometrical configuration of α clusters.The results are shown in TableIV.
c D , and the results are shown in Table For the pure neutron matter energy we use the number of neutrons from 14 to 80 and various box sizes from 6.58 fm to 13.2 fm, and for the symmetric nuclear matter energy we use the number of nucleons from 12 to 160 and a periodic box of length 9.21 fm.For comparison we show the results from variational (APR)

3
He 4 He 6 He 6 Li 7 Li 8 Be 9 Be 10 Be 10 B 11 B 12 C 0 + 13 C 14 C 14 N 15 N 16 O 17 O 18 O 20 O 22 O 24 O 18 F 20 Ne 24 Mg 28 Si 32 S 36 Ar 40 Ca 3H FIG.7.The plots of deviations between the model outputs and experimental data per nucleon at history matching iterations from WAVE 1 to WAVE 5.

TABLE II .
Results for the first three non-zero maximum implausibility measures at WAVE 3. We run the implausibility measure over outputs 3 H,4He,8Be, 12 C, 13 C, 14 C,16O,17O, and 18 O.

TABLE IV .
Results for the RMSD (root mean square deviation) over all calculated nuclear binding energies by using the simplest 3N interactions V D to fit the nuclear binding energies for 4 He,8Be,12C, and16O.All energies are measured in MeV.

TABLE V .
Results for the RMSD over all calculated nuclear binding energies by using one additional 3N interaction as well as V to fit the nuclear binding energies for 4 He, 8 Be,12C, and16O.All energies are measured in MeV.all calculated nuclei is computed roughly with 1.83 % errors.The RMSD analysis and the results shown in TableVII.Now we fit the nuclear binding energies for 3 H, 4 He,7Li, 8 Be,9Be, 10 Be, 10 B, 11 B, 12 C, 13 C, 14 C, 14 N, 15 N,16O,17O, and 18 O using four and five 3N interactions in addition to V , and the results are shown in TableVIII and Table

TABLE VI .
Results for the RMSD over all calculated nuclear binding energies by using two additional 3N interactions as well as V D to fit the nuclear binding energies for 4 He, 8 Be,12C, and16O.All energies are measured in MeV. c

TABLE VII .
Results for the RMSD over all calculated nuclear binding energies by using three additional 3N interactions as well as V D to fit the nuclear binding energies for 3 H, 4 He,7Li, 8 Be,9Be, 10 Be, 10 B, 11 B, 12 C, 13 C, 14 C, 14 N, 15 N,16O,17O, and 18 O.All energies are measured in MeV. c

TABLE VIII .
Results for the RMSD over all calculated nuclear binding energies by using four additional 3N interactions as well as V D to fit the nuclear binding energies for 3 H, 4 He,7Li, 8 Be,9Be, 10 Be, 10 B, 11 B, 12 C, 13 C, 14 C, 14 N, 15 N,16O,17O, and 18 O.All energies are measured in MeV. c

TABLE IX .
Results for the RMSD over all calculated nuclear binding energies by using five additional 3N interactions as well as V to fit the nuclear binding energies for 3 H, 4 He,7Li, 8 Be,9Be, 10 Be, 10 B, 11 B, 12 C, 13 C, 14 C, 14 N, 15 N,16O,17O, and 18 O}.All energies are measured in MeV.

TABLE XII .
Results of calculated pure neutron matter energy using various number of neutrons and box sizes.

TABLE XIII .
Results of calculated symmetric nuclear matter energy using various number of nucleons.