Ab initio predictions link the neutron skin of 208Pb to nuclear forces

Heavy atomic nuclei have an excess of neutrons over protons, which leads to the formation of a neutron skin whose thickness is sensitive to details of the nuclear force. This links atomic nuclei to properties of neutron stars, thereby relating objects that differ in size by orders of magnitude. The nucleus 208Pb is of particular interest because it exhibits a simple structure and is experimentally accessible. However, computing such a heavy nucleus has been out of reach for ab initio theory. By combining advances in quantum many-body methods, statistical tools and emulator technology, we make quantitative predictions for the properties of 208Pb starting from nuclear forces that are consistent with symmetries of low-energy quantum chromodynamics. We explore 109 different nuclear force parameterizations via history matching, confront them with data in select light nuclei and arrive at an importance-weighted ensemble of interactions. We accurately reproduce bulk properties of 208Pb and determine the neutron skin thickness, which is smaller and more precise than a recent extraction from parity-violating electron scattering but in agreement with other experimental probes. This work demonstrates how realistic two- and three-nucleon forces act in a heavy nucleus and allows us to make quantitative predictions across the nuclear landscape.

Heavy atomic nuclei have an excess of neutrons over protons, which leads to the formation of a neutron skin whose thickness is sensitive to details of the nuclear force.This links atomic nuclei to properties of neutron stars, thereby relating objects that differ in size by orders of magnitude.The nucleus 208 Pb is of particular interest because it exhibits a simple structure and is experimentally accessible.However, computing such a heavy nucleus has been out of reach for ab initio theory.By combining advances in quantum many-body methods, statistical tools, and emulator technology, we make quantitative predictions for the properties of 208 Pb starting from nuclear forces that are consistent with symmetries of low-energy quantum chromodynamics.We explore 10 9 different nuclear-force parameterisations via history matching, confront them with data in select light nuclei, and arrive at an importance-weighted ensemble of interactions.We accurately reproduce bulk properties of 208 Pb and determine the neutron skin thickness, which is smaller and more precise than a recent extraction from parity-violating electron scattering but in agreement with other experimental probes.This work demonstrates how realistic two-and three-nucleon forces act in a heavy nucleus and allows us to make quantitative predictions across the nuclear landscape.
Neutron stars are extreme astrophysical objects whose interiors may contain exotic new forms of matter.The structure and size of neutron stars are linked to the thickness of neutron skins in atomic nuclei via the neutronmatter equation of state [1][2][3] .The nucleus 208 Pb is an attractive target for exploring this link in both experimental 4,5 and theoretical 2,6,7 studies, due to the large excess of neutrons and its simple structure.Mean-field calculations predict a wide range for R skin ( 208 Pb) because the isovector parts of nuclear energy density functionals are not well constrained by binding energies and charge radii 2,[7][8][9] .Additional constraints may be obtained 10 by including the electric dipole polarisability of 208 Pb, though this comes with a model dependence 11 which is difficult to quantify.In general, estimation of systematic theoretical uncertainties is a challenge for mean-field theory.
In contrast, precise ab initio computations, which provide a path to comprehensive uncertainty estimation, have been accomplished for the neutron-matter equation of state [12][13][14] and the neutron skin in the medium-mass nucleus 48 Ca 15 .But up to now treating 208 Pb within the same framework was out of reach.Due to breakthrough developments in quantum many-body methods, such computations are now becoming feasible for heavy nuclei [16][17][18][19] .The ab initio computation of 208 Pb we report here represents a significant step in mass number from the previously computed tin isotopes 16,17 , as illustrated in Figure 1.The complementary statistical analysis in this work is enabled by emulators (for mass number A ≤ 16) which mimic the outputs of many-body solvers, but are orders of magnitude faster.
In this paper we develop a unified ab initio framework to link the physics of nucleon-nucleon scattering and few-nucleon systems to properties of medium-and heavy-mass nuclei up to 208 Pb, and ultimately to the nuclear matter equation of state near saturation density.

Linking models to reality
Our approach to constructing nuclear interactions is based on chiral effective field theory (EFT) [21][22][23] .In this theory the long-range part of the strong nuclear force is known and stems from pion exchanges, while the unknown short-range contributions are represented as contact interactions; we also include the ∆ isobar degree of freedom 24 .At next-to-next-to leading order in Weinberg's power counting, the four pion-nucleon low-energy constants (LECs) are tightly fixed from pion-nucleon scattering data 25 .The 13 additional LECs in the nuclear potential must be constrained from data.
We use history matching 26,27 to explore the modeling capabilities of ab initio methods by identifying a nonimplausible region in the vast parameter space of LECs, for which the model output yields acceptable agreement with selected low-energy experimental data-here denoted history-matching observables.The key to efficiently analyze this high-dimensional parameter space arXiv:2112.01125v2[nucl-th] 22 Aug 2022 2015 2015 2016 78 Ni 2016 100 Sn 2018 132  The bars highlight years of first realistic computations of doubly magic nuclei.The height of each bar corresponds to the mass number A divided by the logarithm of the total compute power R TOP500 (in flop/second) of the pertinent TOP500 list 20 .This ratio would be approximately constant if progress were solely due to exponentially increasing computing power.However, algorithms which instead scale polynomially in A have greatly increased the reach.
is the use of emulators based on eigenvector continuation [28][29][30] that accurately mimic the outputs of the ab initio methods at several orders of magnitude lower computational cost.We consider the following historymatching observables: nucleon-nucleon scattering phase shifts up to an energy of 200 MeV; the energy, radius, and quadrupole moment of 2 H; and the energies and radii of 3 Here, experimental observations, z, are related to emulated ab initio predictions M (θ) via random variables ε exp , ε em , ε method , ε model that represent experimental uncertainties, emulator precision, method approximation errors, and the model discrepancy due to the EFT truncation at next-to-next-to leading order, respectively.The parameter vector θ corresponds to the 17 LECs at this order.The method error represents, e.g., model-space truncations and other approximations in the employed ab initio many-body solvers.The model discrepancy ε model can be probabilistically specified since we assume to op-erate with an order-by-order improvable EFT description of the nuclear interaction (see Methods for details).
The final result of the five history-matching waves is a set of 34 non-implausible samples in the 17-dimensional parameter space of the LECs.We then perform ab initio calculations for nuclear observables in 48 Ca and 208 Pb, as well as for properties of infinite nuclear matter.

