Ab initio limits of atomic nuclei

TL;DR

This paper uses advanced ab initio methods with chiral interactions to predict the limits of nuclear existence across a wide range of elements, providing theoretical uncertainties and new predictions for neutron-rich isotopes.

Contribution

It extends ab initio nuclear calculations to medium-mass nuclei up to iron, predicting drip lines with quantified uncertainties and offering new testable predictions.

Findings

Predicted neutron and proton drip lines across light to medium-mass nuclei.

Calculated ground-state and separation energies for nearly 700 isotopes.

Predictions for neutron-rich isotopes to guide future experiments.

Abstract

We predict the limits of existence of atomic nuclei, the proton and neutron drip lines, from the light through medium-mass regions. Starting from a chiral two- and three-nucleon interaction with good saturation properties, we use the valence-space in-medium similarity renormalization group to calculate ground-state and separation energies from helium to iron, nearly 700 isotopes in total. We use the available experimental data to quantify the theoretical uncertainties for our ab initio calculations towards the drip lines. Where the drip lines are known experimentally, our predictions are consistent within the estimated uncertainty. For the neutron-rich sodium to chromium isotopes, we provide predictions to be tested at rare-isotope beam facilities.

Equations

$$\delta S=f(N,Z,S^{\mathrm{th}},\dots )+ϵ({\sigma }^2),$$

$$\delta S\approx -2\frac{\delta R}{R}S+const.,$$

$$\delta S=\mathcal{A}{S}^{\mathrm{th}}+\mathcal{B}+ϵ({\sigma }^2)$$

$${\mathcal{P}}_{\mathrm{bound}}=\prod _{\alpha }{\int }_0^{\infty }d{\tilde{S}}_{\alpha }^{\exp }p({\tilde{S}}_{\alpha }^{\exp }|{\tilde{S}}^{\mathrm{th}},S^{\mathrm{th}},S^{\exp }),$$

$$\delta S_{i,\alpha }={\mathcal{A}}_{\alpha }{S}_{i,\alpha }^{\mathrm{th}}+{\mathcal{B}}_{\alpha }+{ϵ}_{i,\alpha }({\sigma }_{\alpha }^2).$$

$${\mathbf{y}}_{\alpha }={\left(\begin{array}{cccc}\delta S_{1,\alpha }& \delta S_{2,\alpha }& \dots & \delta S_{m,\alpha }\end{array}\right)}^T,$$

$$X_{\alpha }={\left(\begin{array}{cccc}{S}_{1,\alpha }^{\mathrm{th}}& S_{2,\alpha }^{\mathrm{th}}& \dots & S_{m,\alpha }^{\mathrm{th}}\\ 1& 1& \dots & 1\end{array}\right)}^T.$$

$${\mathbf{y}}_{\alpha }=X_{\alpha }{\mathbf{\beta }}_{\alpha }+{\mathbf{ϵ}}_{\alpha }({\sigma }_{\alpha }^2),$$

$$p(\mathbf{\beta },{\sigma }^2)=p({\sigma }^2)p(\mathbf{\beta }|{\sigma }^2),$$

$$p({\sigma }^2)=\frac{b_0^{a_0}}{\Gamma (a_0)}({\sigma }^2{)}^{-a_0-1}\exp (-b_0/{\sigma }^2),$$

$$p(\mathbf{\beta }|{\sigma }^2)\propto \frac{1}{{\sigma }^2}\exp \left[-\frac{(\mathbf{\beta }-{\mathbf{\mu }}_0{)}^T{\Lambda }_0(\mathbf{\beta }-{\mathbf{\mu }}_0)}{2{\sigma }^2}\right].$$

$$p(\mathbf{y}|\mathbf{\beta },{\sigma }^2,X)=$$

$${\left(\frac{1}{2\pi {\sigma }^2}\right)}^{m/2}$$

$$p(\mathbf{\beta },{\sigma }^2|\mathbf{y},X)\propto p(\mathbf{\beta },{\sigma }^2)p(\mathbf{y}|\mathbf{\beta },{\sigma }^2,X).$$

