Quantifying uncertainties and correlations in the nuclear-matter equation of state

We perform statistically rigorous uncertainty quantification (UQ) for chiral effective field theory ($\chi \mathrm{EFT}$) applied to infinite nuclear matter up to twice nuclear saturation density. The equation of state (EOS) is based on high-order many-body perturbation theory calculations with nucleon-nucleon and three-nucleon interactions up to fourth order in the $\chi \mathrm{EFT}$ expansion. From these calculations our newly developed Bayesian machine-learning approach extracts the size and smoothness properties of the correlated EFT truncation error. We then propose a novel extension that uses multitask machine learning to reveal correlations between the EOS at different proton fractions. The inferred in-medium $\chi \mathrm{EFT}$ breakdown scale in pure neutron matter and symmetric nuclear matter is consistent with that from free-space nucleon-nucleon scattering. These significant advances allow us to provide posterior distributions for the nuclear saturation point and propagate theoretical uncertainties to derived quantities: the pressure and incompressibility of symmetric nuclear matter, the nuclear symmetry energy, and its derivative. Our results, which are validated by statistical diagnostics, demonstrate that an understanding of truncation-error correlations between different densities and different observables is crucial for reliable UQ. The methods developed here are publicly available as annotated Jupyter notebooks.

Figure 1

The three columns show three sets of random functions. For each column the three rows show the output of GPs that have the same (arbitrary) mean and marginal variance, but differ in their radial basis function (RBF) kernel: a different length scale is used in each of the three rows. The length scales are indicated by the double-headed arrow at the bottom of each panel Each panel then contains four functions, i.e., four draws from the GP. The $x$-axis could represent an independent variable such as energy, angle, or density and the $y$-axis can be thought of as the order-by-order EFT coefficients for an observable of interest. The dark (light) shaded bands represent one (two) standard deviations from the mean.

Figure 2

Observable coefficients, $c_n$, for $E/N\left(n\right)$ up to ${\mathrm{N}}^3\mathrm{LO}$, as a function of density $n$, obtained using the $\mathrm{\Lambda }=500\text{MeV}$ interactions in Table 1. Markers indicate training points, gray lines indicate $2\overline{c}$, and colored bands are 68% credible intervals of the interpolating GPs. The estimated hyperparameters are given by $\overline{c}=1.0$ and $\ell =0.97{\text{fm}}^{-1}$. Note that the secondary $x$ axis (at the top of the figure) is not linear in $k_{\text{F}}^{\text{PNM}}$.

Figure 3

Observable coefficients $c_n$ for $E/A\left(n\right)$ up to ${\mathrm{N}}^3\mathrm{LO}$, as a function of density $n$, using the $\mathrm{\Lambda }=500\text{MeV}$ interactions in Table 1. See Fig. 2 for the notation. The estimated hyperparameters are given by $\overline{c}=2.9$ and $\ell =0.48{\text{fm}}^{-1}$. Note that the secondary $x$ axis (at the top of the figure) is not linear in $k_{\text{F}}^{\text{SNM}}$.

Figure 4

The posteriors for the EFT breakdown scale ${\mathrm{\Lambda }}_b$ using orders through ${\mathrm{N}}^2\mathrm{LO}$ (blue bands) and ${\mathrm{N}}^3\mathrm{LO}$ (red bands) corresponding to Figs. 2 and 3. The upper pair of posteriors comes from analyzing $E/N$, the middle pair from $E/A$, and the bottom from a combined analysis. In all cases a Gaussian prior centered at ${\mathrm{\Lambda }}_b=600\pm 150\text{MeV}$ is used. The combined ${\mathrm{N}}^3\mathrm{LO}$ posterior is consistent with the ${\mathrm{\Lambda }}_b\approx 600\text{MeV}$ found when considering free-space NN scattering observables [35].

Figure 5