Ab initio computations of 208 Pb
We employ the coupled-cluster (CC) 12,31,32 , the inmedium similarity renormalization group (IMSRG) 33 and many-body perturbation theory (MBPT) methods to approximately solve the Schrödinger equation and obtain the ground-state energy and nucleon densities of 48 Ca and 208 Pb.We analyze the model-space convergence and use the differences between CC, IMSRG and MBPT results to estimate the method approximation errors, see Methods and Extended Data Figures 3 and 4. The computational cost of these methods scales (only) polynomially with increasing numbers of nucleons and single-particle orbitals.The main challenge in computing 208 Pb is the vast number of matrix elements of the three-nucleon force which must be handled.We overcome this limitation by using a recently introduced storage scheme in which we only store linear combinations of matrix elements directly entering the normal-ordered two-body approximation 19 (see Methods for details).
Our ab initio predictions for finite nuclei are summarized in Figure 2. The statistical approach that leads to these results is composed of three stages.First, history matching identified a set of 34 non-implausible interaction parametrizations.Second, model calibration is performed by weighting these parametrizations-serving as prior samples-using a likelihood measure according to the principles of sampling/importance resampling 37 .This yields 34 weighted samples from the LEC posterior probability density function, see Extended Data Figure 5. Specifically we assume independent EFT and many-body method errors and construct a normally distributed datalikelihood encompassing the ground-state energy per nucleon E/A and point-proton radius R p for 48 Ca, and the energy E 2 + of its first excited 2 + state.Our final predictions are therefore conditional on this calibration data.
We have tested the sensitivity of final results to the likelihood definition by repeating the calibration with a non-diagonal covariance matrix or a Student-t distribution with heavier tails, finding small (∼ 1%) differences in the predicted credible regions.The EFT truncation errors are quantified by studying ab initio predictions at different orders in the power counting for 48 Ca and infinite nuclear matter.We validate our ab initio model and error assignments by computing the posterior predictive distributions, including all relevant sources of uncertainty, for both the replicated calibration data (blue colour) and the history-matching observables (green colour), see Figure 2. The percentage ratios σ tot /z of the (theory dominated) total uncertainty to the  1 for numerical specification of experimental data (z), errors (σ i ), medians (white circle) and 68% credibility regions (thick bar).The prediction for R skin ( 208 Pb) in the bottom panel is shown in an absolute scale and compared to experimental results using electroweak 5 (purple), hadronic 34,35 (red), electromagnetic 4 (green), and gravitational waves 36 (blue) probes (from top to bottom; see Extended Data Figure 7b for details).experimental value are given in the right margin.
Finally, having built confidence in our ab initio model and underlying assumptions, we predict R skin ( 208 Pb), E/A and R p for 208 Pb, α D for 48 Ca and 208 Pb as well as nuclear matter properties, by employing importance resampling 37 .The corresponding posterior predictive distributions are shown in the lower panels of Figure 2 (pink colour).Our prediction R skin ( 208 Pb) = 0.14 − 0.20 fm exhibits a mild tension with the value extracted from the recent parity-violating electron scattering experiment PREX 5 but is consistent with the skin thickness extracted from elastic proton scattering 35 , antiprotonic atoms 34 and coherent pion photoproduction 4 as well as constraints from gravitational waves from merging neutron stars 36 .
We also compute the weak form factor F w (Q 2 ) at momentum transfer Q PREX = 0.3978(16) fm −1 , which is more directly related to the parity-violating asymmetry measured in the PREX experiment.We observe a strong correlation with the more precisely measured electric charge form factor F ch (Q 2 ), as shown in Figure 3b.While we have not quantified the EFT and method errors for these observables, we find a small variance among the 34 non-implausible predictions for the difference F w (Q 2 ) − F ch (Q 2 ) for both 48 Ca and 208 Pb as shown in Figure 3c.

Ab initio computations of infinite nuclear matter
We also make predictions for nuclear matter properties by employing the CC method on a momentum-space lattice 38 with a Bayesian machine-learning error model to quantify the uncertainties from the EFT truncation 14 and the CC method (see Methods and Extended Data Figure 6 for details).The observables we compute are the saturation density ρ 0 , the energy per nucleon of symmetric nuclear matter E 0 /A, its compressibility K, the symmetry energy S (i.e. the difference between the energy per nucleon of neutron matter and symmetric nuclear matter), and its slope L. The posterior predictive distributions for these observables are shown in Figure 3a.These distributions include samples from the relevant method and model error terms.Overall, we reveal relevant correlations among observables, previously indicated in mean-field models, and find good agreement with empirical bounds 39 .The last row shows the resulting correlations with R skin ( 208 Pb) in our ab initio framework.In particular, we find essentially the same correlation between R skin ( 208 Pb) and L as observed in mean-field models (See Extended Data Figure 7b).

