Neutron drip line in the Ca region from Bayesian model averaging

TL;DR

This paper uses Bayesian model averaging and machine learning to predict the neutron drip line in calcium isotopes, providing quantified certainty levels and consistent predictions across models.

Contribution

It introduces a Bayesian approach with Gaussian processes to assess the bound nature of neutron-rich nuclei, improving predictions of the calcium drip line with uncertainty quantification.

Findings

$^{68}$Ca has a 76% probability to be bound to two-neutron emission.

$^{61}$Ca likely decays by neutron emission with 46% probability.

Predictions are consistent across different global mass models despite raw prediction variations.

Abstract

The region of heavy calcium isotopes forms the frontier of experimental and theoretical nuclear structure research where the basic concepts of nuclear physics are put to stringent test. The recent discovery of the extremely neutron-rich nuclei around $^{60}$Ca [Tarasov, 2018] and the experimental determination of masses for $^{55-57}$Ca (Michimasa, 2018] provide unique information about the binding energy surface in this region. To assess the impact of these experimental discoveries on the nuclear landscape's extent, we use global mass models and statistical machine learning to make predictions, with quantified levels of certainty, for bound nuclides between Si and Ti. Using a Bayesian model averaging analysis based on Gaussian-process-based extrapolations we introduce the posterior probability $p_{ex}$ for each nucleus to be bound to neutron emission. We find that extrapolations for drip-line locations, at which the nuclear binding ends, are consistent across the global mass models used, in spite of significant variations between their raw predictions. In particular, considering the current experimental information and current global mass models, we predict that $^{68}$Ca has an average posterior probability ${p_{ex}\approx76}$% to be bound to two-neutron emission while the nucleus $^{61}$Ca is likely to decay by emitting a neutron (${p_{ex}\approx 46}$ %).

Equations

$$w_k:=p\left({\mathcal{M}}_k|{}^{52}\text{Cl},{}^{53}\text{Ar},{}^{49}\text{S}\text{exist}\right)$$

Full text

Introduction – How many protons and neutrons can form a bound atomic nucleus? Out of about 3,200 isotopes known Thoennessen (2016) only 286 primordial nuclides have existed in their current form since before Earth was formed. They form the valley of stability on the nuclear landscape. Moving away from the region of stable isotopes by adding neutrons or protons, one enters the regime of short-lived radioactive nuclei, which are beta unstable. Nuclear existence ends at the “drip lines”, where the last nucleons are no longer attached to the nucleus by the strong interaction and drip off. According to current theoretical estimates Erler et al. (2012); Agbemava et al. (2014) the number of bound nuclides with atomic number $Z$ between 2 and 120 is around 7,000.

The particle stability of a nuclide is determined by its separation energy, i.e., the energy required to remove from it a single nucleon or a pair of like nucleons. If the separation energy is positive, the nucleus is bound to nucleon decay; if the separation energy is negative, the nucleus is particle-unstable. In this Letter, we study the one-neutron ($S_{1n}$) and two-neutron ($S_{2n}$) separation energies of neutron-rich nuclei. The drip line is reached when the separation energy reaches zero; hence, one can talk about the one-neutron drip line when $S_{1n}=0$ and the two-neutron drip line when $S_{2n}=0$. Very weakly bound, or unbound, nuclei that lie in the immediate vicinity of drip lines are referred to as threshold systems. The separation energies and drip-line positions are strongly affected by nucleonic pairing, or nuclear superfluidity Brink and Broglia (2005). Since it costs energy to break a nucleonic pair, nuclei with even numbers of nucleons are more bound than their odd-nucleon-number neighbors. As a result, the one-nucleon drip line is reached earlier than the two-nucleon drip line, which results in a highly irregular pattern of nuclear existence that meanders between odd- and even-particle systems.

The territory of neutron-rich nuclei is arguably the most fertile ground for breakthroughs in nuclear structure research and the Ca region is of particular interest. The heaviest Ca isotope discovered to-date is 60Ca Tarasov et al. (2018). This nucleus, having $Z=20$ protons and $N=40$ neutrons, i.e., containing 12 more neutrons than the heaviest stable calcium isotope, was found recently together with seven other neutron-rich nuclei: 47P, 49S, 52Cl, 54Ar, 57K, 59Ca, and 62Sc. In addition, one event consistent with 59K was registered Tarasov et al. (2018). This discovery extends the range of known nuclei in this region, previously established in Refs. Tarasov et al. (2009, 2013). In separate experimental studies, the atomic masses of 55-57Ca were determined Michimasa et al. (2018) and the uncertainties of the 52-55Ti mass values were significantly reduced Leistenschneider et al. (2018).