$${\Lambda }_{\ast }={\Lambda }_0+X^{T}X$$

$${\mathbf{\mu }}_{\ast }={\Lambda }_{\ast }^{-1}({\Lambda }_0{\mathbf{\mu }}_0+X^{T}y)$$

$$a_{\ast }=a_0+\frac{m}{2}$$

$$b_{\ast }=b_0+\frac{1}{2}({{\mathbf{\mu }}_0}^T{\Lambda }_0{\mathbf{\mu }}_0+y^{T}y-{\mu }_{\ast }^T{\Lambda }_{\ast }{\mu }_{\ast }).$$

$$p(\tilde{\mathbf{y}}|\mathbf{y},X,\tilde{X})=\int d\mathbf{\beta }d{\sigma }^{2}p(\tilde{\mathbf{y}}|\mathbf{\beta },{\sigma }^2,\tilde{X})p(\mathbf{\beta },{\sigma }^2|\mathbf{y},X).$$

$$p(\tilde{\mathbf{y}}|\mathbf{y},X,\tilde{X})\approx \mathcal{N}(\tilde{\mathbf{M}},\tilde{\Sigma })$$

$$\tilde{\mathbf{M}}=\tilde{X}{\mathbf{\mu }}_{\ast }$$

$$\tilde{\Sigma }=\frac{b_{\ast }}{a_{\ast }}(\mathbb{I}+\tilde{X}{\Lambda }_{\ast }^{-1}{\tilde{X}}^T).$$

Full text

Atomic nuclei, which form the basis for known matter in the Universe, cannot be made from arbitrary numbers of protons and neutrons. For a given element (i.e., proton number $Z$) a nucleus can support only so many neutrons, $N$, and vice versa. The point at which nucleons no longer form a bound system is referred to as the drip line. Specifically, at the drip line one- or two-nucleon separation energies become negative, and nuclei decay via nucleon emission. The proton drip line is known experimentally to the medium-mass region, but to date, the neutron drip line is established only up to neon ($Z=10$) Thoennessen (2016); Ahn et al. (2019). Pinning down the neutron drip line to calcium and beyond is a flagship scientific motivation for next-generation rare-isotope beam facilities Gade and Sherrill (2016); Durante et al. (2019). Indeed several neutron-rich isotopes, including 60Ca, were recently discovered in this region Tarasov et al. (2018). Furthermore, knowledge of the neutron drip line is important for $r$-process simulations modeling the synthesis of heavy elements Mumpower et al. (2016); Martin et al. (2016) that occurs in neutron-star mergers Abbott et al. (2017).

Predicting the location of the drip lines poses a substantial theoretical challenge, particularly because many nuclei far from known data must be calculated systematically. In a pioneering study, the nuclear landscape was predicted from extrapolations of state-of-the-art nuclear density functional theory, and approximately 7000 nuclei were estimated to exist in nature Erler et al. (2012). Since this work, tremendous progress has been made in statistical analyses of nuclear models Neufcourt et al. (2019); Yoshida et al. (2018) as well as in ab initio nuclear theory. Developments in chiral effective field theory Epelbaum et al. (2009); Machleidt and Entem (2011); Hammer et al. (2013) and similarity renormalization group Bogner et al. (2007, 2010) are pushing nuclear forces to new levels of accuracy and ranges of applicability. Though a robust and systematic theoretical framework has not yet been fully achieved, particular nuclear Hamiltonians have been constructed which reproduce ground-state energies up to the tin region Hebeler et al. (2011); Ekström et al. (2015); Morris et al. (2018). Three-nucleon (3N) forces play a key role for understanding the drip lines Otsuka et al. (2010); Hagen et al. (2012); Cipollone et al. (2013); Hergert et al. (2013); Holt et al. (2013); Hebeler et al. (2015). Moreover, many-body theories Barrett et al. (2013); Hagen et al. (2014); Carlson et al. (2015); Hergert et al. (2016); Stroberg et al. (2019); Barbieri and Carbone (2017) have advanced to treat medium-mass open-shell systems Hergert et al. (2013); Bogner et al. (2014); Jansen et al. (2014); Somà et al. (2014); Stroberg et al. (2017), with the primary limitation being computational resources needed to obtain convergence with respect to basis size, laying the groundwork for a new era of ab initio theory.