Discussion
The predicted range of the 208 Pb neutron skin thickness (see Extended Data Table 2) is consistent with several extractions 4,40,41 , each of which involves some modeldependence, and in mild tension (approximately 1.5σ) with the recent PREX result 5 .Ab initio computations yield a thin skin and a narrow range because the isovector physics is constrained by scattering data 8,13,42 .A thin skin was also predicted in 48 Ca 15 .
We find that both R skin ( 208 Pb) = 0.14 − 0.20 fm and the slope parameter L = 37 − 66 MeV are strongly correlated with scattering in the 1 S 0 partial wave for laboratory energies around 50 MeV (the strongest two-neutron channel allowed by the Pauli principle, with the energy naively corresponding to the Fermi energy of neutron matter at 0.8ρ 0 ), see Extended Data Figure 7a.It is possible, analogous to findings in mean-field theory 1,43 , to increase L beyond the range predicted in this work by tuning a contact in the 1 S 0 partial wave and simultaneously readjusting the three-body contact to maintain realistic nuclear saturation.But this large slope L and increased R skin come at the cost of degraded 1 S 0 scattering phase shifts, well beyond the expected corrections from higher-order terms (see Extended Data Figure 8).The large range of L and R skin obtained in mean-field theory is a consequence of scattering data not being incorporated.It will be important to confront our predictions with more precise experimental measurements 44,45 .
If the tension between scattering data and neutron skins persists, it will represent a serious challenge to our ab initio description of nuclear physics.
Our work demonstrates that ab initio approaches using nuclear forces from chiral EFT can consistently describe data from nucleon-nucleon scattering, few-body systems, and heavy nuclei within the estimated theoretical uncertainties.Information contained in nucleon-nucleon scattering significantly constrains the properties of neutron matter; this same information constrains neutron skins, which provide a non-trivial empirical check on the reliability of ab initio predictions for the neutron matter equation of state.Moving forward, it will be important to extend these calculations to higher orders in the effective field theory, both to further validate the error model and to improve precision, and to push the cutoff to higher values to confirm regulator independence.The framework presented in this work will enable predictions with quantified uncertainties across the nuclear chart, advancing toward the goal of a single unified framework for describing low energy nuclear physics.

