Characterizing the nuclear models informed by PREX and CREX: a view from Bayesian inference

Abstract

New measurements of the weak charge density distributions of ⁴⁸Ca and ²⁰⁸Pb challenge existing nuclear models. In the post-PREX-CREX era, it is unclear if current models can simultaneously describe weak charge distributions along with accurate measurements of binding energy and charge radii. In this letter, we explore the parameter space of relativistic and non-relativistic models to study the differences between the form factors of the electric and weak charge distributions, $\Delta F=F_{ch}-F_{w}$, in ⁴⁸Ca and ²⁰⁸Pb. We show, for the first time, the parts of the mean-field models which are the most important in determining the relative magnitude of the neutron skin in lead and calcium nuclei. We carefully disentangle the tension between the PREX/CREX constraints and the ability of the RMF and Skyrme models to accurately describe binding energies and charge radii. We find that the nuclear symmetry energy coefficient $S_{\mathrm{V}}$ and the isovector spin-orbit coefficient $b^{\prime}_{4}$ play different roles in determining $\Delta F$ of ⁴⁸Ca and ²⁰⁸Pb. Consequently, adjusting $S_{\mathrm{V}}$ or $b^{\prime}_{4}$ shifts predicted $\Delta F$ values toward or away from PREX/CREX measurements. Additionally, $S_{\mathrm{V}}$ and the slope L are marginally correlated given the constraints of our Bayesian inference, allowing us to infer them separately from PREX and CREX data.

Nuclear models describing finite nuclei and nuclear matter – One of the goals in nuclear physics is to construct a unified theory describing the basic properties of a variety of nuclei and neutron star observables. Historically, the binding energy, the charge radii and the response properties of nuclei have been used to calibrate nuclear models. The requirement of simultaneously describing the charge distributions of a variety of isotopes has spurred the improvement of both relativistic mean field (RMF) and non-relativistic Skyrme models, but has remained to be an outstanding challenge for the state-of-the-art nuclear structure theories Lalazissis et al. (1997a); König et al. (2023); Sommer et al. (2022); Malbrunot-Ettenauer et al. (2022); Reinhard and Nazarewicz (2022). It has also been recognized that the neutron distributions in finite nuclei are less sensitive to experimentally measurable quantities like binding energy and charge radii Horowitz and Piekarewicz (2001), leading to less-constrained isovector interactions in traditional nuclear models. Nonetheless, isospin-dependent interactions play a crucial role in describing neutron star properties. Recently, information about differences between the form factors of electric and weak charge distributions $\Delta F=F_{ch}-F_{w}$ of two neutron-rich magic nuclei, ²⁰⁸Pb and ⁴⁸Ca, has been model-independently measured via parity-violating electron scattering experiments Adhikari et al. (2021a, 2022a).Can existing forms of relativistic RMF models and non-relativistic Skyrme models simultaneously describe the weak charge distributions inferred by PREX and CREX, while maintaining consistency with the basic ground state properties of finite nuclei? If so, what would be the implications for the isospin-dependent interactions of such models? In this letter, we study the parameterizations of existing nuclear models informed by the PREX and the CREX experiments. The linear relationship between the isovector-dependent parts of the nuclear models and the neutron skin thickness of two nuclei (²⁰⁸Pb and ⁴⁸Ca) discovered in this work for the first time gives the plausible region of $S_{\mathrm{V}}$, $L$ and $b^{\prime}_{4}$ that are mainly constrained by the PREX and CREX, disentangling the information of iso-vector interactions given by neutron skin thickness measurements from those given by the other nuclei property measurements. This linear relationship together with the region of $S_{\mathrm{V}}$, $L$ and $b^{\prime}_{4}$ enhance our understanding of both the nuclei and neutron star properties.

The picture of $S_{\mathrm{v}}-L$ inferred by joint PREX and CREX measurements – For densities of relevance to nuclei, the energy of asymmetric nuclear matter can be expressed as the sum of symmetric nuclear matter (SNM) and symmetry energies as

$$E(n_{p},n_{n})=E(n_{p}=n_{n}=n/2)+S(n_{p},n_{n})~{},$$ (1)

where $n_{n}$ and $n_{p}$ are the neutron and proton densities, and $n$ is their sum. The SNM energy $E$ can be expanded around the saturation density $n_{0}$:

$$E(n_{p}=n_{n}=n/2)=-B+\frac{K_{0}}{18n_{0}^{2}}(n-n_{0})^{2}+...~{}.$$ (2)

Above, $B$ is the binding energy, and $K_{0}$ is the SNM incompressibility. To good approximation, the symmetry energy $S$ is given by the quadratic term in the expansion (asymmetry) parameter $(n_{n}-n_{p})/n$:

$$S(n_{p},n_{n})=S(n)\left(\frac{n_{n}-n_{p}}{n}\right)^{2}+...~{},$$ (3)

where $S(n)$ is the density-dependent symmetry energy. An alternative definition based on pure neutron matter (PNM), $S(n)=S(n_{n}=n,n_{p}=0)$, has also been used Gandolfi et al. (2014). The deviation of $S(n)$ between these two definitions represents the breaking of the quadratic approximation, which is usually small for most Skyrme models Lee et al. (1998); Zuo et al. (1999); Drischler et al. (2014). However, such a deviation can be larger in some other models Cai and Chen (2019). In RMF models, the inclusion of the $\delta$-meson breaks the quadratic approximation significantly as shown in Salinas and Piekarewicz (2024). Despite the differences in definitions, $S(n)$ can be expanded around the saturation density as

$$S(n) = S_{\mathrm{V}}+\frac{L}{3n_{0}}(n-n_{0})+\frac{K_{sym}}{18n_{0}}(n-n_{0})+...,$$ (4)

where $S_{\mathrm{V}}$ is the bulk symmetry energy, $L$ refers to the slope of the symmetry energy slope, and $K_{sym}$ to the curvature of the symmetry energy.

It has long been recognized that neutron skin thicknesses of neutron-rich nuclei are primarily determined by $S(n)$. The strong correlation between $L$ and the neutron skin thickeness has been extensively studied using Skyrme and RMF models of nuclei Horowitz et al. (2001); Reed et al. (2021). The neutron skin can also be understood using a simple droplet model in terms of the competition between the bulk symmetry term $S_{\mathrm{V}}\cdot I$ and the surface term $\Delta R_{np}\cdot Q$  Myers and Swiatecki (1980); Centelles et al. (2009); Warda et al. (2009). Here, $I=(N-Z)/A$ with $N$ and $Z$ being the neutron and proton numbers in a nucleus and $A$ their sum, $\Delta R_{np}$ is the neutron-proton radius difference, and $Q$ is the surface stiffness parameter. The first of these competing terms is closely related to the properties of bulk matter in the central regions of nuclei and the second relates also to the spin-orbit interaction.

It is not certain that the aforementioned models and approximations work equivalently well to explain the neutron skin thicknesses of all light-to-heavy nuclei. Given the possible sensitivity of neutron skin thickness to both $S_{\mathrm{V}}$ and $L$ (as suggested by previous works) and the weak correlation between them observed in both RMF and Skyrme models, in this work we study the dependence of $\Delta F^{\mathrm{Ca}48}$ and $\Delta F^{\mathrm{Pb}208}$ on $S_{\mathrm{V}}$ and $L$ by extensively surveying the parameter space of both RMF and Skyrme models. Comparing to previous work studying the information content of PREX and CREX Reinhard et al. (2022a); Zhang and Chen (2023); Salinas and Piekarewicz (2023), we apply minimum constraints that are sensitive to isovector interactions except for the PREX and CREX data. This “minimum setting” allows us to characterize the isospin-dependent features of nuclear models mainly informed by PREX and CREX data, rather than the existing extensive measurements of charge radii, binding energy and electric dipole polarizability of finite nuclei.

Figure 1: Symmetry energy parameters $S_{\mathrm{V}}$ and $L$ constrained by parity violating asymmetry of ⁴⁸Ca (CREX) and ²⁰⁸Pb (PREX). Dashed (dotted) line and the corresponding band represent mean and 1-$\sigma$ band of PREX (CREX) alone based on the correlations in Eq. (9) and Eq. (8). Dash-dotted ellipse represent 1-$\sigma$ contour of joint PREX and CREX constraints. The PREX-2 analysis Reed et al. (2021) as well as constraints from the unitary gas conjecture Tews et al. (2017) and dipole polarizability Piekarewicz et al. (2012) are shown for comparison. The horizontally shaded ellipse show the 1-$\sigma$ contour of the posterior with basic nuclei constraints as the yellow distribution for basic nuclei in Fig. 2, whereas the other two shaded ellipses show the posterior with the +PREX+CREX experiment for Skyrme (vertical shade) and RMF (diagonal shade) models, respectively. The green ellipse shows the UNDEF correlation considering properties of more nuclei Kortelainen et al. (2010). The small black ellipse shows the correlation from chiral EFT calculation of PNMDrischler et al. (2024).

In Fig. 1, we present the constraints of $S_{\mathrm{V}}$ and $L$ given the information of $\Delta F^{\mathrm{Ca}48}$ and/or $\Delta F^{\mathrm{Pb}208}$ without assuming the correlation between them. Given the large number of three-dimensional points of $(\Delta F^{\mathrm{Ca}48},\Delta F^{\mathrm{Pb}208},S_{\mathrm{V}})$ and $(\Delta F^{\mathrm{Ca}48},\Delta F^{\mathrm{Pb}208},L)$ sampled from the posterior distribution of RMF/Skyrme models, we find that $S_{\mathrm{V}}$ and $L$ can be well approximated by a linear combination of $\Delta F^{\mathrm{Ca}48}$ and $\Delta F^{\mathrm{Pb}208}$ as

$$S_{\mathrm{V}} = a\Delta F^{Ca48}+b\Delta F^{Pb208}+c$$ (5)

$$L = a^{\prime}\Delta F^{Ca48}+b^{\prime}\Delta F^{Pb208}+c^{\prime}~{}.$$ (6)

	$a$	$b$	$c$	$a^{\prime}$	$b^{\prime}$	$c^{\prime}$	$a^{\prime\prime}$	$b^{\prime\prime}$	$c^{\prime\prime}$
