Introduction

The search for the origin of heavy elements in the Universe remains one of the most formidable challenges in scientific inquiry [1, 2]. While it is widely accepted that approximately half of the elements heavier than iron are produced through the rapid neutron capture process (r-process), the realistic astrophysical sites responsible for r-process are still under debate [3-5]. Calculations of r-process nucleosynthesis involve the properties of thousands of nuclei, spanning from the line of β stability to the neutron dripline. Among these properties, the nuclear mass is arguably the most crucial. On one hand, they affect the energy generation and the reaction path; on the other hand, they influence various properties of nuclei that enter the astrophysical calculations, such as β-decay rates, neutron-capture rates, and β-delayed neutron emission probabilities [6]. Sensitivity study underscores that mass uncertainties of 0.5 MeV significantly impact r-process abundance predictions [7]. A recent study explored how statistical uncertainties in mass models propagate to simulations of r-process nucleosynthesis [8].

Despite considerable efforts to measure the masses of nuclei near the neutron dripline [9, 10], many of which along the r-process path remain inaccessible even with the most advanced radioactive ion beam facilities [11]. Theoretical mass models, although originating from different approaches, are thus the last resort for obtaining mass information on unknown nuclei [12-14], forming the basis for nucleosynthesis calculations. Although most mass models fit well with experimental mass data, they diverge dramatically when extrapolating towards unknown regions [12]. As pointed out in Ref. [15], the accuracy and predictive power of mass models do not always correlate well. Furthermore, the uncertainty in mass models is difficult to quantify, hindering a comprehensive evaluation of their performance in extrapolation. Quantifying uncertainties to assess the predictive capability of mass models is highly demanding [16, 17]. In the 1990s, it was demonstrated that neural networks could be employed to predict the masses and other properties of nuclei [18, 19]. Along this line, neural networks have started to play an important role in predicting nuclear masses [20-25]. Currently, different machine learning approaches equipped with Bayesian methods have gained much attention and are widely used in the nuclear physics community (see Refs. [26-28] and references therein). Other machine learning algorithms, such as kernel ridge regression, have also begun to play an important role in the description of nuclear masses and separation energies [29, 30].

Bayesian-based methods can largely avoid overfitting by incorporating priors over parameters [31]. Bayesian analysis presents its results as a posterior probability distribution over parameters and variables, allowing the mean and variance to be calculated using statistical approaches. Bayesian Neural Networks (BNNs) have been successfully applied to refine model predictions of nuclear masses [32-35], α-decay energies [36], and β-decay half-lives [37]. However, BNNs are often not easy to use in practice because they require many decisions regarding activation functions and learning rates. Moreover, they involve hundreds of parameters that are difficult to interpret. An alternative is to use Gaussian Processes (GPs) [38], which have only a few parameters. GPs have been used to improve model predictions of separation energies [39], nuclear masses [40], and to estimate the limits of the nuclear landscape [41]. Beyond mass models, Bayesian Monte Carlo frameworks have been employed for nuclear data evaluation [42].

Recently, GPs have been extensively employed to refine nuclear mass predictions [43-45]. While these studies demonstrated the power of GPs for mass modeling, our work extends beyond single-model refinement by systematically improving eight commonly used theoretical mass models through their residuals. This ensemble approach combines the strengths of each model while compensating for their individual limitations, which is a crucial advancement for astrophysical applications where extrapolation reliability is essential. To rigorously test the predictive power of our GP, we introduce a novel theory-to-theory validation framework using FRDM12 as “pseudo-experimental” data. This allows us to evaluate the extrapolation performance far beyond the known regions. Furthermore, we provide complete mass predictions for all unmeasured nuclei listed in AME2020 for individual GP-refined mass models, which can serve as a comprehensive source for both astrophysical modeling and future experimental research.