In this Letter we calculate properties of essentially all nuclei from helium to iron ($Z=2-26$), close to 700 in total, to provide a global ab initio survey of ground-state energies and predict the nuclear drip lines. Using two-nucleon (NN) and 3N interactions constrained by only few-body data, we solve the many-body problem with the valence-space formulation of the in-medium similarity renormalization group (VS-IMSRG) Tsukiyama et al. (2012); Hergert et al. (2016); Bogner et al. (2014); Stroberg et al. (2016, 2017, 2019). Our results yield an overall root-mean-square (rms) deviation of 3.3 MeV from absolute experimental energies and 0.7-1.4 MeV from separation energies. In comparison, state-of-the-art energy-density functionals obtain rms devations in the range 0.6–2.0 MeV for energies and 0.4–1.25 MeV for separation energies Goriely et al. (2009); Bulgac et al. (2018); Navarro Pérez et al. (2018); Neufcourt et al. (2018) (note, however, that the density functional rms values are obtained over a much larger range of masses).

While the drip line signature is unambiguous experimentally, from a theoretical perspective an error of a few tens of keV – well beyond current levels of precision – can make the difference between an isotope being bound or unbound. Therefore, an assessment of theoretical uncertainty is mandatory for any meaningful drip line prediction. Ab initio methods present an appealing framework for uncertainty quantification: one begins with the most general Lagrangian compatible with the applicable symmetries, organized by a systematically improvable power counting, then solves the nuclear many-body problem within a controlled and systematically improvable approximation scheme, propagating all uncertainties. Such a prescription has not yet been achieved in practice, so for the present we use a comparison with known data to calibrate a physically motivated model for the error. Recent work in a similar spirit has applied Bayesian machine learning algorithms to global mass models Neufcourt et al. (2018, 2019, 2020). The main advantages of our current approach are (i) the predictions should not be biased towards measured data, because they were not fit to any data beyond helium and (ii) the predictions can be benchmarked where the proton and neutron drip lines are known experimentally (mass models are typically applied to $Z\gtrsim 8$).

In the VS-IMSRG, a valence-space Hamiltonian of tractable dimension is decoupled from the larger Hilbert space via an approximate unitary transformation. We begin in a harmonic-oscillator basis of 15 major shells (i.e., $e=2n+l\leqslant e_{\mathrm{max}}=14$) with an imposed cut of $e_{1}+e_{2}+e_{3}\leqslant E_{\mathrm{3Max}}=16$ for 3N matrix elements. The resulting ground-state energies are converged to better than a few hundred keV with respect to these truncations, and we perform extrapolations in $e_{\rm max}$ to obtain infrared convergence Furnstahl et al. (2012, 2015). Transforming to the Hartree-Fock basis, we capture effects of 3N interactions between valence nucleons via the ensemble normal-ordering of Ref. Stroberg et al. (2017). We then use the Magnus formulation of the IMSRG Morris et al. (2015); Hergert et al. (2016), truncating all operators at the normal-ordered two-body level—the IMSRG(2) approximation—to generate approximate unitary transformations that decouple the core energy and valence-space Hamiltonian for each nucleus to be calculated.

By default, we employ a so-called $0\hbar\omega$ valence space, where valence nucleons occupy the appropriate single major harmonic-oscillator shell (e.g., for $8<N(Z)<20$ the $sd$ shell, $20<N(Z)<40$ the $pf$ shell, etc.). At $N(Z)=2,8,20,40$, we do not decouple a neutron (proton) valence space, and no explicit neutron (proton) excitations are allowed in the calculation. We discuss exceptions to this below. Finally the resulting valence-space Hamiltonians are diagonalized with the NuShellX@MSU shell-model code Brown and Rae (2014) (with the exception of a few of the heaviest Ca, Sc and Ti isotopes, which were computed with the m-scheme code K-shell Shimizu et al. (2019)).