RMF	$-575.2\pm 5.1$	$916.3\pm 4.6$	$32.2\pm 3.7$	$2938.7\pm 43.5$	$2420.6\pm 33.9$	$-149.8\pm 25.6$	$-31.8\pm 0.2$	$23.5\pm 0.1$	$0.94\pm 0.11$
Skyrme	$-503.2\pm 7.8$	$945.2\pm 5.5$	$31.9\pm 2.9$	$1791.2\pm 27.2$	$2652.0\pm 19.0$	$-91.5\pm 10.1$	$-52.0\pm 0.4$	$27.8\pm 0.3$	$1.67\pm 0.15$

Table 1: Parameters for the fitting formulae of $S_{\mathrm{V}}$ (Eq. 5) and $L$ (Eq. 6) are in units of MeV, and for $b^{\prime}_{4}$ (Eq. 14), they are in units of fm⁴.

The slope coefficients $a$, $b$, $a^{\prime}$, and $b^{\prime}$, as well as the intercept parameters $c$ and $c^{\prime}$, are determined by fitting them to the posterior distribution of $(\Delta F^{Ca48},\Delta F^{Pb208},S_{\mathrm{V}})$ and $(\Delta F^{Ca48},\Delta F^{Pb208},L)$, respectively. These parameters may reflect the inherent structure of RMF and Skyrme models and could be insensitive to the choice of model parameters and the experimental constraints applied in the Bayesian inference. All the sampled points from the posterior that show clear linear relationships are presented in the supplemental material. The fitted slope and intercept parameters are summarized in Table 1. The precision of the fitting formula in Eq. (5) remains within 10% (1-$\sigma$) for both Skyrme and RMF models. Given the mean and standard deviation of $\Delta F^{Pb208}$ and $\Delta F^{Ca48}$ measured by PREX and CREX and the linear relationship in Eq. (5) and Eq. (6), we obtain the mean and covariance of $S_{\mathrm{V}}$ and $L$ to be

$$(\bar{S}_{\mathrm{V}},\bar{L}) = (56.7,~{}~{}66.8)_{Skyrme},~{}~{}~{}(53.8,~{}~{}30.9)_{RMF}$$ (7)

$$\sqrt{\mathrm{cov}} = \begin{pmatrix}19.6&31.2\\ 31.2&56.5\end{pmatrix}_{Skyrme},~{}~{}\begin{pmatrix}19.5&31.0\\ 31.0&66.3\end{pmatrix}_{RMF}$$

in units of MeV for Skyrme and RMF models, respectively. See the red and blue dashdotted ellipses in Fig. 1. The covariance is dominated by the experimental error of $\Delta F^{Ca48}$ or $\Delta F^{Pb208}$. The precision of the fitting formulas in Eq. (5) and Eq. (6) serves as a minor modification.

Based on the linear relationship inferred from the posterior distribution of the chosen RMF and Skyrme models, we further examine the $L$ and $S_{\mathrm{V}}$ values inferred from the PREX and CREX experiments, respectively. Equations (5) and (6) can be rewritten as

$$L = \frac{a^{\prime}}{a}S_{\mathrm{V}}-\frac{ba^{\prime}-b^{\prime}a}{a}\Delta F^{\text{Pb208}}-\frac{a^{\prime}c-c^{\prime}a}{a}$$ (8)

$$L = \frac{b^{\prime}}{b}S_{\mathrm{V}}-\frac{b^{\prime}a-ba^{\prime}}{b}\Delta F^{\text{Ca48}}-\frac{b^{\prime}c-c^{\prime}b}{b}$$ (9)

Given the $\Delta F^{\text{Pb208}}$ and $\Delta F^{\text{Ca48}}$ values measured by PREX and CREX, respectively, Eqs. (8) and (9) define a band in the $S_{\mathrm{V}}$-$L$ plane, as shown in Fig. 1. Note that the slope coefficients $a^{\prime}/a$ and $b^{\prime}/b$ are very different. These differing slopes allow us to constrain $S_{\mathrm{V}}$ and $L$ independently using either $\Delta F^{\text{Ca48}}$ or $\Delta F^{\text{Pb208}}$. Our PREX constraint band on $S_{V}-L$ is superior to that of previous PREX-2 analysisReed et al. (2021) as shown in Fig.1, since we disentangle the PREX constraint from the $S_{V}-L$ correlation due to binding energies and charge radii. The overlapping region of these two bands represents the range of ${S_{\mathrm{V}},L}$ inferred by the 1-$\sigma$ experimental measurements from both PREX and CREX, which is quantitatively described by Eq. (7).

The constraints on $S_{\mathrm{V}}$-$L$ correlation outlined above are derived from a combination of the experimental measurements of the neutron skin thickness and the quasi-linear relationships between $\Delta F$ and ${S_{\mathrm{V}},L}$. Notably, these constraints remain independent of the constraints of $S_{\mathrm{V}}$-$L$ based on measurements of nuclei masses and radii, such as those depicted in the UNDEF constraint Kortelainen et al. (2010) shown in Fig. 1. Consequently, one should not anticipate a significant proportion of the sampled points from the posterior distribution of our Bayesian inference to fall within this overlapped region. This is because the Bayesian inference incorporates constraints not only from the PREX and CREX measurements, but also from other fundamental properties of finite nuclei, such as th binding energies and charge radii. Indeed, as illustrated in Fig. 1, the $S_{\mathrm{V}}$ values inferred from measurements of nuclei masses typically range from 30 to 32 MeV, which is in tension with those favored by the combined PREX and CREX measurements.

Figure 2: The posterior distribution of bulk symmetry energy $S_{\mathrm{V}}$ and symmetry energy slope $L$ at saturation density constrained by basic nuclei constraints (dashed yellow), +PREX experiment (dashdotted blue), +CREX experiment (dashdotted red), +PREX+CREX experiment (solid black). Dotted black curve includes additional constraint from dipole Polarizability Piekarewicz et al. (2012). Details of prior and likelihood are discussed in Supplemental Material.

In Fig. 2, we analyze the evolution of the posterior distributions for $S_{\mathrm{V}}$ and $L$ under five different constraint sets: basic nuclear constraints as the ”minimum setting”; inclusion of either PREX or CREX experimental data; inclusion of both PREX and CREX experimental data; and inclusion of additional dipole polarizability constraint. In the upper two panels for the Skyrme model, we observe that the range of $S_{\mathrm{V}}$ is approximately [20,50] MeV, whereas the range of $L$ spans from [0,150] MeV across all Bayesian inferences. The black dot-dashed curve in Fig. 1 illustrates the joint distribution. In the upper right (Skyrme) and lower right (RMF) panels, the posteriors of $L$ exhibit similar trends. The posterior of $L$ constrained by CREX (PREX) measurements has a mean smaller (larger) than that of the posterior constrained by only basic properties of finite nuclei. Moreover, the posterior of $L$ constrained by both PREX and CREX measurements has a mean lower than the one without any constraints from neutron skin measurements, yet slightly higher than the one including only the CREX constraint. This trend is consistent across both our RMF and Skyrme models and aligns with previous studies using Skyrme models Reinhard et al. (2021); Reinhard and Nazarewicz (2022); Zhang and Chen (2023) and RMF models Reed et al. (2021); Salinas and Piekarewicz (2023); Yüksel and Paar (2023). The overlapping integral suggests a probability of inconsistency of 93% (1.8-$\sigma$) for Skyrme models and 76% (1.2-$\sigma$) for RMF models, significantly lower than other analyses showing over 2-$\sigma$ inconsistency Zhang and Chen (2023); Salinas and Piekarewicz (2023); Yüksel and Paar (2023). Further details of our Bayesian inference and Monte Carlo sampling are provided in the supplemental material.

Impact of spin-orbit interactions – Unlike the Hamiltonian of infinite dense matter, the Hamiltonian of finite nuclei incorporates terms proportional to the derivative of isoscalar (vector) densities, primarily arising from spin-orbit interactions. The impact of spin-orbit interactions on single-particle spectra and the charge distributions of various neutron-rich isotopes has been acknowledged and discussed in prior studies Reinhard and Flocard (1995a); Horowitz and Piekarewicz (2012a). However, the isospin-dependent spin-orbit interactions are not well constrained, in comparison with other interactions constrained by isobaric analog states (IAS) Roca-Maza et al. (2018) and mirror nuclei Sagawa et al. (2024). Considering that the values of $F_{\mathrm{ch}}$ and $F_{\mathrm{W}}$ measured in both PREX and CREX are approximately the same, leading to a $\Delta F$ that is roughly one order of magnitude smaller, the influence of spin-orbit terms may play crucial roles in determining the distribution of points of ${\Delta F^{Ca48},\Delta F^{Pb208}}$. For the first time (as far as we are aware), we investigate the influence of spin-orbit interactions on both $\Delta F^{Ca48}$ and $\Delta F^{Pb208}$ across a very large parameter space encompassing both RMF and Skyrme parameterizations. This investigation is based on the posterior distribution incorporating constraints from both PREX and CREX measurements.

Figure 3: Charge and weak from factor difference of ⁴⁸Ca and ²⁰⁸Pb for samples of RMF (left) and Skyrme (right) model colored by the isovector spin-orbit parameter $b^{\prime}_{4}$ in unit fm⁴. Colorful lines corresponds to the minimum $\chi^{2}$ fitting of linear function same as Eq.(5) and Eq.(6) with $b^{\prime}_{4}=[-0.2,0,0.2,0.4,0.6,0.8]$. Constraints from unitary gas conjectureTews et al. (2017) and nuclei propertiesKortelainen et al. (2010) calculated from corresponding curves in Fig.1 based on the correlation Eq.(5) and Eq.(6). Blue dashed line shows 1-,2- and 3-$\sigma$ confidence region of joint PREX and CREX results. ‘x’ marker on left panels are RMF models listed in Supplementary material. Reference and $b^{\prime}_{4}$ for Skyrme models with marker ‘x’ is discussed in the text. ‘+’ marker represent the recently calibrated Skyrme modelZhang and Chen (2023), the relativistic Hartree-Bogoliubov modelYüksel and Paar (2023) and RMF model with tensor interactionSalinas and Piekarewicz (2024) including PREX and CREX experiment which take different form of interactions.

The spin-orbit interactions in Skyrme models are written as,

$$\begin{split}H_{\mathrm{SO}}=b_{4}\textbf{J}\cdot\triangledown n+b^{\prime}_{4}(\textbf{J}_{n}\cdot\triangledown n_{n}+\textbf{J}_{p}\cdot\triangledown n_{p}),\end{split}$$ (10)

where the $b_{4}$ and $b^{\prime}_{4}$ coefficients are corresponding to the isoscalar and isovector constributions of the spin-orbit potentials. To compare the influence of spin-orbit interactions in both relativistic and non-relativistic models in a parallel manner, we recast the Dirac equation in Schrödinger form using the Foldy-Wouthuysen transformation and the spin-orbit contributions in RMF models can be approximated as:

