Propagation of Statistical Uncertainties of Skyrme Mass Models to Simulations of $r$-Process Nucleosynthesis

TL;DR

This paper investigates how uncertainties in nuclear mass models, specifically Skyrme energy density functionals, affect r-process nucleosynthesis simulations, and identifies key parameters influencing abundance predictions, highlighting the importance of improved mass measurements.

Contribution

It introduces a Bayesian framework to propagate uncertainties from Skyrme model parameters to r-process abundance predictions and identifies parameters correlated with abundance features.

Findings

Uncertainty reduction in nuclear masses significantly decreases r-process abundance uncertainties.

Certain Skyrme parameters are strongly correlated with specific abundance pattern features.

Using new mass measurements can substantially improve the accuracy of nucleosynthesis simulations.

Abstract

Uncertainties in nuclear models have a major impact on simulations that aim at understanding the origin of heavy elements in the universe through the rapid neutron capture process ($r$ process) of nucleosynthesis. Within the framework of the nuclear density functional theory, we use results of Bayesian statistical analysis to propagate uncertainties in the parameters of energy density functionals to the predicted $r$-process abundance pattern, by way not only of the nuclear masses but also through the influence of the masses on $\beta$-decay and neutron capture rates. We additionally make the first identifications of specific parameters of Skyrme-like energy density functionals which are correlated with particular aspects of the $r$-process abundance pattern. While previous studies have explored the reduction in the abundance pattern uncertainties due to anticipated new measurements of neutron-rich nuclei, here we point out that an even larger reduction will occur when these new measurements are used to reduce the uncertainty of model predictions of masses, which are then propagated through to the abundance pattern. We make a quantitative prediction for how large this reduction will be.

Figures

Full text

The heaviest elements owe their origins to rapid neutron capture, or $r$-process, nucleosynthesis. In the $r$-process, heavy elements are built up via a sequence of rapid neutron captures and $\beta$-decays that populate nuclei far to the neutron-rich side of stability Burbidge+57; Cameron57. The astrophysical source of the intense neutron flux was initially suspected to be within core-collapse supernovae Meyer+92; Woosley+94, though decades of careful study have shown the required conditions are unlikely to be obtained in this environment Arcones+07; Fischer+10; Hudepohl+10; Roberts+12. Recent evidence Hirai+2015; Ji+2016, including the discovery of GW170817/GRB170817a/SSS17a AbbottGW170817; Cowperthwaite2017, increasingly points to neutron star mergers as the likely $r$-process site. However, many open questions remain. For example, what specific environments within neutron star merger events are responsible for $r$-process production, and what are their properties? Can neutron star mergers account for all galactic $r$-process production, or are there additional astrophysical sites?

The $r$-process astrophysical conditions could in principle be identified by comparing simulations of abundance patterns of elements and observations in the solar system and in old stars. However, analysis of individual environments is complicated by large uncertainties in the astrophysics and nuclear physics Arnould+07. Here we consider the latter. Simulations of the $r$-process are dependent upon nuclear data, including masses, neutron capture rates, and $\beta$-decay and fission properties, for thousands of neutron-rich nuclei Mumpower+16. In spite of a concerted effort at radioactive beam facilities worldwide to measure these properties directly or indirectly, the vast majority of them are as of yet inaccessible and we must rely on theoretical estimates.

Nuclear density functional theory (dft) is currently the only approach that can provide all of these properties in a consistent yet microscopic framework bender2003. Most energy density functionals (edf) are typically characterized by approximately a dozen parameters that are fitted on a small set of nuclear properties. The choices made in selecting the form of the edf and the set of experimental data to fit its parameters lead to both systematic and statistical uncertainties that have an impact on all applications schunck2015-c.

Ideally, one would like to consider simultaneously all sources of uncertainties (systematic, statistical and numerical) and propagate them to all observables (separation energies, $\alpha$-, $\beta$- and $\gamma$-decay rates, fission rates, neutron capture rates) relevant to astrophysical simulations. Such an approach is currently not feasible, partly because of its formidable computational cost, partly because there are still gaps in our understanding of, e.g., $\alpha$-decay, neutron capture or fission. However, we can exploit recent work in determining estimates of theoretical uncertainties to quantify the variations in simulated $r$-process abundances that result from nuclear mass uncertainties alone. Past work in this area has either considered abundance pattern comparisons between distinct mass models, e.g. Martin+2016, or ranges of patterns that result from random, uncorrelated mass variations Mumpower+16; Surman+16.

