Improved determination of the oscillator parameters in nuclei

Our first method follows the conventional framework [2, 1], which employs experimental data on charge radii to constrain the oscillator parameter, while including the five corrective terms discussed in the previous section. The main contributor, $\braket{\bm{r}^{2}_{\pi}}_{pt}$ in Eq. (3) is evaluated using harmonic oscillator radial wave functions, namely

where the sum is taken over all occupied states of protons. The single-particle matrix element,

has been derived using the well-known virial theorem with $n$ and $l$ denoting the radial and orbital angular momentum quantum numbers, respectively. The proton occupation numbers $N_{nl}^{\pi}$ are fixed within the so-called equal-filling approximation. Substituting these expressions into Eq. (3), we obtain

where $\tilde{Z}=\sum_{nl}N_{nl}^{\pi}(2n+l+3/2)$ and $Z/\tilde{Z}\leq 1$. The factor $f=[1-3Z/(2A\tilde{Z})]$ accounts for the center of mass motion and has an effect of enlarging $b_{\pi}^{2}$ especially in the light mass region. The relativistic spin-orbit contribution ($\Delta_{ls}$) is neglected for simplicity.

The experimental charge radius data for 797 nuclei from $A=2$ (excluding ¹H) to $248$ have been taken from the latest compilation of Angeli and Marinova [3]. The updated data in Refs. [4, 5] are also considered. Our results for $b_{\pi}^{2}$ agree very well with those of Kirson as one can see from Fig. 1. It should be noted, however, that our calculation in this section includes both open-shell and closed-shell nuclei, therefore the obtained $b_{\pi}^{2}$ values are significantly scattered off their trend line in the regions where deformation and correlation are dominant. On the other hand, the large scattered points near the coordinate origin might be due to the breakdown of the mean field theory towards the $A\to 1$ limit. It is also seen that our calculation slightly underestimates the result of Blomqvist and Molinari [1] in the light-mass region, but overestimates it for nuclei with a mass number starting around $A=70$. Within this method, we obtain the following expression for protons:

with $\chi^{2}/\nu=48.869$. The number of degrees of freedom, $\nu$ is obtained by subtracting the number of model parameters (4 parameters) from the sample size (1068 nuclei). Beyond the conventional form, which includes only a constant and a term proportional to $A^{\frac{1}{3}}$, we introduce linear and quadratic terms in $I$ to distinguish between isobars where $I=(N-Z)/A$. Without this extension, our fit would yield $\chi^{2}/\nu=54.197$. A large value of $\chi^{2}/\nu$ ($\chi^{2}/\nu\gg 1$) indicates that the uncertainties in the data sample are effectively smaller than the distances of individual data points from their trend line. To account for this inconsistency, we scale the obtained uncertainties on the model parameters with $\sqrt{\chi^{2}/\nu}$ as suggested by the Particle Data Group [6]. The small discrepancy between our results and those of Kirson arise mainly from the difference in the determination of $N_{nl}^{\pi}$ (see Ref. [2] for details).

In order to extract the oscillator length parameter for neutrons, we follow Kirson [2], using the relation

where $\Delta r_{np}$ is the neutron skin thickness. The squared radius of point-like neutron distribution $\braket{\bm{r}^{2}_{\pi}}_{pt}$ can be expressed as Eq. (6) and Eq. (7) for protons except that the proton occupation number $N_{nl}^{\pi}$ and the parameter $b_{\pi}$ must be replaced with those of neutrons. Inserting the expression of $\braket{\bm{r}^{2}_{\nu}}_{pt}$ and of $\braket{\bm{r}^{2}_{\pi}}_{pt}$ into Eq. (10), we obtain

where $\tilde{N}=\sum_{nl}N_{nl}^{\nu}(2n+l+3/2)$ with the sum running over all occupied states of neutrons. Therefore $N/\tilde{N}\leq 1$. The neutron occupation number $N_{nl}^{\nu}$ is determined with the same method as $N_{nl}^{\pi}$. The first term on the right-hand-side (r.h.s) of Eq. (11) is induced by the neutron skin thickness whereas the following term is influenced by the difference between $N$ and $Z$. Unlike the work of Kirson which employed the empirical formula of $\Delta r_{np}$ extracted from antiproton interaction with nuclei [7] or hadronic atom and hadron scattering data [8], we calculate this quantity within the Hartree-Fock-Bogoliubov method using effective Skyrme interaction. We include all even-even nuclei whose charge radius data are available from the above-mentioned compilations. A least squares fit to the results of these calculations yields