We thus calculate ground (and excited) states of all nuclei from helium to iron, except those for which the shell-model diagonalization is beyond our computational limits. For the input NN+3N interaction, we use the 1.8/2.0 (EM) potential of Refs. Hebeler et al. (2011); Simonis et al. (2016), where the 3N couplings were fit to the 3H binding energy and the 4He charge radius. This interaction reproduces experimental ground-state energies of light- to medium-mass nuclei remarkably well Simonis et al. (2017); Morris et al. (2018). Studies of nuclear matter Hebeler et al. (2011); Drischler et al. (2019) have shown that this interaction saturates with slightly too much binding and at at somewhat too high density, leading to too small radii for finite nuclei Simonis et al. (2017). We use this observation below to model our systematic error in the separation energies. In the Supplemental Material, we provide results for absolute and separation energies for all nuclei calculated. In the inset of Fig. 1, we plot the range of agreement with experiment and find an overall rms deviation of 3.3 MeV. The experimental binding and separation energies are taken from the Atomic Mass Evaluation Wang et al. (2017), with additional recent data from Refs. Webb et al. (2019); Leblond et al. (2018); Mukha et al. (2018); Michimasa et al. (2018).

As no experimental input beyond 4He is used in the current calculations, our results should not be biased toward known data111One might argue that the selection of one interaction from the family of interactions used in Simonis et al. (2017) constitutes incorporation of information beyond $A=4$. However, in the cases we have checked, other interactions in the family yield similar separation energies.. Therefore our approach is to use measured data to assess how well separation energies are reproduced in general, then assume our calculations will behave similarly for separation energies which have not yet been measured. This neglects phenomena which may emerge in the neutron-rich region, such as halo structures or island-of-inversion physics Warburton et al. (1990); Caurier et al. (2014), but our results suggest that the impact of these effects on separation energies tends to lie within our estimated uncertainties.

To characterize the quality of the reproduction of experimental data, we assess the residual, $\delta S_{\alpha}=S_{\alpha}^{\rm{th}}-S_{\alpha}^{\rm{exp}}$, for the separation energy in channel $\alpha\in\{n,p,2n,2p\}$. We model the residual as

[TABLE]

where $f$ is a function characterizing the systematic error (see below) and $\epsilon$ is a random number drawn from a Gaussian distribution of mean 0 and variance $\sigma^{2}$.

The main source of many-body error in our calculations is due to the IMSRG(2) approximation. For soft input interactions such as the one used here, this approximation is accurate for binding energies at the level of a few percent Stroberg et al. (2017). Given that binding energies in this region are a few hundred MeV, this would naively suggest an error of several MeV on the separation energies. However, the errors made for neighboring nuclei are strongly correlated (see Supplementary Material) and largely cancel, improving the accuracy of the separation energies (the uncorrelated part will contribute to the random error $\epsilon$). This correlation is deteriorated in the case where a different valence space is used for the two binding energies entering into the separation energy. Wherever possible, we use a consistent valence space to compute separation energies, but this is not always possible. For example, to compute the $S_{2n}$ of 38Cl, we require the ground-state energy of 38Cl ($N$=21), which has a valence neutron in the $pf$ shell, and 36Cl ($N$=19), which has a valence hole in the $sd$ shell. This leads to an increased error for these special cases, and we treat these separately.

We must also consider errors due to deficiencies of the input Hamiltonian. As mentioned above, a notable deficiency of this Hamiltonian is that predicted radii are systematically smaller than experiment by $\approx 4$% for the mass range considered Simonis et al. (2017). This should have an impact on the computed separation energies. In an infinite potential well, decreasing the width of the potential well will spread out the spectrum of single-particle energies, which depend on the radius as $1/R^{2}$, and we expect

[TABLE]

where the radius shift $\delta R/R\approx-0.04$, and the constant shift depends on the potential well depth and any shifts thereof. We therefore anticipate that the residual $\delta S$ will exhibit a linear dependence on the separation energy $S$ with a slope of approximately $0.08$. Continuum and other effects will of course modify this simple picture. However, the angular momentum and Coulomb barriers, as well as the expansion on a harmonic oscillator basis keep deviations from this linear behavior within the error band, even down to slightly negative separation energies (see discussion below and Fig. 2).