METHODS
Hamiltonian and model space.The many-body approaches used in this work [CC, IMSRG, and manybody perturbation theory (MBPT)] start from the intrinsic Hamiltonian Here T kin is the kinetic energy, T CoM is the kinetic energy of the center of mass, V NN is the nucleon-nucleon, and V 3N is the three-nucleon interaction.In order to facilitate the convergence of heavy nuclei, the interactions employed in this work used a non-local regulator with a cutoff Λ = 394 MeV/c.Specifically, the Results should be independent of this choice, up to higherorder corrections, provided renormalization-group invariance of the EFT.However, increasing the momentum scale of the cutoff leads to harder interactions, considerably enlarging the required computational effort.We represent the 34 non-implausible interactions that resulted from the history-matching analysis in the Hartree-Fock basis in a model-space of up to 15 major harmonic oscillator shells (e = 2n + l ≤ e max = 14 where n and l denote the radial and orbital angular momentum quantum numbers, respectively) with oscillator frequency ω = 10 MeV.Due to storage limitations, the three-nucleon force had an additional energy cut given by e 1 + e 2 + e 3 ≤ E 3max = 28.After obtaining the Hartree-Fock basis for each of the 34 non-implausible interactions, we capture 3N force effects via the normal-ordered two-body approximation before proceeding with the CC, IMSRG and MBPT calculations 46,47 .The convergence behaviour in e max and E 3max is illustrated in Extended Data Figure 3.In that figure, we use an interaction with a high likelihood that generates a large correlation energy.Thus, its convergence behaviour represents the worst case among the 34 non-implausible interactions.The model-space converged results are investigated with E 3max → 3e max and e max → ∞ extrapolations.The functions 2(e max + 7/2)b (b is the harmonic-oscillator length; c i s and d i s are the fitting parameters) are used as the asymptotic forms for E 3max 19 and e max 48,49 , respectively.Through the extrapolations, the ground-state energies computed with e max = 14 and E 3max = 28 is shifted by −75 ± 60 MeV.Likewise, the extrapolations of proton and neutron radii with the functional form given in Refs. 19,48,49yield a small +0.005 ± 0.010 fm shift of the neutron skin thickness.
In-medium similarity renormalization group calculations.The IMSRG calculations 33,50 were performed at the IMSRG(2) level, using the Magnus formulation 51 .Operators for the point-proton and pointneutron radii, form factors, and the electric dipole operator were consistently transformed.The dipole po-larizablility α D was computed using the equations-ofmotion method truncated at the 2-particle-2-hole level (the EOM-IMSRG(2,2) approximation 52 ) and the Lanczos continued fraction method 53 .We compute the weak and charge form factors using the parameterization presented in Ref. 54 , though the form given in Ref. 55 yields nearly identical results.
Many-body perturbation theory calculations.MBPT theory calculations for 208 Pb were performed in the Hartree-Fock basis to third order for the energies, and to second order for radii.
Coupled-cluster calculations.The CC calculations of 208 Pb were truncated at the singles-and-doubles excitation level, known as the CCSD approximation 12,31,32 .We estimated the contribution from triples excitations to the ground-state energy of 208 Pb as 10% of the CCSD correlation energy (which is a reliable estimate for closedshell systems 32 ).
Extended Data Figure 4 compares the different manybody approaches used in this work, i.e.CC, IMSRG, MBPT, and allows us to estimate the uncertainties related to our many-body approach in computing the ground-state observables for 208 Pb.The point proton and neutron radii are computed as ground-state expectation values (see e.g.Ref. 15 for details).For 48 Ca we used a Hartree-Fock basis consisting of 15 major oscillator shells with an oscillator spacing of ω = 16 MeV, while for 3N forces we used E 3max = 16, which is sufficiently large to obtain converged results in this mass region.Here we computed the ground-state energy using the Λ-CCSD(T) approximation 56 which include perturbative triples corrections.The 2 + excited state in 48 Ca was computed using the equation-of-motion CCSD approach 57 , and we estimated a −1 MeV shift from triples excitations based on EOM-CCSD(T) calculations of 48 Ca and 78 Ni using similar interactions 58 .
For the history-matching analysis we used an emulator for the 16 O ground-state energy and charge radius that was constructed using the recently developed subspace coupled-cluster method 30 .For higher precision in the emulator we went beyond the SP-CCSD approximation used in Ref. 30 and included leading-order triples excitations via the CCSDT-3 method 59 .The CCSDT-3 ground-state training vectors for 16 O were obtained starting from the Hartree-Fock basis of the recently developed chiral interaction ∆NNLO GO (394) of Ref. 60 in a model-space consisting of 11 major harmonic oscillator shells with the oscillator frequency ω = 16 MeV, and E 3max = 14.The emulator used in the history matching was constructed by selecting 68 different training points in the 17-dimensional space of LECs using a spacefilling Latin hypercube design in a 10% variation around the ∆NNLO GO (394) LECs.At each training point we then performed a CCSDT-3 calculation in order to obtain the training vectors for which we then construct the sub-space projected norm and Hamiltonian matrices.Once the SP-CCSDT-3 matrices are constructed we may obtain the ground-state energy and charge radii for any target values of the LECs by diagonalizing a 68 by 68 generalized eigenvalue problem (see Ref. 30 for more details).We checked the accuracy of the emulator by cross-validation against full-space CCSDT-3 calculations as demonstrated in Extended Data Figure 4a and found a relative error that was smaller than 0.2%.
The nuclear matter equation of state and saturation properties are computed with the CCD(T) approximation which includes doubles excitations and perturbative triples corrections.The three-nucleon forces are considered beyond the normal-ordered two-body approximation by including the residual three-nucleon force contribution in the triples.The calculations are performed on a cubic lattice in momentum space with periodic boundary conditions.The model space is constructed with (2n max +1) 3 momentum points, and we use n max = 4(3) for pure neutron matter (symmetric nuclear matter) and obtain converged results.We perform calculations for systems of 66 neutrons (132 nucleons) for pure neutron matter (symmetric nuclear matter) since results obtained with those particle numbers exhibit small finite size effects 38 .
Iterative history matching.In this work we use an iterative approach known as history matching 26,27 in which the model, solved at different fidelities, is confronted with experimental data z using Eq. ( 1).Obviously, we do not know the exact values of the errors in Eq. ( 1), hence we represent them as random variables and specify reasonable forms for their statistical distributions, in alignment with the Bayesian paradigm.
For many-body systems we employ quantified method and (A = 16) emulator errors as discussed above and summarized in Extended Data Table 1.For A ≤ 4 nuclei we use the no-core shell model in Jacobi coordinates 61 and eigenvector continuation emulators 29 .The associated method and emulator errors are very small.Probabilistic attributes of the model discrepancy terms are assigned based on the expected EFT convergence pattern 62,63 .For the history-matching observables considered here we use point estimates of model errors from Ref. 64 .
The aim of history matching is to estimate the set Q(z) of parameterizations θ, for which the evaluation of a model M (θ) yields an acceptable-or at least not implausible-match to a set of observations z.History matching has been employed in various studies involving complex computer models [65][66][67][68] ranging, e.g., from effects of climate modeling 69,70 to systems biology 71 .
We introduce the individual implausibility measure which is a function over the input parameter space and quantifies the (mis-)match between our (emulated) model output M i (θ) and the observation z i for an observable in the target set Z. We mainly employ a maximum implausibility measure as the restricting quantity.Specifically, we consider a particular value for θ as implausible if with c I ≡ 3.0 appealing to Pukelheim's three-sigma rule 72 .In accordance with the assumptions leading to Eq. ( 1), the variance in the denominator of Eq. ( 3) is a sum of independent squared errors.Generalizations of these assumptions are straightforward if additional information on error covariances or possible inaccuracies in our error model would become available.An important strength of the history matching is that we can proceed iteratively, excluding regions of input space by imposing cutoffs on implausibility measures that can include additional observables z i and corresponding model outputs M i with possibly refined emulators as the parameter volume is reduced.The history matching process is designed to be independent of the order in which observables are included, as is discussed in 67 .This is an important feature as it allows for efficient choices regarding such orderings.The iterative history matching proceeds in waves according to a straightforward strategy that can be summarized as follows: 1.At wave j: Evaluate a set of model runs over the current NI volume Q j using a space-filling design of sample values for the parameter inputs θ.Choose a rejection strategy based on implausibility measures for a set Z j of informative observables.
2. Construct or refine emulators for the model predictions across Q j .
3. The implausibility measures are then calculated over Q j using the emulators, and implausibility cutoffs are imposed.This defines a new, smaller non-implausible volume Q j+1 which should satisfy 4. Unless (a) computational resources are exhausted, or (b) all considered points in the parameter space are deemed implausible, we may include additional informative observables in the considered set Z j+1 , and return to step 1.

