Emulating ab initio computations of infinite nucleonic matter

We construct efficient emulators for the ab initio computation of the infinite nuclear matter equation of state. These emulators are based on the subspace-projected coupled-cluster method for which we here develop a new algorithm called small-batch voting to eliminate spurious states that might appear when emulating quantum many-body methods based on a non-Hermitian Hamiltonian. The efficiency and accuracy of these emulators facilitate a rigorous statistical analysis within which we explore nuclear matter predictions for > 10 6 different parametrizations of a chiral interaction model with explicit ∆-isobars at next-to-next-to leading order. Constrained by nucleon-nucleon scattering phase shifts and bound-state observables of light nuclei up to 4 He, we use history matching to identify non-implausible domains for the low-energy coupling constants of the chiral interaction. Within these domains we perform a Bayesian analysis using sampling/importance resampling with different likelihood calibrations and study correlations between interaction parameters, calibration observables in light nuclei, and nuclear matter saturation properties.


I. INTRODUCTION
Infinite nuclear matter is an idealized system of strongly interacting nucleons that holds translational invariance without Coulomb and surface effects.Studies of its equation-of-state (EOS) at nuclear densities allow to explore properties of the microscopic interactions between constituent nucleons and to understand bulk properties of finite nuclei.It also provides an important anchor point for extrapolations to higher densities which is needed for the description of neutron stars and their mergers [1][2][3][4][5][6][7].While previous theoretical studies are mostly based on fixed parametrizations of certain interaction models [8][9][10][11][12][13][14][15][16][17][18][19][20][21][22][23] a rigorous statistical analysis will require to incorporate all relevant sources of uncertainty including the parametric one.In this work we introduce several key developments that are needed to perform such an analysis and we study in particular the multi-dimensional parameter domain of ∆-full chiral effective field theory (χEFT) at next-to-next-to-leading order (NNLO).The sensitivity of nuclear matter predictions to the calibration observables is explored in a companion paper [24].
Recent calculations have shown that an accurate description of bulk properties of finite nuclei and nuclear matter involves fine-tuning of the underlying nuclear Hamiltonian [11,20,21,50,52,53].The fact that these nuclear Hamiltonians give similar results for scattering phase shifts and few-body observables, but differ for many-body systems.indicates the challenge we face when constructing nuclear forces.In principle, microscopic approaches based on nucleonic degrees of freedom and chiral interactions should be able to describe both the two-and few-body sectors, as well as infinite nuclear matter.It is not clear whether going to higher orders in the effective field theory (EFT) expansion will resolve the fine-tuning, or whether the free parameters of the chiral Hamiltonian, i.e., the low-energy constants (LECs), are not sufficiently constrained by standard fitting observables in the few-nucleon sector.Therefore, the construction of nuclear forces to meet the required pre-cision and accuracy is still an ongoing and complex task.Systematic studies, such as the present work, are needed to better understand how details of the interaction influence properties of nucleonic matter, how predictions for different observables may be correlated, and also to rigorously quantify uncertainties.
Although significant progress has been made towards quantifying uncertainties of ab initio nuclear matter predictions and identifying possible correlations with observables in finite nuclei [11, 21-23, 44, 54-59], a full exploration of the sensitivity to variations of the LECs has been considered too difficult due to the immense computational cost.
In this work we address this problem by developing accurate emulators and adapting a robust statistical approach known as history matching [60][61][62] to explore the entire LEC parameter space.History matching is specifically designed to aid the analysis and calibration of highdimensional computationally expensive physical models.By monitoring predictions of scattering phase shifts and few-body observables we can iteratively identify the (nonimplausible) region of the parameter space that gives results consistent with a set of data.This approach then provides a finite domain in which the probability of finding accurate interaction parametrizations is higher and that makes it feasible to perform rigorous statistical studies even with limited computational resources.
The history matching procedure and the statistical analysis require a significant amount of computations.Therefore, emulators-that mimic the outputs of the exact calculations at a fraction of the computational costare required to bring this kind of statistical approach into practical use.Recently, model reduction methods [63] such as eigenvector continuation (EC) [44,[64][65][66][67][68] has proved to be an efficient and accurate approach to emulate the predictions of ab initio many-body methods.In this paper, we generalize the subspace-projected coupled cluster (SPCC) [65] to construct emulators for the energy per particle of PNM and SNM at different densities that work for a wide range of LECs.We will demonstrate that these emulators provide fast and accurate approximations to CC calculations of nuclear matter properties [15].While the SPCC method has been successfully applied to a global sensitivity analysis of bulk properties of 16 O [65], the challenge of the present work is that the bi-variational CC energy functional [69,70], combined with the increased level density in infinite nucleonic matter and large number of LECs, may give rise to spurious states when diagonalizing the subspace-projected non-Hermitian CC Hamiltonian.We therefore introduce a new algorithm called small-batch voting to efficiently locate the physical ground state and to identify and eliminate spurious states.These nuclear matter emulators, along with other emulators for light nuclear systems, are then applied to 1.7 × 10 6 different chiral interactions acquired by a history matching procedure.This allows us to carry out a comprehensive study of correlations between nuclear matter properties and observables in finite nuclei.A subsequent Bayesian analysis is performed by establishing error models for EFT truncations, method uncertainties, and using sampling/importance resampling [71,72] to obtain probabilistic distributions of both LEC parameters and posterior predictions.Details of nuclear matter posterior predictive distributions calibrated by different sets of observables are elucidated in a companion paper [24].

II. METHOD
In this work we consider ∆-full χEFT at NNLO [43,[73][74][75][76][77][78].The explicit inclusion of the ∆-isobar degree of freedom is beneficial since it increases the breakdown scale of the χEFT and gives a better description of nuclear matter properties [43,44,78].We use standard nonlocal regulators: for the NN and 3NF interactions respectively, with n = 4 and a fixed cutoff Λ = 394 MeV.The chiral Hamiltonian of ∆NNLO is parametrized with 17 LECs which are here represented by the vector ⃗ α.Following Refs.[65,66] it can be written as: with h 0 = t kin + V 0 .Here t kin is the kinetic energy and V 0 represents the constant potential term.One of the most important conclusions of EC is that the trajectory of the eigenvectors as a function of the smoothly varying control parameters of the Hamitlonian (in our case the LECs) can be well described by a finite-dimensional manifold [64].With this statement, the ground-state eigenvector of any target Hamiltonian H(⃗ α ⊚ ) can be well approximated by some linear combinations of ground-state eigenvectors of a finite set of training Hamiltonians H(⃗ α 1 ), In practice, to create a subspace emulator we therefore need to use a many-body solver to generate N sub different ground-state eigenvectors for a set of training points ⃗ α i .Any target Hamiltonian (H(⃗ α ⊚ )) is then projected onto this subspace, and the approximate ground state is obtained by diagonalizing a generalized eigenvalue problem.By choosing an optimal set of training points, the ground-state eigenvalue and eigenvector from this subspace rapidly converge to the full-space solution as the number of training points is increased [67,79].In addition, since the target Hamiltonian H(⃗ α ⊚ ) is expressed by a finite number of terms that depend linearly on the LECs, one can project each term h i onto the subspace such that the resulting matrices can be stored and used to quickly construct the projection of any target Hamiltonian H(⃗ α ⊚ ) onto the subspace of training vectors.
The proof of convergence of EC emulators as outlined in Ref. [67] has not been generalized to the case of non-Hermitian Hamiltonians as encountered in the CC method.However, the studies of Refs.[44,65] showed a rapid convergence to the full-space solution also in this case.In this work we employ no-core shell model (NCSM)-based EC emulators from Hu et al. [44] for 2,3 H and 4 He observables.In addition, we construct a new emulator for 6 Li using JupiterNCSM [80,81] to perform model checking.For 16 O and nuclear matter we use the CC method [70,[82][83][84][85][86][87][88][89][90][91][92][93] as the many-body solver and employ the corresponding SPCC method [65] to construct emulators.For 16 O observables that we use in our statistical analysis we adopt the SPCC-emulators based on the CC method with singles-doubles and leading-order triples excitations (CCSDT-3) as described in Ref. [44].The novel small-batch voting algorithm, that will be presented below, is exclusively used for the nuclear matter SPCC emulators.

A. Subspace-projected coupled-cluster
The essence of the CC method is its similarity transformed Hamiltonian: where T (⃗ α) is the cluster operator that induces all possible particle-hole excitations.For the nuclear matter emulators in this work we adopt the CC method with doubles (CCD) approximation, so that the cluster operator is truncated at two-particle-two-hole excitations, i.e.T (⃗ α) = T 2 (⃗ α).We note that there are no oneparticle-one-hole excitations in infinite nuclear matter, i.e., T 1 (⃗ α) = 0, due to momentum conservation.Note that the transformation in Eq.( 2) is non-unitary and the resulting similarity transformed Hamiltonian H(⃗ α) is non-Hermitian.The direct consequence is that the CC method is non-variational, but it follows a bi-variational principle [69] by parametrizing the left and right CC ground states as: where Λ(⃗ α) = Λ 2 (⃗ α) is a two-hole-two-particle deexcitation operator.We note here that the parametrization of the left state is the first-order approximation to Arponen's extended CC method [94] where the left state is parametrized in a more symmetric way by writing α) .In this work we determined Λ(⃗ α) amplitudes by solving the eigenvalue problem for the left ground state [90].The bi-orthonormality is assured by ⟨ Ψ|Ψ⟩ = 1 when both left and right ground-states are acquired by the same ⃗ α.The reference state |Φ 0 ⟩ is chosen to be the closed-shell configuration on a discrete lattice in momentum space with periodic boundary conditions.The model-space for which we solve the CCD equations has (2n max + 1) 3 momentum points and we set n max = 4, which is sufficiently large to obtain converged results [15].To minimize finite-size effects we use 132 nucleons for SNM and 66 neutrons for PNM, respectively [10,15].
In order to construct the subspace projected target Hamiltonian H(⃗ α ⊚ ) we solve for the left and right CC ground-states for a set of N sub training Hamiltonians H(⃗ α 1 ), • • • , H(⃗ α N sub ), and subsequently project H(⃗ α ⊚ ) and the identity matrix onto this subspace giving, where e X = e −T ′ +T and we indicate quantities related to different training points, ⃗ α and ⃗ α ′ , by unprimed and primed symbols.With Eqs. ( 4) and ( 5) one can easily acquire the ground-state energy for the nuclear matter system by solving a N sub × N sub generalized eigenvalue problem.Note that the N sub subspace vectors should not be linearly dependent as it would induce numerical instability when solving the generalized eigenvalue problem.
Another important aspect of the SPCC method is the selection of an appropriate set of training points to construct the subspace.To ensure that the selected training points will lead to accurate emulators in the relevant parameter domain, we first apply history matching to restrict the LEC ranges.Furthermore, to maximize the worth of full computations we select the non-implausible samples with highest likelihood in the final Bayesian analysis as training vectors.More details about the history matching and Bayesian analysis can be found in Sec.II E and Sec.IV.The training points used in this work are shown in Fig. 1 and, as can be seen, they cover a very broad LEC range.

B. Emulators for a single LEC
We first consider an illustrative example with simpler nuclear matter emulators.For this purpose we use an interaction model with a single free parameter and perform nuclear-matter modeling with a smaller number of particles.Specifically, we employ the ∆NNLO GO (394) interaction [43] and allow the single C 1S0 LEC to vary and we use 14 (28) neutrons (nucleons) for PNM and SNM, respectively.Fig. 2 shows the calculated energy per neutron (E/N ) and energy per nucleon (E/A) for PNM and SNM, respectively.The SPCC predictions using three or five subspace vectors are compared with full space CCD results for a wide range of the low-energy constant C 1S0 (the remaining LECs are kept fixed).As we can see, using N sub = 5 training points chosen in a small region, the SPCC method already accurately reproduces the full space CCD calculations over a large range for the C 1S0 LEC.As expected, if we reduce the number of training points to N sub = 3, the SPCC predictions of SNM start to deviate more from the exact solutions in the case of large exptrapolations.However, the predictions for PNM remain accurate over the whole range considered.The choice of training points in Fig. 2

C. Small-batch voting
When building SPCC emulators for systems with 132 nucleons for SNM one suffers from a persistent spurious state problem.We find that there can be several eigenstates of the N sub ×N sub matrix that have lower eigenvalues than the corresponding full-space CCD result.The exact reason for the appearance of these states is not yet fully understood, but it could be a consequence of several factors: (i) the SPCC Hamiltonian is by construction non-Hermitian and the variational theorem does not apply; (ii) for increasing number of nucleons (132 nucleons and 66 neutrons in our case) the level density increases which more easily leads to the occurence of these states, and (iii) increasing correlation energies associated with less perturbative interactions might produce more spurious states.
As we consider these spurious states to be unphysical we seek a method to identify them a priori in order to remove them from the spectrum.Recall that CC theory fulfills a bivariational theorem and the physical solution is a stationary point with respect to variations of the CC amplitudes.Whether the bivariational property of CC theory also holds for the SPCC remains to be shown, but it is reasonable to assume that it holds as long as the subspace is sufficiently large.In this section we will show   394) [43].The red diamonds correspond to exact CCD calculations.The C1S0 value of ∆NNLOGO( 394) is indicated with a dashed vertical line.This one-dimensional emulator demonstration is obtained for a model with N = 14 (A = 28) for PNM (SNM).
how we can use the bivariational property to efficiently identify the physical solution within the SPCC spectrum using a method we call small-batch voting.
Fig. 3(a)(b) illustrates the relative errors (E SPCC − E CCD )/|E CCD | between emulator predictions and exact CCD calculations for PNM (SNM) modeled with 66 (132) neutrons (nucleons).The emulator predictions are chosen as the SPCC solution with the lowest (real) energy.The validation points correspond to 50 random parametrizations of the ∆NNLO(394) interaction with the energy per particle computed for five densities: ρ ∈ {0.12, 0.14, 0.16, 0.18, 0.20} fm −3 .Thus, there are 250 points in total.In these calculations, N sub = 64 subspace vectors are used to construct the SPCC emulators.Note that the validation interactions are selected randomly within a constrained parameter domain (resulting from history matching).For SNM, most errors are negative which indicates that the emulator predictions correspond, in fact, to spurious states that give lower energies than the corresponding full-space CC calculations.Based on this, it seems the occurrence of spurious states seriously hampers the predictive power of the SPCC method for SNM.The fact that the SPCC method works better for PNM is probably due to the smaller correlation energies for this system, i.e., it is more perturbative.However, when comparing the full eigenspectrum of the subspace problem to the full-space CC result, we find that the physical ground-state indeed is contained therein.The challenge then is to efficiently identify this state without actually doing any full-space CC calculations.Assuming that the subspace is large enough to accurately describe the exact CC ground-state solution, the eigenvector should also fulfill the bivariational principle.This means we can write the exact solution as, here ⃗ c ⋆ is the SPCC eigenvector that corresponds to the physical ground-state (denoted by ⋆).The physical solution converges rapidly with increasing number of training points.By utilizing the bivariational property, it is then reasonable to assume that the physical state should remain stable (displaying small variations in the corresponding eigenenergy) when removing a small portion ∆N of the subspace, while the spurious states and their energies should change significantly.Here we introduce batches, N batch = N sub − ∆N , that should be large enough to still provide an accurate representation of the physical ground state.
Based on this argument we develop a new algorithm called small-batch voting to efficiently identify the physical ground state and counter the spurious-state problem.The procedure of small-batch voting can be summarized as follows: 1. Solve the target Hamiltonian H(⃗ α ⊚ ) using the SPCC method in a subspace with relatively large N sub to ensure that the physical ground-state is well established in the spectrum.The eigenvalues of H(⃗ α ⊚ ) are stored as E i=1,...N sub 2. Construct k different small batches by randomly picking subsets of N batch (< N sub ) vectors from the original subspace.
3. For each batch, solve the generalized N batch × N batch eigenvalue problem.Compare the eigenvalues e r=1,...N batch with the eigenvalues of the original subspace E i=1,...N sub .If the same eigenvalue occurs in both spectra (within a specified relative tolerance 4. Repeat step 3 for all k small batches. 5. The eigenvector with the highest number of votes v * = max(v i ) is assumed to correspond to the physical ground state (with the lowest energy as a deciding vote if there is a draw) and its eigenvalue E * is used as the emulator prediction.
In this work, we set N sub = 64, k = 100, N batch = 30, and the relative tolerance ε = 0.02.We checked that N batch = 30 is sufficiently large to reproduce the fullspace CC solution to within 1% when we know exactly which is the physical ground-state in the spectrum.To summarize, the essential idea of the algorithm is that by varying the composition of the subspace, the eigenvalues of the spurious states will be shifted dramatically since they are not stationary solutions, while the physical ground-state remains relatively unchanged.Fig. 3 summarizes our results with and without smallbatch voting.Panels (a) and (b) show the relative error between the lowest SPCC energies (without small-batch voting) and the corresponding exact CCD results, while panel (c) show the results of SNM with small-batch voting.For the latter we apply a 1% mean shift up in energy since the voting procedure favors the lowest state that is found within the 2% tolerance window.The comparison demonstrates that the small-batch voting algorithm successfully removes most of the spurious states and that the emulator predictions are much improved.Note that for PNM the SPCC predictions are already extremely accurate thus we do not apply small batch voting for the PNM emulator.The standard deviation of the relative error is σ = 0.002% (0.8%) for PNM (SNM).We note that there are still a few spurious states that remain after applying the small-batch voting algorithm for the more complex non-perturbative SNM case.The total computational cost of the SPCC emulators for SNM with smallbatch voting is six orders of magnitude smaller than the corresponding full space CC calculations.For PNM emulators, where we don't need small-batch voting, we have a computational speedup of more than eight orders of magnitude.

D. Nuclear matter saturation observables
At this point we are able to construct emulators for PNM and SNM and predict the energy per particle at a given density using the SPCC method with small-batch voting.In order to study nuclear matter saturation properties, e.g., the saturation density ρ 0 , the saturation energy E 0 /A, the symmetry energy S, the symmetry energy slope parameter L and the incompressibility K, one needs to acquire the EOS for both pure neutron and symmetric nuclear matter around the saturation point.Ideally one would like to include the density ρ parameter in the eigenvector continuation scheme and build an emulator that works for different LECs and at arbitrary densities.However, changing the density leads to different discretizations of the momentum space lattice and one would therefore need to work out matrix elements connecting different reference states and lattices.
Fortunately, we are not completely ignorant about the properties of the EOS of nuclear matter since the energy per particle should be represented by continuous smooth functions of ρ.This smoothness implies that we do not need many density points to obtain sufficient information about the EOS.In this work, we construct SPCC emulators for both PNM and SNM at five different densities: ρ = 0.12, 014, 0.16, 0.18, 0.20 fm −3 .We choose to study this density region simply because the empirical saturation density is around 0.16 fm −3 [11,95].The nuclear matter EOS is then interpolated within this range by using Gaussian processes (GPs) [96] as the interpolation method.We choose the radial basis function (RBF) as the correlation function to ensure the smoothness of the EOS.The hyperparameter (correlation length l) of the GP is learned from a validation data set which contains 50 interaction samples that are generated by the same history matching process mentioned in Sec.II C. The PNM and SNM correlation lengths optimized from the validation set are 0.297 fm −3 and 0.259 fm −3 , respectively.In the end, we take the more conservative value l = 0.25 fm −3 for both PNM and SNM so that we do not overestimate the correlation length.
The major advantage of using GPs is that they are infinitely differentiable under the RBF kernel and that derivative properties such as L and K can easily be obtained.For any given interaction with its saturation point ρ 0 ∈ [0.12, 0.20] fm −3 we can therefore extract all saturation properties from the corresponding GPs and its (first and second) derivatives.
Fig. 4 shows the performance of GP interpolation with different density spacing ρ gap from 0.02 to 0.0025 fm −3 .It can be seen that L and K, which correspond to first and second derivatives of the EOS, have larger deviation when the density spacing is increased.In principle, we expect better accuracy with smaller density spacing.However, the values of saturation properties at ρ gap = 0.02 fm −3 differ by less than 0.1% compared with ρ gap = 0.0025 fm −3 .This difference is rather small compared to other sources of uncertainty.In practice  we therefore use ρ gap = 0.02 fm −3 and ignore the GP interpolation error.

E. History matching
In this work we use an iterative history matching approach [44,[60][61][62] with selected experimental fewnucleon data to study and reduce the huge parameter space of our χEFT interaction model.For each wave of history matching we need to establish a quantitative criterion that determines if a parametrization ⃗ α yields acceptable (or at least not implausible) model predictions when confronted with the selected set of observations Z.We first introduce the individual implausibility measure which includes the squared difference between the model prediction M i (⃗ α) and the observation z i for observable i from the target set Z. The total variance in the denominator of Eq. ( 7) is here constructed under the assumption of independent errors.It is therefore a sum of variance terms that in our case include experimental, model, method, and emulator errors.Unless differently specified we use the maximum of the individual measures to define the implausibility constraint where the default choice is c I ≡ 3.0 inspired by Pukelheim's three-sigma rule [97].
History matching proceeds by reducing the parameter space iteratively.In each wave one removes regions that are deemed implausible by failing the constraint in Eq. (8) .A visualization of this process is shown in Fig. 5.We first use a space-filling Latin hypercube design [98] to generate well-spaced interaction samples in the input parameter domain.Then we use fast modeling or emulation to compute the implausibility measures and apply the maximum implausibility constraint.The remaining nonimplausible interaction samples are kept and define the non-implausible region for the next wave.Note that we are using parallelograms to define the bivariate surfaces of the non-implausible volume which allows to incorporate parameter correlations.In this work the iterative history matching is carried out in five waves as shown in Table I.Initial waves comprise selected groups of observables and subsets of active input parameters.This enables reaching sufficiently high resolution in the spacefilling design of interaction samples.Figs. 6 and 7 show the non-implausible volume for each wave including the final one.In waves 1 and 2 we constrain all relevant LECs (except c D and c E ) grouped by partial waves ( 1 S 0 , 3 S 1 , 1 P 1 , 3 P 0 , 3 P 1 , 3 P 2 ) using neutron-proton scattering phase shifts at six energies (T lab = 1, 5, 25, 50, 100, 200).
In wave 3 we include only the deuteron ground-state energy, point-proton radius and quadrupole moment as target observables and consider C3S1 , C 3S1 , C E1 and c 1,2,3,4 as active parameters.In wave 4, we add the 3 H binding energy and the 4 He binding energy and point-proton radius to the set of target observables.Here we consider also the three-nucleon force parameters c D , c E along with the other LECs.In this wave, however, we fixed the four ) to the values from the ∆NNLO GO (394) interaction [43], since the selected target observables are not very sensitive to these parameters.We added an additional method uncertainty to the denominator of the implausibility metric (7) to capture the reduced precision of the model with fixed P -wave parameters.Following a sensitivity study we set the standard deviation of this additional error to 100 and 400 keV for the 3 H and 4 He binding energy, respectively, and to 0.03 fm 2 for the squared point-proton radius of 4 He.Finally in wave 5, all 17 model parameters are active and the non-implausible domain is explored by 1 × 10 9 space-filling design samples.In the end we find that 1.7 × 10 6 of them pass the implausibility constraint with the same set of few-nucleon target observables as in wave 4. It is worth to mention that the order at which observables are considered in history matching waves is irrelevant as long as one uses the maximum implausibility as the constraint.With the maximum implausibility measure, the final non-implausible region is the intersection of constrained parameter regions from different waves and is therefore unrelated to wave ordering.

F. Bayesian machine learning error model
In order to connect the emulator predictions with actual nuclear matter properties one needs to incorporate errors from different sources.A statistical model for the density-dependent energy per particle, that accounts for the most relevant sources of uncertainty, can be written as (9) where y k (ρ) is the nuclear matter emulator prediction using our EFT model truncated at order k at given density ρ, while the stochastic terms ε k (ρ), ε method (ρ), and ε emu (ρ) correspond to the EFT truncation error (model discrepancy), the method error, and the emulator error, respectively.To quantify the distributions of these stochastic variables we apply and extend a Bayesian machine learning error model that was originally proposed by Drischler et al. [22].Within this error model we construct multitask GPs [100] to estimate both the variance and the covariance of target errors as a function of density and proton-to-neutron fraction (SNM and PNM) from given prior information.
Let us first consider the EFT truncation error.We follow Refs.[101,102] and write the EFT expansion for an observable, truncated at order k, as Here y ref is a reference scale for the observable y, c n are dimensionless expansion coefficients (with c 1 = 0 in Weinberg power counting), while the expansion parameter Q = k F /Λ b is the ratio of the Fermi momentum k F and the breakdown momentum Λ b .Here we use Λ b = 600 MeV.Note that for simplicity we fixed Λ b and used the Fermi momentum as the energy scale for the EOS.In future studies one could attempt to infer Λ b or Q directly when making inferences to data that is more sensitive to higher energies.With a simple extension of Eq. 10 to infinite order, the truncation error ε k at order k can be expressed as We will infer the expansion coefficients c n (ρ) given EFT convergence assumptions and choose y ref (ρ) = y LO (ρ).We further assume that c i (ρ) and c j (ρ), at different orders i and j, should be independent and identically distributed random functions of natural size.Thus, the error model assumes that they can be described by a single underlying GP where we use a (Gaussian) radial basis correlation function r(ρ, ρ ′ ; l) with c2 and l the GP hyperparameters corresponding to the variance and the correlation length, respectively.Note that in these equations ρ and l are measured in fm −1 as we translate from density to the corresponding Fermi momentum for each type of nuclear matter.The mean function of the GP is taken to be zero since the correction at each order can be positive or negative.With Eqs.(11) and ( 12) one can easily derive the EFT truncation error as a geometric sum of independent normally distributed variables.Its distribution then follows with We note that k = 3 for ∆NNLO.
Having defined the GP that describes the truncation error, the hyperparameters c2 and l can be inferred from order-by-order EFT predictions and expert elicitation [98].Using data D, corresponding to order-by-order predictions at several densities ρ, and the incorporation of EFT expectations via a prior pr(c 2 , l), the posterior for c2 and l becomes Here we use a scaled inverse-chi-squared distribution [103] as the prior for c2 and a uniform prior for l (with l ∈ [0, 10]).In practice, we first train two GPs separately with order-by-order predictions of the EOS for PNM and SNM (at discrete density points ρ = 0.06, 0.14, ...0.38 fm −3 ) using leading order (LO), next-to-leading order (NLO) and NNLO ∆-full interactions [104] that were optimized using the protocol described in Ref. [78].See the Supplemental Material [105] for numerical values of the LECs that define these convergence-study interactions.For simplicity we then used the maximum a posteriori (MAP) value as a point estimate for the hyperparameters.The hyperparameters learned in this way from the training data are c1 = 0.99 and l 1 = 0.88 fm −1 for PNM and c2 = 1.66 and l 2 = 0.45 fm −1 for SNM.These two GPs (separately trained for PNM and SNM) describe the correlation of truncation errors as a function of density for either system individually.As discussed in Refs.[22,23] it is crucial to also account for the cross correlation between PNM and SNM truncation errors.This is important to avoid overestimating the total uncertainty for observables such as the symmetry energy S that corresponds to the difference between E/N and E/A.Applying a multitask GP model, the truncation errors of PNM and SNM become jointly distributed where K ii is the covariance matrix generated by the kernel function c2 i R ε k (ρ, ρ ′ ; l i ) of PNM (i = 1) and SNM (i = 2) respectively, as described above, while K 12 = K T 21 is the cross-covariance that we describe with the kernel function ρ 12 c1 c2 R ε k (ρ, ρ ′ ; l 12 ).Following Ref. [23] we set (Color online) Non-implausible parameter domains for (C1P 1, C3P 0, C3P 1, C3P 2, cD, cE and c1,2,3,4) at the ends of the five history matching waves.The initial parameter domain is represented by the axes limits for all panels (except c1,2,3,4).The volume of the non-implausible domain is iteratively reduced in waves 2, 3, 4, and 5 (shown by green/dash-dotted, blue/dashed, black/dotted, red/solid rectangles, respectively).The non-implausible samples in the final wave are shown as two-dimensional histograms (purple).Note that the sampled volume for c1,2,3,4 (illustrated by red/solid and dashed contour lines denoting 68% and 90% credible regions, respectively) remain the same in all waves.In practice, for these LECs, we use a four-dimensional hypercube mapped onto the multivariate Gaussian probability density function resulting from a Roy-Steiner analysis of pionnucleon scattering data [99].The contour lines (purple/solid and dashed) for the non-implausible samples identified in the final wave are shown for comparison.
In this work we extend the statistical error model by also incorporating the method error and its correlation structure into the analysis.Specifically we consider our two main sources of method uncertainty, namely the truncation of the cluster operator and the finite-size effect of the cubic momentum lattice.Assuming independence, the total method error can then be written as ε method = ε cc + ε fs .We will again use the GP error model such that where the subscript "κ" can be either the cluster operator truncation "cc" or finite-size effect "fs", and As before, we will use the corresponding Fermi momentum (in fm −1 ) as the independent variable rather than the density.For the cluster operator truncation we estimate the density-dependent mean error and covariance using results from a previous convergence study with 34 ∆-full interactions at NNLO [44] (with the same 394 MeV momentum cutoff as here).In that study, computations were performed at CCD(T) level which is a more accurate, but computationally heavier, CC approximation that includes doubles excitations and perturbative triples corrections [15].In particular, triples correla-tion energies-the difference between CCD and CCD(T) results-were extracted for the 34 interactions.We use the average (per density) triples correlation energy as our mean errors, resulting in µ cc,PNM (ρ) = 0.16k 2 F − 0.50k F + 0.32 MeV/nucleon for PNM and µ cc,SNM (ρ) = −1.28kF +0.80 MeV/nucleon for SNM.At ρ = 0.16 fm −3 these mean shifts are −0.04MeV/nucleon for PNM and −0.91 MeV/nucleon for SNM.Furthermore, we also use these mean triples correlation energies as our reference scale y cc,ref (ρ) which implies that ccc should be interpreted as the ratio between the cluster truncation error and the triples correlation energy.Following previous CC convergence studies [106,107] and expert elicitation we conservatively assign ccc = 0.1 for both PNM and SNM corresponding to ±20% of the triples correlation energy as a 95% degree-of-belief error estimate of corrections beyond the triples approximation.We also assume that the observed density-dependence of the energy differences between the CCD and CCD(T) results for the 34 interactions in Ref. [44] provides a relevant measure of the correlation structure of the CC truncation method error.This data is therefore used to train the GPs and the inferred correlation lengths (within the density range ρ ∈ [0.12, 0.20] fm −3 ) are l cc,PNM = 0.50 fm −1 and l cc,SNM = 0.58 fm −1 .
For the finite-size effect, we use the CCD groundstate energy as the reference scale and set correlation length l cc,PNM = 0.50 fm −1 and l cc,SNM = 0.58 fm −1 .Following the study in Ref [15], we take ±0.5%(±4%) for PNM(SNM) as estimates of the 95% credible intervals of finite-size errors for each system.This gives cfs,PNM = 0.0025 (c fs,SNM = 0.02) for the nuclear-matter calculations in this work.
Finally, we also use the GP error model described in Eqs. ( 17) and ( 18) to incorporate the emulator error ε emu (ρ) in Eq. ( 9).These GPs were trained by the differences between emulator predictions and CCD results with the latter then used as reference scale.The training data are taken from the CCD computations and emulator predictions of the 34 interactions in Ref. [44].For the SNM emulator error we found l emu,SNM = 0.38 fm −1 and we use µ emu,SNM = 0 and cemu,SNM = 0.01.We ignore the small emulator error for PNM due to the high accuracy achieved by the PNM emulator predictions.
Following these assignments, the full posterior predictive distribution (PPD) for nuclear matter observables, incorporating all relevant sources of uncertainty, can be sampled according to Eq. (9).In particular, it becomes straightforward to sample the error terms from the corresponding covariance matrices once the multitask GPs are determined.In practice, this task is efficiently performed using with L being the Cholesky decomposition of the cross covariance matrix K (K = LL T ) and x a standard normal random vector.Note that we emulate results at five densities for both PNM and SNM.Thus ε is a 10dimensional vector and the cross covariance is a 10 × 10 matrix.This sampling procedure is crucial for generating the PPD of nuclear matter properties.The emulator predictions for the nuclear matter EOS and the corresponding 2σ (95%) credible interval for errors are illustrated in Fig.  this figure it is also clear that the method error for PNM is quite small.This can be understood since emulator errors, finite-size effects and CC correlation energies are all rather small for PNM.

III. HISTORY MATCHING ANALYSIS
The tremendous computational speed-up offered by our novel nuclear-matter emulators allows to perform a detailed statistical analysis of observable predictions using the χEFT model.The results shown in this section represent general outcomes of the interaction model described in Sec.II used within ab initio computations.
As a first step of this analysis we apply the history matching procedure as described in Sec.II E with five waves of global parameter search to iteratively reduce the LEC domain.The history-matching is performed using neutron-proton phase shifts in S-and P -waves plus fewbody (A = 2 − 4) bound-state observables in the target sets.The experimental values and the error assignments for the nuclear bound-state observables can be found in Table II.The experimental targets are from Refs.[108][109][110].Note that the target point-proton radii were transformed from experimental charge radii using the same relation as in Ref. [50].For the deuteron quadrupole moment we use the theoretical result obtained by the CD-Bonn [110] model with a 4% error bar.The theoretical model ε model (method ε method ) errors are estimated from the EFT (CC) convergence pattern as in Ref. [44] while emulator errors ε em are estimated from cross validation.This selection of target data is representative of what could have been considered when seeking an optimal interaction model.However, the aim of our approach is fundamentally different.Rather than seeking a single optimum, we consider all non-implausible parametrizations in order to make a comprehensive study of the behaviour of our model.Furthermore, we consider much simpler linearised probability distributions, with just mean values and variances as specifiers, to identify the interesting parameter domain.Finally, we just divide the parameter space into implausible or non implausible.All samples from the latter domain are included in this part of the analysis without any probability weighting.
In the final wave, we explored 1 × 10 9 samples from a space-filling design in the non-implausible domain that was established at the end of wave 4. We confronted the model predictions for the six A = 2 − 4 observables and found 1.7 × 10 6 non-implausible interaction parametrizations.At this point, we did not see any need to proceed with another wave since there were no signs of further reduction of the parameter domain.The 1.7 × 10 6 samples constitute a good representation of all non-implausible interactions.
Predictions for different few-body observables are shown in Fig. 9.Here we compare model predictions made with the 1 × 10 9 random samples generated at the start of wave 5 (hashed histograms) with the results ob-tained with the 1.7 × 10 6 samples that survive the implausibility constraint.As shown, the predictions with the random samples are characterized by very large variances.Clearly, the 17-dimensional LEC domain is still quite large even after the history matching waves.As for the 1.7 × 10 6 non-implausible samples, all of them give results within ±3σ error regions for the A = 2 − 4 observables (as indicated by the red band) since these observables were included as target data in the history matching procedure and we used c I = 3 for the implausibility constraint (8).
The prediction for the 6 Li ground state energy (last panel of Fig. 9) serves as model validation since it was not included in the history matching.We can see that the mode of the E( 6 Li) histogram is reasonable, and clearly within the 3σ region, which indicates a very reasonable model performance for light nuclei.We then consider model predictions for the infinite nuclear matter systems using the SPCC emulators from 0.12 0.14 0.16 0.18 0.20 Sec.II A with small-batch voting and GP interpolation as described in Secs.II C and II D. In this particular analysis we don't perform a full sampling of the error model outlined in Sec.II F but only include the mean shift of the EOS for SNM that is expected from triples corrections.This shift is applied for all results shown in Figs. 10, 11, and 12. Saturation properties for the non-implausible interaction samples are shown in Fig. 10.All results are obtained using the nuclear matter emulator outputs as described in Sec.II D. Interactions that give a saturation density outside of the interval ρ ∈ [0.12, 0.20] fm −3 (about 27% of all non-implausible interactions) are not shown since our emulators are only constructed within this density interval.It is quite clear from Fig. 10 that the modes for saturation density, saturation energy and the symmetry energy deviate from the empirical region and that here is a very large variance (in particular for the saturation density).We hypothesize that this is a consequence of the large extrapolation from the history matching observables in light nuclei (A = 2 − 4) to properties of inifinite nuclear mater, and to the limited correlation (see below) between few-body observables and properties of nuclear matter.
The large number of non-implausible interaction samples together with access to fast and accurate emulators enable an extensive study of correlations between properties of finite nuclei and infinite nuclear matter obatined  ), saturation energy E0/A, symmetry energy S, slope L, and incompressibility K (all in MeV) and selected observables of finite nuclei: ground-state energies (in MeV) and radii (in fm) of 4 He, 6 Li, 16 O.All results are obtained with the 1.7 × 10 6 non-implausible interactions from the fifth wave of history matching.The Pearson correlation coefficient r is indicated in each panel.Note that 4 He observables (thick, blue panel axes) are included in the history matching procedure while the other observables are pure predictions.
with chiral forces.The results of such a correlation study are shown in Fig. 11.Note that the 6 Li and 16 O observables are model predictions while the 4 He ones are part of the history-matching procedure.We observe a positive correlation between ground-state energies of finite nuclei and the saturation energy E 0 /A.This correlation is getting stronger from 4 He via 6 Li to 16 O (r = 0.29, 0.85, and 0.98, respectively).This is reasonable since the central density of heavier system is closer to the density of nuclear matter at saturation.On the other hand we find an anti-correlation between ground-state energies and the saturation density ρ 0 .This negative correlation is also getting more prominent in heavier system.As for the 16 O radius we observe a similar correlation structure as for the energy meaning a positive correlation with E 0 /A and anti-correlation with ρ 0 .We stress that these correlations are general results of the design of the ∆NNLO interaction model and that they are characteristic features of the corresponding Hamiltonian.
The correlation between selected LECs and nuclear matter properties are shown in Fig. 12.Even though 3NFs should be important for an accurate description of the saturation of nuclear matter, we find that the correlation between the 3NF short-range contacts c D and c E and nuclear matter properties is weak (as indicated by small Pearson correlation coefficients).This observation is consistent with the fact that c D and c E are not well constrained by the present history matching observables and that their contribution to the 3NFs are not exclusive since other LECs such as c 1,3,4 also play an important role.It is interesting to note that the singlet S-wave contact C 1S0 gives the strongest correlation with the slope L of the symmetry energy, which is also known to be strongly correlated with the neutron skin thickness of 208 Pb [34,44].This result indicates that this particular LEC serves as a bridge between neutron-proton scattering in the 1 S 0 partial wave and the thickness of neutron skins in finite nuclei.A more detailed discussion of this constraint on the allowable range of the 48 Ca and 208 Pb neutron skin thicknesses can be found in Hu et al. [44].12. (Color online) Correlation structure between nuclear matter saturation properties: saturation density ρ0 (in fm −3 ), saturation energy E0/A, symmetry energy S, slope L, and incompressibility K (all in MeV) and selected LECs (cD, cE, C1S0, and C3P 0).The parameters ci and C1S0 are in units of GeV −1 and 10 4 GeV −4 , respectively.Results are obtained with the 1.7 × 10 6 non-implausible interactions from the fifth wave of history matching.The Pearson correlation coefficient r is indicated in each panel.

IV. BAYESIAN ANALYSIS: POSTERIOR PREDICTIVE DISTRIBUTIONS
The history matching results do not offer a probabilistic interpretation since no actual probability distributions were invoked.For the target data we only considered an implausibility criterion rather than a fully specified likelihood.The non-implausible samples do, however, offer an excellent starting point for a Bayesian analysis.In order to acquire PPDs for nuclear matter and finite nuclei observables we proceed as follows: First, since the final wave of history matching procedure does not include phase shifts, we confront all 1.7 × 10 6 nonimplausible samples with the phase shift targets from the first wave (including S and P partial waves up to T lab = 200 MeV) and apply the implausibility constraint (8).We also examine whether the samples give an np bound state in the 1 S 0 channel as a sanity check.Just a few samples failed this test.Taken together, these constraints reduce the number of non-implausible samples to 8,218.Second, we use the method of sampling/importance resampling [71,72] to extract an approximate posterior probability density function (PDF) of the LECs via Bayes' theorem We assume a uniform prior probability distribution, pr(⃗ α), for all LECs except c 1,2,3,4 for which the prior is described by a multivariate normal distribution originating in the Roy-Steiner analysis of πN scattering data performed in Ref. [99].The history matching procedure provides a set of samples from this prior.Although the full data likelihood is not involved, the incorporation of implausibility constraints guarantees that samples with a negligible contribution to the posterior PDF are removed.Operating with the remaining large set of prior samples {⃗ α i } n i=1 we now specify a data likelihood and evaluate ω i ≡ L (D cal | ⃗ α i ) and so called importance weights q i ≡ ω i / n j=1 ω j .Finally, we resample a set {⃗ α * i } N i=1 from the discrete distribution {⃗ α i } n i=1 according to the importance weights q i .This resampled set will then be approximately distributed according to the target distribution pr (⃗ α | D cal ) ∝ L (D cal | ⃗ α) pr(⃗ α).See Ref. [72] for a recent importance resampling review with a nuclear theory perspective.We have studied the convergence of the posterior and found that a resampling set of size N = 10, 000 is sufficient.We also found that a rather large subset of > 400 samples from {⃗ α i } n i=1 provide 95% of the posterior PDF samples.
In order to examine how the choice of calibration data in the likelihood L (D cal | ⃗ α) affects the χEFT prediction we considered two different versions: (i) D cal = D A=2,3,4 encompassing binding energies and radii of 2,3 H and 4 He including the quadrupole moment of the deuteron, and (ii) D cal = D A=2,3,4,16 where we also include the energy and radius of 16 O.The default choice for the functional form of the likelihood is a normal distribution with independent errors (as summarized in Table II).The sensitivity to this specification of uncorrelated, Gaussian errors was tested using two alternative likelihood forms, namely a non-correlated Student-t distribution (with ν = 5 degrees of freedom implying heavier tails) and a multivariate normal distribution with positive correlation (ρ = 0.7) between the ground-state energy and radius of the same nucleus and between ground-state energies of different nuclei.In the end, we found no significant impact using the alternative distributions and therefore only show results obtained with the default, uncorrelated Gaussian likelihoods.See the Supplemental material [105] for numerical LEC values of the nonimplausible samples, the corresponding observable predictions using the emulators described in the text, and the likelihood PDF values.
The model PPD for an observable can be written as the set of model predictions evaluated for samples drawn from the parameter posterior for which we use the resampled set {⃗ α * i } N i=1 .From Eq. 21 it is clear that the predictive distribution is conditional on the selected calibration data D cal .
Fig. 13 shows the predicted distribution of nuclear matter properties calibrated by either D A=2,3,4 or D A=2,3,4,16 .Here we collect samples from the full PPD for which we also sample different sources of uncertainty as discussed in Sec.II F. For PPD A=2,3,4 (PPD A=2,3,4,16 ) we find that 15%(60%) of the samples are drawn from the training set and therefore have no emulator error.The PPD A=2,3,4 is shown in the left column of Fig. 13.The modes of the marginal distributions for saturation density, saturation energy and symmetry energy still deviate from the empirical values.We note that the saturation density PPD is asymmetric and our result almost indicates a bimodal distribution.The PPD A=2,3,4,16 is shown in the right column of Fig. 13 and provides an improved prediction of nuclear matter saturation with better precision.The saturation energy is slightly lower (more binding) compared with the empirical range and the mode of the incompressibility K is shifted to larger values.The predictions for S and L are not significantly affected by the addition of 16 O to the calibration data.The comparison of PPD A=2,3,4 and PPD A=2,3,4,16 reveals that the description of nuclear matter properties is quite sensitive to the choice of calibration observables.The reason is quite clear from the correlation results shown in Fig. 11.It is clear that the 16 O ground-state energy and radii are strongly correlated with nuclear matter properties.Thus a likelihood that contains these observables provides a more precise nuclear matter prediction.The full PPD for the energy per particle of PNM (top panel) and SNM (bottom panel) as a function of density is shown in Fig. 14.
Given the multi-stage analysis one can ask whether the final result is sensitive to the randomness of nonimplausible samples that results from the space-filling designs used in history matching.For this reason we performed two independent runs, from the start of the history matching to the final Bayesian analysis, with each one producing ≈ 5,000 non-implausible samples.It is the sum of those two runs that is presented in the upper half of each panel in Fig. 13.The lower half displays the smaller statistics result that is obtained with just the first run.The similarity of the upper and lower halves indicates the robustness of the approach and the fact that the convergence of the sampling/importance resampling 0.12 0.14 0.16 0.18 0 Run 1+2 Run 1 0.12 0.14 0.16 0.18 step is sufficient to accurately represent the target distribution.Finally, in Fig. 15 the LEC parameter PDF is shown conditional on the two calibration data sets.The marginal distribution of c D , c E , C 3P 0 , C 3P 1 and C 3P 2 are the most sensitive ones with respect to the choice of calibration data, while the differences in the other marginal distributions are barely distinguishable.This finding suggests that these terms in the chiral Hamiltonian are the most important ones to describe nuclear matter and heavier mass nuclei-yet remains poorly constrained by observables in the few-nucleon sector.

V. SUMMARY
We have constructed nuclear matter emulators using the SPCC method that works for a large 17-dimensional LEC hyperspace of ∆-full χEFT at NNLO.In particular, we have developed a small-batch voting algorithm to handle the spurious-state problem that can occur when emulating quantum many-body methods employing a non-Hermitian Hamiltonian.These nuclear matter emulators are then applied to 1.7 × 10 6 non-implausible interaction samples generated via five waves of history matching with A = 2 − 4 observables.This allows to study properties of the χEFT model including the correlation structure between nuclear matter saturation properties and observables of finite nuclei without bias from a specific optimization scheme.In particular we find an increasing correlation between saturation energy/density and the ground-state energy/radius of finite nuclei as the mass number of nuclei increases.In addition, a positive correlation between C 1S0 and the symmetry energy slope L is observed.
Starting from the history matching samples we performed a Bayesian analysis including relevant sources of uncertainty and using a correlated error model for the nuclear EOS.We applied the method of sampling/importance resampling method to obtain approximate samples of two parameter posterior PDFs with two different calibration data sets, D A=2,3,4 and D A=2,3,4,16 .The corresponding nuclear matter predictions (given by PPD A=2,3,4 and PPD A=2,3,4,16 ) illustrate the sensitivity to the calibration data.We found that predictions of nuclear matter saturation is more precise when incorporating the 16 O energy and radius in the likelihood calibration.We conclude that observables from 16 O are informative, but we note that they are not the only choice.We have seen that predictions for S and L were not significantly affected by the addition of 16 O to the calibration data.It will therefore be interesting to explore the information content of observables from neutron-rich systems.Furthermore, one should also consider other few-nucleon observables, such as the 3 H β decay rate for which emulators are also available [113], to monitor how they constrain the chiral interaction model and to explore whether this can lead to a satisfactory description of nucleonic matter ranging from light nuclei to infinite nuclear matter.

FIG. 2 .
FIG. 2. (Color online) Demonstration of SPCC predictions of PNM and SNM at ρ = 0.16 fm −3 , using three or five training points (circle markers), for different values of the low-energy constant C1S0.The other LECs are fixed at their values in ∆NNLOGO(394)[43].The red diamonds correspond to exact CCD calculations.The C1S0 value of ∆NNLOGO(394) is indicated with a dashed vertical line.This one-dimensional emulator demonstration is obtained for a model with N = 14 (A = 28) for PNM (SNM).

FIG. 3 .
FIG. 3. (Color online) Relative errors between SPCC predictions and exact CCD calculations for PNM (top) and SNM (bottom).Panels (a) and (b) show validation results without small-batch voting.For SNM we use small-batch voting for the final emulator.The significantly improved validation results are shown in panel (c).

FIG. 4 .
FIG. 4. (Color online)The convergence of saturation properties extracted by GP interpolation with different density spacings.The hyperparameter of the RBF kernel used in the GP is l = 0.25 fm −3 .In this convergence study, the energy per particle at discrete density points is obtained with the CC method for a single interaction.

FIG. 5 .
FIG. 5. (Color online) Visualization of the history matching procedure in a trivariate (C1S0, cD, cE)-subspace.The LECs cD and cE are dimensionless while C1S0 is in units of 10 4 GeV −4 .The non-implausible interaction samples are shown as dark dots.These dots are projected on different LEC surfaces and the outlines of the bounding regions are represented by parallelograms.The purple box outlines the final nonimplausible LEC domain.
FIG.6.(Color online) Non-implausible parameter domains for (C1P 1, C3P 0, C3P 1, C3P 2, cD, cE and c1,2,3,4) at the ends of the five history matching waves.The initial parameter domain is represented by the axes limits for all panels (except c1,2,3,4).The volume of the non-implausible domain is iteratively reduced in waves 2, 3, 4, and 5 (shown by green/dash-dotted, blue/dashed, black/dotted, red/solid rectangles, respectively).The non-implausible samples in the final wave are shown as two-dimensional histograms (purple).Note that the sampled volume for c1,2,3,4 (illustrated by red/solid and dashed contour lines denoting 68% and 90% credible regions, respectively) remain the same in all waves.In practice, for these LECs, we use a four-dimensional hypercube mapped onto the multivariate Gaussian probability density function resulting from a Roy-Steiner analysis of pionnucleon scattering data[99].The contour lines (purple/solid and dashed) for the non-implausible samples identified in the final wave are shown for comparison.

FIG. 7 .
FIG. 7. (Color online)Non-implausible parameter domains for ( C1S0,np, C1S0,nn, C1S0,pp, C3S1, C1S0, C3S1, CE1) at the ends of the five history matching waves.Note that the input volume for the isospin symmetry breaking LECs ( C1S0,np, C1S0,nn, C1S0,pp) are strongly correlated.The initial parameter domain is displayed by the black/solid quadrilaterals.The volume of the nonimplausible domain is iteratively reduced in waves 2, 3, 4, and 5 as shown by the green/dash-dotted, blue/dashed, black/dotted, red/solid quadrilaterals.The non-implausible samples in the final wave are illustrated as two-dimensional histograms (purple).Enlargements of the most relevant region of the LEC pairs of ( C3S1, C1S0,np, C1S0,nn, C1S0,pp) are shown in the top right panels.

FIG. 8 .
FIG. 8. (Color online)The EOS for PNM (top) and SNM (bottom) calculated for one representative interaction with the nuclear matter emulators (open squares) plus the mean value of the method error (solid squares).The bands indicate two standard deviations of the truncation error (green), method error (blue) and emulator error (pink) from the GP error models described in the text.The errors at different densities are correlated as illustrated by three random samples shown by dashed curves.Correlations extend between PNM and SNM (sampled error curves in the same colour).

FIG. 9 .
FIG. 9. (Color online) Histograms of A = 2 − 6 few-body observables predicted during history matching.The hashed histograms represent results obtained with 1 × 10 9 random samples from the wave 5 input LEC domain.The solid histograms correspond to model predictions with the final set of 1.7 × 10 6 non-implausible samples.Energies in (MeV), square of point-proton radii in (fm 2 ) and the deuteron quadrupole moment in (e 2 fm 2 ).The red dashed lines denote the experimental values and the red bands indicate the ±3σ error region.

55 FIG. 11
FIG. 11. (Color online)Correlation structure between nuclear matter saturation properties: saturation density ρ0 (in fm −3 ), saturation energy E0/A, symmetry energy S, slope L, and incompressibility K (all in MeV) and selected observables of finite nuclei: ground-state energies (in MeV) and radii (in fm) of4 He,6 Li,16 O.All results are obtained with the 1.7 × 10 6 non-implausible interactions from the fifth wave of history matching.The Pearson correlation coefficient r is indicated in each panel.Note that4 He observables (thick, blue panel axes) are included in the history matching procedure while the other observables are pure predictions.

30 S 75 L
FIG.12.(Color online) Correlation structure between nuclear matter saturation properties: saturation density ρ0 (in fm −3 ), saturation energy E0/A, symmetry energy S, slope L, and incompressibility K (all in MeV) and selected LECs (cD, cE, C1S0, and C3P 0).The parameters ci and C1S0 are in units of GeV −1 and 10 4 GeV −4 , respectively.Results are obtained with the 1.7 × 10 6 non-implausible interactions from the fifth wave of history matching.The Pearson correlation coefficient r is indicated in each panel.
FIG. 14. (Color online)The PPD for the EOS around saturation density, pr (E(ρ)/N, E(ρ)/A | DA=2,3,4,16).The sampling of the PPD includes all relevant errors as described in Sec.II F as well as the parametric uncertainty.

TABLE I .
Properties and summary statistics of the five waves of history matching performed in this work.See the text for details on the few-nucleon observables that are included in the target sets.The number of active inputs correspond to the dimensionality of the LEC subspace that is being explored in each wave.The second to last column indicates the fraction of input samples that passed the implausibility test, while the last column shows how large proportion of the initial volume that remains.

TABLE II .
Experimental values and error assignments for observables used in the fifth wave of the iterative history matching (history-matching observables) and for observables used in the model validation/calibration (predicted observables).Energies E in (MeV), point-proton radii rp in (fm), and the deuteron quadrupole moment Q in (e 2 fm 2 ).See the text for details.