{'fraction_non_alphanumeric': 0.03616160318581797, 'fraction_numerical': 0.024439591495921384, 'mean_word_length': 4.646237533998187, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Beyond the proton drip line: Bayesian analysis of proton-emitting nuclei

Léo Neufcourt, 2 Yuchen Cao (曹宇晨), 3 Samuel Giuliani, 3 Witold Nazarewicz, 4 Erik Olsen, and Oleg B. Tarasov Department of Statistics and Probability, Michigan State University, East Lansing, Michigan 48824, USA Facility for Rare Isotope Beams, Michigan State University, East Lansing, Michigan 48824, USA National Superconducting Cyclotron Laboratory, Michigan State University, East Lansing, Michigan 48824, USA Department of Physics and Astronomy, Michigan State University, East Lansing, Michigan 48824, USA Institut d’Astronomie et d’Astrophysique, Université Libre de Bruxelles, 1050 Brussels, Belgium (Dated: October 29, 2019)

Much of the nuclear landscape of bound nuclides remains unexplored [1,2]. The one-and two-proton drip lines lie relatively close to the line of beta stability due to the presence of the Coulomb barrier that has a confining effect on the proton density. As a result, relatively long-lived, proton-unstable nuclei can exist beyond the drip line [3][4][5][6][7]. The vast territory of proton-unstable nuclides contains rich and unique information on nuclear structure and dynamics in the presence of the low-lying proton continuum.
The positions of particle decay thresholds are determined by the nuclear binding energy through measured masses. In the regions where experimental mass data are absent, nuclear models must be deployed to provide the missing information about the topography of the mass surface. In this context, the quality of theoretical mass predictions can be significantly improved when aided by the current experimental information through machine learning techniques [29][30][31][32][33][34][35][36][37][38][39]. Recently, we developed the statistical framework of Bayesian Gaussian process techniques to quantify patterns of systematic deviations between theory and experiment by providing statistical corrections to average prediction values, and develop full uncertainty quantification on predictions through credibility intervals [40]. The quantified predictions of individual models enabled us to carry out Bayesian model averaging (BMA) analysis of nuclear masses [41]. In this way, the "collective wisdom" of several relevant models could be maximized by providing the best prediction rooted in the most current experimental information.
In this paper, using several variants of BMA, we quantify the predicted binding-energy surface in the region of 2p radioactivity. To this end, we employ several global mass models to determine the posterior probability for each proton-rich nucleus to exhibit ground-state proton or 2p emission. We find that extrapolations for proton drip-line locations are fairly consistent across the mass models used, in spite of significant variations between their raw predictions. In this respect, this study should be considered as an extension of the previous work [42,43], in which the most promising candidates arXiv:1910.12624v1 [nucl-th] 28 Oct 2019 for 2p emitters were identified by considering several mass models. Here, we limit our investigations to nuclei with Z ≤ 82 as it is predicted [42,43] that above lead the α-decay mode dominates and no measurable candidates for 2p emission are expected.
The paper is organized as follows. Section II presents the global nuclear mass models used in our study. The statistical methodology is described in Sec. III. The results obtained in this study are discussed in Sec. IV. Finally, Sec. V contains a summary and conclusions.

II. NUCLEAR MASS MODELS
Our global mass calculations are based on nuclear density functional theory (DFT) with several energy density functionals (EDFs). As in the previous studies [40,41], we considered the Skyrme functionals SkM * [44], SkP [45], SLy4 [46], SV-min [47], UNEDF0 [48], UNEDF1 [49], and UNEDF2 [50]. In this work, we have enriched the set of nuclear models with two additional EDFs: D1M and BCPM. The functional D1M [51] is a modern parametrization of the finite-range Gogny interaction, optimized to 2149 measured masses from the 2003 mass evaluation (AME2003) [52], charge radii, and nuclear matter properties. In the functional BCPM [53], the bulk part of the functional is given by a fit to the microscopic EOS in both neutron and symmetric nuclear matter. This formulation of the functional results in a relatively small number of free parameters that are adjusted to reproduce the experimental binding energies of 579 even-even nuclei of AME2003.
For each Skyrme EDF, the mass table of even-even nuclei was computed self-consistently by solving the Hartree-Fock-Bogoliubov (HFB) equations. 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. In this respect, this work follows closely the methodology described in Refs. [1,40,41,43]. For D1M and BCPM, the binding energies of odd-A and odd-odd nuclei are computed by solving the HFB equations for one-and two-quasiparticle configurations with the appropriate constraint on particle number [54]. The above set of DFT models was augmented by two mass models commonly used in nuclear astrophysics studies: FRDM-2012 [55] and HFB-24 [56].
It is to be noted that while the proton chemical potential is positive for proton unbound nuclei, the HFB results obtained with the discretized continuum are very stable in the considered range of binding energies. This is because the Coulomb barrier tends to confine the proton density in the nuclear interior and effectively pushes the continuum up in energy [57,58] on the proton-rich side. Consequently, the proton separation energies S 1p and S 2p and the corresponding Q-values can be obtained safely from the predicted binding energies.
The candidates for the true 2p decay were selected according to the energy criterion used in the global survey [42]: This condition corresponds to the simultaneous emission of two protons as the sequential emission of two protons is energetically forbidden.

III. STATISTICAL METHODS
Our Bayesian methodology for building Gaussian process (GP) emulators to produce quantified extrapolations of theoretical nuclear model predictions beyond the experimental data range has been extensively developed in our previous work [40,41]. Here, we incorporate two statistical innovations, a non-zero GP mean parameter and a new Bayesian calculation of model mixing weights.

A. Gaussian process
The Bayesian statistical model for the separationenergy residuals (i.e., differences between experimental and theoretical values) y i = y exp (x i ) − y th (x i ) can be written as: where the function f (x, θ) represents the systematic deviations and σ the statistical uncertainty propagated from the uncertainty on the statistical model parameters.
Quantified extrapolations y * are obtained from the posterior distribution p(y * |y) using a stationary Markov chain. Similarly to our previous studies, we model independently S 1p and S 2p on the subsets of nuclei defined by particle-number parities (even-even, even-odd, etc.). By doing this we are ignoring some (slight) correlations between systematic uncertainties.
For the function f we take a Gaussian process, i.e., a Gaussian functional on the two-dimensional nuclear domain indexed by x = (Z, N ), characterized by its mean µ and covariance k: We have found in a previous study [40] that Gaussian processes overall outperform Bayesian neural networks, achieving similar root-mean-squared (rms) deviations with a more faithful uncertainty quantification and considerably fewer parameters. In our previous studies, we took for simplicity the GP mean to be uniformly zero.
Here we take it as an additional scalar parameter. The results below show that this can improve the rms deviation by an additional 15%, compared to the initial 25% refinement brought by the GP. In order to model the "spatial" dependence of nearby nuclei, we use an exponential quadratic covariance kernel: where the parameters θ ≡ {η, ρ Z , ρ N } have a straightforward interpretation: η defines the scale and ρ Z and ρ N are characteristic correlation ranges in the proton and neutron directions, respectively. Consequently, our statistical model has four parameters θ := (µ, η, ρ Z , ρ N ). Samples from posterior distributions were obtained from 50,000 iterations of Monte Carlo Markov Chains (MCMC), after the stationary state was reached (50,000 samples are previously burnt-in), which were used in turn to generate 10,000 mass tables. Each series of simulations required 16,128 core hours distributed on 96 machines on computing clusters.

B. Datasets
Our dataset combines all experimental masses from AME2003 [52] and AME2016+, which contains the AME2016 dataset [59] augmented by masses from Refs. [60][61][62][63][64][65][66]. For nuclei where experiments have been repeated, we keep the most recent value. For testing purposes we split this dataset into a training set (AME2003) and a testing set (AME16-03: all masses in AME2016+ that are not in AME2003). For prediction purposes, we use the full dataset AME2016+ for training, and carry out calculations based on a large set of nuclei for which raw theoretical separation energies are not too negative; this includes all proton-bound nuclei. Nuclei with negative experimental separation energies have not been used for training.

C. Bayesian model averaging
When combining several models M k , the classical literature uses Bayesian posterior weights conditional to the data y given by [67][68][69] where π(M k ) are the prior weights and p(y|M k ) are evidence (integrals) obtained by integrating the likelihood over the parameter space. In our study, we shall use uniform priors. First, we carry out model mixing calculations with the prior average of the models, i.e., with uniform weights. In the absence of additional information and costly posterior computations, the choice of uniform weights is essentially optimal [69]. This variant is denoted as BMA-0.
Similarly as in Ref. [41], in the context of extrapolations, we want to select the models with the best predictive power, and avoid overfitting. To this end, we evaluate the evidence integrals on new independent data that are not included in the training of any individual statistical model. (In contrast, stricto sensu BMA would compute the evidence based on the whole same dataset as used in the training of each model.) In a first variant of BMA calculations (denoted BMA-I), in the spirit of [41] we consider simplified Bayesian weights where the evidence is replaced by the posterior probability that each model accounts correctly for the signs of the experimental Q 2p and S 1p values of the five 2p emitters x 2p,known according to (1). That is, in this variant, the weights are computed based on the ability of model M k to predict the set x 2p,known as true 2p emitters. Here, the resulting weights are: In the second variant of BMA calculations (denoted BMA-II), the evidence p(y * |y, M k ) is defined based on the ability of model M k to predict known Q 2p values of five experimentally known 2p emitters (this set of nuclei is called x 2p,known ≡{ 19 Mg, 45 Fe, 48 Ni, 54 Zn, 67 Kr} in the following). The resulting weights are: Finally, we also consider BMA-III, a trivial version of BMA-II, consisting in the Gaussian likelihood of x 2p,known evaluated at the posterior mean and posterior variance values, assuming that these quantities are independent. Considering the Q 2p -residuals y i of the five known 2p emitters, the corresponding weights are: Compared to BMA-II, this approximation neglects all correlation effects as well as Gaussianity defects of the posterior predictions at the five locations that we meticulously added to the Gaussian process.

D. Model weights computation
The evidence integrals p(y * |y, M k ) are obtained by "recycling" Monte-Carlo samples using the estimate where θ (i) k is the i th parameter sample of M k . This is justified by the formula of total probability, p(y * |y, M k ) = p(y * |y, M k , θ k )dP (θ k |M k , y). (10) For each model, this calculation can be efficiently performed in two steps. In the first step, we compute q In the second step, the weights are obtained as The testing dataset used to compute the weights overlaps the sets of even-N and odd-N nuclei. In our final analysis, we have assumed that the statistical models are independent of N -parity, meaning that these datasets can be divided as y * := (y * e , y * o ). From the independence of y * e |y e , M k , θ k and y * o |y e , M k , θ k for each model it follows: (11) and the final result is a sum i . Note that the assumption of independence is an important caveat, and that this calculation would be enhanced by better estimates of the correlations between the systematic errors of S 1p and S 2p .

IV. RESULTS
A. Model performance and assessment Figure 1 shows the comparison between theoretical predictions and experimental data for 19 Mg, 45 Fe, 48 Ni, 54 Zn, and 67 Kr. We find an overall good agreement with experimental values, up to error bars, for most statistically corrected models. It is seen that individual models, e.g., SkM * , HFB-24, and D1M, are behaving singularly for some nuclei; this shall impact the model weights used in the following BMA analysis. The worst performer is the traditional SkM * model, which practically misses all experimental data points. The best performers are the UNEDF1, UNEDF2, and SLy4 models, which provide the lowest rms deviation from experiment. Table I illustrates the global performance of individual models for the measured proton separation energies contained in the AME16-03 testing dataset. It is seen that the statistical-model correction brings a significant improvement to the predictions. Namely, in the GP (µ = 0) variant, the rms deviation is reduced by about 25%, similarly to our previous studies [40,41], while the µ = 0 calculation brings another 15% reduction. The improvement is most significant on the S 2p values, and for the models with the largest raw (i.e., statistically-uncorrected) rms deviations. The improvement for FRDM-2012 and HFB-24, optimized to the experimental mass table, is minor, and both variants of GP calculations yield practically identical results. As in the previous work [40], the rms deviations are quite similar across models following the statistical treatment, which suggests that most of the systematic trends have been captured by our statistical models.

B. Gaussian process parameters
To gain some insight into our statistical model, Fig. 2 shows the univariate and bivariate distributions of the posterior samples of µ, η, ρ Z , and ρ N computed for SLy4 and UNEDF1. Kernel Density Estimates (KDEs) are classically obtained using a Gaussian kernel. Bivariate Highest Density Regions (HDR)s are taken as level lines of the continuous KDE [70].
We can see that the parameters are both well constrained and uncorrelated. The only exception is the pair (η, ρ N ) for which the correlation coefficient reaches 0.59 for SLy4, which indicates a fairly low correlation. The dependence of GP parameters on the nuclear model used is also weak, with UNEDF1 producing slightly more localized KDEs.
The posterior mean and standard deviation of the GP parameters for the two-proton separation energies of even-even nuclei are listed in Table II. We observe an overall stability of the parameters ρ Z and ρ N , with correlation effects occurring within the range of ±(2 − 3) particle numbers. Symmetrically to what we have earlier noticed for neutron separation energies [40,41], we see that the correlations effects are here substantially smaller in the proton direction than the neutron direction, consistently with variations in the proton separation energies, which are stronger overall in the Z direction than along N .
There are more disparities on the parameters µ and η, which are directly related to the scale of the statistical corrections. The parameter µ is maximal at more than 1 MeV for the model SkM * and within one standard deviation of zero for the more phenomenological models HFB-24 and FRDM-2012. The values of η are also significantly higher for SkM * and lower for the two phenomenological models. This confirms the commonsense expectation that there is not much left for a GP to capture when the nuclear model has already exploited enough of the data structure.
The scale of the mean and scale parameters µ and η appear consistent with the rms improvement in the second row of Table I. As stated above, their largest values were obtained for SkM * , for which the rms improvement is also the best. In this case only 7 masses of spherical nuclei (in addition to other observables) were used in its optimization process. Interestingly, the performance of SkP and SLy4 is outstanding, considering their limited mass input datasets. In contrast these parameter values are particularly low for HFB-24 and FRDM-2012, which were optimized to very large sets of masses. This confirms that the statistical correction has a balanced autoscaling, and is itself robust to overfitting.
A slightly different argument applies to D1M and BCPM, for which the values of µ and η are relatively large in spite of these functionals being optimized to a large set of nuclei. It is to be noted, however, that the rotational corrections were added atop the HFB binding energies during the optimization process of D1M and BCPM. Since such corrections result in unphysical staggering of separation energies in transitional nuclei [71], they have been neglected in this work. This resulted in increased values of the mean and scale GP parameters. The error bars produced for HFB-24 and FRDM-2012, which are scaled by η 2 , must be taken with caution, and are certainly undervalued in the propagation of the training error. Indeed, the rms deviations are among the highest on the test data (including points more recent than the models) for HFB-24, both before and after statistical correction.

C. Model mixing performance and assessment
The model weights obtained in the BMA variants used are listed in Table III. The weights of BMA-I are rather uniformly distributed, with the highest values for UN-EDF1, UNEDF2, D1M, FRDM-2012, and HFB-24; the weight 0 for SkM * is due to the fact that it misses completely the sign of the second and third testing points (see Fig. 1), which is eliminatory. In the case of the likelihood-based BMA-II, only the Skyrme models UN-EDF1 and UNEDF2 practically contribute. This can be explained by the concavity of the evidence with respect to the data, which heavily penalizes large deviations at single locations. Figure 1 shows the results of BMA for the five true 2p emitters considered. It is encouraging to see that the BMA predictions agree well with experimental values, i.e., their standard error bars overlap. As expected theoretically [72], the BMA results should achieve the lowest rms deviations among all models. This is confirmed in Table I, which shows that the rms deviations for S 1p and S 2p obtained in all BMA calculations are indeed below those of individual nuclear mass models. We see in Table I that all BMA variants perform very similarly overall: the rms deviations are around 0.38 keV for S 1p and 0.36 MeV for S 2p . It may be thus tempting to associate these values with the precision limit that current EDFs can achieve in the description of proton separation energies. Interestingly, the simple models BMA-0 and BMA-III yield rms deviations comparable, if not lower, to the other elaborate variants. This is not entirely surprising: the posterior weights favor by construction the statistically corrected models best fitted at locations y * . Indeed, some of the nuclear models employed, such as HFB-24 and FRDM-2012 are less robust to out-of-sample testing data, which can possibly be linked to some overfitting. Figure 1 shows that the posterior means of BMA-III are -by construction -close to the experimental Q 2p values. The disadvantage of this approximation to BMA-II is that it neglects all correlation effects and long range dependencies meticulously added to the GP. Consequently, the differences between the BMA-II and BMA-III weights can be large for some models, cf. the weights for SLy4 and UNEDF1 in Table III. This result can be traced back to non-zero covariances between the five locations and Gaussianity defaults. First, the (posterior) covariances between the different Q 2p values are non-zero in the actual samples, but assumed uniformly zero in BMA-III. Second, the evidence calculation in BMA-II includes an integration over the parameter space, which is absent in BMA-III. For instance, the parameter ρ N is 27% larger for SLy4 than in UNEDF1, and this means that one more neighbor in the neutron direction is roughly included in the GP calculations based on SLy4. While the parameter difference is not dramatic, it does significantly impact both the predictions and model weight estimates.

D. Predictions of BMA calculations: quantified landscape of proton-rich nuclei
The quantified nuclear binding-energy landscape for proton rich nuclei, predicted in BMA-I and BMA-II, is displayed in Fig. 3. To facilitate the presentation, the information for each isotope is given relative to the neutron number N 0 of the lightest proton-bound isotope for which an experimental one-or two-proton separation energy value is available. In analogy with Ref. [41] we show the probability p ex that a given isotope is proton-bound, i.e., that S 2p > 0 for even-Z nuclei and S 1p > 0 for odd-Z nuclei. Formally, in the Bayesian language, this quantity can be defined as: The proton drip line corresponds to p ex = 0.5. The reference values of N 0 (Z) are listed in Tables IV (for even-Z   Tables IV and V, is the neutron number of the lightest proton-bound isotope for which an experimental one-or two-proton separation energy value is available. The domain of nuclei that have been experimentally observed is marked by stars; those within FRIB's experimental reach are marked by dots. See text for details. nuclei) and V (for odd-Z nuclei), together with the range of observed nuclei and proton drip-line nuclei. Figure 3 and Tables IV and V also show the isotopes that will be accessible at the future Facility for Rare Isotope Beams (FRIB) [73,74]. FRIB will accelerate ion species up to 238 U with energies of no less than 200 MeV/u (beam power up to 400 kW). The FRIB production rates have been calculated with the LISE ++ code [75] designed to predict intensity and purity for future experiments using rare beams with in-flight separators. The EPAX2.15 cross-section systematics [76] and the LISE ++ 3EER Abrasion-Fission model [77,78] were used to calculate production cross-sections for projectile fragmentation and fission reactions correspondingly. The multi-step reactions in thick targets were taken into ac-count. In this process, the projectile undergoes a series of successive reactions until the fragment of interest is produced. FRIB Rates and details of their calculations are available online [79]. In our estimates, we assumed the experimental limit for the confirmation of existence of an isotope to be 1 event/2.5 days.
In general, the drip-line predictions of BMA-I and BMA-II are very similar; some differences can be see for the elements just below Pb, with BMA-I calculating slightly higher values of p ex . In this region of nuclei, FRIB is expected to have a particularly high discovery potential: by exploring the vast territory of protonunstable isotopes, it will extend the domain of known nuclei considerably.
The territory of the true 2p emitters predicted in our BMA calculations is shown in Fig. 4. Here, the quantity of interest is the posterior probability p 2p that S 2p < 0 and S 1p > 0: see Eq. (1). Again, BMA-I and BMA-II predictions are close. The isotopes that are potential candidates for 2p radioactivity lie in a band corresponding to p 2p ≥ 0.5. Most of those isotopes are within the range of FRIB.

E. Two-proton radioactivity: lifetime considerations
In the previous sections, the discussion of proton radioactivity was solely based on energy arguments. However, the energy relations do not tell the full story. Indeed, the proton decays corresponding to very large Q 1p/2p -values are going to be too fast to be observed. If the Q 1p/2p -values are very low, the proton-decay rate is going to be negligible compared to other decay modes, such as β or α decays. When it comes to 2p decay, the practical range of lifetimes is [42]: The lower bound of 100 ns corresponds to the typical sensitivity limit of in-flight, projectile-fragmentation techniques. The upper bound of 100 ms ensures that the 2p decay will not be dominated by β decay.
To get an order-of-magnitude estimate of 2p lifetimes, we used the simple semiclassical Wentzel-Kramers-Brillouin (WKB) approximation, and assumed a diproton decay with = 0. For details of the WKB calculations, see [42,80]. The value of the proton overlap O 2 has been fitted to match the measured lifetimes of 19 Mg, 45 Fe, 48 Ni, 54 Zn, yielding O 2 = 0.0011. Our WKB approach agrees very well with the semiclassical effective liquid drop model analysis of Ref. [81]. Figure 5 shows the Q 2p values predicted in BMA-I together with the lifetime range (14). It is important to note that the uncertainties on the predicted values of Q 2p usually correspond to several decades of the 2p-decay lifetime. As seen in Fig. 5, the known 2p emitters 45 Fe, 48 Ni, 54 Zn, and 67 Kr consistently fall within the lifetime range (14). The most promising other candidates for the true 2p radioactivity are: 30 Ar, 34 Ca, 39 Ti, 42 Cr, 58 Ge, 62 Se, 66 Kr, 70 Sr, 74 Zr, 78 Mo, 82 Ru, 86 Pd, 90 Cd, and 103 Te.
For nuclei with Z ≥ 54 that are within the lifetime  It is also to be noted that many of the extremely proton-rich nuclei in Fig. 5 with small p 2p values, such as 131,132 Dy, 134,135 Er and 144,145 Hf, are excellent candidates for the sequential emission of two protons (pp) [43].
Our BMA findings are fairly consistent with predictions of other papers. The nuclei 39 Ti (p 2p = 0.74) and 42 Cr (p 2p = 0.60) are expected to be excellent 2p-decay candidates according to the phenomenological analysis based on the modified Kelson-Garvey mass relations [82] and shell model analysis [83]. As discussed in [13], 39 Ti primarily decays by beta disintegration. This is not inconsistent with the low Q 2p value predicted in BMA-I. Other 2p-decay candidates predicted by BMA-I discussed in the literature include: 26 [84]. The nucleus 103 Te, which has been predicted [43] to exhibit a competition between alpha decay and 2p radioactivity, is expected to have T 2p > 0.1 s. In 145 Hf, alpha decay is predicted [43] to compete with sequential pp emission.

V. CONCLUSIONS
In this paper, we employed the Bayesian Model Averaging framework to quantify the proton-stability of the nucleus. To this end, we introduced the probability p ex , which is the Bayesian posterior probability that the oneor two-proton separation energy of a nucleus is positive. We also evaluated the posterior probability p 2p that a nucleus is a true 2p emitter.
We demonstrate that the statistical-model correction improves predictions significantly. Overall, for the testing dataset AME16-03, the rms deviation from experimental S 2p values is in the 400-600 keV range in the GP (µ = 0) variant across theoretical models employed in this study. This is consistent with the previous analysis of S 2n values [40] and indicates that our GP model captures a significant part of the systematics and brings a sound refinement to the nuclear theory models.
Following the model averaging, the rms deviations from experiment for proton separation energies are surprisingly similar for all the BMA variants used: they are around 360 keV for S 2p and 380 keV for S 1p . This result suggests that the further rms reduction cannot probably be expected without dramatic improvements of fidelity of nuclear mass models.  NFRIB  17  Cl  14  11  14  14  19  K  16  16  16  16  21  Sc  19  18  19  18  23  V  20  20  20  19  25  Mn  22  21  22  21  27  Co  24  23  24  23  29  Cu  27  26  27  25  31  Ga  30  29  29  28  33 As 33 31 In general, our results are fairly consistent with the current experimental data on the proton drip line position and the appearance of 2p radioactivity. Our calculations suggest that no true 2p emission is expected in the lifetime range (14) above Z = 54.
In contrast to an increasing number of studies applying directly statistical models to nuclear physics observables, the sound nuclear physics model underneath our statistical emulators, trained on the residuals, guarantees that our predictions are globally consistent with known physics. The GP parameters are well constrained with relatively weak posterior variances and correlations, and the additional mean parameter brings another reduction of rms deviation. The GP could potentially be refined in future studies by including additional degrees of freedom to describe different regions of the nuclear domain.
We proposed Bayesian average models obtained by using several variants of the BMA weights; the weights are calibrated on independent test data, and thus directly related to the extrapolative power of each model. While we observe significant variations in the weights obtained for different BMA variants, all average models achieve a lower rms deviation than individual model constituents. This validates empirically the essence of the recent BMA analysis [72] where it has been established that the BMA estimator achieves the lowest posterior variance among all models and all model combinations. This result suggests that even the simplest uniform model mixing carried out in several previous studies involving different quantities [1,2,41,85,86] can provide a very valuable information.
It can appear statistically disappointing that the BMA-0, BMA-I, and BMA-III variants achieve the better testing rms performance as compared to the most sophisticated BMA-II method. Surely the comparison of the BMA weights (BMA-II) with their counterparts obtained from model conditional likelihoods taken at posterior mean values (BMA-III) highlights the Gaussianity defects and correlations, hidden in a standard analysis of rms deviations and error bands. The philosophy of using model mixing to attenuate the individual defects of individual models makes it desirable to include a greater diversity of models, and our results suggest that difficulties are to be expected when none of the proposed statistical models gives an accurate description of the covariance structure between data points.
The extrapolation outcomes discussed in this study will be tested by experimental data from rare-isotope facilities. As illuminated by our Bayesian analysis, experimental discoveries of new proton-rich nuclides will be crucial for delineating the detailed behavior of the nuclear mass surface in the vast unexplored region of the nuclear landscape beyond the proton drip lines.