$$V^{p/n}_{SO} = -\frac{2b_{4}}{r}\frac{d\rho_{B}}{dr}-\frac{b^{\prime}_{4}}{r}\frac{d}{dr}\left[\rho_{B}\pm(\rho_{p}-\rho_{n})\right]~{},$$ (11)

with

$$b_{4} \approx \frac{1}{4m^{2}}\left(\frac{g_{\sigma}^{2}}{m^{2}_{\sigma}}+\frac{g_{\omega}^{2}}{m^{2}_{\omega}}-\frac{g_{\delta}^{2}}{4m^{2}_{\delta}}-\frac{g_{\rho}^{2}}{4m^{2}_{\rho}}\right),$$ (12)

$$b^{\prime}_{4} \approx \frac{1}{8m^{2}}\left(\frac{g_{\delta}^{2}}{m^{2}_{\delta}}+\frac{g_{\rho}^{2}}{m^{2}_{\rho}}\right).$$ (13)

The detailed derivation of spin-orbit contributions in RMF models is provided in the supplemental material.

In Fig. 3, we present the distributions of $\{\Delta F^{Ca48},\Delta F^{Pb208}\}$ colored by $b^{\prime}_{4}$ in both RMF and Skyrme models. Interestingly, we observe that an increase of $b^{\prime}_{4}$ (which represents the isovector contribution to the spin-orbit interactions) leads to a larger $\Delta F^{Pb208}$ while resulting in a smaller $\Delta F^{Ca48}$. This causes the points of ${\Delta F^{Ca48},\Delta F^{Pb208}}$ to collectively move towards the mean of joint PREX and CREX measurements. Neutron skin thickness $\Delta R_{np}$ follows the same trend as $\Delta F$ when varying $b^{\prime}_{4}$. The detailed proton and neutron density profiles are presented in the supplemental material. We can fit the parameter $b^{\prime}_{4}$ as:

$$b^{\prime}_{4}=a^{\prime\prime}\Delta F^{Ca48}+b^{\prime\prime}\Delta F^{Pb208}+c^{\prime\prime},$$ (14)

where the corresponding fitting parameters $a^{\prime\prime}$, $b^{\prime\prime}$, and $c^{\prime\prime}$ for Skyrme and RMF models are shown in Table 1. The accuracy of the fitting is captured in $c^{\prime\prime}=1.69\pm 0.15$ ($0.93\pm 0.11$) in Skyrme (RMF), which is well within the experimental error of PREX and CREX. Based on this correlation, $b^{\prime}_{4}=1.37\pm 0.49$ fm⁴ for Skyrme models can be estimated from PREX and CREX measurements of $\Delta F$, corresponding to a 90% lower bound of $b^{\prime}_{4}>0.744$ fm⁴. The widely used SLy4 model has $W_{0}\in 123$ MeV fm⁵Chabanat et al. (1998), equivalent to $b_{4}=b^{\prime}_{4}=0.312$ fm⁴, which lies on the 2-$\sigma$ contour of the PREX and CREX measurement, as seen in Fig. 3. The first Skyrme models, SkI3 and SkI4, with $b_{4}\neq b^{\prime}_{4}$, are introduced with $b^{\prime}_{4}\leq 0<b_{4}$, placing them further away from the PREX and CREX measurement Reinhard and Flocard (1995b). SV classes of Skyrme model prefer $0<b^{\prime}_{4}/b_{4}<0.5$, see Fig. 9 in Klüpfel et al. (2009). TOV-min take a relatively large $b^{\prime}_{4}=0.426$Erler et al. (2013), placing it closer towards the right bottom corner in Fig. 3. More recent calibrated models UNDEF0 with $b_{4}=0.634$, $b^{\prime}_{4}=-0.463$ and UNDEF1 with $b_{4}=0.194$, $b^{\prime}_{4}=0.362$ have $b^{\prime}_{4}$ of opposite signKortelainen et al. (2012). Therefore, UNDEF1 is much closer to the PREX and CREX measurements compared to UNDEF0. In the case of RMF models, the effective parameter is $b^{\prime}_{4}=1.02\pm 0.37$, corresponding to a 90% lower bound of $b^{\prime}_{4}>0.544$ fm⁴ or $g_{\delta}^{2}+1.65g_{\rho}^{2}>2433$ after incorporating the masses of nucleons and mesons. Widely used RMF models without delta mesons, such as FSUGold2 Chen and Piekarewicz (2014) with $g_{\delta}^{2}=80.4656$ and FSUGarnet Chen and Piekarewicz (2015) with $g_{\delta}^{2}=192.9274$, typically exhibit a small effective $b^{\prime}_{4}<0.1$. Consequently, such RMF models generally fall short of reaching the 2-$\sigma$ confidence region of PREX and CREX, as demonstrated by recent efforts like HPUA-CSharma et al. (2023). New RMF models featuring $\delta$ mesons with small isovector Yukawa coupling, such as BSRV with $g_{\delta}^{2}<20$ and $g_{\rho}^{2}<265$, corresponding to $b^{\prime}4\lesssim 0.1$, also do not show significant improvement Kumar et al. (2023). However, models like DINOa-c, which display large isovector Yukawa couplings ($g_{\delta}^{2}>1100$ and $g_{\rho}^{2}>800$), all surpass the 90% lower bound of $b^{\prime}_{4}>0.544$ Reed et al. (2024). Consequently, DINOa-c lies close to the 1-$\sigma$ confidence region of PREX and CREX, as depicted in Figure 3. Exploring the parameter space with large $g_{\delta}$ and $g_{\rho}$, we identify a new parameter set, CPREX1 ($b^{\prime}_{4}=0.652$), which resides on the 1-$\sigma$ upper (lower) edge of $\Delta F^{Ca48}$ ($\Delta F^{Pb208}$), and well within the 1-$\sigma$ joint contour. CPREX1 also yields reasonable nuclear properties for${}^{48}Ca$, ${}^{90}Zr$, and ${}^{208}Pb$, along with decent neutron star properties, featuring a maximum mass of M${}_{max}=$ 2.04 M_⊙ and a dimensionless tidal deformability of $\Lambda_{1.4}=584$ within the GW170817 constraintAbbott et al. (2018). Additionally, we introduce CPREX2 with an intermediate $b^{\prime}_{4}=0.386$, closely aligned with Skyrme models like UNDEF1 and SV. Model parameters and other properties of CPREX are tabulated in the supplementary material. Their corresponding neutron star EOSs with crust calculated from the liquid droplet model are available from the public GitHub repositoryZhao .

Conclusion and discussion – In this letter, we first demonstrated that by using the existing forms of nuclear models, it’s possible to find a set of RMF parameterizations that falls within the 1-$\sigma$ of the PREX and CREX respectively, while maintaining good agreement ($<1\%$ deviation) with the binding energy, the charge radii of ${}^{208}\mathrm{Pb}$, ${}^{90}\mathrm{Zr}$ and ${}^{48}\mathrm{Ca}$ as well as tidal deformability and maximum mass of neutron stars.

Based on large samples of RMF and Skyrme models from our Bayesian inference, we further identified a linear equation describing the relationship between the symmetry energy parameters ($S_{\mathrm{V}}$ and $L$) and the $\Delta F$ values of ²⁰⁸Pb and ⁴⁸Ca. This linear relationship exhibits a weak dependency on specific choices of model parameterizations. Unlike $L$, which is positively related to the $\Delta F$ values of both the ²⁰⁸Pb and ⁴⁸Ca, $S_{\mathrm{V}}$ impacts ²⁰⁸Pb and ⁴⁸Ca differently. This suggests that $S_{\mathrm{V}}\gtrsim 40$ MeV is favored by PREX and CREX measurements. The large $S_{\textrm{V}}$ compared with empirical value (around 32 MeV) demonstrates the tension between the current nuclear models and joint PREX and CREX experiments.

Besides $S_{\mathrm{V}}$, we found that the isovector parts of the spin-orbit interactions corresponding to the $b^{\prime}_{4}$ term in the Skyrme model also influence the neutron skin of ²⁰⁸Pb and ⁴⁸Ca differently. We derived the equivalent density-independent $b^{\prime}_{4}$ parameter in terms of $g_{\delta}$ and $g_{\rho}$ in the RMF model and observed a similar trend. Although the correlation between $b^{\prime}_{4}$ and $\Delta F^{\mathrm{Ca}48/\mathrm{Pb}208}$ may not be strong, as also observed in Reinhard et al. (2022b), $b^{\prime}_{4}$ can be quite sensitive to a linear combination of $\Delta F^{\mathrm{Ca}48}$ and $\Delta F^{\mathrm{Pb}208}$, such as $a^{\prime\prime}\Delta F^{\mathrm{Ca}48}+b^{\prime\prime}\Delta F^{\mathrm{Pb}208}$, and can effectively control the discrepancy between the model predictions of ${\Delta F^{\mathrm{Pb208}},\Delta F^{\mathrm{Ca48}}}$ and the combined PREX and CREX measurements. We also realized that in RMF models, the increase of $b^{\prime}_{4}$ usually comes with the amplification of the (isovector) nucleon density fluctuation in the nuclei.

In RMF models, the large $b^{\prime}_{4}$ is achieved by large $g_{\delta}$ and $g_{\rho}$, which break the common quadratic dependence of symmetry energy on asymmetry. For example, the CPREX1 RMF model, calibrated to reach the $1-\sigma$ bound of PREX and CREX, ended up with $S_{\mathrm{V}}=32.9$ MeV around SNM, but $S_{\mathrm{V}}=54.3$ MeV defined with PNM. Although the slope $L$ is close to zero, the neutron star EOS in $\beta$-equilibrium remains stable due to the large L defined with PNM. CPREX1 and CPREX2 could be interesting for both the nuclear and astrophysics communities, as the neutron star properties satisfy maximum mass and tidal deformability constraints.

Follow-up neutron skin experiments, like the Mainz Radius EXperiment (MREX), will continue to refine the possible uncertainty of the PREX experiment Sfienti . The linear relations between the isovector parameters ($S_{\mathrm{V}}$, $L$, and $b^{\prime}_{4}$) and the $\Delta F$ values of ²⁰⁸Pb and ⁴⁸Ca may be extended to other double magic nuclei, serving as a simple “emulator” for rapidly and reliably approximating $S_{\mathrm{V}}$, $L$, and $b^{\prime}_{4}$ in the investigations of future neutron skin thickness experiments.