This paper is organized as follows: In Sect. 2, we discuss the current mass models that are widely used in nuclear physics and their comparisons with experimental data. In Sect. 3, we outline the concept of GP regression and its implementation in mass models. In Sect. 4, we illustrate the performance of the GP models. In Sect. 5, we fine-tune our GP model and test its predictive capability. Finally, we conclude our work in Sect. 6 and offer future perspectives.

Nuclear Mass Model

Current global mass models can be categorized into three classes [12]: macroscopic-microscopic, microscopic, and phenomenological. The first nuclear mass model, the Bethe-Weizsäcker formula, treats the nucleus as a charged liquid drop [46, 47]. While it can predict the binding energies for most nuclei, it fails to describe those near the magic numbers. Based on the images of liquid drops, macroscopic-microscopic mass models incorporate shell corrections and pairing effects. Two widely used models in this category are the Finite-Range Droplet Model (FRDM) [48] and the Weizsäcker-Skyrme (WS) mass model [49]. The microscopic mass model is primarily based on the mean-field theory using energy density functionals. For example, self-consistent Hartree-Fock-Bogoliubov (HFB) calculations with Skyrme and Gogny interactions [50] provide the most accurate results, achieving a root-mean-square (rms) deviation of 0.5 MeV [51]. A more complex microscopic model, which further improves the description of fission barriers and large deformations, provides global descriptions of both ground-state and excited-state properties of nuclei, albeit with a modest accuracy of 1-4 MeV [52]. Phenomenological mass models, on the other hand, use a set of parameters fitted to experimental mass data, achieving an rms deviation of approximately 0.4 MeV [53]. Improving the accuracy of these mass models is crucial for meeting the requirements of astrophysical calculations.

Figure 1 displays the residuals between different mass models and the experimental data in the Atomic Mass Evaluation AME2016 [54]: the macroscopic-microscopic mass models FRDM95 [55], FRDM12 [48], and WS3 [56], microscopic mass models HFB27 [51], UNEDF1 [52], and RMF [57], and the phenomenological models DZ28 [53] and KTUY05 [58]. Most models generally described the mass data well, particularly DZ28, WS3, and HFB27. However, the microscopic mass models RMF and UNEDF1 exhibit prominent systematic patterns in the nuclear chart. For UNEDF1, there appears to be a separation where the mass-excess values of proton-rich nuclei are underestimated, whereas those of neutron-rich nuclei are overestimated. This suggests that some important physics may be missing or not considered in their effective interactions. The models DZ28, FRDM95, FRDM12, and WS3 appear to underpredict the mass values for nuclei with Z>100. The predictive power of all mass models diminishes as they approach the driplines.

Fig. 1Full-screen View

(Color online) Residuals (in MeV) between theoretical mass models and experimental data from AME2016

The residuals between the mass models and experimental data not only quantify the model accuracy but also provide valuable insights into the model deficiencies. A detailed study of these residuals can help improve the mass models. GPs are a powerful tool for capturing systematic patterns and underlying structures in the data. By carefully designing the framework for GPs, the model accuracy can be significantly improved.

Gaussian-Process Regression

The Gaussian-Process regression is a non-parametric Bayesian technique that provides exceptional flexibility in modeling complex functions of arbitrary forms [38]. By employing a kernel function, it can effectively capture the underlying trends and quantify uncertainties, establishing correlations between any two points in the input space. In this work, we model the residuals between theoretical mass predictions () and experimental data () using GP regression. The GP-refined predictions of mass models can be easily recovered by using . The mathematical derivations of Gaussian-Process regression are detailed in the Supplementary Material (SM) [59]. Our GP model employs two kernel functions to describe spatial dependencies between nuclei: the Matérn (MAT32) and the RBF kernel (see Eqs. 4-5 in SM for definitions [59]). The input space includes proton and neutron numbers (Z, N), the parity number of proton and neutron , the distance to the nearest proton magic number and neutron magic number , where , 8, 20, 28, 50, 82 and , 8, 20, 28, 50, 82, 126, the neutron-excess number , and the nuclear collectivity indicator [60]. The GP model with the normal two-dimensional input space (Z, N) serves as a baseline in our work and is denoted as GP-2D, and the model with the full input space () is denoted as GP-8D.