The Ca region is arguably the most critical one to look at from a theory perspective, because it provides an exciting opportunity to bridge the refined methods based on realistic interactions, in which all $A$ nucleons are considered as elementary degrees of freedom, with nuclear density functional theory (DFT) employing energy density functionals (EDFs) expressed in terms of proton and neutron local densities and currents. During recent years, the $A$-body approaches reached selected medium-mass nuclei and provided predictions for their global properties and spectroscopy Hagen et al. (2014); Hergert et al. (2016). Nuclear DFT offers a more coarse-grained picture of nuclei than $A$-body approaches, but can be applied globally across the nuclear chart from light to superheavy nuclides Erler et al. (2012); Agbemava et al. (2014). The associated EDFs are primarily constrained by global nuclear properties such as binding energies and radii Bender et al. (2003). By considering symmetry-breaking effects, nuclear DFT can describe on the same footing spherical nuclei close to magic shells and deformed open-shell systems. A well-controlled link between $A$-body methods and DFT is essential if one aims to understand nuclei and nucleonic matter from a bottom-up perspective Forssén et al. (2013); Nazarewicz (2016).

When it comes to the Ca chain itself, $A$-body methods provide an excellent description of binding energies, charge radii, and spectroscopy up to 54Ca, depending on the interaction used Hagen et al. (2012); Holt et al. (2014); Somà et al. (2014); Hergert et al. (2014); Hagen et al. (2015); Garcia Ruiz et al. (2016); Simonis et al. (2017); Stroberg et al. (2017). Likewise, DFT approaches with globally-optimized EDFs reproduce measured global properties. However, there is no consensus when it comes to extrapolations towards the neutron drip line. Namely, $A$-body methods with two- and three-body interactions predict the two-neutron drip-line around 60CaHagen et al. (2013); Stroberg et al. (2017) while the DFT approaches locate it around 70CaErler et al. (2012); Forssén et al. (2013).

In this Letter, we investigate what global nuclear mass models, aided by Bayesian machine learning, can tell us about the topography of the mass surface and neutron drip-lines in the Ca region. Our methodology roughly follows the recent paper Neufcourt et al. (2018). Since Bayesian machine learning requires a sufficient number of data points in order to make extrapolations with reasonable certainty, one must work with models which are mostly global. To this end, we consider global models based on nuclear DFT with realistic Skyrme EDFs as well as the more phenomenological mass models FRDM-2012 and HFB-24 rooted in the mean-field theory.

Density functional theory calculations – We used the DFT mass predictions based on SkM∗ Bartel et al. (1982), SkP Dobaczewski et al. (1984), SLy4 Chabanat et al. (1995), SV-min Klüpfel et al. (2009), UNEDF0 Kortelainen et al. (2010), and UNEDF1 Kortelainen et al. (2012) EDFs stored in the theoretical database MassExplorer mas . The UNEDF2 Kortelainen et al. (2014) mass table has been computed exactly in the same way as in Ref. Erler et al. (2012). The DFT predictions are compared to the results of the global mass models FRDM-2012 Möller et al. (2016) and HFB-24 Goriely et al. (2013).

DFT calculations were carried out for even-even nuclei as we want to avoid additional complications and uncertainties related to the choice and treatment of quasi-particle configurations in odd-$A$ and odd-odd systems Bonneau et al. (2007); Schunck et al. (2010); Afanasjev (2015). Binding energies of odd-$A$ and odd-odd nuclei were obtained from the binding energy values and average pairing gaps computed for even-even neighbors. The associated error on $S_{1n}$ is expected to be 200-300 keV Schunck et al. (2010). For completeness, in our analysis we also considered even-even nuclei predicted to lie just beyond the two-neutron drip line, i.e., those with a slightly positive neutron chemical potential. Those results should be considered as rough estimates as the HFB theory does not guarantee that the nucleonic densities and fields are localized in this regime Dobaczewski et al. (1984); Dobaczewski and Nazarewicz (2013).