In this work, we perform the first rigorous propagation of statistical uncertainties of nuclear mass models based on dft. We generate fifty different edfs by sampling the Bayesian posterior distribution of the unedf1 edf. For each sample, we compute a full nuclear chart and update neutron capture rates and $\beta$-decay properties to be consistent with each table. We implement these sets of nuclear data in $r$-process simulations to place “error bars” due to nuclear masses on $r$-process abundances and to identify correlations between theoretical model parameters and abundance pattern features. Such correlations could possibly lead to additional constraints on $r$-process conditions or, e.g., the unedf1 parameters themselves. Finally, we provide a quantitative estimate of the improvements to $r$-process pattern uncertainties expected from anticipated mass measurements at current and upcoming facilities and concurrent advancements in theoretical models.

We begin by computing atomic mass tables within the nuclear dft approach to nuclear structure with Skyrme edfs. Our starting point is the unedf1 parametrization, in which the coupling constants were optimized globally on select experimental nuclear masses, radii, deformations and excitation energies of fission isomers in the actinides Kortelainen:2011ft. While the r.m.s. deviation on nuclear binding energies of unedf1 is only 1.8 MeV, it goes down to 0.45 MeV for 2-neutron separation energies. Bayesian inference methods were later used to compute the posterior distribution of the unedf1 parameters Higdon:2014tva and propagate theoretical statistical uncertainties in predictions of nuclear masses, two-neutron drip line, and fission barriers McDonnell:2015sja. Here, we sample the same posterior distribution within the $90\%$ confidence region to generate fifty different parameter sets for the Skyrme edf.

For each sample, we compute the nuclear ground-state binding energy of all even-even nuclei from Hydrogen to $Z=120$ by solving the Hartree-Fock-Bogoliubov (HFB) equation. The limits of nuclear stability (proton and neutron drip lines) are reached when the value of the two neutron (proton) separation energy becomes negative. Compared to alternative options based, e.g., on the value of the Fermi energy, this criterion offers the advantage of being model-independent since binding energies are true observables. With this criterion, each mass table contains approximately 2,000 even-even nuclei. For each even-even nucleus, the ground-state is determined by exploring locally the potential energy surface of the nucleus for a range of eleven axial quadrupole deformations $\beta_{2}$ between -0.5 and +0.5. The configuration with the lowest energy defines the ground-state. Details of the exploration of the even-even nuclear landscape with the numerical solver hfbtho can be found in perez2017. With this procedure, computing 50 mass tables requires of the order of 1 million HFB calculations.

Although odd-even and odd-odd binding energies could be computed with the blocking procedure, see, e.g. schunck2010, this would require about an order of magnitude more HFB calculations. Instead we adopt a standard approximation for the binding energy of odd nuclei that combines information about binding energies and HFB pairing gaps in neighboring isotopes/isotones; see Supplemental Material of erler2012. This procedure yields an excellent approximation of, in particular, one-particle separation energies.

For each of the fifty mass tables thus described, we calculate a self-consistent set of all nuclear data inputs required for *r-*process calculations. We calculate neutron capture and neutron-induced fission rates using the Los Alamos Hauser-Feshbach code CoH Kawano+16 and $\beta$-decay half-lives with probabilities for delayed emission of one or more neutrons using the QRPA+HF framework of MumpQRPA+HF and unmodified strength data from Moller+97. We repeat these calculations using the masses given in the 2016 Atomic Mass Evaluation (AME2016) wang2017; where possible, these results are taken to replace those based on the unedf1 mass tables. The decay properties of the Nubase 2016 compilation NUBASE2016 are further taken to replace any calculated values based on either AME2016 or unedf1 nuclear masses. For all fissioning nuclei, we use a symmetric, two-particle product distribution.

We implement each of these datasets into the nuclear reaction network code PRISM Mumpower+17b; Mumpower+18; Holmbeck+18 to simulate nucleosynthesis for three distinct types of astrophysical conditions where *r-*process nucleosynthesis may occur: (1) a supernova-type high-entropy wind, with entropy $s/k=300$, dynamical timescale $\tau=80$ ms, and electron fraction $Y_{e}=0.30$, (2) a parameterized merger accretion disk wind with $s/k=30$, $\tau=80$ ms, and $Y_{e}=0.21$, and (3) fission-recycling outflow from a neutron star merger Mendoza+15. For each simulation, we dynamically update the evolution of temperature with respect to the release of energy from nuclear reactions, decays, and fission, with an assumed thermalization efficiency of $10\%$ for all energy released.