The training dataset consists of 2272 nuclei from AME2016 with uncertainties keV and . To validate the extrapolation of our GP model, we tested the model on new masses from AME2020 [61] that were not included in AME2016. This non-random partitioning ensured a rigorous evaluation of the predictive capability. Our framework was implemented using the GPy software developed by the Sheffield Machine Learning Group [62] and applied to refine the residuals of the mass models.

Results and Validation

In this study, the mean and variance of our estimator δ^GP are calculated analytically from the posterior distribution under the assumption of a Gaussian-distributed likelihood function. The approach yields results that are approximately equal to those obtained via MCMC sampling but with significantly reduced computation costs [39].

Figure 2 displays the rms deviations of the mass models for the training dataset. The GP-2D regression with only the two-dimensional input space (Z, N) markedly enhanced the accuracy of all mass models. σ^GP is consistently lower for the RBF kernel than that for the MAT32 kernel in most cases, indicating that it might produce more accurate mass predictions. The greatest enhancements were observed in the two microscopic mass models, RMF and UNEDF1, where the rms deviations dropped from the original values of approximately 2 MeV to approximately 0.1 MeV. For the other models, the accuracy improved by factors of 2–6 compared to their original values.

Fig. 2Full-screen View

(Color online) Root-mean-square (rms) deviations (in MeV) of mass models before and after Gaussian-Process refinement using (Z,N) input (GP-2D). The results for both the RBF and MAT32 kernels are shown, with a gray dashed line indicating 100 keV accuracy

In Bayesian analysis, the credibility interval (CI) is the interval within which an estimate is expected to fall with a particular probability. As pointed out in Ref. [39], the rms deviation is not the sole quantity used to address the performance of a mass model. Therefore, we evaluate two empirical coverage probabilities for CIs for our estimator δ^GP, which correspond to 68%(1CI) and 95%(2CI) posterior probabilities of lying for the MAT32 and RBF kernels used in this study. The results are summarized in Table 1.

The credibility intervals for both the MAT32 and RBF kernels fell below their nominal values, suggesting a potential underestimation of the uncertainties in these GP regressions. Ideally, we expect a model whose credibility interval matches the expected values, although this is not always the case in practice. The 1CI and 2CI empirical coverage probabilities show that the MAT32 kernel often performs closer to the expected 68% and 95% ranges, which may indicate better uncertainty quantification in some instances. Furthermore, the RBF kernel might oversmooth the mass surface, while the Matérn kernel can capture sharp jumps in the residuals and avoid the “concentration of measure" problem in high-dimensional spaces [63]. Regarding error estimation, the GP with the MAT32 kernel performed better than the RBF kernel. Given its more conservative uncertainty estimates and better empirical coverage, we selected the MAT32 kernel for subsequent GP model refinement.

Figure 3 illustrates the performance of our MAT32 kernel regression for the Sn isotopic chain. After GP-2D refinement, all mass models, regardless of their theoretical foundations, showed excellent agreement with the experimental data. For example, the two bulges from the HFB27 model are flattened, and the downsloping trends of the RMF and UNEDF1 predictions towards the neutron-deficient side of Sn are significantly corrected. This demonstrates that the GP can capture the systematic patterns inherent in each model’s residuals, ultimately improving the model’s accuracy (see Fig. 1 in SM for a full residual map [59]). More importantly, the inclusion of statistical uncertainties with each prediction provides valuable insight into the model’s reliability, which is particularly crucial for estimating the masses of exotic nuclei far from stability, which are key to astrophysical studies.

Fig. 3Full-screen View

(Color online) Residuals for the Sn isotopic chain before (lines) and after (symbols with error bars) the application of the MAT32 kernel in the (Z,N) input space (GP-2D) in the training dataset

Fine-tuning and Predictive power