I Supplemental Mateiral

I.1 Relativistic mean field model and Skyrme model

		Experiment	NL3	FSU2	IOPB-I	IUFSU	BigApple	HPUC	BSRV	DINOa	DINOb	DINOc	CPREX1	CPREX2
	$B/A$ [MeV]	7.87	7.88	7.87	7.86	7.88	7.85	7.85	7.84	7.87	7.87	7.87	7.84	7.86
	$R_{ch}$ [fm]	5.50	5.51	5.49	5.52	5.49	5.50	5.56	5.53	5.51	5.51	5.51	5.49	5.49
208Pb	$\Delta R_{np}$ [fm]	0.159$\pm$ 0.017	0.2797	0.2862	0.2195	0.1618	0.1508	0.1196	0.2595	0.1746	0.1993	0.2235	0.1905	0.1525
	$F_{ch}$ []	0.409	0.4067	0.4094	0.4052	0.4106	0.4080	0.3992	0.4043	0.4074	0.4075	0.4073	0.4100	0.4092
	$\Delta F$ []	0.041$\pm$0.013	0.0414	0.0423	0.0319	0.0233	0.0214	0.0168	0.0378	0.0262	0.0303	0.0342	0.0282	0.0222
	$B/A$ [MeV]	8.67	8.65	8.62	8.64	8.53	8.52	8.65	8.66	8.67	8.67	8.67	8.64	8.66
	$R_{ch}$ [fm]	3.48	3.45	3.43	3.45	3.44	3.46	3.46	3.44	3.47	3.47	3.47	3.48	3.46
48Ca	$\Delta R_{np}$ [fm]	0.137$\pm$0.015	0.2255	0.2318	0.1995	0.1736	0.1690	0.1479	0.2196	0.0994	0.1054	0.1141	0.1252	0.1357
	$F_{ch}$ []	0.158	0.1604	0.1665	0.1616	0.1647	0.1582	0.1577	0.1621	0.1591	0.1589	0.1585	0.1537	0.1571
	$\Delta F$ []	0.0277$\pm$0.0055	0.0551	0.0564	0.0490	0.0435	0.0413	0.0391	0.0527	0.0330	0.0345	0.0364	0.0335	0.0362

Table 2: Experimental data for the binding energy per nucleonWang et al. (2012), charge radiiAngeli and Marinova (2013), neutron skins (excluding PREX and CREX)Lattimer (2023), charge from factor and form factor difference from PREXAdhikari et al. (2021b) for ²⁰⁸Pb and CREXAdhikari et al. (2022b) for ⁴⁸Ca. Also displayed are the theoretical results obtained with NL3Lalazissis et al. (1997b), FSUGold2Chen and Piekarewicz (2014), IOPB-IKumar et al. (2018), IUFSUFattoyev et al. (2010), BigAppleFattoyev et al. (2020), HPUCSharma et al. (2023), BSRVKumar et al. (2023), DINOa-cReed et al. (2024) and the two new parameterizations, CPREX1 and CPREX2.

		NL3	FSU2	IOPB-I	IUFSU	BigApple	HPUC	BSRV	DINOa	DINOb	DINOc	CPREX1	CPREX2
	ns [fm-3]	0.1483	0.1504	0.1495	0.1546	0.1546	0.1490	0.1480	0.1522	0.1525	0.1519	0.1516	0.1518
	M∗ [MeV]	558.7	557.0	557.2	572.1	572.8	572.9	565.3	587.4	593.0	593.9	692.8	648.1
	B [MeV]	16.24	16.26	16.10	16.40	16.34	15.98	16.10	16.16	16.21	16.21	16.29	16.14
SNM	K [MeV]	271.6	237.7	222.6	231.3	227.0	220.2	227.2	210.0	207.0	206.0	223.8	223.5
	SV [MeV]	37.3	37.6	33.3	31.3	31.3	28.4	36.1	31.4	33.1	34.6	32.9	29.8
	L [MeV]	118.2	112.7	63.6	47.2	39.8	41.6	84.6	50.0	70.0	90.0	-3.5	0.4
	Ksym [MeV]	101.0	25.4	-37.0	28.5	87.5	81.1	-73.2	506.0	609.1	714.8	-418.4	-239.8
	M${}^{*}_{n}$ [MeV]	569.2	566.0	566.7	580.5	582.8	581.4	573.3	352.1	333.0	320.5	377.4	465.6
	M${}^{*}_{p}$ [MeV]	569.2	566.0	566.7	580.5	582.8	581.4	574.8	908.8	948.2	969.1	1062.5	870.1
PNM	SV [MeV]	38.3	38.6	34.7	32.9	33.1	29.9	37.2	46.5	50.6	53.4	54.3	38.4
	L [MeV]	121.2	115.9	67.7	49.5	40.6	42.7	88.7	172.1	216.4	247.8	211.2	75.9
	Ksym [MeV]	100.3	27.2	-45.5	23.1	74.3	89.2	-70.6	726.7	907.2	1021.2	801.8	76.4
	Mmax [M⊙]	2.77	2.07	2.15	1.94	2.60	2.05	2.04	2.17	2.15	2.15	2.04	2.12
	R1.0 [km]	14.4	14.1	13.2	12.6	12.8	12.6	13.6	14.4	14.8	15.1	13.9	12.9
NS	R1.4 [km]	14.5	13.9	13.2	12.6	13.1	12.8	13.4	14.4	14.6	14.9	13.4	12.9
	$\Lambda_{1.0}$ []	7797	6473	4347	3384	3918	3752	4903	6623	7572	8579	4543	3544
	$\Lambda_{1.4}$ []	1275	876	687	500	719	593	689	1065	1150	1256	584	570

Table 3: Saturation properties and neutron star properties of RMF models listed in Table 2. Saturation properties for SNM and PNM are defined in the letter. Neutron star properties are calculated with the crust EOSs constructed with the compressible liquid droptlet model repectively for various RMF models with fixed surface tension parameters $\sigma_{s}=1.2$ MeV fm^-2, $S_{S}=48$ MeVLattimer and Swesty (1991).

Table 4: $\sigma$ meson mass and 7 coupling constant of RMF models.

	$m_{\sigma}$ (MeV)	$g_{\sigma}^{2}$	$g_{\delta}^{2}$	$g_{\omega}^{2}$	$g_{\rho}^{2}$	$\kappa$ (MeV)	$\lambda$	$\Lambda_{\omega\rho}$	$\zeta$
CPREX1	452.852	60.4489	1348.64	96.8552	948.802	12.1881	-0.0342058	0.00227575	0.00587231
CPREX2	484.272	77.601	754.613	120.463	589.842	7.6226	-0.0211199	0.00335068	0.0116056

In this section, we provide an overview of the theoretical framework utilized. For a more comprehensive understanding, interested readers can refer to the detailed discussions available in various sources such as Kumar et al. (2017, 2018); Yang and Piekarewicz (2020) along with their cited references. The starting point for our relativistic calculation of nuclear response is based on the covariant model presented in Ref. Chen and Piekarewicz (2014). Additionally, we incorporate the $\delta$-meson to accommodate the CREX experiment, which prefers a softer symmetry energy Adhikari et al. (2022b). The interacting Lagrangian density can be expressed as:

$${\mathscr{L}}_{\rm int} = \bar{\psi}\left[g_{\sigma}\sigma+\frac{g_{\delta}}{2}{\mbox{\boldmath$\tau$}}\cdot{\mbox{\boldmath$\delta$}}\!-\!\left(g_{\omega}\omega_{\mu}\!+\!\frac{g_{\rho}}{2}{\mbox{\boldmath$\tau$}}\cdot{\mbox{\boldmath$\rho$}}_{\mu}\!+\!\frac{e}{2}(1\!+\!\tau_{3})A_{\mu}\right)\gamma^{\mu}\right]\psi$$ (15)

$$- \frac{\kappa}{3!}(g_{\sigma}\sigma)^{3}\!-\!\frac{\lambda}{4!}(g_{\sigma}\sigma)^{4}\!+\!\frac{\zeta_{\omega}}{4!}g_{\omega}^{4}(\omega_{\mu}\omega^{\mu})^{2}+\Lambda_{\omega\rho}\Big{(}g_{\rho}^{2}\,{\mbox{\boldmath$\rho$}}_{\mu}\cdot{\mbox{\boldmath$\rho$}}^{\mu}\Big{)}\Big{(}g_{\omega}^{2}\omega_{\nu}\omega^{\nu}\Big{)}\;.$$

The Lagrangian density used in our study includes the isodoublet nucleon field $\psi$ as the basic degree of freedom. These nucleons interact via photon ($A_{\mu}$) exchange and through the exchange of four massive mesons: a scalar-isoscalar ($\sigma$), scalar-isovector ($\delta$), a vector-isoscalar ($\omega^{\mu}$), and a vector-isovector ($\rho^{\mu}$). This formulation is based on the works of Walecka (1974); Serot and Walecka (1986, 1997). The coupling constant $\Lambda_{\omega\rho}$ for the isoscalar-isovector interaction plays a crucial role in determining the density dependence of the symmetry energy, particularly its slope at saturation density $L$ as discussed in Horowitz and Piekarewicz (2001). It is noteworthy that the model can be calibrated using an analytic one-to-one correspondence between bulk parameters of infinite nuclear matter and the various coupling constants. Specifically, the values of the symmetry energy $S_{\mathrm{V}}$, its slope $L$, and curvature $K_{sym}$ at saturation density are correlated to the isovector parameters ($g_{\delta}$, $g_{\rho}$, and $\Lambda_{\omega\rho}$) as explained in Chen and Piekarewicz (2014); Singh et al. (2014). Table 4 shows the $\sigma$ meson mass and 8 coupling constants for CPREX1 and CPREX2 proposed in this letter. The masses of other massive mesons and the free nucleon are $m_{\delta}=980$ MeV, $m_{\omega}=782.5$ MeV, $m_{\rho}=763$ MeV and $m=939$ MeV. The nuclear properties of ⁴⁸Ca and ²⁰⁸Pb are shown in Table 2 along with experimental measurement and our Hatree calculations for various models for comparison. Table 3 shows the infinite nuclear matter properties at saturation density and the neutron star properties calculated from EOSs in $\beta$-equilibrium. The crust EOSs are constructed with the compressible liquid droplet model with fixed surface tension parameters $\sigma_{s}=1.2$ MeV fm^-2, $S_{S}=48$ MeVLattimer and Swesty (1991). The tabulated complete EOSs of these RMF models can be found in the GitHub repository [].