The range in final abundance patterns across these fifty calculations is shown in Fig. 1 for each set of astrophysical conditions we consider. In each case, the shaded band represents the propagation of statistical uncertainties from the unedf1 nuclear mass model to the corresponding $r$-process simulation. The influence of masses on reaction and decay rates contribute to the width of the band, via mechanisms described in, e.g., Mumpower+16 and references therein. In addition, the location of the neutron drip line is of key importance for some types of astrophysical conditions. In the fifty mass tables considered here, the location of the one-neutron drip line varies by more than ten neutron numbers. Notably, the band is the widest for the astrophysical conditions in which the $r$-process path is the closest to the drip line, the fission recycling example of Fig. 1(c).

With the wealth of data available from these $r$-process simulations, we can search for correlations between unedf1 functional parameters and the formation of abundance pattern features. Here we demonstrate how such analyses might proceed.

Several of the unedf1 functional parameters are poorly constrained by data near stability. One such parameter is the isovector surface coupling constant $C_{1}^{\rho\Delta\rho}$, with unedf1 value $-145.382\pm 52.169$; see Table II in Kortelainen:2011ft. For a low-entropy hot wind $r$-process environment, this parameter is correlated with the formation of the rare earth peak, the small feature around $A\sim 160$ in the solar $r$-process isotopic pattern. The top panel of Fig. 2 shows $C_{1}^{\rho\Delta\rho}$ versus the abundance-weighted average $A$ of the rare earth peak for the $r$-process simulations from the middle panel of Fig. 1. The placement of the rare earth peak is calculated from the solar $r$-process abundances of Arnould+07 and Sneden+08 and is given by the shaded vertical band. Correlations between $C_{1}^{\rho\Delta\rho}$ and rare earth peak placement are weaker for other types of $r$-process environments. The $r$-process path in the high-entropy wind case is not so neutron-rich and thus not as sensitive to $C_{1}^{\rho\Delta\rho}$. The fission recycling example has a distinct rare earth peak formation mechanism Mumpower+17 that is not particularly active with the unedf1 masses, resulting in a comparatively weaker correlation. However, recent studies mclaughlin+2005; surman+2006; Just+2015; Martin+2015; Wanajo+2014; Siegel+2017 favor $r$-process conditions that are most similar to those of our low-entropy wind where this correlation is strongest, suggesting that the $r$-process abundance pattern may provide an important additional constraint on the value of $C_{1}^{\rho\Delta\rho}.$

In all of the astrophysical environments considered, we found the proton pairing strength $V_{0}^{p}$ to be correlated with the ratio between the summed abundances of the rare earth and $A\sim 195$ peak regions, as illustrated in the bottom panel of Fig. 2, where the solar values are given by the shaded band. The correlations in each case are distinct, with different astrophysical conditions picking out different preferred values of $V_{0}^{p}$. Only the least negative values of $V_{0}^{p}$ considered reproduce solar values for the high-entropy conditions, while values of $V_{0}^{p}$ that tend towards the center of the distribution reproduce solar values for the low-entropy conditions. Within the range of values we consider, the fission recycling conditions fail to reproduce solar values, with the correlation suggesting an even more negative value of $V_{0}^{p}.$ Thus, if $V_{0}^{p}$ could be more tightly constrained, the simulated ratio of the rare earth and $A\sim 195$ peak regions could be used as a diagnostic of $r$-process conditions.