If 4(a
) is true we generate a number of acceptable runs from the final non-implausible volume Q final , sampled according to scientific need.
The ab initio model for the observables we consider comprises at most 17 parameters; four subleading pionnucleon couplings, 11 nucleon-nucleon contact couplings, and two short-ranged three-nucleon couplings.To identify a set of non-implausible parameter samples we performed iterative history matching in four waves using observables and implausibility measures as summarized in Extended Data Figure 1b.For each wave we employ a sufficiently dense Latin hypercube set of several million candidate parameter samples.For the model evaluations we utilized fast computations of neutron-proton scattering phase shifts and efficient emulators for the few-and many-body history-matching observables.See Extended Data Table 1 and Extended Data Figure 2 for the list of history-matching observables and information on the errors that enter the implausibility measure (3).The input volume for wave 1 incorporates the naturalness expectation for LECs, but still includes large ranges for the relevant parameters as indicated by the panel ranges in Extended Data Figure 1a.In all four waves the input volume for c 1,2,3,4 is a four-dimensional hypercube mapped onto the multivariate Gaussian probability density function (PDF) resulting from a Roy-Steiner analysis of pion-nucleon scattering data 73 .In wave 1 and wave 2 we sampled all relevant parameter directions for the set of included two-nucleon observables.In wave 3, the 3 H and 4 He observables were added such that the three-nucleon force parameters c D and c E can also be constrained.Since these observables are known to be rather insensitive to the four model parameters acting solely in P −waves, we ignored this subset of the inputs and compensated by slightly enlarging the corresponding method errors.This is a well known emulation procedure called inactive parameter identification 26 .For wave 4 we considered all 17 model parameters and added the ground-state energy and radius of 16 O to the set Z 4 and emulated the model outputs for 5 × 10 8 parameter samples.By including oxygen data we explore the modeling capabilities of our ab initio approach.Extended Data Figure 1a summarizes the sequential non-implausible volume reduction, wave-by-wave, and indicates the set of 4,337 non-implausible samples after the fourth wave.We note that the use of history matching would in principle allow a detailed study of the information content of various observables in heavy-mass nuclei.Such an analysis, however, requires an extensive set of reliable emulators and is beyond the scope of the present work.The volume reduction is determined by the maximum implausibility cutoff (4) with additional confirmation from the optical depths (which indicate the density of non-implausible samples; see Eqs. ( 25) and ( 26) in Ref. 71 ).The nonimplausible samples summarise the parameter region of interest, and can directly aid insight regarding interdependencies between parameters induced by the match to observed data.This region is also where we would expect the posterior distribution to reside and we note that our history-matching procedure has allowed us to reduce its size by more than seven orders of magnitude compared to the prior volume (see Extended Data Figure 1b).
As a final step, we confront the set of non-implausible samples from wave 4 with neutron-proton scattering phase shifts such that our final set of non-implausible samples has been matched with all history-matching observables.For this final implausibility check we employ a slightly less strict cutoff and allow the first, second and third maxima of I i (θ) (for z i ∈ Z final ) to be 5.0, 4.0, and 3.0, respectively, accommodating the more extreme max-ima we may anticipate when considering a significantly larger number of observables.The end result is a set of 34 non-implausible samples that we use for predicting 48 Ca and 208 Pb observables, as well as the equation of state of both symmetric nuclear matter and pure neutron matter.
Posterior predictive distributions.The 34 nonimplausible samples from the final history matching wave are used to compute energies, radii of proton and neutron distributions, and electric dipole polarizabilities (α D )for 48 Ca and 208 Pb.They are also used to compute the electric and weak charge form factors for the same nuclei at a relevant momentum transfer, and the energy per particle of infinite nuclear matter at various densities to extract key properties of the nuclear equation of state (see below).These results are shown as blue circles in Figure 3.
In order to make quantitative predictions, with a statistical interpretation, for R skin ( 208 Pb) and other observables we use the same 34 parameter sets to extract representative samples from the posterior PDF p(θ|D cal ).Bulk properties (energies and charge radii) of 48 Ca together with the structure-sensitive 2 + excited-state energy of 48 Ca are used to define the calibration data set D cal .The IMSRG and CC convergence studies make it possible to quantify the method errors.These are summarized in Extended Data Table 1.The EFT truncation errors are quantified by adopting the EFT convergence model 74,75 for observable y with observable coefficients c i that are expected to be of natural size, and the expansion parameter Q = 0.42 following our Bayesian error model for nuclear matter at the relevant density (see below).The first sum in the parenthesis is the model prediction y k (θ) of observable y at truncation order k in the chiral expansion.The second sum than represents the model error as it includes the terms that are not explicitly included.We can quantify the magnitude of these terms by learning about the distribution for c i which we will assume is described by a single normal distribution per observable type with zero mean and a variance parameter c2 .We employ the nuclear matter error analysis for the energy per particle of symmetric nuclear matter (described below) to provide the model error for E/A in 48 Ca and 208 Pb.For radii and electric dipole polarizabilities we employ the next-to leading order and next-to-next-to leading order interactions of Ref. 60 and compute these observables at both orders for various Ca, Ni, and Sn isotopes.The reference values y ref are set to r 0 • A 1/3 for radii and to the experimental value for α D .From this data we extract c2 and perform the geometric sum of the second term in Eq. ( 5).The resulting standard deviations for model errors are summarized in Extended Data Table 1.
At this stage we can approximately extract samples from the parameter posterior p(θ|D cal ) by employing the established method of sampling/importance resampling 37,76 .We assume a uniform prior probability for the non-implausible samples and we introduce a normallydistributed likelihood, L(D cal |θ), assuming independent experimental, method, and model errors.The prior for c 1,2,3,4 is the multivariate Gaussian resulting from a Roy-Steiner analysis of πN scattering data 73 .Defining importance weights we draw samples θ * from the discrete distribution {θ 1 , . . ., θ n } with probability mass q i on θ i .These samples are then approximately distributed according to the parameter posterior that we are seeking 37,76 .
Although we are operating with a finite number of n = 34 representative samples from the parameter PDF, it is reassuring that about half of them are within a factor two from the most probable one in terms of the importance weight, see Extended Data Figure 5. Consequently, our final predictions will not be dominated by a very small number of interactions.In addition, as we do not anticipate the parameter PDF to be of a particularly complex shape, based on the results of the history match, consideration of the various error structures in the analysis, and on the posterior predictive distributions (PPDs) shown in Figure 3, and as we are mainly interested in examining such lower 1-or 2-dimensional PPDs, this sample size was deemed sufficient and the corresponding sampling error assumed subdominant.We use these samples to draw corresponding samples from This PPD is the set of all model predictions computed over likely values of the parameters, i.e., drawing from the posterior PDF for θ.The full PPD is then defined, in analogy with Eq. ( 7), as the set evaluation of y which is the sum where we assume method and model errors to be independent of the parameters.In practice, we produce 10 4 samples from this full PPD for y by resampling the 34 samples of the model PPD (7) according to their importance weights, and adding samples from the error terms in (8).We perform model checking by comparing this final PPD with the data used in the iterative historymatching step, and in the likelihood calibration.In addition, we find that our predictions for the measured electric dipole polarizabilities of 48 Ca and 208 Pb as well as bulk properties of 208 Pb serve as a validation of the reliability of our analysis and assigned errors.See Figure 2 and Extended Data Table 1.
In addition, we explored the sensitivity of our results to modifications of the likelihood definition.Specifically, we used a student-t distribution (ν = 5) to see the effects of allowing heavier tails, and we introduced an error covariance matrix to study the effect of possible correlations (with ρ ≈ 0.7) between the errors in binding energy and radius of 48 Ca.In the end, the differences in the extracted credibility regions was ∼ 1% and we therefore present only results obtained with the uncorrelated, multivariate normal distribution.
Our final predictions for R skin ( 208 Pb), R skin ( 48 Ca) and for nuclear matter properties are presented in Figure 3 and Extended Data Table 2.For these observables we use the Bayesian machine learning error model described below to assign relevant correlations between equationof-state observables.For model errors in R skin ( 208 Pb) and L we use a correlation coefficient of ρ = 0.9 as motivated by the strong correlation between the observables computed with the 34 non-implausible samples.It should be noted that S, L, and K are computed at the specific saturation density of the corresponding non-implausible interaction.
Bayesian machine learning error model.Similar to Eq. ( 1) the predicted nuclear matter observables can be written as: where y k (ρ) is the CC prediction using our EFT model truncated at order k, ε k (ρ) is the EFT truncation (model) error, and ε method (ρ) is the CC method error.In this work we apply a Bayesian machine learning error model 14 to quantify the density dependence of both method and truncation errors.The error model is based on multitask Gaussian processes that learn both the size and the correlations of the target errors from given prior information.Following a physically-motivated Gaussian process (GP) model 14 , the EFT truncation errors ε k at given density ρ are distributed as: with Here k = 3 for the ∆NNLO(394) EFT model used in this work, while c2 , l and Q are hyperparameters corresponding to the variance, the correlation length, and the expansion parameter.Finally, we choose the reference scale y ref to be the EFT leading-order prediction.The mean of the Gaussian process is set to be zero since the order-by-order truncation error can either be positive or negative and the correlation function r(ρ, ρ ; l) in ( 11) is the Gaussian radial basis function.
We employ Bayesian inference to optimize the Gaussian process hyperparameters using order-by-order predictions of the equation of state for both pure neutron matter and symmetric nuclear matter with the ∆-full interactions from Ref. 64 .In this work, we find cPNM = 1.00 and l PNM = 0.92 fm −1 for pure neutron matter and cSNM = 1.55 and l SNM = 0.48 fm −1 for symmetric nuclear matter.
The above Gaussian processes only describe the correlated structure of truncation errors for one type of nucleonic matter.In addition, the correlation between pure neutron matter and symmetric nuclear matter is crucial for correctly assigning errors to observables that involve both E/N and E/A (such as the symmetry energy S).For this purpose we use a multitask Gaussian process that simultaneously describes truncation errors of pure neutron matter and symmetric nuclear matter according to: where K 11 and K 22 are the covariance matrices generated from the kernel function c2 R ε k (ρ, ρ ; l) for pure neutron matter and symmetric nuclear matter, respectively, while K 12 (K 21 ) is the cross-covariance as in Ref. 77 .
Regarding the CC method error, different sources of uncertainty should be considered.The truncation of the cluster operator and the finite-size effect are the main ones and the total method error is then ε method = ε cc + ε fs .Following the Bayesian error model we have a general expression for the method error: with R me (ρ, ρ ; l me ) = y me,ref (ρ)y me,ref (ρ )r(ρ, ρ ; l me ).( 14) Here the subscript "me" stands for either the cluster operator truncation "cc" or finite-size effect "fs" method error.For the cluster operator truncation errors ε cc the reference scale y me,ref is taken to be the CCD(T) correlation energy.The Gaussian processes are then optimized with data from different interactions by assuming that the energy difference between CCD and CCD(T) can be used as an approximation of the cluster operator truncation error.The correlation lengths learned from the training data are l me,PNM = 0.81 fm −1 for pure neutron matter and l me,SNM = 0.34 fm −1 for symmetric nuclear matter.Based on the convergence study we take ±10% of the correlation energy as the 95% credible interval which gives cme = 0.05 for ε cc .As for the finite-size effect ε fs , the reference scale is taken to be the CCD(T) ground-state energy.Then following Ref. 38, we use ±0.5% (±4%) of the ground-state energy of the pure neutron matter (the symmetric nuclear matter) as a conservative estimation of the finite-size effect (95% credible interval) when using periodic boundary conditions with 66 neutrons (132 nucleons) around the saturation point.This leads to cme,PNM = 0.0025 and cme,SNM = 0.02 for ε fs .The finite-size effects of different densities are clearly correlated while there are insufficient data to learn its correlation structure.Here we simply used 0.81 fm −1 (0.34 fm −1 ) as the correlation length for pure neutron matter (symmetric nuclear matter) and assume zero correlation between pure neutron matter and symmetric nuclear matter.
Once the model and method errors are determined, it is straightforward to sample these errors from the corresponding covariance matrix and produce the equationof-state predictions using Eq. ( 9) for any given interaction.This sampling procedure is crucial for generating the posterior predictive distribution of nuclear matter observables shown in Figure 3a.CCD(T) calculations for nuclear-matter equation of state and the corresponding 2σ credible interval for method and model errors are illustrated in Extended Data Figure 6.The sampling procedure is made explicit with three randomly sampled equation-of-state predictions.Note that even though the sampled errors for one given density appear to be random, the multitask Gaussian processes will guarantee that the sampled equation of state of nuclear matter are smooth and properly correlated with each other.Finally, the proportion of the parameter space deemed non implausible is listed in the last column.Note that no additional reduction of the non-implausible domain is achieved in the fourth and final waves, in which 16 O observables are included, but that parameter correlations are enhanced.Histogram of importance weights for the 34 non-implausible interaction samples.These are obtained from likelihood calibration as defined in Eq. (6).shown along with the corresponding method error (blue shade) and EFT truncation error (green shade) for one representative interaction.Errors are correlated as a function of density ρ and the dashed orange, green and purple curves illustrate predictions with randomly sampled method and model errors drawn from the respective multitask Gaussian processes.Correlations extend between pure neutron matter (E/N ) and symmetric nuclear matter (E/A) energies per particle which is represented here by curves in the same colour., p34 , and GW170817 36 .All these results involve modeling input as the neutron skin thickness cannot be measured directly.The quoted experimental error bars include statistical and some systematical uncertainties except for Ref. 34 that is statistical only and the GW170817 constraint which is a 90 % upper bound from relativistic mean-field modeling of the tidal polarizability extracted in Ref. 36 .
Extended Data Table 2 | Predictions for the nuclear equation of state at saturation density and for neutron skins.Medians and 68%, 90% credible regions (CR) for the final PPD including samples from the error models (see also Figure 3 and text for details).The saturation density, ρ 0 , is in (fm −3 ), the neutron skin thickness, R skin ( 208 Pb) and R skin ( 48 Ca), in (fm), while the saturation energy per particle (E 0 /A), the symmetry energy (S), its slope (L), and incompressibility (K) at saturation density are all in (MeV).Empirical regions shown in Figure 3 are E 0 /A = −16.0± 0.5, ρ 0 = 0.16 ± 0.01, S = 31 ± 1, L = 50 ± 10 and K = 240 ± 20 from Refs. 39,88,89 4), Extended Data Table 1 and Extended Data Figure 2. b, Illustration of the freedom in Skyrme parametrizations to adjust L while preserving ρ 0 and E 0 /A.The parameters x 0 , t 0 , x 3 , t 3 correspond to the functional form given in e.g. 90 .The black circles correspond to different parameter sets, while the red line indicates the result of starting with the SKX interaction and modifying the x 0 , x 3 parameters while maintaining the binding energy per nucleon E/A of 208 Pb.The right column also shows the 208 Pb point-proton and point-neutron radii (R p and R n , respectively) and neutron skin thickness R skin for different parametrisations.The gray bands indicate a linear fit to the black points with r 2 the coefficient of determination.Skyrme parameter sets included are SKX, SKXCSB 9 , SKI, SKII 91 , SKIII-VI 92 , SKa, SKb 93 , SKI2, SKI5 94 , SKT4, SKT6 95 , SKP 96 , SGI, SGII 97 , MSKA 98 , SKO 99 , SKM * 100 .