We therefore perform a Bayesian linear regression for the model

[TABLE]

to obtain a posterior $p(\mathcal{A},\mathcal{B},\sigma^{2}|S^{\rm th},S^{\rm exp})$ for the unknown parameters $\mathcal{A},\mathcal{B},\sigma^{2}$, given the theoretical and experimental data. We then marginalize over the posterior to obtain a posterior predictive distribution (PPD) $p(\tilde{S}^{\rm exp}|\tilde{S}^{\rm th},S^{\rm th},S^{\rm exp})$ for a not-yet-measured separation energy $\tilde{S}^{\rm exp}$, given the theoretical value, and all the known data. Details can be found in the Supplemental Material. The result is shown in Fig. 2. As expected based on the discussion of Eq. (2), we obtain a slope of approximately 0.08 for each separation energy. Importantly, if we did not account for this systematic effect our uncertainty assessment would be biased by the relative abundance of data on well-bound isotopes.

The cases requiring inconsistent valence spaces (which were not used in the regression) are marked with crosses in Fig. 2. In these cases the error due to incomplete cancellation of induced many-body effects is more difficult to model and so we are more conservative. For each valence-space boundary at $N,Z=8,20,40$ we apply an additional $S$-independent shift to the PPD—equal to the mean deviation from the regression line—and inflate the standard deviation by the size of the shift.

The probability $\mathcal{P}_{1n}$ that an isotope is bound with respect to one-neutron emission is given by the fraction of the PPD for which $S_{n}$ is positive. The total probability to be bound is given by the fraction of the joint PPD for which all four separation energies are positive,

[TABLE]

where as above $\alpha\in\{n,p,2n,2p\}$. As an illustration, the calculated separation energies, with the 68% uncertainty band, are shown in Fig. 3, for chlorine isotopes. Analogous figures for all isotopic chains studied are included in the Supplemental Material, and a complete data table is provided as a Supplemental File.

We translate this large-scale analysis into the main result of this Letter in Fig. 1. For each calculated nuclide from helium to iron, we assign a probability that it is bound with respect to one- or two-nucleon emission. Every nuclide calculated is represented by a box in the plot, where its height and color denotes this probability: a full box is bound with probability 1, and an empty gray box is bound with probability 0. For ease of interpretation, we employ a color code with divisions at $\mathcal{P}_{\rm bound}=(0.05,0.32,0.68,0.95)$. We denote experimentally known drip lines with a filled symbol and the heaviest (lightest) observed isotopes in the neutron-rich (deficient) regions with an open symbol.

Qualitatively, the location of the known drip lines appear to be reproduced well in Fig. 1, both on the proton-rich side, and on the neutron-rich side where it is known up to $Z=10$. (Further quantitative validation of the approach is presented in Supplemental Material). Even the well known halo systems 11Li and 22C are predicted to be either bound or marginal, implying that physics of threshold systems lie within our estimated error bands. For all isotopic chains from sodium to chromium, our calculations indicate the likely existence of at least one isotope beyond the current known limits.

In calcium, earlier ab initio calculations have generally found that 62Ca Hagen et al. (2013); Holt et al. (2014) is the heaviest bound isotope (see also Ref. Hergert (2020)), with a very flat trend in binding energies beyond, leaving the location of the drip line ambiguous. The present analysis reflects that ambiguity; similar to oxygen the final bound nucleus could be closer to stability, but there is a reasonable probability that the drip line extends beyond 70Ca, as predicted in the statistical analysis of Ref. Neufcourt et al. (2019). We note the remarkable similarities of the latter results to our ab initio predictions, which thus provides a consistent picture of the neutron drip line up to calcium from independent theoretical approaches.