Statistical analysis – We first compute the so-called separation energy residuals, i.e., the differences $\delta_{1n/2n}(Z,N):=S_{1n/2n}^{\mathrm{exp}}(Z,N)-S_{1n/2n}^{\mathrm{th}}(Z,N)$, between experimental values and model predictions of $S_{1n/2n}$, based on the training datasets AME2003 Audi et al. (2003) and AME2016* consisting of AME2016 masses Wang et al. (2017) supplemented by the recently updated 52-55Ti masses Leistenschneider et al. (2018) (the subscripts $1n/2n$ are used to indicate either one-neutron or two-neutron separation energies). Using the values of $\delta_{1n/2n}(Z,N)$ for those training nuclei $(Z,N)$, we construct emulators $\delta_{1n/2n}^{\mathrm{stat}}(Z,N)$ using a Bayesian machine learning analysis of extrapolations via Gaussian processes (GP) following the methodology previously developed in Ref. Neufcourt et al. (2018). Our likelihood, the GP model, is a popular way Rasmussen and Williams (2006) of interpolating or extrapolating quantities from neighboring ones. It is strongly based on the assumption of a local spatial structure in the data, and contains required uncertainty modeling. We took the GP model in a form of a mean-zero Gaussian random field with a quadratic exponential spatial covariance kernel MacKay (2005) featuring three parameters: its scale $\eta$, which represents a noise intensity, and two characteristic spatial lengths $\rho$: one in the proton direction and one in the neutron direction; see Supplemental Material (SM) SM for details. We performed the statistical analysis independently on the sets of $S_{1n}$ and $S_{2n}$, respectively for odd-$N$ and even-$N$ nuclei, and independently for odd-$Z$ and even-$Z$ nuclei.

Posterior samples are obtained via 100,000 iterations of the Metropolis Hastings (1970) algorithm, from which the posterior mean value provides our predictions while Bayesian credibility intervals (CIs) are built using the corresponding posterior quantiles, symmetric around the mean, at all uncertainty levels. We also evaluate the performance of the prediction via a comparison of the rms deviation before and after statistical refinement. It is worth noting that our statistical CI estimates take into account all possible sources of uncertainty including statistical, numerical, and systematic uncertainty, including model approximations and modeling uncertainty within the DFT framework.

The unknown separation energies are predicted statistically by combining the theoretical predictions and the credibility intervals for the estimated residuals. For instance, the estimated prediction at $(Z,N)$ for the one-neutron (resp. two-neutron) separation energy is $S_{1n/2n}^{\mathrm{est}}(Z,N)=S_{1n/2n}^{\mathrm{th}}(Z,N,)+\delta_{1n/2n}^{\mathrm{stat}}(Z,N)$ where the last term is the Bayesian posterior mean prediction for the residual. A similar strategy of correcting model predictions outside the training domain by estimated residuals has recently been applied in Ref. Utama and Piekarewicz (2018). (Applications of statistical methods to mass predictions have been described in several papers Athanassopoulos et al. (2004); Bayram and Akkoyun (2017); Yuan (2016); Utama et al. (2016); Utama and Piekarewicz (2017); Bertsch and Bingham (2017); Zhang et al. (2017); Niu and Liang (2018), primarily in the context of interpolations.) The new aspect of our work lies in that we apply the Bayesian method to provide a full quantification of the uncertainty surrounding the point estimate. For more details we refer the reader to SM.

Results – GP’s superior predictive power was assessed in Ref. Neufcourt et al. (2018) for the $S_{2n}$ of even-even nuclei. The present work achieves comparable performances for odd-$Z$ nuclei and for $S_{1n}$ values, with prediction improvements ranging from 20 % to 40 % for most models (see SM). To further assess the performance of our approach, we apply it to the recently measured masses of 55-57Ca Michimasa et al. (2018). As seen in Fig. 1, the predicted $S_{1n}$ values for 55,57Ca are consistent with experiment for most models while the $S_{2n}$ of 56Ca is slightly overestimated. The impact of newer mass measurements beyond AME2003 on our predictions is minor; this is because very few datapoints that can impact our local GP model were added in the Ca region. The large deviation in the $S_{1n}$ of 55Ca in HFB-24 is noteworthy. As illustrated in SM and Ref. Neufcourt et al. (2018), neutron separation energies predicted by this model often exhibit irregular behavior.

Figure 2 shows extrapolated separation energies for the Ca isotopic chain for three global mass models corrected with the GP emulator. (Here and in the following we shall use the notation “model+GP” (e.g., UNEDF0+GP) to emphasize that the statistical corrections are done with the GP emulator.) The models are consistent overall once the statistical correction and uncertainty are taken into account. According to the computed empirical coverage probabilities Gneiting and Raftery (2007); Gneiting et al. (2007), our credibility intervals are slightly conservative for large credibility levels (see Sec. I.C of SM for more discussion).

