Bayesian model averaging for nuclear symmetry energy from effective proton-neutron chemical potential difference of neutron-rich nuclei

The data-driven Bayesian model averaging is a rigorous statistical approach to combining multiple models for a unified prediction. Compared with the individual model, it provides more reliable information, especially for problems involving apparent model dependence. In this work, within both the non-relativistic Skyrme energy density functional and the nonlinear relativistic mean field model, the effective proton-neutron chemical potential difference $\Delta \mu^*_{\rm{pn}}$ of neutron-rich nuclei is found to be strongly sensitive to the symmetry energy $E_{\rm{sym}}(\rho)$ around $2\rho_0/3$, with $\rho_0$ being the nuclear saturation density. Given discrepancies on the $\Delta \mu^*_{\rm{pn}}$-$E_{\rm{sym}}(2\rho_0/3)$ correlations between the two models, we carry out a Bayesian model averaging analysis based on Gaussian process emulators to extract the symmetry energy around $2\rho_0/3$ from the measured $\Delta \mu^*_{\rm{pn}}$ of 5 doubly magic nuclei $^{48}$Ca, $^{68}$Ni, $^{88}$Sr, $^{132}$Sn and $^{208}$Pb. Specifically, the $E_{\mathrm{sym}}(2\rho_0/3)$ is inferred to be $E_{\mathrm{sym}}(2\rho_0/3) = 25.6_{-1.3}^{+1.4}\,\mathrm{MeV}$ at $1\sigma$ confidence level. The obtained constraints on the $E_{\mathrm{sym}}(\rho)$ around $2\rho_0/3$ agree well with microscopic predictions and results from other isovector indicators.