To gain a complementary understanding of both PREX and CREX measurements, we also study the Skyrme interactions inferred from those neutron skin thickness measurements. The Skyrme Hamiltonian for infinite nuclear matter is written as,

$$H_{\mathrm{Sk}}=\frac{k_{Fn}^{5}}{10\pi^{2}m_{n}^{*}}+\frac{k_{Fp}^{5}}{10\pi^{2}m_{p}^{*}}+H_{\mathrm{pot}}(n_{n},n_{p}),$$ (16)

where the first two terms are the kinetic terms for the neutrons and the protons with effective mass contribution and $H_{\mathrm{pot}}$ is the potential term expressed as:

$$\begin{split}H_{\mathrm{pot}}=&\frac{1}{2}n^{2}t_{0}\left[1+\frac{x_{0}}{2}\right)-\frac{1}{2}(n_{n}^{2}+n_{p}^{2})t_{0}\left(\frac{1}{2}+x_{0}\right)\\ +&\frac{1}{24}n^{\gamma}t_{3}\left(n^{2}(2+x_{3})-(n_{n}^{2}+n_{p}^{2})(1+2x_{3})\right].\end{split}$$ (17)

The neutron and proton effective masses are functions of densities as,

$$\begin{split}\frac{m_{n/p}}{m^{*}_{n/p}}=&1+\frac{m_{n/p}}{4}\Big{\{}n\big{[}t_{1}(2+x_{1})+t_{2}(2+x_{2})\big{]}\\ +&n_{n/p}\big{[}-t_{1}(1+2x_{1})+t_{2}(1+2x_{2})\big{]}\Big{\}}\,.\end{split}$$ (18)

To describe a finite nucleus, two additional terms are included in the Skyrme Hamiltonian Reinhard and Flocard (1995a); Chabanat et al. (1998). These terms are:

$$\begin{split}H_{\mathrm{SO}}=b_{4}\textbf{J}\cdot\triangledown n+b^{\prime}_{4}(\textbf{J}_{n}\cdot\triangledown n_{n}+\textbf{J}_{p}\cdot\triangledown n_{p}),\end{split}$$ (19)

and

$$\begin{split}H_{\mathrm{J}}=-\frac{1}{16}(t_{1}x_{1}+t_{2}x_{2})\textbf{J}^{2}+\frac{1}{16}(t_{1}-t_{2})(\textbf{J}^{2}_{n}+\textbf{J}^{2}_{p}),\end{split}$$ (20)

where $H_{\mathrm{SO}}$ are spin-orbit interactions and $H_{\mathrm{J}}$ are the central tensor terms. The $\bf{J}_{n/p}$ is the spin-orbit density $\bf{J}_{n/p}=\psi^{\dagger}_{n/p}\overrightarrow{\sigma}\times\overrightarrow{\triangledown}\psi_{n/p}$ and $\bf{J}=\bf{J}_{n}+\bf{J}_{p}$.

I.2 The approximated spin-orbit potential in RMF models

The large component of Dirac equation follows,

$$\left[\vec{\sigma}\cdot\vec{p}\frac{1}{2M+\epsilon-S-V}\vec{\sigma}\cdot\vec{p}-S+V\right]g(r) = \epsilon g(r)$$ (21)

where $S=g_{\sigma}\ \sigma_{0}\pm g_{\delta}\delta_{0}$, $V=g_{\omega}\ \omega\pm g_{\rho}\rho_{0}$, are mean field potential from classical expectation value of meson fields. By reducing the spherically symmetric Dirac equation with scalar and vector potentials to the non-relativistic Schrödinger equation to order $(p/m)^{2}$Horowitz and Serot (1981); Reinhard (1989), the effective single nucleon classical Hamiltonian is,

$$H = \vec{p}\frac{1}{2M^{*}}\vec{p}-S(r)+V(r)+V_{SO}\vec{s}\cdot\vec{l}$$ (22)

$$V_{SO} = -\frac{1}{r}\frac{d}{dr}\left[\frac{1}{M^{*}}\right]$$ (23)

where $\vec{l}=\vec{r}\times\vec{p}$, $\vec{s}=\vec{\sigma}/2$ and $M^{*}=m-S(r)$ assuming that the relativistic single particle energy $\epsilon+S(r)-V(r)<<2(m-S)$. However, it’s possible to rearrange terms of higher order in $(p/m)^{2}$ to obtain an equivalent form, $M^{*}=m-V(r)$, as shown in the Appendix of Reinhard (1989). Although these two expressions are equivalent to order $(p/m)^{2}$, their ratio can be of the order of 10%. Therefore, it’s reasonable to consider the average of these two choices,

$$M^{*} = m-\frac{S(r)+V(r)}{2},$$ (24)

which can also be justified by taking $\epsilon<<2m-S-V$ Ring (1996); Ebran et al. (2016). Since scalar and vector potentials are mainly determined by the corresponding densities of the two-form current, modified by Coulomb interaction and the self and mixing couplings of mesons, we can simplify by considering only Yukawa coupling to generate the potentials. Furthermore, vector densities are strongly correlated with scalar densities, allowing us to approximate scalar densities with vector densities:

$$S(r) = \frac{g_{\sigma}^{2}}{m^{2}_{\sigma}}n_{B}\pm\frac{g_{\delta}^{2}}{4m^{2}_{\delta}}(n_{p}-n_{n})$$ (25)

$$V(r) = \frac{g_{\omega}^{2}}{m^{2}_{\omega}}n_{B}\pm\frac{g_{\rho}^{2}}{4m^{2}_{\rho}}(n_{p}-n_{n})$$ (26)

where the upper sign is for protons and the lower sign for neutrons. The difference of isospin potential between proton and neutron is,

$$V^{p}_{SO}+V^{n}_{SO} = -\frac{d}{rdr}\left[\frac{M^{*}_{p}+M^{*}_{n}}{M^{*}_{p}M^{*}_{n}}\right]$$ (27)

$$V^{p}_{SO}-V^{n}_{SO} = \frac{d}{rdr}\left[\frac{M^{*}_{p}-M^{*}_{n}}{M^{*}_{p}M^{*}_{n}}\right]$$ (28)

$$M^{*}_{p}-M^{*}_{n} = -(\frac{g_{\delta}^{2}}{m^{2}_{\delta}}+\frac{g_{\rho}^{2}}{m^{2}_{\rho}})\frac{n_{p}-n_{n}}{4}$$ (29)

$$M^{*}_{p}+M^{*}_{n} = 2m-(\frac{g_{\sigma}^{2}}{m^{2}_{\sigma}}+\frac{g_{\omega}^{2}}{m^{2}_{\omega}})n_{B}$$ (30)

since $(M^{*}_{p}+M^{*}_{n})^{2}>>(M^{*}_{p}-M^{*}_{n})^{2}$, we can take $M^{*}_{p}M^{*}_{n}\approx(M^{*}_{p}+M^{*}_{n})^{2}/4$, keeping only linear term in vector densities,

$$V^{p}_{SO}+V^{n}_{SO} \approx -\frac{\frac{g_{\sigma}^{2}}{m^{2}_{\sigma}}+\frac{g_{\omega}^{2}}{m^{2}_{\omega}}}{rm^{2}}\frac{dn_{B}}{dr}$$ (31)

$$V^{p}_{SO}-V^{n}_{SO} \approx -\frac{\frac{g_{\delta}^{2}}{m^{2}_{\delta}}+\frac{g_{\rho}^{2}}{m^{2}_{\rho}}}{4rm^{2}}\frac{d(n_{p}-n_{n})}{dr}.$$ (32)

A similar form without the delta meson has been used to study the isospin dependence of the spin-orbit interaction, contrasting with the spin-orbit force in non-relativistic Skyrme models Sharma et al. (1995); Lalazissis et al. (1998); Pearson (2001). Isospin dependence of spin-orbit interaction is introduce to Skyrme models as:

$$V^{p/n}_{SO} = -\frac{2b_{4}}{r}\frac{dn_{B}}{dr}-\frac{b^{\prime}_{4}}{r}\frac{d}{dr}\left[n_{B}\pm(n_{p}-n_{n})\right]$$ (33)

to address the problem of isotope shifts in the Pb regionReinhard and Flocard (1995b). Note that older models with only one spin-orbit potential parameter, $W_{0}=2b_{4}=2b^{\prime}_{4}$ is about 100 MeV$\cdot$fm⁵Chabanat et al. (1997); Pearson (2001).

$$V^{p}_{SO}+V^{n}_{SO} = -\frac{4b_{4}+2b^{\prime}_{4}}{r}\frac{dn_{B}}{dr}$$ (34)

$$V^{p}_{SO}-V^{n}_{SO} = -\frac{2b^{\prime}_{4}}{r}\frac{d}{dr}\left[n_{p}-n_{n}\right]$$ (35)

By comparing Eq. (34-35) with Eq. (31-32), the model parameter $b_{4}$ and $b_{4}^{\prime}$ can expressed as:

$$b_{4} \approx \frac{1}{4m^{2}}\left(\frac{g_{\sigma}^{2}}{m^{2}_{\sigma}}+\frac{g_{\omega}^{2}}{m^{2}_{\omega}}-\frac{g_{\delta}^{2}}{4m^{2}_{\delta}}-\frac{g_{\rho}^{2}}{4m^{2}_{\rho}}\right),$$ (36)

$$b^{\prime}_{4} \approx \frac{1}{8m^{2}}\left(\frac{g_{\delta}^{2}}{m^{2}_{\delta}}+\frac{g_{\rho}^{2}}{m^{2}_{\rho}}\right).$$ (37)

In RMF models, we obtain $b_{4}$ and $b_{4}^{\prime}$ similar to those in non-relativistic models, as described in the above equations. In principle, $b_{4}$ and $b_{4}^{\prime}$ are density dependent in RMF models. We take the leading order term so that $b_{4}$ and $b_{4}^{\prime}$ are essentially a linear combination of coupling Yukawa coupling strength.

I.3 Weak Charge form factors with spin-orbit contributions

Form factors are key observable in electron scattering experiments. To calculate the weak charge form factor measured in PREX and CREX. We need to perform finite nuclei calculations to obtain the proton and neutron profiles then perform Fourier transform on the profile to get the corresponding form factor. In this work, we consider the impact of spin-orbit currents on the electromagnetic and weak charge form factors of ${}^{48}\mathrm{Ca}$ and ${}^{208}\mathrm{Pb}$. It has been shown the contributions from spin-orbit currents are comparable to the present CREX experimental error barsHorowitz and Piekarewicz (2012b), and therefore should be necessary to incorporate spin-orbit currents in the calculation of weak form factors as:

$$\begin{split}Z\mathrm{F}_{\mathrm{ch}}(q)&=\sum_{i=p,n}\{\mathrm{G}_{\mathrm{E}}^{i}(q^{2})\mathrm{F}_{\mathrm{V}}^{i}(q)+\frac{\mathrm{G}_{\mathrm{M}}^{i}(q^{2})-\mathrm{G}_{\mathrm{E}}^{i}(q^{2})}{1+\tau}\\ &\times[\tau\mathrm{F}_{\mathrm{V}}^{i}(q)+\frac{q}{2\mathrm{m}}\mathrm{F}_{\mathrm{T}}^{i}(q)]\}~{},\end{split}$$

and

$$\begin{split}Q_{\mathrm{W}}\mathrm{F}_{\mathrm{W}}(q)&=\sum_{i=p,n}\{\tilde{\mathrm{G}}_{\mathrm{E}}^{i}(q^{2})\mathrm{F}_{\mathrm{V}}^{i}(q)+\frac{\tilde{\mathrm{G}}_{\mathrm{M}}^{i}(q^{2})-\tilde{\mathrm{G}}_{\mathrm{E}}^{i}(q^{2})}{1+\tau}\\ &\times[\tau\mathrm{F}_{\mathrm{V}}^{i}(q)+\frac{q}{2\mathrm{m}}\mathrm{F}_{\mathrm{T}}^{i}(q)]\}~{},\end{split}$$

where the $\mathrm{G}^{i}_{\mathrm{E,M}}$ are electric and magnetic single nucleon form factors and the $\tilde{\mathrm{G}}^{i}_{\mathrm{E,M}}$ are electric and magnetic form factors for the weak-neutral current. Similar as Horowitz and Piekarewicz (2012b), we adopted a simple dipole parametrization of these single-nucleon form factors.

The vector and tensor form factors $\mathrm{F}_{\mathrm{v}}(q)$ and $\mathrm{F}_{\mathrm{T}}(q)$ are obtained by:

$$\mathrm{F}_{\mathrm{V}}(q)=\sum_{n\kappa}(2j+1)\int_{0}^{\infty}[g^{2}_{n\kappa}(r)+f^{2}_{n\kappa}(r)]j_{0}(qr)dr~{},$$

and

$$\mathrm{F}_{\mathrm{T}}(q)=\sum_{n\kappa}2(2j+1)\int_{0}^{\infty}g_{n\kappa}(r)f_{n\kappa}(r)j_{1}(qr)dr~{},$$

where $g_{n\kappa}(r)$ and $f_{n\kappa}$ are upper and lower component of Dirac spinor in RMF models. Wave functions in non-relativistic Skyrme model corresponds to the larger upper component $g_{n\kappa}(r)$. The lower component $f_{n\kappa}$ is obtained by assuming a free space relation Horowitz and Piekarewicz (2012b); Adhikari et al. (2022b).

$$f_{n\kappa}=\frac{1}{2\mathrm{m}}(\frac{d}{dr}+\frac{\kappa}{r})g_{n\kappa}(r).$$

Table 5: Prior distribution of Skyrme and RMF parameters. $n_{s}$, $B$, $M^{*}$, $K$ and $S_{V}$ map to $t_{0}$, $t_{1}$, $t_{2}$, $t_{3}$, $\gamma$ for Skyrme models and $g_{\sigma}$, $g_{\omega}$, $g_{\rho}$, $\kappa$, $\lambda$ for RMF models.

	parameter	prior
	$n_{s}$ [MeV]	[0.14,0.165]
	$B$ [MeV]	[-15.5,-16.5]