In summary, we have calculated ground-state energies of essentially all nuclei from helium to iron in the ab initio VS-IMSRG starting from NN and 3N interactions fit to few-body systems only. Using available experimental data to quantify our theoretical error, we provide drip line predictions in the neutron-rich region above neon to guide ongoing and future efforts at rare-isotope beam facilities worldwide. This work also advances ab initio theory to global calculations, highlighting the rapidly increasing scope of the field, and the potential to provide predictions beyond where data exists with uncertainty estimates. In principle, we might have expected a global survey to uncover deficiencies in the 1.8/2.0 (EM) interaction which were not apparent based on prior calculations of closed-shell nuclei or selected isotopic chains; however, we globally find an impressive agreement. While significant challenges remain in improving the rigor of theoretical error estimates from nuclear forces and many-body methods, the approach presented here indicates a path – once current computational limitations can be overcome – for ab initio input for nucleosynthesis calculations probing the r-process region of extreme neutron-rich nuclei.

I Supplemental Material

I.1 Bayesian linear regression

Our assessment of the residuals $\delta S$ employs a Bayesian linear regression with conjugate priors (see, e.g., Refs. Gelman et al. (2014); Fahrmeir et al. (2009)). Our linear model is

[TABLE]

The index $i$ labels a specific nuclide, and the index $\alpha\in\{n,p,2n,2p\}$ labels a particular separation energy channel. The coefficients $\mathcal{A}_{\alpha}$, $\mathcal{B}_{\alpha}$, and the variances $\sigma^{2}_{\alpha}$ are unknown parameters with distributions to be determined. The term $\epsilon_{i,\alpha}(\sigma^{2}_{\alpha})$ denotes a random number drawn from a Gaussian distribution of mean zero and variance $\sigma^{2}_{\alpha}$. We collect $\mathcal{A}_{\alpha}$ and $\mathcal{B}_{\alpha}$ into a two-dimensional column vector $\bm{\beta}_{\alpha}=\begin{pmatrix}\mathcal{A}_{\alpha}&\mathcal{B}_{\alpha}\end{pmatrix}^{T}$. Given $m$ data points, the residuals for a given channel $\alpha$ form an $m$-dimensional column vector

[TABLE]

and the ab initio separation energies $S^{\rm th}_{i,\alpha}$ are collected in an $m\times 2$ matrix

[TABLE]

Our linear model becomes

[TABLE]

where the bold symbol $\bm{\epsilon}_{\alpha}$ indicates a column vector of $m$ random numbers drawn from a Gaussian of variance $\sigma^{2}_{\alpha}$.

In the following, we focus on a single channel and suppress the label $\alpha$ in order to avoid cluttering the notation. Our priors for $\bm{\beta},\sigma^{2}$ are taken to be a normal-inverse-gamma distribution

[TABLE]

with $p(\sigma^{2})$ given by an inverse-gamma distribution

[TABLE]

where $\Gamma$ is the gamma function, and $p(\beta|\sigma^{2})$ a multivariate normal distribution

[TABLE]

The parameters $a_{0},b_{0},\bm{\mu_{0}},\Lambda_{0}$ characterize our priors. We choose relatively noninformative values of $a_{0}=b_{0}=1$, $\bm{\mu_{0}}=\begin{pmatrix}0&0\end{pmatrix}^{T}$, and $\Lambda_{0}=\begin{pmatrix}0&0\\ 0&0\end{pmatrix}$. The likelihood function for $\bm{y}$ given $\bm{\beta},\sigma^{2},X$ is a multivariate Gaussian

[TABLE]

The posterior distribution for the parameters $\bm{\beta},\sigma^{2}$ is obtained with Bayes’ theorem

[TABLE]

Due to our choice of conjugate priors, the posterior distribution has the same functional form as our priors, but with updated parameters Fahrmeir et al. (2009), which we indicate with a subscript $*$:

[TABLE]

The posterior predictive distribution (PPD) for new data $\bm{\tilde{y}}$ given the training data $\bm{y},X$ and corresponding ab initio calculations $\tilde{X}$ is

[TABLE]

This integral can be evaluated analytically, resulting in a multivariate Student-t distribution with $2a_{*}$ degrees of freedom. For many degrees of freedom, the Student-t distribution approaches a normal distribution, and given that in the present case $m>100$, we take the PPD to be a multivariate normal. The PPD for new data $\bm{\tilde{y}}$, given the training data $\bm{y},X$, and additional ab initio calculations $\tilde{X}$ is therefore