with $\chi^{2}/\nu=1.616$. Note that we also exclude cases with $A<10$ from our mean field calculations; consequently, our sample size for fitting Eq. (12) is reduced to 317. Then, the number of degrees of freedom in this process becomes $317-3=314$. This small $\chi^{2}/\nu$ value indicates that the model’s errors are consistent with the uncertainties in the data sample. Several well-established Skyrme force parameterizations are considered, namely SLY4/SLY5 [9], SKM* [10], SGII [11], SII/SIII/SIV [12], and UNEDF0/UNEDF1/UNEDF2 [13, 14, 15]. Besides the spread of the results among the selected Skyrme parameterizations, we account for an uncertainty of $\pm 0.05$ fm for $\Delta r_{np}$ propagated from the symmetry energy slope ($L$) for ²⁰⁸Pb, which was estimated to be $64\pm 39$ MeV [16]. The pairing correlation and deformation are found to be significant for certain individual cases; however, their impact on the global trend of neutron skin thickness is generally negligible. By substituting the expression (12) into Eq. (11) and then performing a least squares fit with a model similar to Eq. (9) for $b_{\pi}^{2}$, we obtain

The resulting $\chi^{2}/\nu$ value is 0.861. We notice a somewhat diminished uncertainty in $b_{\nu}^{2}$ compared to $b_{\pi}^{2}$. This reduction is attributed to the factors $N/\tilde{N}$ and $(N/\tilde{N})\times(\tilde{Z}/Z)$ in Eq. (11). Given that $N/\tilde{N}<1$ and $Z/\tilde{Z}<1$ as noted above, the factor $(N/\tilde{N})\times(\tilde{Z}/Z)$ is also effectively less than 1, owing to the predominance of neutron-rich species among the majority of nuclei. Again, our result for $b_{\nu}^{2}$ is not very far from that of Kirson in Eq. (5), even though we use a different method for the determination of neutron skin thicknesses.

Despite the unrealistic nature of the oscillator potential itself, the single-particle Hamiltonian (1) does not account for the repulsive Coulomb force among protons and the symmetry effect induced by the difference between neutron and proton numbers. This means that replacing the $\bm{r}^{2}$ operator in Eq. (6) by unity, we will get a probability density distribution of uncharged particles of an $N=Z$ nucleus. To go beyond this conventional picture, we add the following Coulomb term derived from an assumption of a uniformly charged sphere, to the oscillator Hamiltonian:

where $e$ is the elementary charge and $t_{z}$ is the isospin projection of the nucleon with the convention of $t_{z}=\frac{1}{2}$ for neutrons and $-\frac{1}{2}$ for protons. The following convention for isospin must be applied when $t_{z}$ is used as an index: $t_{z}=\pi$ (for protons) or $\nu$ (for neutrons). In fact, Eq. (14) is correct only for $r<R_{C}$. However, the significance of this is negligible because the potential, $m\omega^{2}\bm{r}^{2}/2$ diverges rapidly as distances increase. For simplicity, the Coulomb exchange term is neglected. The Coulomb radius, $R_{C}$ is parameterized in literature as $R_{C}\approx 1.26~{}\text{fm}\times A^{\frac{1}{3}}$ [17]. For the present work, we fix this parameter with the measured mean squared charge radii through the following formula [18],

where $\braket{\bm{r}^{2}_{\pi}}_{pt}$ can be converted into $\braket{\bm{r}^{2}}_{ch}$ using the relation Eq. (3).

The first term on the r.h.s. of Eq. (14) does not depend on $\bm{r}$, so it has no effect on the eigenfunctions. Meanwhile, the second term is proportional to $\bm{r}^{2}$, so it can be merged with the oscillator potential. Similarly, the symmetry effect in an $N\neq Z$ system can be accounted for by introducing a factor $[1-2t_{z}\kappa(N-Z)/A]$ to the oscillator potential, such that

where $\kappa$ is a free parameter characterizing the strength of the symmetry term. In principle, the spin-orbit term of Thomas [19] can also be added in a straightforward fashion. However, within the oscillator potential the radial form of this spin-orbit term is reduced to a constant, so it has no effect on the radial component of wave functions, regardless of energy splits between spin-up and spin-down states.

