Nuclear-matter saturation and symmetry energy within $\Delta$--full chiral effective field theory

Nuclear saturation and the symmetry energy are key properties of low-energy nuclear physics that depend on fine details of the nuclear interaction. The equation-of-state around saturation is also an important anchor for extrapolations to higher densities and studies of neutron stars. Here we develop a unified statistical framework that uses realistic nuclear forces to link the theoretical modeling of finite nuclei and infinite nuclear matter. We construct fast and accurate emulators for nuclear-matter observables and employ an iterative history-matching approach to explore and reduce the enormous parameter domain of $\Delta$-full chiral interactions. We perform rigorous uncertainty quantification and find that model calibration including \nuc{16}{O} observables gives saturation predictions that are more precise than those that only use few-body data.

Introduction-A key question in low-energy nuclear physics is whether it is possible to successfully describe all systems from finite nuclei to infinite nuclear matter using nucleons as effective degrees of freedom.Realistic interaction models based on chiral effective field theory (χEFT) used in combination with ab initio methods, that solve the many-body problem with controlled approximations, have the potential to deliver on that research program [1][2][3][4][5][6][7][8][9].However, using a complicated interaction model and calibration data with limited independence implies a risk for overfitting unless relevant theoretical uncertainties are accounted for.Furthermore, the extrapolation from well-studied few-nucleon systems to heavier nuclei and infinite nuclear matter will lead to increasing variances in model predictions.Therefore it becomes important to study the precision of these predictions and to explore the sensitivity to different choices of calibration data.The emergence of nuclear saturation-represented by a minimum in the equation-of-state (EOS) for infinite symmetric nuclear matter (SNM)-is particularly important since it affects bulk properties (such as binding energies and radii) of finite atomic nuclei.Furthermore, the density dependence of the EOS for pure neutron matter (PNM) is central for the physics of neutron stars [10][11][12].
In this Letter we develop fast and accurate emulators of coupled-cluster computations of PNM and SNM [29,30] starting from ∆-full χEFT interaction models at NNLO [5,[31][32][33][34][35][36].Our construction of emulators is based on the subspace-projected coupled cluster (SPCC) method [25], which we here extend with small batch voting [37]-a new on-the-fly validation approach that addresses the possible appearance of spurious states when diagonalizing a Hamiltonian in a subspace of biorthogonal coupled-cluster solutions [38].Using these emulators we can reproduce full-space coupled-cluster computations of infinite nuclear matter with high precision at a tiny fraction of the computational cost.
After validating our approach, we use history matching to identify a region of the 17-dimensional parameter space of the low-energy constants (LECs) at NNLO that give acceptable results when confronted with few-body data.This iterative parameter search is enabled by employing emulators of few-body systems up to 4 He and allows us to collect 1.7 × 10 6 non-implausible interaction samples.Finally, we calibrate our ab initio model with two alternative LEC posterior probability density functions (PDFs) and use the principle of importance resampling [39,40] to quantify the sensitivity of nuclear matter predictions to calibration data.Here we focus on the results of this analysis, and the emergence of nuclear saturation, while details of the small-batch voting scheme and the history match are presented in a companion paper [37].
Method-The chiral Hamiltonian of ∆NNLO is parametrized with 17 LECs [5], and following Refs.[25][26][27] it can be written as Here, h 0 = t kin + v 0 where t kin is the kinetic energy and arXiv:2212.13203v2[nucl-th] 13 Jun 2024 v 0 represents the constant potential term without LECdependence.Note that ⃗ α is a vector that denotes all LECs.In this work we use non-local regulators with a cutoff Λ = 394 MeV/c.
Recently, model reduction methods [41] such as eigenvector continuation (EC) [42,43], has proven to be both efficient and accurate for emulating the output of ab initio computations of both scattering [44,45] and boundstate observables [25][26][27][28]46].These methods employ the fact that the eigenvector trajectory generated by smooth changes of the Hamiltonian matrix is well approximated by a very low-dimensional manifold [42].Following Refs.[25,26] we can therefore obtain good approximations of the ground-state of a target Hamiltonian H(⃗ α ⊚ ) by projecting it on a subspace of N sub different ground-state training eigenvectors and solving the corresponding N sub ×N sub -dimensional generalized eigenvalue problem.In this work we employ emulators for 2,3 H and 4 He for history matching, and the SPCC method with singles-, doubles-, and leading order triples excitations (CCSDT-3) for subsequent model calibration on 16 O (see Ref. [7] for details).We also construct a new emulator for 6 Li based on Ref. [28] for model checking.Here we used N sub = 32 training points reaching a relative accuracy of 10 −3 for the non-implausible parametrizations [47].
In this Letter we use the coupled-cluster method [38,[48][49][50][51][52][53][54][55][56][57][58][59] within the doubles approximation (CCD) and solve for pure neutron and symmetric nuclear matter on a cubic lattice in momentum space using periodic boundary conditions [30].The model-space has (2n max + 1) 3 momentum points with n max = 4.We use 132 nucleons for SNM and 66 neutrons for PNM which allows to minimize finite-size effects [30,60].We then use the EC-inspired SPCC method [25] to construct subspace emulators at five different nuclear densities.Since the subspace projected coupled-cluster Hamiltionian is non-Hermitian, the variational theorem does not hold [61].This might lead to the appearance of "spurious" states in the subspace spectrum that are lower in energy than the corresponding full-space coupled-cluster solution for certain combinations of the LECs.These states appear due to the non-symmetric treatment of the left and right bi-orthogonal coupled-cluster states [56].Therefore, we have developed a new algorithm-called small-batch voting-which allows to identify the physical ground state using a set of different subspace projections.This algorithm employs the strong sensitivity of spurious states on the specific choice of basis, and the contrasting stability of physical states.Details of this algorithm are given in the companion paper [37].
Results-We used history matching to collect a large number of 1.7 × 10 6 model parametrizations that exhibit acceptable (or at least non-implausible) fits to the history-matching observables.The latter include a total of 36 neutron-proton phase shifts in S and P waves up to T lab = 200 MeV scattering energy and six boundstate observables for A = 2, 3, 4 systems.This set is relevant for model calibration, allows fast model simulation or emulation, and permits the construction of simple implausibility measures.In the final wave we defined the non-implausible volume using a rotated hyperrectangle that captured parameter correlations and allowed to significantly increase the number of collected non-implausible samples.Further details of the history match are presented in the companion paper [37].We then strategically selected 64 of the most accurate nonimplausible samples for which we performed full CCD computations for SNM and PNM at five different densities ρ ∈ {0.12, 0.14, 0.16, 0.18, 0.20} fm −3 .Together with small-batch voting this allowed to create SPCC emulators for E/A and E/N at these densities.Finally, the nuclear matter EOS around saturation could be obtained by interpolation using Gaussian processes [62] as described in more detail in Ref. [37].
Cross validation of emulator performance is shown in Fig. 1.The 50 validation interactions were randomly selected from the non-implausible volume in a spacefilling manner.We conclude that our emulators predict the energy per particle for SNM (PNM) with < 10 −2 (< 10 −4 ) relative precision at a computational cost that is six (eight) orders of magnitude smaller than the full solution.Furthermore, Gaussian-process interpolation of the emulated EOS allows to extract empirical nuclearmatter properties with ≈ 1% precision (increasing to 3% and 10% for derivative quantities L and K, respectively).
Having access to fast and accurate emulators, and having identified a non-implausible region in the parameter space of our chiral interaction model, we can proceed to study the general behaviour of nuclear matter model predictions and possible correlations between different properties of these systems.Here we are not interested in the usual optimization approach that results in a single ("best fit") model parametrization.Instead, we consider all non-implausible samples from the history match.About 73% of these samples predict saturation within the studied density region and are kept for further consideration.The outputs from the nuclear matter emulators for these 1.2 × 10 6 interaction samples are shown in the upper triangle of Fig. 2.Here we have applied a density-dependent energy shift to approximately account for triples corrections [37].
We observe a strong anti-correlation between the saturation energy E 0 /A and the saturation density ρ 0 (Pearson correlation coefficient r = −0.92).This finding is in agreement with the Coester line [63] for nucleon-nucleon interactions.Similarly, the symmetry energy S and its slope L show a positive correlation (r = 0.60).These correlations have been seen in DFT calculations [64,65] and have been indicated with small families of EFTinspired [15] or (∆-less) chiral interactions [16].The com- parison with empirical nuclear matter properties reveals that the model predictions from the history match are clustered in a region with too small |E 0 /A|, ρ 0 , and S.This finding is consistent with previous results from various few-body optimized interactions [20,36,66].
The non-implausible predictions of nuclear matter properties that is shown in the upper triangle of Fig. 2 is a form of Bayes linear forecasting.By just considering expectation values such as means and variances-thereby avoiding the full probabilistic specification of uncertain quantities-we made the analysis simpler and more technically straightforward.While this allowed us to identify the parameter region of interest and to explore the model's forecasting capabilities, we now make an effort to extract a posterior predictive distribution (PPD) with which we can make probabilistic statements.
First, we revisit all remaining interaction parametrizations and keep only those that survive an extra implausibility analysis involving the selected S and P wave phase shifts and that give an unbound 1 S 0 ground-state in the neutron-proton system.This step results in a significant reduction to 8218 acceptable samples.Second, we introduce a set of calibration observables D cal and a normally-distributed likelihood, L (D cal | ⃗ α), assuming independent experimental, method, emulator, and model errors [37].At this stage we employ the established method of sampling/importance resampling [39,40] to approximately extract samples from the parameter posterior pr (⃗ α | D cal , I).In this step we assume a uniform prior probability for all non-implausible samples [71].
In this work we especially considered two different sets of model calibration data: (i) D cal = D A=2,3,4 encompassing binding energies and point-proton radii of 2,3 H and 4 He plus the quadrupole moment of the deuteron, and (ii) D cal = D A=2,3,4,16 where we complemented the above set of observables with the energy and radius of 16 O.Note that our choice of history-matching observables allows for both PPD resampling analyses to be performed from the same set of non-implausible prior samples.
Having access to approximate parameter posteriors we first extract samples from the model PPD defined as calibration data, respectively.Both predictions are consistent with the experimental value S d = 1.474MeV [72], but the latter one is significantly more precise.
We then predict nuclear matter properties.For these predictions we collect samples from the full PPD, which includes the EFT truncation error, CC method errors, and the SPCC emulator error (estimated from cross validation shown in Fig. 1).The density dependence and relevant cross correlation of these errors are described by a Bayesian machine learning error model [19,20,37].Our strategic emulator construction-using training samples with high importance weight-implies that up to 60% of resampled interactions have zero emulator error.Using the two different sets of calibration data, D A=2,3,4 and D A=2,3,4,16 , we compare the resulting PPDs that we label PPD A=2,3,4 and PPD A=2,3,4,16 , respectively.These are shown in the lower triangle of Fig. 2.
The marginal distributions on the diagonal reveal that ρ 0 , E 0 /A and S for PPD A=2,3,4 are characterized by low precision and mild tension with the empirical region-as previously observed with the non-importance-weighted samples.In other words, even by enforcing a higher accuracy for two-and few-body systems via the data likelihood, the general description of nuclear matter properties is not improved.On the other hand, the predictions become more precise and accurate when we include the 16 O observables.In particular, the saturation point is more precisely predicted while its mode shifts to larger saturation density and binding energy.We actually find that PPD A=2,3,4 displays a significant asymmetry in some dimensions, e.g., ρ 0 , which hints to a possible bimodality.Samples from the tail region of PPD A=2,3,4 correspond largely to the mode of PPD A=2,3,4,16 .This finding shows   that the emergence of saturation, represented by the position and shape of this mode, depends on the choice of calibration observables.
Our predictions of the nuclear EOS around saturation, with quantification of relevant sources of uncertainty, can serve as an important anchor for extrapolations to higher densities and studies of neutron star physics.To illustrate, we show a simple extrapolation based on PPD A=2,3,4,16 -which gives the most precise predictions-for empirical parameters in Fig. 3.It higlights the fact that an uncertainty band is obtained, which is key for rigorous error estimates at higher densities.A multivariate Gaussian approximation of our PPD for nuclear-matter parameters is provided as Supplemen-tal Material [67].The interaction samples with importance weights are provided in the companion paper [37].Summary and outlook-In this work we constructed emulators that accurately reproduce full-space coupledcluster computations of nucleonic matter starting from ∆-full χEFT at NNLO while reducing the computational cost by several orders of magnitude.The small-batch voting algorithm was developed to remove spurious states that occured in the SPCC method.Using these tools, together with emulators of light nuclei, we employed history matching to identify more than one million acceptable interaction samples, and could reveal correlations between different properties of infinite nuclear matter.We then employed importance resampling to obtain PPDs, and studied the sensitivity to the choice of calibration data.Broad, asymmetric marginal distributions for the saturation energy and density were observed when making predictions conditional on few-nucleon data only, while the predictions become more precise and accurate when adding the 16 O energy and radius to the calibration data set.Binding energies and radii of finite nuclei are obviously useful for model calibration, but are not the only choice.Other observables such as the 3 H beta-decay rate and three-nucleon scattering cross sections should also be considered within a statistical framework involving many-body predictions.Ongoing work shows that also these observables can be efficiently computed [74,75] or emulated [27,45].Our rigorous error bands for the nuclear matter EOS around saturation will be important for advances in studies of dense, neutron-rich matter and for the interpretation of nascent observations from multimessenger astonomy [10][11][12].
The introduction of small-batch voting demonstrates how emulators can be successfully constructed for challenging many-body observables.Further developments within various many-body computational frameworks are expected to follow [76] in both nuclear physics and beyond.Future work using history matching, Bayesian PPD sampling, and many-body emulators, as exemplified by this study, will help to elucidate the information content of various low-energy observables, the order-byorder convergence of χEFT, and the predictive power of ab initio modeling across the nuclear landscape.

FIG. 1 .
FIG. 1. (Color online)Cross-validation of the SPCC emulator with exact CCD calculations using 50 interaction samples.Shown here from top to bottom is the nuclear matter saturation density (ρ0), saturation energy (E0/A), and symmetry energy (S).Histograms of relative errors are shown in the right column.

FIG. 2 .
FIG. 2. (Color online)Upper triangle: nuclear matter emulator output (saturation density ρ0, saturation energy E0/A, symmetry energy S, symmetry energy slope L and incompressibility K) for the non-implausible interactions from the fifth wave of history matching.The axes limits are the same as in the corresponding panels in the lower triangle.Lower triangle and diagonal: PPDs for nuclear matter properties using two different PDFs for the LECs plus error sampling.These predictions are based either on few-body (A = 2, 3, 4) calibration data (orange PPD) or the addition of16 O to the calibration data set (blue PPD).See the Supplemental Material [67] for a Gaussian approximation of PPDA=2,3,4,16.The contour lines on the bivariate distributions denote 68% and 90% credible regions.The red bands indicate empirical ranges for E0/A = −16.0± 0.5 MeV, ρ0 = 0.16 ± 0.01 fm −3 , S = 31 ± 1, L = 50 ± 10 and K = 240 ± 20 MeV from Refs.[68][69][70].
FIG. 3. (Color online) Dark (light) blue bands show the PPDA=2,3,4,16 for the energy per particle of PNM (upper panel) and SNM (lower panel) computed in the ρ ∈ [0.12, 0.20] fm −3 range.The corresponding purple bands represent simple extrapolations based on the empirical nuclear-matter parameters.For each sample of the PPD the EOS is generated by using the expansion of Ref.[73] and retaining up to quadratic terms in x = (ρ − ρ0)/3ρ0.Ten such samples are shown with thin, gray lines.The extrapolation uncertainty is not incorporated in the error bands.