While the GP-2D model has already demonstrated significant improvements, further enhancement of its predictive capability requires the incorporation of additional nuclear structure features. Thus, we conducted GP regression by integrating the expanded set of physics-informed inputs, as mentioned above. As shown in Table 2, the full 8D model, GP-8D, outperforms the baseline model GP-2D. The rms deviations in GP-8D for RMF, FRDM95, HFB27, and WS3 decreased by 72.8%, 50.3%, 64.1%, and 63.3%, respectively. These enhancements demonstrate that our GP framework further captures subtle nuclear interactions that are often missed by conventional mass models. The accuracy achieved in this study is comparable to the mass predictions using the 12 Garvey-Kelson relations [64] and the mass relations of mirror nuclei [65]. One of us (H.Z.L.) used BNN to learn the mass surface of even-even nuclei and the correlation energies, reaching a similar accuracy of 84 keV [35]. The predictions for the masses of the Sn isotopic chain in the full input space are displayed in Fig. 2 in SM [59].

The exceptional accuracy of our refined GP models highlights the importance of examining the length-scale parameters in the GP framework (the definition of length scales is seen in Eqs. 4–5 in SM [59]). These parameters govern the decay of correlations between nuclei across the (Z, N) plane, directly affecting the model’s ability to extrapolate to regions far from the experimental data. As shown in Table 2 (last two columns), the length scales for Z and N exhibited significant changes when transitioning from GP-2D to GP-8D. For models that utilize only (Z, N), these parameters typically fall within the range of 2 to 14. However, upon the implementation of more physical inputs, these length-scale values increase significantly. For example, the length scales for Z and N for FRDM12 increase from 4 and 5 to 29 and 24, respectively. In GP regression, larger length scales indicate a slower decay in the correlations between data points as their spatial separation increases. This behavior aligns with the nature of the GP model, reflecting stronger interrelations among the expanded set of input features. We favor larger length-scale values in the GP model because of their robustness against outliers. A shorter length scale could overly rely on neighboring nuclei for predictions, potentially compromising accuracy if the last few nuclei are measured inaccurately. Although it is not the scope of this study to discuss the impact of the length scale in detail, its critical role in model performance merits further investigation.

To evaluate the predictive power of our refined GP models, we tested them against the 74 new masses adopted in the AME2020 (see Fig. 4a). The residuals of each mass model before and after the GP-8D refinements are compared in Fig. 4(c-j). The original models exhibited large deviations, whereas the GP-8D refinement substantially reduced the residuals across all cases. Specifically, the refined models demonstrated a closer alignment with the experimental data, as evidenced by the narrowing of the residual distributions and reduction in the rms deviations. If we constrain the test set for nuclei with uncertainties smaller than 100 keV, the rms deviations are further reduced. For example, the rms deviations for DZ28 and WS3 were reduced to 149 keV, crossing the 150 keV level. It should be noted that the obtained 150 keV rms deviation represents only short-range (local) extrapolation capability, which cannot be extrapolated to the r-process region without additional validation. We quantify the ensemble performance of our GP models by calculating the mean and standard deviation across all GP-refined models, which are defined as and (m=8 and n=74), respectively. As shown in Fig. 4b, the mean predictions match the accuracy of the best-performing models (DZ28 and WS3).

Fig. 4Full-screen View

(Color online) (a) Nuclear chart showing the new masses from AME2020 (colored squares) and AME2016 training data (gray). The color shows the residuals between the simple average of the GP-8D estimates and the new mass values (see text for details). (b) Histogram of residuals between the model average and new masses in AME2020. (c-j) Histograms of the residuals for each mass model. The numbers in the histograms show the calculations of the standard deviations of the original mass models(σ_th), GP-8D refined models(σ_GP), and GP-8D refined models compared to masses with uncertainties smaller than 0.1 MeV ()