For a given isotopic chain and nuclear model, one obtains an upper bound on the location of the first isotope at which the binding energy becomes negative, depending on the choice of credibility level. For instance, the posterior mean values (full lines) of the UNEDF0+GP model place the $2n$ drip line for Ca around $N=54$, while considering the lower bound of the one-sigma credibility intervals provides that it is placed beyond $N=46$ with probability 84%. This very wide interval suggests that the posterior distribution of the separation energies is perhaps not the most appropriate quantity to consider. To this end, for each model, we consider the probability $p_{ex}(Z,N)$ of the predicted separation energy $S_{1n/2n}^{*}(Z,N)$ to be positive under the posterior probability distribution conditioned on the experimental masses available. In the Bayesian paradigm, this probability is $p_{ex}(Z,N):=p(S_{1n/2n}^{*}(Z,N)>0|S_{1n/2n})$. The insert in Fig. 2 shows $p_{ex}$ for the Ca chain. The model-averaged existence probabilities for the Ca region are shown in Fig. 3(a) assuming uniform prior weights. (For the values of $p_{ex}$ for individual models, see Sec. III.C of SM.) As noticed in Ref.Tarasov et al. (2018), the $N=35$ isotones 52Cl and 53Ar, as well as 49S represent a challenge for nuclear mass models. Our results in Fig. 3(a) confirm this finding through the low calculated prior-average $p_{ex}$ values for these nuclei. Indeed, with the exception of SV-min, UNEDF0, and FRDM-2012, other models calculate them to be either marginally bound or to lie outside the one-neutron drip line. Since 49S, 52Cl and 53Ar do exist Tarasov et al. (2009, 2018), this prior knowledge can inform the model averaging process Hoeting et al. (1999); Wasserman (2000); Bernardo and Smith (1994) through posterior weights:

[TABLE]

(see additional discussion in SM). The weight $w_{k}$ reflects the ability of the model $\mathcal{M}_{k}$ to predict the existence of nuclei in the Ca region. In this respect UNEDF0+GP is superior, see Table S1 in SM. We emphasize that conditioning with respect to these three nuclei corresponds actually to conditioning over the observed nuclei in the whole Ca region, since other experimentally-observed isotopes are predicted to be bound by the global models considered. The values of $p_{ex}$ obtained in this way are shown in Fig. 3(b).

As shown in Figs. 2 and 3, the nucleus 68Ca is expected to be bound. However, as seen in Fig. 2, $S_{2n}$ approaches zero very gradually; this results in a spread of predictions of individual models. According to the average $p_{ex}$, 61Ca and 71Ti are expected to be $1n$-unstable while the $2n$ drip line extends all the way to 72Ca and 78Ti. The nucleus 59K – for which one event was registered in Ref. Tarasov et al. (2018) – is expected to be firmly neutron-bound. By comparing Figs. 3(a) and (b) one can immediately assess the impact of the discovery of 52Cl, 53Ar, and 49S on drip-line predictions: the $2n$ drip line obtained with posterior weights generally extends by two neutron numbers for odd-$Z$ chains.

Conclusions – In summary, in this Letter we quantified the neutron-stability of the nucleus in terms of its existence probability $p_{ex}$, i.e., the Bayesian posterior probability that the neutron separation energy is positive. Our results are fairly consistent with recent experimental findings Tarasov et al. (2018): 60Ca is expected to be well bound with $S_{2n}\approx 5$ MeV while 49S, 52Cl, and 53Ar are marginally-bound threshold systems.

We emphasize that the nuclear model itself is not capable of gauging the likelihood of existence. To overcome this problem, we introduce a machine learning algorithm, with a stochastic exploration part and a deterministic modeling part, which, when combined, result in Bayesian statistical machine learning. One could say this is supervised learning, with the nuclear modeling and the choice of priors representing two aspects of the supervision.

The Bayesian model averaging employed in this Letter is based on global DFT/mean-field models. Therefore the computed probabilities of existence are conditional on the correctness of the DFT framework. Currently, many $A$-body methods based on realistic inter-nucleon interactions calculate the two-neutron drip line at 60Ca. Since Bayesian machine learning requires a sufficient number of data points to extrapolate with reasonable certainty, $A$-body models are not yet amenable to statistical analysis as the corresponding global mass tables are difficult to compute. It will be extremely valuable to apply a Bayesian uncertainty quantification analysis to $A$-body mass tables when those become available.

The extrapolation outcomes discussed in this Letter will be tested by experimental data from rare-isotope facilities. New mass measurements on neutron-rich nuclei will help to develop increasingly more quantitative models of the atomic nucleus and also allow for a higher-fidelity statistical analysis. As illuminated by our Bayesian analysis of 49S, 52Cl, and 53Ar, experimental discoveries of new nuclides will also be crucial for delineating the detailed behavior of the nuclear mass surface, including the placement of particle drip lines.