Figure 1 |
Figure 1 | Trend of realistic ab initio computations for the nuclear A-body problem.The bars highlight years of first realistic computations of doubly magic nuclei.The height of each bar corresponds to the mass number A divided by the logarithm of the total compute power R TOP500 (in flop/second) of the pertinent TOP500 list20 .This ratio would be approximately constant if progress were solely due to exponentially increasing computing power.However, algorithms which instead scale polynomially in A have greatly increased the reach.

Figure 3 |
Figure3| Posterior predictive distribution for R skin ( 208 Pb) and nuclear matter at saturation density.a, Predictions for the saturation energy per particle E 0 /A and density ρ 0 of symmetric nuclear matter, its compressibility K, the symmetry energy S, and its slope L are correlated with the those for R skin ( 208 Pb).The bivariate distribution include 68% and 90% credible regions (black lines) and a scatter plot of the predictions with the 34 non-implausible samples before error sampling.Empirical nuclear-matter properties are indicated by purple bands (see Extended Data Table2).b, Predictions with the 34 non-implausible samples for the electric F ch versus weak F w charge form factors for 208 Pb at the momentum transfer considered in the PREX experiment5 .c, Difference between electric and weak charge form factors for48 Ca and 208 Pb at the momentum transfers Q CREX = 0.873 fm −1 and Q PREX = 0.3978 fm −1 that are relevant for the CREX and PREX experiments, respectively.Experimental data (purple bands) in panels b and c are from Ref.5 , the size of the markers indicate the importance weight, and blue lines correspond to weighted means.