Measurements of the masses of increasingly neutron-rich nuclei are the focus of a number of experimental efforts worldwide, for example at the Canadian Penning Trap at CARIBU Hirsh+2016; Orford+2018, JYFLTRAP at Jyväskylä Kankainen+2013; Vilen+2018, ISOLTRAP at CERN Lunney+2017, TITAN at TRIUMF Lascar+2017, and storage rings at GSI in Germany, IMP in China, and RIKEN in Japan Zhang+2016. Next-generation radioactive ion facilities, such as the Facility for Rare Isotope Beams (FRIB), will have unprecedented access to isotopes far from stability FRIB. New mass measurements improve the reliability of $r$-process simulations in two ways: directly, by dramatically reducing the uncertainty in the masses of newly measured nuclei, and indirectly, by enabling improvements to mass modeling. Theoretical mass models are all calibrated to known data, so known masses tend to be well reproduced by theory. Outside the known region, theoretical predictions tend to diverge. The variations among our fifty unedf1 mass tables, shown for the tin isotopes in Fig. 3, clearly demonstrate this behavior. For this element, unedf1 fits known masses to about $\sigma_{rms}\sim 1$ MeV, and variations increase sharply past the $N=82$ closed shell. Additional measurements increase the available data with which to constrain theory and thus hold the potential to reduce uncertainties outside the measured region. We simulate this effect by generating an adjusted set of fifty mass tables, in which the variations from the mean are reduced to match the experimentally-known region for nuclei within the FRIB range and increase with the same slope outside this range. The variations of the simulated set of tables are shown in the dark shaded region of Fig. 3.

Our two sets of unedf1 mass tables can be used to quantify the reductions in $r$-process abundance pattern uncertainties that have already been achieved by measurements to date and that are anticipated from future mass measurements. We rerun the example $r$-process simulations from Fig. 1 using three different sets of nuclear data. The first set is a theory-only set, with all quantities derived exclusively from our fifty unedf1 tables. The second set is that used in Fig. 1, where experimental nuclear data is additionally incorporated. The third set is constructed to mimic the influence of anticipated mass measurements. Experimental values or values derived from the nominal unedf1 mass tables are held fixed for all nuclei within the FRIB reach; elsewhere we use theory values derived from our set of fifty simulated unedf1 mass tables.

Fig. 4 shows the abundance pattern variations normalized by the mean for the three astrophysical trajectories in Fig. 1, each calculated with the three data sets described above. For the high entropy wind, Fig. 4(a), many of the relevant nuclear properties have already been measured, so there is significant improvement realized between the theory-only (lightest shaded band) calculations and those that include current experimental values (medium shaded band). Looking forward to FRIB, the majority of nuclei along the equilibrium $r$-process path in the $N=82$ and rare earth regions will be within reach. Thus systematic measurement campaigns at FRIB have the potential to essentially remove mass as a source of uncertainty in simulated $r$-process abundances below $A\sim 170$ for high entropy winds.

However, in the currently-favored potential $r-$process astrophysical site of neutron star/neutron star-black hole mergers, the environments are likely lower entropy, $s/k\sim 5-50$, and more neutron-rich, similar to the conditions used for the middle and bottom panels of Figs. 1 and 4. The $r$-process equilibrium paths are farther from stability in these cases, thus the current reach of experimental data results in more modest improvements, as indicated when comparing the light- and medium-shaded bands in Fig. 4(b) and (c). Prospects for the future, however, are encouraging. For the low-entropy wind example of Fig. 4(b), FRIB can reach the majority of the key nuclei and the remaining uncertainty band should be similar to the high-entropy wind case. In particular the excellent precision anticipated for abundances $140<A<170$ can facilitate the use of the rare earth peak as a key $r$-process diagnostic Mumpower+12; Orford+2018.

For the fission recycling example, uncertainties in the location of the drip line and in the fission properties of heavy nuclei near the drip line dominate the uncertainty bands. Even with FRIB at full power these uncertainties are unlikely to be resolved with direct measurements. Here nuclear theory will play a critical role. The simulated improvements to theory anticipated in our approach do result in a narrowing of the uncertainty band, as seen in a comparison between the medium- and dark-shaded bands of Fig. 4(c). Further potential improvements to nuclear EDF theory, e.g. Neufcourt+2018, and its full application to the problem of fission, e.g. Giuliani+2018; Bulgac+18, are not captured in our approach. Therefore, there remains the possibility for more significant improvements to the uncertainty band associated with fission-recycling conditions with concurrent advances in experiment and theory.

The origins of the heaviest elements have remained mysterious for decades. Thanks to concerted efforts in astrophysical modeling, spectroscopic observations, neutrino and nuclear experiment and theory, and, now, gravitational wave astronomy, a detailed understanding of $r$-process nucleosynthesis finally seems within reach. Still, further advances are needed in each of these areas. Here we have highlighted how careful quantification of nuclear physics uncertainties has the potential to provide crucial insight into $r$-process astrophysical conditions and the nuclear models themselves.