Therefore, the inclusion of the Coulomb and symmetry terms simply results in an effective oscillator frequency, namely

with $b_{t_{z}}^{2}=\hbar/(m\omega_{t_{z}})$ and $b^{2}=\hbar/(m\omega)$. Since the Coulomb and symmetry terms are separated out, the parameter $b$ or $\omega$ in Eq. (17) should be isospin-invariant. Subsequently, the following relation can be derived:

These effective oscillator parameters correspond to those extracted from the experimental data within the method I in the previous subsection. Therefore, with a known $b_{\pi}$ and $b_{\nu}$, Eq. (18) may be regarded as an alternative tool for the determination of the symmetry parameter $\kappa$. Nevertheless, we found that the $\kappa$ values extracted from Eq. (18), i.e. $\kappa=0.423(441)$, vary strongly from nucleus to nucleus and their average underestimates the values used in literature considerably. For the present calculations, we adopt $\kappa=0.75(12)$ which is taken from the global fits of the Woods-Saxon parameter sets [17, 20, 21, 22, 23]. The resulting isospin-invariant $b^{2}$ values are represented by the expression

with $\chi^{2}/\nu=6.853$ and $\nu=1064$. Note that $b$ would correspond to the length parameter obtained for neutrons in subsection II.1 if only nuclei with $N=Z$ (excluding the symmetry effect) are considered. Additionally, within Eq. (18), $b^{2}$ can be extracted from either $b^{2}_{\pi}$ or $b^{2}_{\nu}$. For this study, we adopt the average of these two values and treat their difference as an uncertainty source.

As a general feature, the inclusion of the Coulomb term leads to a reduction in the effective oscillator frequency or an increase in the effective length parameter, because it reduces the proton binding energies. The symmetry term has an opposite effect on protons in neutron-rich nuclei, contributing about 50 % to the squared length parameter, as depicted in Fig. 1, compared to the Coulomb contribution. Consequently, a substantial cancellation occurs between the Coulomb and symmetry contributions almost everywhere throughout the mass range under consideration. The symmetry effect on neutrons in neutron-rich nuclei is repulsive, thus leading to a larger effective length parameter for neutrons. The impact of the Coulomb and symmetry terms on the effective length parameter for protons is illustrated in Fig. 1.

These results suggest that a distinction of the effective oscillator parameters between protons and neutrons, $b_{\pi}\neq b_{\nu}$, should be made for a calculation in which isospin-symmetry breaking is taken into account. This distinction would be particularly important for the shell model description of isospin mixing, where the configuration space is extremely limited. Conversely, a single oscillator parameter, $b$ should be used instead when isospin symmetry is assumed, such as in the conventional shell model calculations [24, 25, 26].

One may further argue that, in addition to its unrealistic form in the space coordinate, the oscillator Hamiltonian (1) also misses an appropriate spin-orbit term. The presence of a spin-orbit term breaks the degeneracy of the oscillator states, leading to a difference in single-particle configurations, which then yields different expectation values for the mean squared radii. Furthermore, the experimental data of charge radii may contain a considerable contribution of deformation and other effects beyond the spherical mean field description. We therefore apply in this subsection an alternative method for determining the oscillator length parameter without the utilization of the experimental data. The idea is to directly fit the harmonic oscillator radial wave functions to the eigenfunctions of spherical SHF mean field by varying the length parameter while the Coulomb and nuclear charge-dependent forces are turned off. It is expected that the influence of the spin-orbit coupling, as well as spurious isospin mixing and beyond spherical mean field effects, can be avoided by considering only $s_{1/2}$ states of the doubly magic $N=Z$ nuclei, namely ¹⁶O, ⁴⁰Ca, ⁵⁶Ni and ¹⁰⁰Sn. Although the sample size of this method would be too small for a prediction of the trend of the squared length parameter, it may provide a meaningful test of the previous methods through a point-by-point comparison. The fit of this method is carried out by maximizing the overlap integral defined below

where $\psi_{nl}^{t_{z}}(r)$ and $\phi_{nlj}^{t_{z}}(r)$ are, respectively, the harmonic oscillator and SHF radial wave functions. Note that the indices $l$ and $j$ in Eq. (20) can be omitted if only $s_{1/2}$ states are considered. Obviously, $\Omega_{nlj}^{t_{z}}$ reduces to the normalization integral if $\psi_{nl}(r)$ and $\phi_{nlj}^{t_{z}}(r)$ are perfectly coincided.