1 |
History-matching waves.a, The initial parameter domain used at the start of history-matching wave 1 is represented by the axes limits for all panels.This domain is iteratively reduced and the input volumes of waves 2, 3, and 4 are indicated by green/dash-dotted, blue/dashed, black/solid rectangles.The logarithm of the optical depths log 10 ρ (indicating the density of non-implausible samples in the final wave) are shown in red with darker regions corresponding to a denser distribution of non-implausible samples.b, Four waves of history matching were used in this work plus a fifth one to refine the final set of non-implausible samples.The neutron-proton scattering targets correspond to phase shifts at six energies (T lab = 1, 5, 25, 50, 100, 200 MeV) per partial wave: 1 S 0 , 3 S 1 , 1 P 1 , 3 P 0 , 3 P 1 , 3 P 2 .The A = 2 observables are E( 2 H), R p ( 2 H), Q( 2 H), while A = 3, 4 are E( 3 H), E( 4 He), R p ( 4 He).Finally, A = 16 targets are E( 16 O), R p ( 16 O).The number of active input parameters is indicated in the fourth column.The number of inputs sets being explored, and the fraction of non-implausible samples that survive the imposed implausibility cutoff(s) are shown in the fifth and sixth columns, respectively.

Extended Data Figure 2 | 14 Extended Data Figure 3 |
Neutron-proton scattering phase shifts.34 interaction samples survive the final implausibility cutoff with respect to neutron-proton phase shifts δ in S and P waves up to 200 MeV.The red circles are from the Granada phase shift analysis79 , while the 2σ error bars are dominated by the estimated EFT truncation errors64 .Convergence of energy and radius observables of 208 Pb with the e max and E 3max truncations.a, Ground state energy as a function of E 3max .The dashed lines indicate a Gaussian fit.b, Ground state energy (extrapolated in E 3max as a function of e max .The smaller error bar on the adopted value indicate the error due to model space extrapolation, and the larger error bar also includes the method uncertainty.c, Neutron skin as a function of E 3max .d, Neutron radius as a function of oscillator basis frequency ω for a series of e max cuts.

4 |Extended Data Figure 5 |
Precision of sub-space coupled-cluster emulator and many-body solvers.a, Cross-validation of the SP-CCSDT-3 emulator for the ground-state energy of 16 O.Results from full computations using CCSDT-3 are compared with emulator predictions for 50 samples from the 17 dimensional space of LECs.The standard deviation for the residuals ∆E SPCC−CC is 0.19 MeV.b,c,d, Differences between IMSRG and CC results versus differences between MBPT and CC results for the ground-state energy per nucleon ∆E/A (panel b), the point-proton radius ∆R p (panel c), and the neutron-skin ∆R skin (panel d) of 208 Pb using the 34 non-implausible interactions obtained from history matching (see text for more details).The CC results for the ground-state energy include approximate triples corrections.Importance weights.

7 |
Correlation of R skin ( 208 Pb) with scattering data and L. a Correlation of computed R skin ( 208 Pb) with the proton-neutron 1 S 0 phase shift δ( 1 S 0 ) at a laboratory energy of 50 MeV, shown in blue.The error bars represent method and model (EFT) uncertainties.The green band indicates the experimental phase shift 79 , while the purple line (band) indicate the mean result (one-sigma error) of the PREX experiment 5 .The dashed line indicates the linear trend of the ab initio points with r 2 the coefficient of determination.b Correlation of neutron skin R skin ( 208 Pb) vs slope of the symmetry energy L. Relativistic and non-relativistic mean-field calculations are indicated with open symbols 87 , while ab initio results using the 34 non-implausible samples are indicated with filled circles.Experimental extractions of R skin ( 208 Pb) shown in the figure are from PREX 5 , MAMI 4 , RCNP

64 Extended Data Figure 8 |
Parameter sensitivities in ab initio models and Skyrme parametrizations.a, Tuning the C 1S0 LEC in our ab initio model to adjust the energy slope parameter L while compensating with the three-nucleon contact c E to maintain the saturation density ρ 0 and energy per nucleon E 0 /A of symmetric nuclear matter.The green pentagons correspond to results with one of the 34 interaction samples while the black squares indicate the results after tuning the C 1S0 and c E of that interaction.The right column shows the scattering phase shift δ in the 1 S 0 channel at 50 MeV, the ground-state energies in 3 H and16 O and the point-proton radius R p in16 O.The red diamonds and the dashed lines indicate the experimental values of target observables and the red bands indicate the corresponding c I = 3 non-implausible regions, see Eq. ( H, 4 He, and 16 O.We perform five waves of this global parameter search-see Extended Data Figures 1 and 2-sequentially ruling out implausible LECs that yield model predictions too far from experimental data.For this purpose we use an implausibility measure (see Methods) that links our model predictions and experimental observations as