[TABLE]

with mean vector

[TABLE]

and covariance matrix

[TABLE]

The above formulation does not account for experimental error bars and, particularly near the drip lines, experimental errors can be comparable to the theoretical uncertainty. To preserve the advantages of conjugate priors, we incorporate the experimental uncertainties by Monte Carlo sampling the experimental data, interpreting the reported error bars as indicating one standard deviation of a normal distribution. The parameters $\Lambda_{*}$ and $a_{*}$ are unaffected by this procedure, since they do not depend on $\bm{y}$, but we obtain distributions for $b_{*}$ and $\bm{\mu_{*}}$. These distributions are unimodal and relatively sharply peaked. We have confirmed that taking the mean values of these distributions yields a PPD indistinguishable from the PPD obtained by marginalizing over their distributions. (The resulting PPD is, however, slightly broader than what would be obtained by neglecting the experimental uncertainties altogether).

I.2 Correlation of residuals

between neighboring nuclides

We have claimed that the error due to the IMSRG(2) truncation leads to a many-body error which is correlated between neighboring nuclei, thus leading to a greater precision in the separation energies than would be naively estimated. To support this claim, we plot in Fig. 4 the ground-state energy residuals $\delta E=E^{\rm th}-E^{\rm exp}$ for neighboring nuclides. We find that the errors in neighboring nuclides are indeed strongly correlated, with a correlation coefficient $\rho_{xy}=\frac{\langle xy\rangle-\langle x\rangle\langle y\rangle}{\sigma_{x}\sigma_{y}}\approx 0.9$. (The angle brackets $\langle\,\cdots\rangle$ denote the mean value).

I.3 Correlation of residuals for a given nuclide

For a given nuclide, the residuals $\delta S_{n}$ and $\delta S_{2n}$ are correlated, and this should be accounted for when estimating the probability that a given isotope is bound with respect to both $1\mathrm{n}$ and $2\mathrm{n}$ emission. This discussion also applies to the proton emission case. The proton and neutron separation energy residuals are also correlated (e.g., $\delta S_{n}$ is correlated with $\delta S_{p}$) but this has no impact on the drip line. As a result, in the actual calculation, the four-dimensional integral in  (4) can be factored into a product of two-dimensional integrals. The impact of correlations in the residuals reaches at most 0.1 (i.e., changing the probability by 10%) and is typically below 0.01.

I.4 Validation

To validate the reliability of our approach, we perform an empirical coverage test. We partition our data into “training” and “validation” sets, apply our analysis to the training set, and then bin the validation data according to the predicted probability to be bound. In each bin, we evaluate what fraction of the data points are experimentally bound. If our predictions are reliable, the fraction bound should coincide with the mean value of the bin.

We perform the partitioning into training and validation sets in two different ways. In the first way, we take as training data all nuclides for which all theoretical separation energies are above a cutoff $S_{\rm cut}$, and validate using the remaining cases. The first three panels of Fig. 5 show the results for $S_{\rm cut}=1,2,5$ MeV. The second way we randomly partition the data into either set with equal probability. The result of this test is shown in the fourth panel of Fig. 5, labeled “random”. The vertical error bars represent a 68% confidence interval, due to finite sample size, estimated using a Wilson score with continuity correction (method 4 from Ref. Newcombe (1998)). The error bars do not get noticeably smaller with increasing $S^{\rm cut}$ because nearly all the additional validation points end up in the upper-right corner with a probability to be bound of approximately 1. Importantly, the $S_{\rm cut}=5$ MeV figure demonstrates that we can reliably use well bound nuclei to predict the drip line.

I.5 Separation energy plots

We provide plots of the separation energies, equivalent to Fig. 3 in the main text, for all the elements studied here ($2\leqslant Z\leqslant 26$). In addition, we plot in Fig. 6 the range of agreement with the experimental ground-state energies as a function of $Z$ to complement the information given in the inset to Fig. 1 as a function of $N$.