For the SHF calculations, we employ the same sets of Skyrme parameterizations as considered in subsection II.1. Since the Coulomb repulsion is not included, the obtained oscillator length parameters must be compared with the isospin-invariant one given in Eq. (19) or the one obtained for neutrons within Method I. Note that these correspondences are valid only for self-conjugate nuclei, where the symmetry effect is absent. The results for the four above-mentioned doubly magic $N=Z$ nuclei are listed in Table 1. For comparison, some other $N\neq Z$ species including closed-shell and closed-subshell nuclei are also listed. Unfortunately, it is found that the results are highly sensitive to the fine details of the SHF mean field due to the missing of the centrifugal barrier. Furthermore, the obtained length parameter values for a given doubly magic $N=Z$ nucleus vary remarkably depending on the radial quantum number, even when solely focusing on $s_{1/2}$ states. In principle, the state that generates the largest overlap integral, $\Omega_{nlj}$ should be selected, typically the $0s_{1/2}$ state, as evidenced by Table 1. In particular, this method consistently yields $b$ values that tend to exceed those obtained with method I, except for light nuclei, typically below ⁴⁰Ca. Generally, the largest $b$ value is obtained for the lowest state, specifically the $0s_{1/2}$, followed by a gradual decrease with increasing energy or radial quantum number. This observation suggests that, considering specific higher single particle energy occupied states with $l\neq 0$ for particular nuclei could potentially enhance the agreement with method I, despite the lack of clear theoretical justification. In our opinion, the dependence on $n$, or more generally, on single-particle orbits, may be attributed to the weakly bound effect, which cannot be seen within the harmonic oscillator potential which rises to $\infty$ rapidly with increasing of the spatial distance. In addition, the eigenfunctions of $s_{1/2}$ states may still be influenced indirectly by the spin-orbit potential of the other occupied states with $l\neq 0$ due to the self-consistency of the SHF equation.

In order to gain a better insight, we perform a further analysis with the inclusion of open-shell nuclei. We employ the standard SHF procedure with the presence of the Coulomb interaction, for which the Coulomb exchange contribution is evaluated using the Slater approximation. The obtained results are given in Fig. 2. The $b_{\nu}$ values illustrated on the top-left panel of Fig. 2 are obtained by fitting to neutron eigenfunctions whereas the $b_{\pi}$ values on the top-right panel of the same figure are obtained by fitting to proton eigenfunctions. The corresponding values of the overlap integral are shown on the bottom panels of the same Figure. It is interesting to remark that the fitted values of $b_{\nu}$ and $b_{\pi}$ in Fig. 2 have a similar pattern except that $b_{\pi}$ is slightly larger than $b_{\nu}$ for the effects discussed in subsection II.2. In general, a lower SHF state yields a larger length parameter value, as noted for closed-shell nuclei in the above paragraph. We notice that a notable dependence on $l$ persists in the absence of the spin-orbit term, indicating a distinct impact of the centrifugal barrier on both the SHF and harmonic oscillator potential. Despite the dependence on $n$ and $l$, it is seen from Fig. 2 that only fitting to the last occupied states achieves reasonable agreement with the results obtained from method I for both neutrons and protons, irrespective of the weakly-bound effect. However, the corresponding overlap integral values do not consistently exhibit a similar pattern. On average, though, these values for the last occupied states tend to approach unity more closely than those for lower-energy states. With the inclusion of the Coulomb interaction and the consideration of opened-shell nuclei, the best-fit length parameter for last occupied states is represented by

for protons, and by

for neutrons. The corresponding $\chi^{2}/\nu$ values are 151.12 and 64.47, respectively. The isospin-invariant oscillator parameter can, in principle, be extracted from these results following the method described in subsection II.2. The resulting expression is,

The parameterizations of the oscillator parameter determined by the present method differ considerably from those discussed in the previous subsections, in particular for neutrons, and the isospin-invariant version. The substantial uncertainties on the coefficients in Eq. (21), Eq. (22), and Eq. (23) reflect a large scatter in the data samples. Nevertheless, these parameterizations are much improved compared to those obtained by fitting to $s_{1/2}$ states alone as listed in Table 1. This improvement occurs despite the ambiguity in the fit quality, which is more pronounced around the Fermi level.