As nuclear symmetry energy E sym (ρ) governs the isospin-dependence of nuclear matter equation of state (EOS), its density dependence has long been a subject of intense study in nuclear physics and astrophysics [1][2][3][4].Thanks to the rapid development in terrestrial nuclear experiments and astrophysical observations, great progress has been achieved recently in exploring the density dependence of the symmetry energy [5][6][7][8].However, various challenges emerge concurrently.For example, very recently, the PREX-2 and CREX collaborations have respectively measured the neutron skin thickness (one of the most important probes of the density dependence of symmetry energy) in 208 Pb and 48 Ca with minimum strong-interaction model dependence [9,10].Unfortunately, the PREX-2 and CREX results exhibit a significant tension, which challenges current theoretical models and our understanding on the density dependence of the symmetry energy [11][12][13][14][15][16].The accurate determination of the density dependence of symmetry energy remains an open challenge.
An important issue in the study of nuclear symmetry energy is the model dependence.Although there are some isovector indicators, like the neutron skin thickness [17][18][19][20] and the electric dipole polarizability [21][22][23][24], whose correlations to the symmetry energy have been proved to be model independent, it is not the case for most of other probes.Analyses with different models on the same observable may result in different and even contradictory conclusions, e.g., the case of charged pion ratio [25][26][27][28][29][30][31][32].Even for the widely used nuclear energy density functional theory (EDF), such as the Skyrme EDF and the nonlinear relativistic mean field (RMF) model, there exist systematic discrepancies among their predictions for the saturation properties of nuclear matter [33][34][35][36][37].The model dependence issue hinders the reliable uncertainty quantification of the density dependence of symmetry energy.
As a powerful statistical inference method, the Bayesian approach has been extensively used to extract physics information from experiments and observations.Within the Bayesian framework, the model dependence issue can be addressed via the Bayesian model averaging(BMA), which mixes predictions of multiple models according to their capabilities of reproducing given data [38][39][40].In nuclear physics studies, the BMA has been successfully used to predict the neutron drip line using various nuclear EDFs [41] and to extract the shear and bulk viscosities of quark-gluon plasma from measurements in heavy-ion collision [42].It has also been proposed to quantify the inter-model uncertainties in some important nuclear physics problems, such as the nuclear matrix elements governing neutrinoless double-beta decay [43].
In this work, within both the non-relativisitic Skyrme EDF and the nonlinear RMF model, we find the effective proton-neutron chemical potential difference ∆µ * pn is strongly correlated with the symmetry energy around a subsaturation density of 2ρ 0 /3, with ρ 0 being the saturation density.However, there exists systematic discrepancy between the correlations predicted by the two models.Based on Gaussian process (GP) emulators trained respectively with a number of Skyrme EDFs and nonlinear RMF parameter sets, the symmetry energy around 2ρ 0 /3 are extracted from Bayesian analyses of the ∆µ * pn in 48 Ca, 68 Ni, 88 Sr, 132 Sn and 208 Pb.Further Bayesian model averaging over the results from the two models arXiv:2312.07031v2[nucl-th] 18 Jan 2024 gives more reliable constraints on the symmetry energy around 2ρ 0 /3 by taking into account the model dependence.
The neutron and proton chemical potential in a nucleus with N neutrons and Z protons can be approximated respectively to be with B(N, Z) being the binding energy of the nucleus.
Although single nucleon separation energy can also be used to approximate the nucleon chemical potential, it suffers from pairing effects, and strong shell effects at magic numbers.Comparing with the single nucleon separation energy, the adopted approximation for µ n(p) is relatively free of pairing and shell effects.
The effective proton-neutron chemical potential difference is then defined by The sensitivity of the ∆µ * pn to the symmetry energy can be expected from the semi-empirical mass formula is the isospin asymmetry, a v , a s , a c and a sym are the coefficients of volume, surface, Coulomb and asymmetry terms, respectively, and the microscopic correction E mic includes pairing and shell effects.Considering 2/A as a small quantity, one can derive ∆µ * pn ≃ −2a c Z A 1/3 + 4a sym I, which suggests that the ∆µ * pn is essentially determined by the symmetry energy of finite nuclei, since the Coulomb coefficient a c ≃ 0.71 MeV [44] is well determined.Given the well-confirmed empirical relationship a sym ≃ E sym (ρ A ) at a subsaturation reference density ρ A in numerous studies [45][46][47], one can expect that there exist strong correlations between ∆µ * pn of neutron-rich doubly magic nuclei and E sym (ρ) at subsaturation densities.We note that the ∆µ * pn is similar to the neutronproton Fermi energy difference proposed in Ref. [47].Nevertheless, the ∆µ * pn can be directly calculated by nuclear energy density functionals and deduced from experimental measured nuclear masses.
We further study the correlation between the effective proton-neutron chemical potential difference and the symmetry energy at subsaturation densities using a number of non-relativistic Skyrme and covariant EDFs.The Skyrme EDF is taken to be the conventional standard form which contains 9 Skyrme parameters of x 0 − x 3 , t 0 − t 3 and α, and a spin-orbit parameter W 0 [48].The 9 Skyrme parameters t 0 − t 3 , x 0 − x 3 and α can be expressed in terms of 9 pseudo-observables in nuclear matter: the nuclear matter saturation density ρ 0 , the energy per particle of symmetric nuclear matter E 0 (ρ 0 ), the incompressibility K 0 , the isoscalar effective mass m * s,0 , the isovector effective mass m * v,0 , the gradient coefficient G S , the symmetry-gradient coefficient G V , and the magnitude E sym (ρ 0 ) and density slope L of the nuclear symmetry energy at ρ 0 [49][50][51].Consequently, the Skyrme EDF has the following 10 parameters: For given θ Sky , the ground-state properties of even-even nuclei are calculated using the Skryme-Hartree-Fock code reported in [52] with the pairing correlation included via the constant-gap approach.
As for the covariant EDF, we take the nonlinear RMF model based on the effective Lagrangian density [54]: where ψ is the nucleon field, A µ is the photon field, and σ, ω µ , ⃗ ρ µ represent the isoscalar-scalar, isoscalar-vector and isovector-vector meson fields, respectively.
the field tensors.⃗ τ is the isospin Pauli matrices, and we take its third component τ 3 = 1(−1) for proton(neutron).In addition, m = 939 MeV is the nucleon bare mass, while m σ , m ω , m ρ are meson masses.The coupling constants of mesons to nucleons are denoted by g σ , g ω and g ρ , while b σ , c σ and c ω describe meson selfcoupling and Λ V is the ρ-ω coupling constant.In this work, the m ω and m ρ are taken to be their experimental values of m ω = 782.5 MeV and m ρ = 763 MeV, respectively.The adopted RMF model then has 8 adjustable parameters of m σ , g σ , g ω , g ρ , b σ , c σ , c ω and Λ V .As in the Skyrme EDF, these parameters can be analytically deduced from pseudo-observables in nuclear matter (see., e.g., [55,56]), and we use following Ref.[56].Here, the m * Dirac is the nucleon Dirac effective mass in symmetric nuclear matter at saturation density.Based on the Lagrangian density given in Eq.( 5), the DIRHF package [57] is modified to calculate the ground-state properties of finite nuclei by solving the stationary relativistic Hartree-Bogoliubov equations.The pairing force is deduced from the D1S Gogny force (see Ref. [57] for details).
Using the sampled 50 Skyrme EDFs and 50 covariant EDFs, we calculate the ∆µ * pn of 5 doubly magic nuclei 48 Ca, 68 Ni, 88 Sr, 132 Sn and 208 Pb as well as the symmetry energy at subsaturation densities.The Pearson correlation coefficients R between the ∆µ * pn and the E sym (ρ) within the Skyrme EDF and the nonlinear RMF model are calculated and shown as functions of nucleon density ρ in Fig. 1 (a) and (b), respectively.As can be seen, although the Pearson correlation coefficient depends on specific model and nucleus, they all reach a maximum larger than 0.96 at subsaturation densities of 1/2 ≤ ρ/ρ 0 ≤ 3/4, which indicates a strong linear correlation between the ∆µ * pn and the symmetry energy at subsaturation densities.
To see more clearly the ∆µ * pn -E sym (2ρ 0 /3) correlations, we plot in Fig. 2 data-to-data relations predicted by the sampled parameter sets for the Skyrme EDF and the nonlinear RMF model by open circles and triangles, respectively.The experimental values [53] for the ∆µ * pn are indicated by the grey dotted lines.One sees that, the predicted values of ∆µ * pn for the chosen 5 doubly magic nuclei are rather sensitive to the E sym (2ρ 0 /3) within both the two models, although there exist model discrepancies on the predicted ∆µ * pn -E sym (2ρ 0 /3) relations.Given the strong ∆µ * pn -E sym (2ρ 0 /3) correlations, we take the predictions of the sampled EDFs as training points to construct emulators using the surmise python package [58] for the ∆µ * pn (N, Z) as functions of E sym (2ρ 0 /3) via Gaussian processes (see, e.g., [59]).Particularly, the effects of the other model parameters except E sym (2ρ 0 /3) are mimicked by noise terms in GPs.The GP predictions from the Skyrme EDF and the nonlinear RMF model are shown as the red and blue lines labeled by GP-Sky and GP-RMF in Fig. 2, respectively.The shaded region represent the 1σ uncertainty bands of GP predictions.It is seen that the GPs well reproduce the dependence of the ∆µ * pn on the E sym (2ρ 0 /3) and reasonably describe the model uncertainties.
Based on the GPs constructed from the Skyrme EDFs and the nonlinear RMF models, we further extract the E sym (2ρ 0 /3) from Bayesian analyses of the ∆µ * pn of 5 doubly magic nuclei, 48 Ca, 68 Ni, 88 Sr, 132 Sn and 208 Pb.The prior of E sym (2ρ 0 /3) is taken to be uniform in the range of 15-35 MeV, and the likelihood is given by where the model parameter θ is E sym (2ρ 0 /3), and y and y exp are the GP predictions and experimental values for the ∆µ * pn of the 5 doubly magic nuclei, respectively.The covariance matrix Σ = Σ GP +diag(σ 2 ) includes the uncertainties and correlations of the emulator predictions via the GP covariance matrix Σ GP , and the adopted error σ which takes into account the deficiency of theoretical models and experimental errors.In prior, the model deficiency is unknown and we assume that σ lies uniformly between 0 and 10 MeV.The posterior distributions of E sym (2ρ 0 /3) and σ are then generated using sequential Monte Carlo algorithm from PyMCv4.0 [60].Based on GP-Sky and GP-RMF, the posterior σ are inferred to be 0.4 +0.4 −0.2 and 0.8 +0.6 −0.3 MeV, respectively.The smaller σ value of GP-Sky indicates the Skyrme EDF is more capable of simultaneously reproducing the ∆µ * pn of the 5 doubly magic nuclei.The posterior distributions of E sym (2ρ 0 /3) from GP-Sky and GP-RMF are shown in Fig. 3 as the red dashed and blue dashed-dotted lines, respectively.Quantitatively, we obtain MeV and E sym (2ρ 0 /3) = 24.9 ± 1.1 MeV at 68.3% confidence level from the Skyrme and RMF models, respectively.To handle with the model dependence of the inferred E sym (2ρ 0 /3), we further carried out Bayesian model averaging [38][39][40] over the results from the two models.Within the BMA approach, the posterior distribution of any quantity O is obtained by mixing the posterior predictive distributions p(O|y, M i ) of each individual model M i with a probability weight p(M i |y), i.e., p (O|y) = i p(O|y, M i )p(M i |y). ( The model weight, reflecting the validity of M i given experimental data y, is given by where the prior model probability π(M i ) reflects our preference on M i before seeing the data, and the evidence p (y | M i ) measures the probability that the model reproduces the experimental data y with given prior distributions of model parameter and adopted errors π(θ, σ|M i ), i.e., In our case, in prior, the Skryme and nonlinear RMF models are equally weighted with π(M i ) = 1/2, and the evidences p (y | M i ) of the two models are estimated with the sequential Monte Carlo method.The obtained ratio of the GP-Sky evidence to that of GP-RMF is estimated to be 3.3, which means the measured ∆µ * pn for the chosen 5 doubly magic nuclei slightly prefers the Skyrme EDF rather than the nonlinear RMF model.Invoking Eqs.( 7) and ( 8), we calculate the BMA posterior distribution of the E sym (2ρ 0 /3) based on those from the Skyrme EDF and nonlinear RMF model, and show it as the purple solid line in Fig. 3. Quantitatively, the E sym (2ρ 0 /3) is inferred to be 25.6 From similar analyses, we extract the symmetry energy from ρ 0 /2 to 3ρ 0 /4 by BMA over the predictions of the Skyrme EDF and the nonlinear RMF model, and show in Fig. 4 the median value and the 1σ confidence region as the purple solid line and shaded band, respectively.For comparison, Fig. 4 also includes the results from fitting various terrestrial and astrophysical constraints with a cubic polynomical form of the potential part of symmetry energy (Lynch & Tsang) [61], and the predictions of many-body perturbation calculations using N3LO chiral interactions with momentum cutoff Λ = 450 and 500 MeV (N3LO-Λ450 and N3LO-Λ500 ) [62].Note that for the N3LO-Λ450 and N3LO-Λ500 results, the saturation density is taken to be their respective predicted values of 0.17 and 0.173 fm −3 .In addition, we also exhibit as symbols in Fig. 4 six constraints on the symmetry energy around 2ρ 0 /3 extracted from the properties of doubly magic nuclei (Brown) [63], and the giant dipole resonance (GDR) in 208 Pb [64], the binding energy difference of heavy isotope pairs (∆B) [65], the neutronproton Fermi energy difference (∆ϵ F ) [47], and a recent Skyrme EDF analysis of PREX-2 and CREX results on the charge-weak form factor difference in 208 Pb and 48 Ca together with properties of doubly magic nuclei (Zhang & Chen) [16].Our results agree well with these constraints and predictions, and provide more statistically reliable predictions.
In conclusion, based on a number of Skyrme energy density functionals and nonlinear relativistic mean field model parameter sets, we extract the symmetry energy E sym (ρ) around 2ρ 0 /3 from a Bayesian model averaging analysis of the effective proton-neutron chemi-cal potential difference ∆µ * pn for 5 doubly magic nuclei, 48 Ca, 68 Ni, 88 Sr, 132 Sn and 208 Pb.Particularly, we infer E sym (2ρ 0 /3) = 25.6 +1.4  −1.3 MeV at 1σ confidence level.The obtained constraints on the E sym (ρ) around 2ρ 0 /3 are consistent with microscopic predictions and results from analyses of other isovector indicators, and are potentially useful in the studies of, e.g., the neutron drip line, r-process path, the inner crust of neutron stars and the core-collapse supernovae [6,66].
Since both the intra-and inter-model uncertainties are taken into account in our BMA analyses, the present results are statistically more reliable.Moreover, allowing further inclusion of more experimental observables and theoretical models, the employed BMA method provides a rigorous statistical approach towards the accurate determination of nuclear symmetry energy.

FIG. 2 .
FIG. 2. Effective proton-neutron potential differences ∆µ * pn of 48 Ca (a), 68 Ni (b), 88 Sr (c), 132 Sn (d) and 208 Pb (e) versus Esym(2ρ0/3) predicted by the sampled Skyrme EDFs (open cycles) and nonlinear RMF model parameter sets (open triangles).The predictions of Gaussian Process emulators trained using the Skyrme and RMF parameter sets are respectively shown as the red and blue lines.The shaded bands present the 1σ error bars of GP predictions.The experimental values [53] for the ∆µ * pn are indicated by the grey dotted lines.

FIG. 3 .
FIG. 3. Posterior probability density distribution (PDF) ofEsym(2ρ0/3) inferred from the Skyrme EDF (red dashed line), the RMF model (blue dashed-dotted line), and the Bayesian model averaging (purple solid line).The prior distribution of Esym(2ρ0/3) is plotted as the black dotted line.The grey shaded area indicates the 68.3% credible interval of posterior from BMA.
Pearson correlation coefficients between the Esym(ρ) and the ∆µ * pn in 48 Ca, 68 Ni, 88 Sr, 132 Sn and 208 Pb as functions of nucleon density ρ predicted by the sampled 50 Skyrme EDFs (left window) and nonlinear RMF parameter sets (right window).The vertical grey dotted line indicates ρ = 2ρ0/3.The green shaded region represents R > 0.96.