Both	$M^{*}$ [MeV]	[0.5,0.8]$\times$ 939
	$K$ [MeV]	[210,250]
	$S_{V}$ ([MeV]	[20,50]
	$m_{\sigma}$ [MeV]	[450,550]
	$m_{\delta}$ [MeV]	980
	$m_{\omega}$ [MeV]	782.5
RMF	$m_{\rho}$ [MeV]	763
	$L$ [MeV]	[$L^{-}$,$L^{+}$]
	$g_{\delta}^{2}[]$	[0,1500]
	$\zeta_{\omega}$ []	[0,0.03]
	$x_{0}$ []	[-1.81,2.15]
	$x_{1}$ []	[-7.53, 3.77]
Skyrme	$x_{2}$ []	[-49.90, 91.93]
	$x_{3}$ []	[-3.41, 3.73]
	$b_{4}$ [$\mathrm{fm}^{4}$]	[-0.36, 0.72]
	$b^{\prime}_{4}$ [$\mathrm{fm}^{4}$]	[-0.36, 0.72]

I.4 Prior

In this and the following two sections, we describe the prior, likelihood and Posterior in our Bayesian analysis. Our RMF model have 12 parameters. We fixed the $\delta$, $\omega$ and $\rho$ meson masses, $m_{\delta}=980$ MeV, $m_{\omega}=782.5$ MeV and $m_{\rho}=763$ MeV as in many FSU type models, leaving us with 9 free parameters. $\sigma$ meson mass is a free parameter with uniform distribution $[450,550]$ MeV. We also take uniform prior in SNM properties $n_{s}$, $BE$, $M^{*}$ and $K$, see Table 5. These 4 bulk nuclear properties can be mapped to 4 isoscalar coupling constant, as in Chen and Piekarewicz (2014), except for $\omega$ meson self-coupling $\zeta_{\omega}$ which take flat distribution $[0,0.03]$. Isovector sector is governed by $g_{\delta}$, $g_{\rho}$ and $\Lambda_{\omega\rho}$ which can be mapped to symmetry energy properties $S_{\mathrm{V}}$ and $L$, assuming a flat distribution of $g_{\delta}\in[0,1500]$. $S_{\mathrm{V}}$ take uniformly prior while $L$ is uniformly distributed in $[L^{-},L^{+}]$. Because $L^{-}$ and $L^{+}$ are not fixed but determined by other parameters. And flat distributions of these parameters finally result in a non-flat prior distribution of L,

In this and the following two sections, we describe the prior, likelihood, and posterior in our Bayesian analysis. Our RMF model has 12 parameters. We fixed the $\delta$, $\omega$, and $\rho$ meson masses at $m_{\delta}=980$ MeV, $m_{\omega}=782.5$ MeV, and $m_{\rho}=763$ MeV, as in many FSU-type models, leaving us with 9 free parameters. The $\sigma$ meson mass is a free parameter with a uniform distribution in the range $[450,550]$ MeV. We also take a uniform prior for the symmetric nuclear matter (SNM) properties $n_{s}$, $BE$, $M^{*}$, and $K$, as shown in Table 5. These 4 bulk nuclear properties can be mapped to 4 isoscalar coupling constants, as in Chen and Piekarewicz (2014), except for the $\omega$ meson self-coupling $\zeta_{\omega}$, which takes a flat distribution in the range $[0,0.03]$. The isovector sector is governed by $g_{\delta}$, $g_{\rho}$, and $\Lambda_{\omega\rho}$, which can be mapped to the symmetry energy properties $S_{\mathrm{V}}$ and $L$, assuming a flat distribution of $g_{\delta}$ in the range $[0,1500]$. $S_{\mathrm{V}}$ takes a uniform prior, while $L$ is uniformly distributed in $[L^{-},L^{+}]$. Because $L^{-}$ and $L^{+}$ are not fixed but determined by other parameters, the flat distributions of these parameters result in a non-flat prior distribution of $L$:

$$P(L)=\int^{\zeta_{\omega}^{+}}_{\zeta_{\omega}^{-}}\dots\int^{m_{\sigma}^{+}}_{m_{\sigma}^{-}}\int_{L^{-}(m_{\sigma},\dots,\zeta_{\omega})}^{L^{+}(m_{\sigma},\dots,\zeta_{\omega})}\delta(L^{\prime}-L)dL^{\prime}dm_{\sigma}\dots d\zeta_{\omega}$$ (38)

$L^{-}$ is a physical bound beyond which the quadratic formula for mapping $L$ to $\Lambda_{\omega\rho}$ has no solution Chen and Piekarewicz (2014). $L^{+}$ can be solved by setting $\Lambda_{\omega\rho}=0$, so that $L<L^{+}$ corresponds to $\Lambda_{\omega\rho}>0$, ensuring that $\beta$-equilibrium matter remains neutron-rich ($n_{p}<n_{n}$) at high density. Positive $\Lambda_{\omega\rho}$ is introduced to reduce the symmetry energy, similar to the positive $\zeta_{\omega}$, which serves to reduce the SNM energy. Thus, we take $\Lambda_{\omega\rho}=0$ as our natural upper bound for the symmetry energy slope in the RMF model.

The prior distribution of Skyrme parameters is chosen similarly to ensure consistency with the prior distribution used for the RMF models. We use flat prior distributions for $n_{s}$, $B$, $M^{*}$, $K$, and $S_{\mathrm{V}}$, as well as for the $x_{i}$, $b_{4}$, and $b^{\prime}_{4}$ parameters with bounds in Table 5. The range of the $x_{i}$, $b_{4}$, and $b_{4}^{\prime}$ parameters in the Skyrme prior distributions are determined as follows. First, we examined 11 widely used Skyrme models (SGII, NRAPR, UNEDF0, UNEDF2, SLy4, SV-min, SkO, SkOp, BSk16, KDE0v, and Gs). We then identified the maximum and minimum $x_{i}$, $b_{4}$, and $b^{\prime}_{4}$ parameters in these 11 models, denoted as $x_{i}^{\mathrm{max}}$, $x_{i}^{\mathrm{min}}$, $b_{4}^{\mathrm{max}}$, $b_{4}^{\mathrm{min}}$, $b_{4}^{\prime\mathrm{max}}$, and $b_{4}^{\prime\mathrm{min}}$. To cover a wide range of possible Skyrme parameters, we determine the prior distribution of $x_{i}$, $b_{4}$, and $b^{\prime}_{4}$ parameters as $[x_{i}^{\mathrm{min}}-\Delta x_{i},x_{i}^{\mathrm{max}}+\Delta x_{i}]$, $[b_{4}^{\mathrm{min}}-\Delta b_{4},b_{4}^{\mathrm{max}}+\Delta b_{4}]$, and $[b_{4}^{\prime\mathrm{min}}-\Delta b_{4}^{\prime},b_{4}^{\prime\mathrm{max}}+\Delta b_{4}^{\prime}]$, where $\Delta x_{i}=x_{i}^{\mathrm{max}}-x_{i}^{\mathrm{min}}$, $\Delta b_{4}=b_{4}^{\mathrm{max}}-b_{4}^{\mathrm{min}}$, and $\Delta b_{4}^{\prime}=b_{4}^{\prime\mathrm{max}}-b_{4}^{\prime\mathrm{min}}$.

Additionally, we constrain our RMF and Skyrme models to ensure that the pressure of matter in $\beta$-equilibrium is always increasing with increasing energy density and to ensure that the proton fraction in $\beta$-equilibrium never vanishes. Note that this doesn’t guarantee that symmetry energy slope $L$ defined around SNM has to be positive.

Table 6: The list of experiment and observation with adopted errors. Charge and weak form factor $F_{ch}$ and $F_{W}$ listed here correspond to momentum transfer $q=0.8733$ fm^-1 for ⁴⁸Ca in CREX or $0.3977$ fm^-1 for in PREX.

	Property	Mean	Standard Deviation
	$R_{ch}^{{}^{48}Ca}$ [fm]	3.48	0.070
	$R_{ch}^{{}^{90}Zr}$ [fm]	4.27	0.085
	$R_{ch}^{{}^{208}Pb}$ [fm]	5.50	0.11
Basic nuclei	BE${}^{{}^{48}Ca}$ [MeV]	8.67	0.173
constraints	BE${}^{{}^{90}Zr}$ [MeV]	8.71	0.174
	BE${}^{{}^{208}Pb}$ [MeV]	7.87	0.157
	$F_{ch}^{{}^{48}Ca}$ []	0.1581	0.001
	$F_{ch}^{{}^{208}Pb}$ []	0.409	0.001
CREX	$F_{W}^{{}^{48}Ca}-F_{ch}^{{}^{48}Ca}$ []	0.0277	0.0055
PREX	$F_{W}^{{}^{208}Pb}-F_{ch}^{{}^{208}Pb}$ []	0.041	0.013

I.5 Likelihood

In the Bayesian analysis of RMF and Skyrme models, we use the same constraints and uncertainties. To account for theoretical uncertainty, we incorporate a 2% uncertainty in the measured value of $R_{ch}$ and a 5% uncertainty in the standard deviation of $B/A$. It is also necessary to constrain $F_{ch}$ accurately to ensure that the electric weak form factor difference, $F_{ch}-F_{W}$, can be interpreted as a measurement in PREX and CREX. The constraints on $B/A$, $R_{ch}$, and $F_{ch}$ together form the basic nuclear constraints in our Bayesian analysis, as listed in Table 6.

Sampling with basic nuclear constraints results in drastically different posteriors for $L$. In the Skyrme model, negative $L$ is suppressed by the constraint that the pressure of matter in $\beta$-equilibrium is always increasing with energy density. However, RMF models with negative $L$ around SNM can have a large symmetry energy slope around PNM due to the breakdown of the quadratic formula of symmetry energy with large $g_{\delta}$ and $g_{\rho}$. Since our goal is not to constrain the nuclear model with basic nuclear constraints, we add an additional likelihood factor to our RMF models to ensure the $S_{\mathrm{V}}$-$L$ distribution is compatible with only basic nuclei constraints, as shown by the black dashed ellipse in Fig. 1 in the letter.

In addition to the basic nuclei constraints, we impose additional PREX and CREX seperately and jointly, forming 4 sets of likelyhoods: basic nuclei, +PREX, +CREX and +PREX+CREX. We take the electric weak form factor difference $F_{ch}-F_{W}$ instead of weak form factor $F_{W}$ as constraints of PREX and CREX experiments as tabulated in Table 6. In principle, weak form factor $F_{W}$ can be determined from the parity violating asymmetry $A_{PV}$ measured in PREX and CREX experiments by the Coulomb distortion calculation assuming a electric charge distribution. As the electric charge distribution has already been determined experimentally with great accuracy, both $F_{W}$ and $F_{ch}-F_{W}$ are determined experimentally from the measured $A_{PV}$ without any model dependence. However, the theoretical models like RMF and Skyrme models optimized with charge radius doesn’t have enough accuracy to ensure consistency with the charge distribution assumed in Coulomb distortion calculation. Therefore, we need to make sure that the value of $F_{ch}$ used in our calculations is close to that used in analysis of PREX and CREX. To achieve this, we take $\sigma_{F_{ch}}=0.001$ which is much smaller than $\sigma_{F_{W}}$. And due to the remaining uncertainty in $F_{ch}$, we consider $F_{ch}-F_{W}$ in the likelihood for its better constraints than $F_{W}$.

In addition to the basic nuclear constraints, we impose additional PREX and CREX constraints separately and jointly, forming four sets of likelihoods: basic nuclei, +PREX, +CREX, and +PREX+CREX. We take the electric weak form factor difference $F_{ch}-F_{W}$ instead of the weak form factor $F_{W}$ as constraints of PREX and CREX experiments, as tabulated in Table 6. In principle, the weak form factor $F_{W}$ can be determined from the parity-violating asymmetry $A_{PV}$ measured in PREX and CREX experiments through Coulomb distortion calculations, assuming an electric charge distribution. Since the electric charge distribution has already been determined experimentally with great accuracy, both $F_{W}$ and $F_{ch}-F_{W}$ are determined experimentally from the measured $A_{PV}$ without any model dependence. However, theoretical models like RMF and Skyrme models, optimized with charge radius, do not have enough accuracy to ensure consistency with the charge distribution assumed in Coulomb distortion calculations. Therefore, we need to make sure that the value of $F_{ch}$ used in our calculations is close to that used in the analysis of PREX and CREX. To achieve this, we take $\sigma_{F_{ch}}=0.001$, which is much smaller than $\sigma_{F_{W}}$. Due to the remaining uncertainty in $F_{ch}$, we consider $F_{ch}-F_{W}$ in the likelihood for its better constraints compared to $F_{W}$.

Finally, we want to comment on additional constraints from giant monopole resonances (GMR) and dipole polarizability, which we do not consider explicitly in the likelihood. Instead, we use a tight prior on the compressibility $K$ within the empirical window $[210,250]$ MeV to account for GMR. Although we don’t include dipole polarizability explicitly, it is possible to consider dipole polarizability as an additional independent constraint on $S_{\mathrm{V}}$ Piekarewicz et al. (2012):

$$S_{\mathrm{V}}=\frac{L+\left[146\pm 1\right]\rm{~{}MeV}}{\left[6.11\pm 0.316\right]}$$ (39)

This is shown in Fig. 1 and applied in Fig. 2 in the letter.

Figure 4: The posterior distributions of model parameters, saturation properties and $\Delta F$ measured in PREX and CREX for Skyrme model (upper) and RMF (lower) model.

Figure 5: The Pearson correlation coefficient between model parameters, saturation properties and $\Delta F$ measured in PREX and CREX for Skyrme model (upper) and RMF (lower) model.

Figure 6: The posterior distributions of basic nuclear properties used in the likelihood as well as neutron skins of ⁴⁸Ca and ²⁰⁸Pb for Skyrme model (upper) and RMF (lower) model.

I.6 Posterior

With the prior and likelihood discussed above, we obtain the posterior distributions of model parameters, saturation properties, and $\Delta F$ measured in PREX and CREX, as shown in 4. To highlight the significant aspects of the posterior distribution, we also plot the Pearson correlation coefficient between any two parameters, as shown in Fig. 5. Some model parameters, such as $t_{0}$, $t_{3}$, and $\alpha=\gamma-2$, are tightly correlated due to the mapping between saturation properties used in our prior. Isoscalar and isovector Yukawa couplings are also tightly correlated due to the priors on binding energy $B$ and effective mass $M^{*}$. For both Skyrme and RMF models, the symmetry energy $S_{\mathrm{V}}$ and its slope $L$ are very sensitive parameters related to the form factor difference measured in PREX and CREX. From the posterior of Bayesian analysis, we obtain the mean and covariance of $S_{\mathrm{V}}$ and $L$ with +PREX+CREX constraints,

$$(\bar{S}_{\mathrm{V}},\bar{L}) = (34.1,~{}~{}40.3)_{Skyrme},~{}~{}~{}(31.4,~{}~{}49.1)_{RMF}$$ (40)

$$\sqrt{\mathrm{cov}} = \begin{pmatrix}6.5&9.1\\ 9.1&21.3\end{pmatrix}_{Skyrme},~{}~{}\begin{pmatrix}5.2&9.6\\ 9.6&27.4\end{pmatrix}_{RMF}$$

compared with the posterior with basic nuclei constraints only,

$$(\bar{S}_{\mathrm{V}},\bar{L}) = (32.8,~{}~{}66.8),~{}~{}~{}\sqrt{\mathrm{cov}}=\begin{pmatrix}5.9&11.7\\ 11.7&31.2\end{pmatrix}$$ (41)

.

Interestingly, despite the weak correlation between $b^{\prime}_{4}$ and $\Delta F^{\mathrm{Ca}48/\mathrm{Pb}208}$ in the Skyrme model, as observed in Reinhard et al. (2022b), $b^{\prime}_{4}$ can be quite sensitive to a linear combination of $\Delta F^{\mathrm{Ca}48}$ and $\Delta F^{\mathrm{Pb}208}$, such as $a^{\prime\prime}\Delta F^{\mathrm{Ca}48}+b^{\prime\prime}\Delta F^{\mathrm{Pb}208}$, significantly influencing the distribution of $\Delta F^{\mathrm{Ca}48},\Delta F^{\mathrm{Pb}208}$. The posterior for $b^{\prime}_{4}$ increases to $b^{\prime}4=0.310^{+0.24}{-0.26}$ with additional PREX and CREX constraints from $b^{\prime}_{4}=0.244^{+0.21}{-0.19}$ with basic nuclei constraints only. A similar trend is also observed in RMF models. The detailed impact of $b^{\prime}_{4}$ on the nucleon distribution is discussed in a separate section below.

Figure 6 shows the finite nuclei properties of ⁴⁸Ca, ⁹⁰Zr, and ²⁰⁸Pb of the Bayesian posterior samples with +PREX+CREX constraints. This distribution is consistent with the likelihood imposed in Table 6, with the exception that the charge posterior is more constrained than the 2% charge radius constraints we impose. This tighter constraint on charge radius is indeed due to the charge form factor being more constraining. We also noticed a much weaker $\Delta R_{np}$-$\Delta F$ correlation in ⁴⁸Ca than in ²⁰⁸Pb, since the momentum transfer of electron elastic scattering is much higher in CREX than in PREX due to the experimental setup. Therefore, CREX should be viewed more as a measurement of form factor difference than the neutron skin of ⁴⁸Ca.

Figure 7: Similar to Fig. 3 in the letter but colored by symmetry energy parameter $S_{V}$. Colorful lines corresponds to the minimum $\chi^{2}$ fitting of Eq.(5) in the letter with $S_{V}=[20,25,30,35,40,45]$ MeV.

Figure 8: Similar to Fig. 3 in the letter but colored by symmetry energy parameter$L$. Colorful lines corresponds to the minimum $\chi^{2}$ fitting of Eq. (6) in the letter with $L=[0,25,50,75,100,125,150]$ MeV.

Figure 9: Pearson covariance coefficient between neutron skin (form factor difference) of ⁴⁸Ca (solid) and ²⁰⁸Pb (dashed), and pressure at various baryon number densities. Blue and Red correspond to RMF models and Skyrme models.

I.7 Linear Correlation Between $S_{\mathrm{V}}$, $L$, $\Delta F^{\mathrm{Ca48}}$, and $\Delta F^{\mathrm{Pb208}}$

Based on all the samples used in the Bayesian analysis, we found that any three of the variables $S_{\mathrm{V}}$, $L$, $\Delta F^{\mathrm{Ca48}}$, and $\Delta F^{\mathrm{Pb208}}$ can form semi-linear relations with each other, resulting in four relations as shown in Eqs. (5-9) in the letter. These pocket formulas could help interpret future parity-violating asymmetry experiments. In this section, we explicitly show how these formulas perform for Skyrme and RMF models.

Figures 7-8 show the distributions of $(\Delta F^{\mathrm{Ca48}},\Delta F^{\mathrm{Pb208}},S_{\mathrm{V}})$ and $(\Delta F^{\mathrm{Ca48}},\Delta F^{\mathrm{Pb208}},L)$, where the linear correlations in Eqs. (5-6) in the letter are observed. The precision of the fitting formula in Eq. (5) remains within 10% (1-$\sigma$) for both Skyrme and RMF models, while Eq. (6) shows about a 10 MeV deviation for the Skyrme model and a 25 MeV deviation in the RMF model. An interesting trend is noted in Fig. 7 for both RMF and Skyrme models: the predicted points of ${\Delta F^{\mathrm{Ca48}},\Delta F^{\mathrm{Pb208}}}$ move collectively from left to right (towards the mean of joint PREX and CREX measurements) with increasing $S_{\mathrm{V}}$. This collective movement results from $S_{\mathrm{V}}$ being sensitive to the neutron skin of ²⁰⁸Pb while being insensitive to the neutron skin of ⁴⁸Ca in both RMF and Skyrme models. The difference in sensitivity of $\Delta F^{\mathrm{Pb208}}$ and $\Delta F^{\mathrm{Ca48}}$ to $S_{\mathrm{V}}$ is reflected by the slope coefficients $a$ and $b$ in Eq. (5) in the letter, where $b$ is much larger than $a$ in both Skyrme and RMF fittings.

In Fig. 8, we present the corresponding $L$ values of our sampled points of ${\Delta F^{\mathrm{Ca48}},\Delta F^{\mathrm{Pb208}}}$. As expected, $L$ is sensitive to both the neutron skin of ⁴⁸Ca and ²⁰⁸Pb, with the points moving collectively from lower left to upper right with increasing $L$. This movement highlights the challenge of finding parameterizations that satisfy both PREX and CREX measurements by only tuning parameters controlling the $L$ values. Notably, for the same $L$, Skyrme models predict slightly smaller $\Delta F$ due to the larger slope of the constant $L$ relation for $\Delta F$ compared to RMF models, as seen in Fig. 8. This can be understood from the Pearson correlation coefficient between $\Delta F$ and $L(n_{B})$ in Fig. 9. The Pearson correlation peaks around two-thirds of the saturation density for ²⁰⁸Pb, where the neutron skin correlates well with the form factor difference $\Delta F$ for ²⁰⁸Pb. The neutron skin is known to correlate best with the symmetry energy slope at two-thirds of the saturation density for the Skyrme model Brown (2000), a correlation confirmed by RMF models as well Typel and Brown (2001). Our Pearson correlation study for $\Delta F^{\mathrm{Ca48}}$ suggests that $\Delta F^{\mathrm{Ca48}}$ measured in CREX is most sensitive to the symmetry energy slope $L$ at higher densities: $n_{B}=0.12$ fm^-3 for the RMF model and $n_{B}=0.16$ fm^-3 for the Skyrme model, around saturation density. The Pearson correlation is equal between the two nuclei for the RMF model, while the Pearson correlation for ²⁰⁸Pb is much higher than for ⁴⁸Ca in Skyrme models, explaining the larger slope of the constant $L$ relation for $\Delta F$.

Figure 10: Vector densities of proton (lower) and neutron (upper) in ⁴⁸Ca, ⁹⁰Zr and ²⁰⁸Pb. Different colored band represent 18%, 50% and 82% percentile under various constraints.

Figure 11: Similar to Fig. 10. Different colored band represent RMF models with various $b^{\prime}_{4}$ defined in Eq.(37).

I.8 Nucleon density distribution of ⁴⁸Ca and ²⁰⁸Pb in RMF model

In this section, we delve into the proton and neutron density distributions of ⁴⁸Ca and ²⁰⁸Pb within RMF models, utilizing samples with four different likelihoods. We computed the baryon number density, denoted as $n(r)$, as a function of radius $r$ and assessed the 18%, 50%, and 82% percentiles of $n(r)$ across all radii.

Figure 10 illustrates the 18%, 50%, and 82% percentiles of proton and neutron densities in ⁴⁸Ca and ²⁰⁸Pb under various likelihood scenarios. The yellow band represents the basic nuclear constraints outlined in Table 6. Interestingly, this band significantly overlaps with the blue band, which incorporates additional constraints from PREX. This overlap is expected due to the abundance of RMF models in our prior with large $S_{\mathrm{V}}$ and $L$, coupled with the relatively larger statistical error of PREX compared to CREX. Similar overlaps are observed between the black and red bands. However, the imposition of CREX constraints notably alters the proton and neutron distributions, as evidenced by the discernible differences between the yellow and red bands (or black and blue bands). In ⁴⁸Ca, the neutron bulge around 2.5 fm increases by approximately 0.005 fm^-3 with the imposition of CREX constraints, while the distributions for $r<1$ fm remain unchanged. Consequently, the neutron distribution around the surface contracts inward by about 0.05 fm, resulting in smaller values for $\Delta R_{np}$ and $\Delta F$. Conversely, in ²⁰⁸Pb, the imposition of CREX constraints enhances neutron bulges at 1.5 fm and 4 fm, while also increasing neutron deficits and proton bulges at the center. Since the central region contributes minimally to the neutron skin, the overall neutron distribution shifts inward, albeit not as significantly as in ⁴⁸Ca.

As highlighted in our letter, the degeneracy (correlation) between $\Delta R_{np}$ (or $\Delta F$) of ⁴⁸Ca and ²⁰⁸Pb may be mitigated by breaking $S_{\mathbf{V}}$-$L$ relations or adjusting $b^{\prime}_{4}$. We categorized the sample into five bins based on $b^{\prime}_{4}$ and assessed the 18%, 50%, and 82% percentiles of $n(r)$ accordingly. Figure 11 elucidates how varying $b^{\prime}_{4}$ modulates the density distribution of protons and neutrons. For ⁴⁸Ca, increasing $b^{\prime}_{4}$ amplifies the neutron bulge and suppresses the proton bulge around 2.5 fm. Given that the bulge around 2.5 fm exerts the greatest influence on the overall nucleon distribution due to the $r^{2}$ factor in the volume weighting, the heightened isovector density in the bulge precipitates a notable decrease in $\Delta R{np}$ and $\Delta F$ for ⁴⁸Ca. The behavior of the intermediate nucleus ⁹⁰Zr appears to mirror that of ⁴⁸Ca, albeit with slightly smaller variations in magnitude. However, for ²⁰⁸Pb, increasing $b^{\prime}_{4}$ deepens the neutron deficit and enhances the proton bulge at the center, while also diminishing the isovector density at 3 fm. Although neutron excesses increase around 1.5 fm and 4 fm, the overall neutron excess remains relatively unchanged. Notably, $\Delta R_{np}$ and $\Delta F$ exhibit insensitivity or slight positive correlation with $b^{\prime}_{4}$. This can also been seen from Fig. 5.

Across all nuclei, a larger $b^{\prime}_{4}$ accentuates the isovector density oscillations throughout the nuclei. However, in ⁴⁸Ca, these oscillations manifest as a neutron excess in the bulk interior, while maintaining a relatively consistent average neutron excess in the interior of ²⁰⁸Pb. Since a substantial neutron excess in the nuclear interior correlates with a smaller neutron skin and a diminished form factor difference, adjusting $b^{\prime}_{4}$ may break the correlation between $\Delta R_{np}$ (or $\Delta F$) of ⁴⁸Ca and ²⁰⁸Pb.

In summary, the CREX experiment necessitates a small $\Delta R_{np}$ and $\Delta F$ for ⁴⁸Ca, while the PREX experiment demands a large $\Delta R_{np}$ and $\Delta F$ for ²⁰⁸Pb. Given the larger error bars in $\Delta F$ for the PREX experiment (three times standard deviation) compared to the CREX experiment, the joint PREX and CREX results are predominantly influenced by the CREX experiment, thereby shifting $\Delta R_{np}$ and $\Delta F$ for both ⁴⁸Ca and ²⁰⁸Pb towards smaller values, consistent with other Bayesian analysesZhang and Chen (2023); Salinas and Piekarewicz (2023). If future experiments, such as MREX, confirm the large $\Delta F^{Pb208}$ by narrowing down the statistical error, one potential solution involves introducing a significant isospin-dependent term of the spin-orbit interaction to amplify the variation of isovector density. Such variation impacts the average interior neutron excess differently between ⁴⁸Ca and ²⁰⁸Pb. In this study, we achieve this by adjusting $b^{\prime}_{4}$, although other interactions such as the isovector tensor interaction, may serve a similar purpose Salinas and Piekarewicz (2024). Accounting for multiple interactions that excite isovector density oscillations may offer insights into resolving the tension between PREX and CREX.