Length-scale posteriors organized similarly to Fig. 4. A scale-invariant prior proportional to $1/\ell$ is used, and each length scale is relative to the $k_{\text{F}}$ of each system. If we were to use a single $k_{\text{F}}$ prescription, the length scales in PNM and SNM, ${\ell }_{\text{PNM}}$ and ${\ell }_{\text{SNM}}$, respectively, would transform just as $k_{\text{F}}$, making the posteriors shift towards agreement. See the discussion of the input space in the main text.

Figure 6

Energy per particle in PNM with truncation errors using the $\mathrm{\Lambda }=500\text{MeV}$ interactions in Table 1. From left to right, top to bottom, the panels show the order-by-order progression of EFT uncertainties as the $\chi \mathrm{EFT}$ order increases. The bands indicate 68% credible intervals.

Figure 7

Similar to Fig. 6 but for SNM. The gray box depicts the empirical saturation point, $n_0=0.164\pm 0.007{\text{fm}}^{-3}$ with $E/A\left(n_0\right)=-15.86\pm 0.57\text{MeV}$, obtained from a set of energy-density functionals [18, 51] (see the main text for details).

Figure 8

Credible-interval diagnostics for the $E/N\left(n\right)$ (left-hand side) and $E/A\left(n\right)$ uncertainty bands (right-hand side) for the $\mathrm{\Lambda }=500\text{MeV}$ interactions in Table 1; for details see Ref. [25]. At each order we construct an uncertainty band for the upcoming correction (not the full truncation error) and test whether the next order is contained within it at a specific credible interval. The expected size of fluctuations due to the finite effective sample size of the curves is depicted using dark (light) gray bands for the 68% (95%) interval. Both bands are quite large, which shows that correlations are crucial to assess whether truncation errors have been properly assigned.

Figure 9

Predicted nuclear saturation point of SNM at ${\mathrm{N}}^2\mathrm{LO}$ (blue) and ${\mathrm{N}}^3\mathrm{LO}$ (red) with 95% ($2\sigma$) credible interval ellipses including correlated truncation errors. Panels (a) and (b) show the $\mathrm{\Lambda }=500\text{MeV}$ and $\mathrm{\Lambda }=450\text{MeV}$ interactions described in Table 1, respectively. The colored crosses depict the minimum of each EOS obtained at that order without accounting for truncation errors. In contrast, each dot is the minimum of a curve sampled from the corresponding $E/A\left(n\right)$ GP in Fig. 7, which are used to fit the ellipses. These should be compared to the $2\sigma$ uncertainty bands in the same color, which are estimated using ${\mathrm{\Lambda }}_b=600\text{MeV}$. The gray box again depicts the empirical saturation point.

Figure 10

Correlations between the observable coefficients $c_n\left(k_{\text{F}}\right)$ from $E/N\left(n\right)$ and $E/A\left(n\right)$. The points on each coefficient curve are at the same density. (a) The coefficients from the $\mathrm{\Lambda }=500\text{MeV}$ interactions. (b) Toy coefficients created using a correlated multioutput GP trained on the coefficients shown in (a). Each dimension of the multi-output GP ellipse takes into account the associated $\overline{c}$ from Figs. 2 and 3, which need not be the same as the empirical variances of the ellipse in (a). For the $\mathrm{\Lambda }=450\text{MeV}$ interactions, the corresponding plots are shown in (c) and (d).

Figure 11

Total correlation matrix of $E/N\left(n\right)$ and $E/A\left(n\right)$ assuming a multitask GP model that was trained to order-by-order results of the $\mathrm{\Lambda }=500\text{MeV}$ interactions. Each submatrix uses the same grid spaced linearly in $k_{\text{F}}$. The diagonal blocks show the autocorrelation, and the off-diagonal block is known as the cross correlation. All use RBF kernels [see Eqs. (6), (C12), and (C16)]. The length scale of the $E/A\left(n\right)$ blocks has been transformed to $k_{\text{F}}^{\text{PNM}}$ as discussed in the text. Hence, points of equal density between $E/N\left(n\right)$ and $E/A\left(n\right)$ lie on the diagonal band of the off-diagonal blocks, making them the most highly correlated, but the $E/A\left(n\right)$ autocorrelation is unchanged. The length scales were determined by fitting to the $E/N\left(n\right)$ and $E/A\left(n\right)$ coefficients independently as shown in Sec. 3b. The correlation $\rho =0.95$ is a prediction (which agrees with the empirical correlation) given these length scales.