To rigorously assess the extrapolation capability of our GP ensemble approach, we performed a theory-to-theory validation using FRDM12 as “pseudo-experimental” data (see SM for details [59]). This validation revealed three characteristic regimes in the trends of the rms deviation as a function of the extrapolation distance r (c.f. Figure 3 in SM [59]). In the short-range regime (r<5), the rms deviations of the averaged GP models increase with r. This behavior reflects the expected uncertainty growth in the GP predictions. In the intermediate regime , the rms deviation first decreases and then reaches a plateau at approximately 200 keV. This behavior indicates that the GP can better capture systematic trends even far from stability. For long-range extrapolation (r>13), the rms deviation starts to increase again and appears to reach a plateau at approximately 300 keV. The theory-to-theory validation provides quantitative confidence intervals for GP predictions in the unknown regions of the nuclear chart.

Figure 5 compares the original and average GP-refined mass models for the Sc and Pr isotopic chains. The new mass values for ^54-60Sc were taken from Refs. [66-69] and for Pr taken from Ref. [70]. Our GP estimates show good agreement with experimental data and the new mass data (within 2σ). The emergence of systematic deviations for neutron-rich ^56-60Sc should be interpreted with caution. These masses come from Bρ-ToF measurements [68, 69] and have larger uncertainties (>100 keV). In AME2020, five of the 15 results in Ref. [69] were excluded. Specifically, AME2020 suggests ⁶⁰Sc should be 1 MeV more bound. For ^54,55Sc, they were measured using high-precision MR-ToF [66, 67] and our GP predictions agree with their results. The deviations for neutron-rich Sc isotopes may reflect the limitations of the Bρ-ToF technique rather than the GP’s extrapolation deficiency. High-precision mass measurements of Sc isotopes are required for more critical validation.

Fig. 5Full-screen View

(Color online) Mass differences between different mass models relative to FRDM95 for Sc (a) and Pr (b) isotopic chains. The black solid points correspond to the experimental data in AME2016, while the red ones are new masses taken from [66-70]. The dark gray (blue) and light gray (blue) bands show the and standard deviations of the simple average of GP-8D (original mass models), respectively

We also note in Fig. 5 that the averaged uncertainty bands for Pr and Sc isotopes do not grow at the same pace—the uncertainties of Pr grow more slowly than those of Sc isotopes. This is mainly due to the fact that the mass data around neutron-rich Pr isotopes are more constrained (with more higher precision experimental data) than that of neutron-rich Sc isotopes. The Pr isotopes are located in the heavier region, where more nuclides are available than those of Sc. Because the GP is completely driven by the data, we obtain a more reliable extrapolation if the surrounding training dataset is more accurate. Another possibility is the way in which each model describes the nuclear structure effects in its respective region. While the uncertainty of the averaged GP-refined models increases steadily as they move away from the known region, this represents a limitation of the current experimental constraints rather than inherent model deficiencies.

Conclusion

The present study demonstrates that Gaussian Processes significantly improve the accuracy of nuclear mass models. By capturing the inherent and systematic biases rooted in each model, GPs enhance the reliability of mass predictions, particularly for nuclei far from stability, which are critical for astrophysical modelings. Through the incorporation of physics-informed features, our refined models achieved a sub-100 keV accuracy for nuclei with , representing a substantial improvement over the traditional mass models. By comparing with high-precision ( keV) mass data from AME2020, we confirmed that our GP-refined mass models achieved an rms deviation below 150 keV for short-range extrapolation, with the Matérn kernel offering particularly robust uncertainty quantification compared to the RBF kernel. While the GP-refined models show robust performance near known regions, their reliability for r-process regions requires further validation. Through theory-to-theory validation, we noticed that the rms deviations first increased as one moved away from the known region. However, the flat area in the intermediate region () indicates that our GP model can better capture systematic patterns even at larger extrapolation distances. A comprehensive set of mass predictions (Table 1 in SM [59]) for unknown masses in AME2020 serves as a valuable benchmark for future mass measurements and astrophysical calculations. Given the inherent capability of GP for interpolation, we can also use our models to identify outliers in the existing mass data. Future studies should explore the impact of length-scale parameters and validate models using high-precision data.