Figure 12

Similar to Fig. 6, but these are the order-by-order predictions of the symmetry energy $S_2\left(n\right)$ for the $\mathrm{\Lambda }=500\text{MeV}$ interactions.

Figure 13

Order-by-order predictions of the pressure $P\left(n\right)$ of SNM, including differentiation and truncation uncertainty, for the $\mathrm{\Lambda }=500\text{MeV}$ interactions. See the main text for details.

Figure 14

Violin plots of the incompressibility $K$ of SNM, shown order by order for the $\mathrm{\Lambda }=500\text{MeV}$ interactions. Curves show the entire smoothed posterior (and its reflection). Each posterior includes differentiation and truncation uncertainty, and are marginalized over all plausible saturation densities $n_0=0.17\pm 0.01{\text{fm}}^{-3}$; see the main text. Dots and bars indicate the mean value, along with the $1\sigma$ and $2\sigma$ uncertainties. The black line and gray band extends the ${\mathrm{N}}^3\mathrm{LO}$ mean and $1\sigma$ uncertainty to more easily compare $\chi \mathrm{EFT}$ orders.

Figure 15

Observable coefficients for (a) $E/A\left(n\right)$ (SNM) and (b) $E/N\left(n\right)$ (PNM) up to ${\mathrm{N}}^3\mathrm{LO}$ using the $\mathrm{\Lambda }=450\text{MeV}$ interactions in Table 1. Markers indicate training points, gray bands indicate $2\overline{c}$ and colored bands are 68% credible intervals of the interpolating GPs. The estimated hyperparameters are given by $\overline{c}=3.0$ and $\ell =0.50{\text{fm}}^{-1}$ for SNM and $\overline{c}=0.96$ and $\ell =0.81{\text{fm}}^{-1}$ for PNM.

Figure 16

Model checking diagnostics applied to (a) SNM and (b) PNM coefficients from the $\mathrm{\Lambda }=450\text{MeV}$ interactions in Fig. 15. The MD computed against the underlying process is shown in the left panel of each subplot. The interior line, box end caps, and whiskers on the box plot show the median, 50% credible intervals, and 95% credible intervals, respectively. The right panel shows the PC diagnostic ${\mathbf{D}}_{\text{PC}}$ vs index, with gray lines that represent its $2\sigma$ error bands. Both diagnostics point to the $c_3\left(k_{\text{F}}\right)$ coefficient as a possible outlier. See Ref. [25] for more details about analyzing these plots.

Figure 17

Similar to Fig. 16 but for the $\mathrm{\Lambda }=500\text{MeV}$ interactions in Figs. 2 and 3.

Figure 18

Observables coefficients and model diagnostics for $E/N\left(n\right)$ (PNM) up to ${\mathrm{N}}^3\mathrm{LO}$ for the $\mathrm{\Lambda }=500\text{MeV}$ interactions in Table 1 using separate reference scales for NN-only (denoted with a superscript “(2)”) and 3N contributions (denoted with a superscript “(3)”). The model for $y_{\mathrm{ref}}\left(k_{\text{F}}\right)$ is explained in the text [see Eq. (A1)].

Figure 19

The posteriors for the EFT breakdown scale ${\mathrm{\Lambda }}_b$ as in Fig. 4 but using the alternative model for $y_{\mathrm{ref}}\left(k_{\text{F}}\right)$ as explained in the text [see Eq. (A1)].

Figure 20

Length-scale posteriors as in Fig. 5 but using the alternative model for $y_{\mathrm{ref}}\left(k_{\text{F}}\right)$ as explained in the text [see Eq. (A1)].

Figure 21

Same as Fig. 18 but for $E/A\left(n\right)$ (SNM).