Realistic evaluation of Coulomb potential in spherical nuclei and test of the traditional approach

The charge form factors for some stable nuclei were accurately measured by electron scattering. These data are usually analyzed within two different model-independent approaches to extract numerical values for the charge density as a function of $r$ namely the Fourier-Bessel (FB) [27] and the sum-of-Gaussians (SOG) [28] expansion. In the former approach, the charge density is expanded in terms of the spherical Bessel function of order zero ($j_{0}$) namely,

where $a_{\nu}$ are the expansion coefficients and $R_{cut}$ is a cut-off radius beyond which the charge density is sufficiently small and is equated to zero. The first $n_{max}$ coefficients of this series expansion are obtained directly from the experimental data [20]. Here $n_{max}=R_{cut}q_{max}/\pi$ where $q_{max}$ is the maximum momentum transfer up to which the charge form factor data are determined. For the normalization, the integral of $\rho_{\text{FB}}(r)$ over all spaces must be equal to the total nuclear charge $+Ze$.

In the sum-of-Gaussians approach, the charge density is expressed as

where $y_{i}(r)=(r+R_{i})/\gamma$ and $x_{i}(r)=(r-R_{i})/\gamma$. The expansion coefficients $A_{i}$ are given by

where $\gamma$ is the width of the Gaussians and is related to the root-mean-square radius (RMS) through $R_{g}=\gamma\sqrt{3/2}$. According to Sick [28], $\gamma$ is chosen equal to the smallest width of the peaks in the nuclear radial wave functions calculated using the HF method. The author reported that the $\gamma$ values extracted from the harmonic oscillator and WS radial wave functions yield almost identical results. The charge fraction $Q_{i}$ must be normalized such that $\sum_{i}Q_{i}=1$.

In practice, only first few terms are included within the sum in Eq. (14) and Eq. (15). For example, a truncation with $n_{max}\leq 17$ and $m_{max}\leq 12$ was imposed for the analyses carried out in Ref. [20]. The parameters of $\rho_{\text{FB}}(r)$ and $\rho_{\text{SOG}}(r)$ deduced from experiments are given in the data compilation [20].

We remark that the charge density obtained from the FB approach has an undesired property–it contains an oscillatory component and, sometimes, has a negative value near the cut-off radius $R_{cut}$ as illustrated in the right panel of Fig. 1. Although these oscillations are extremely small in magnitude, they are greatly amplified by the first and second derivatives, then leading to a large error, if the GGA is employed for treating the Coulomb exchange term. Oscillations in the asymptotics of charge density obtained from the SOG approach are also visible for some nuclei, however their magnitude are generally much smaller (larger width) as shown in the left panel of Fig. 1. This unphysical property originates from the incompleteness in the expansions and too small cut-off radius in the FB analysis. Fundamentally, the tail of charge density should decrease monotonically since the asymptotic radial wave functions decay exponentially.

In order to avoid this problem, we replace both FB and SOG data from a point just before the oscillations occur (denoted as $R_{x}$) towards infinity with the 2pF function (see subsection IV.3 for the definition of 2pF). We select for both data sets $R_{x}=4.5$ fm for $A\leq 58$ and $8.5$ fm for heavier nuclei. Instead of performing a global fit, the 2pF function parameters are here determined by matching it as well as its first and second derivatives with those of the data at $R_{x}$. Therefore, the resulting 2pF function might not be optimal for the whole available data in the domain of $[R_{x},+\infty]$. Nevertheless, it is the only way to ensure the continuity of the charge density at $R_{x}$. This replacement would not cause a serious error since the charge density in this domain is small.

For a sphere of radius $R_{C}$ containing a total charge of $+Ze$ uniformly distributed throughout its volume, the charge density is written as

and the normalization condition implies that

Here we follow Ref. [29] wherein $\rho_{\text{Unif}}(r)$ is regarded as a pointed proton distribution and is normalized to $Z$ instead of $+Ze$. Therefore, its RMS radius (denoted as $R_{\text{Unif}}$) can be calculated analytically, namely

In this uniform distribution, the corrections for the finite size, Darwin-Foldy term and COM motion must be introduced in order to convert $R_{\text{Unif}}$ into the charge radius. Eq. (19) provides an experimental constraint for $R_{C}$, at least, for the cases where the experimental data are available. Otherwise it can be parameterized as usual, namely $R_{C}=r_{C}(A-1)^{\frac{1}{3}}$ where $r_{C}\approx 1.26$ fm [9]. The validity of different methods for determining $R_{C}$ can also be tested within the framework of the present study. More discussions on this point are given in the next section.

Because of its simplicity, the uniform charge distribution is widely used in nuclear physics, especially within the nuclear optical model. However, one should note that it is an assumption of classical electromagnetism and has no quantum mechanical equivalent since wave functions, which are the building-blocks of charge density, are required to be continuous in coordinate space. Furthermore, the two aforementioned approaches to the Coulomb exchange term (Slater and GGA) are evidently inapplicable for this distribution because of its discontinuity and nondifferentiability at the surface of a nucleus.

We also consider a realistic phenomenological model for the charge density distribution namely the two-parameter Fermi function,

Unlike the previous models, the 2pF is continuous and differentiable, and moreover it decays monotonically towards large distances. Although the 2pF function looks similar to the $f_{i}(r)$ functions of the WS potential in Eq. (4), in general $c$ and $z$ are smaller than the length and the surface diffuseness parameter of the WS potential. It is also evident from the Skyrme HF theory that the mean field potential is not a linear function of densities. Thus it is not necessary that their geometrical characterizing parameters a the same with the WS potential parameters. More details on this point are given in Ref. [30].

As its name indicates, the 2pF is determined by two parameters because $\tilde{\rho}_{0}$ is obtained via the normalization condition,

where $F_{2}(c/z)$ is the second order Fermi integral (see Appendix C of Ref. [16]). Similar to the previous subsection, we employ the convention that $\rho_{2pF}(r)$ is normalized to $Z$, therefore the parameter $c$ can be extracted from the charge radius by solving the following equation [16]

As before, the last three terms account for the finite size of proton, Darwin-Foldy term and center-of-mass motion, respectively. $F_{4}(c/z)$ is the fourth order Fermi integral (see Appendix C of Ref. [16]).

Recently, Horiuchi [31] proposed a new method for determining the surface thickness of charge density. Implementing the Taylor expansion of $\rho_{2pF}(r)$ at $r=c$ and retaining up to the first-order term, he derived the following relation,

where $\rho_{2pF}^{\prime}(c)$ denotes the first derivative of $\rho_{2pF}(r)$ at $r=c$. By matching $\rho_{2pF}$ on the right-hand-side of Eq. (23) with the charge density constructed from eigenfunctions of the WS potential whose Coulomb term is, in turn, a function of $\rho_{2pF}$, Eq. (23) can be solved in a self-consistent manner. This offers an alternative means for constraining $z$ when experimental data are unavailable.

Although this distribution is nicely representative for the diffuseness at the nuclear surface region, it is not able to describe the oscillations of the charge density observed in the nucleus interior. It has been shown that this difficulty can be overcome by extending Eq. (20) to the three-parameter Fermi (3pF) function [32]. However, the 3pF is not a good choice for a general application because of lack of experimental data for constraining the third parameter.

Within a microscopic nuclear structure model, the charge density is basically decomposed into three components

where $\rho_{ch}^{p}(r)$/$\rho_{ch}^{n}(r)$ come from the finite charge distribution of the proton/neutron folded with the point-like proton/neutron density, and $\rho_{ch}^{ls}(r)$ is the relativistic electromagnetic correction which depends on the spin-orbit coupling. However, the shape of the charge-density distribution is mainly determined by the shape of the point-like proton density. The contributions $\rho_{ch}^{n}(r)$ and $\rho_{ch}^{ls}(r)$ were found to be negligible [12]. Furthermore, they tend to cancel out each other, thus will not be considered here.

By definition $\rho_{ch}^{p}(r)$ is given by

where $\boldsymbol{r}$ is the position vector in $\mathbb{R}^{3}$. The effective electromagnetic form factor $G_{p}$ is taken as a sum of three Gaussians as described in Ref. [18]. The point-like proton density $\rho_{p}(r)$ can be defined in terms of proton radial wave functions $R_{\alpha}(r)$ namely

where $\alpha$ stands for the spherical quantum numbers $nlj$ and the sum runs over all occupied states. The proton occupation number $n_{\alpha}$ is obtained with the so-called equal-filling approximation. Therefore, for a closed-shell configuration, $n_{\alpha}=(2j+1)$ for the occupied orbits and 0 for the unoccupied orbits.

The radial wave functions $R_{\alpha}(r)$ are taken as the eigenfunctions of the Skyrme HF mean field. For the Skyrme Hamiltonian, we select the SLY5 parameter set [33] which is invariant under rotation in the isospin space, while the Coulomb exchange term is treated using the Slater approximation. The charge-symmetry and charge-independence breaking forces [34] are neglected. We use the HFBRAD program [7] to solve the spherical Skyrme HF equation with which the Coulomb terms are evaluated using the point-like proton density. The finite size effect is corrected for by using Eq. (25) within an external program after the HF variation is terminated.

Before performing a comparative study of Coulomb potential, it is necessary to figure out the impact of the differences between the FB and SOG data for charge density on proton single-particle energies. For simplicity, we neglect the uncertainties on the coefficients of FB and SOG expansions even though they are available from the data compilations. Noted that the quantity which will be considered throughout this section is the proton single-particle energies. The other observables such as the RMS radii and neutron skin thicknesses were found to be much less sensitive to small variation in the Coulomb terms.

Furthermore, instead of a point-by-point comparison, we define the following radial mismatch factor,

as an effective measure for the charge density differences, where $\Omega$ is the overlap integral between the FB and SOG data:

whereas $\bar{\Omega}$ is the integral of their average squared:

The resulting $\Lambda$ values in % are 3.441$\times 10^{-3}$ (¹⁶O), 5.659$\times 10^{-3}$ (²⁸Si), 0.209 (³²S), 2.601$\times 10^{-2}$ (⁴⁰Ca), 4.126$\times 10^{-2}$ (⁴⁸Ca), 3.866$\times 10^{-3}$ (⁵⁸Ni), 3.641$\times 10^{-4}$ (²⁰⁵Tl), 5.048$\times 10^{-4}$ (²⁰⁶Pb) and 1.459$\times 10^{-2}$ (²⁰⁸Pb).

We recall that the exact Coulomb exchange term is nonlocal in coordinate spaces and cannot be written as a function of charge density, it is thus inappropriate for the present study which replaces the self-consistent HF mean field with the phenomenological WS potential. In order to demonstrate whether or not the selection of Coulomb exchange functional is matter to the impact of the FB-SOG charge density differences on proton single-particle energies, we consider three different approximations for Coulomb potential: one of them consists purely of the direct term whereas the other two include also an exchange term employing either the Slater approximation or the GGA. Then together with the nuclear component described in Eq. (3), the single-particle energies and wave functions can be obtained by solving the radial Schrödinger’s equation. The splits in proton energies due to the differences in the charge density data are listed in Table 1. As a convention, the negative sign indicates that the FB data yields a lower proton energy level relative to that yielded by the SOG data i.e., $\Delta E^{\rm FB}_{\rm SOG}=E_{\rm FB}-E_{\rm SOG}$. The meaning of the positive sign is opposite. For completeness, the averaged proton energies between the FB and SOG results calculated using the GGA functional [denoted as $\bar{E}^{\rm FB}_{\rm SOG}=(E_{\rm FB}+E_{\rm SOG})/2$] are graphically illustrated in Fig. 2. It is clearly observed from Table 1 that the splits in proton energies induced by the charge density differences are insensitive to the Coulomb exchange term. This witnesses that any existing Coulomb functionals can be employed for our test of the charge density models which will be discussed in the following subsection.

We remark that although the amount of energy splits is mainly determined by the size of $\Lambda$, it can be greatly enhanced by the shell-structure effect. As a general feature, the solution of the radial Schrödinger equation for protons will be particularly sensitive to the Coulomb terms when the last occupied orbit is highly filled and at the same time, has a low centrifugal barrier (low orbital angular momentum). In the remainder of this paper, this phenomenon will be referred to as the “weakly bound effect”. This is the main reason behind the large energy split of about $-500$ keV for ³²S. We notice that this number is comparable to the energy splits due to the use of the Slater approximation within the Skyrme HF framework studied in Ref. [25]. Therefore, a special attention must be paid to the inclusion of this nucleus for our analyses in the following subsection. The impact for the two Calcium (${}^{40}\rm{Ca}$,${}^{48}\rm{Ca}$) isotopes is also quite large namely about $-100$ keV, while the impact for the other cases is less than 50 keV in magnitude.

Because of its discontinuity and nondifferentiability, the uniform distribution is inappropriate to be used as an input for calculating the Coulomb exchange term within any existing approximations. In order to evaluate the Coulomb potential for this hypothetical distribution as well as the other charge density models listed in section IV on an equal footing, the Coulomb exchange term will be omitted for our study in this subsection. It was shown in the previous subsection that the absence of the Coulomb exchange term does not significantly affect the energy splits induced by the differences in the charge density data ($\Delta E^{\rm FB}_{\rm SOG}$), even for ³²S where $\Lambda$ is as large as 0.209 %.

Here we use the averaged proton energy between FB and SOG values as a reference (denoted as $\bar{E}_{R}$) relative to which the result obtained for a given charge density model is compared. As an exception, only the FB data are used for ²⁰⁹Bi because a SOG analysis has not been carried out for this nucleus. The proton energies evaluated with the uniform distribution, the 2pF function and the microscopic Skyrme-HF charge density are denoted as $E_{\rm Unif}$, $E_{2pF}$ and $E_{HF}$¹¹1$E_{HF}$ denotes the proton single-particle energies evaluated with the WS potential as a nuclear component and the charge density obtained from the microscopic Skyrme HF calculation as an input for the Coulomb potential. It should not be confused with the eigenvalues of the HF mean field., respectively. Our results for the proton energy differences including $\Delta E^{\rm Unif}_{R}=E_{\rm Unif}-\bar{E}_{R}$, $\Delta E^{2pF}_{R}=E_{2pF}-\bar{E}_{R}$ and $\Delta E^{HF}_{R}=E_{HF}-\bar{E}_{R}$ are given in Table 2.

We remark that the energy differences obtained in this subsection for all charge-density models are generally small in magnitude. For most cases, their magnitude is only of a few tens keV larger than those induced by the experimental uncertainties obtained in the previous subsection ($\Delta E^{\rm FB}_{\rm SOG}$). The values obtained for ³²S are somewhat larger compared with the other nuclei and positive for all bound orbitals, because of large uncertainties in the model-independent data for charge density and an additional enhancement due to the weakly bound effect. It is also seen that all the selected charge density models provide a similar accuracy for proton single-particle energies. In particular, the microscopic Skyrme HF model for charge density works adequately regardless of the various deficiencies related to the isospin-symmetry breaking discussed in Ref. [12]. We notice that, in general the calculations using the uniform distribution or the 2pF function are strongly parameter-dependent. For example, if the Coulomb radius is parameterized as $R_{C}\approx 1.26~{}\text{fm}\times(A-1)^{\frac{1}{3}}$ as usual, considerably larger proton single-particle energy differences are obtained in many cases. Similar problem was also observed within the 2pF function. Therefore, it is necessary to constrain their parameters using the relevant experimental data if these distribution functions are selected as a charge-density model.

The present subsection is concerned with the compensation for the omission of Coulomb exchange term which is unavoidable when using the uniform distribution for charge densities. We are interested specifically in the method employed in Eq. (1) by excluding the last proton contribution to Coulomb direct term. In order to verify this traditional idea of self-interaction correction, it is instructive to perform an analytic analysis of its impact before discussing the numerical results.

Since the Coulomb direct term (9) is linear in charge density, it can be written as

where $\rho_{ch}^{Z}$ is the total charge density which can be decomposed as

with $\rho_{ch}^{1}$ being the last proton contribution and $\rho_{ch}^{Z-1}$ the contribution of the remaining $Z-1$ protons. In order to gain a better insight into the functional structure of the Coulomb direct term, we further assume that $\rho_{ch}^{Z}$ and $\rho_{ch}^{Z-1}$ have identical radial form and they differ from each other only by the normalization condition, namely $\rho_{ch}^{Z}$ is normalized to $Z$ whereas $\rho_{ch}^{Z-1}$ is normalized to $Z-1$ (this is equivalent to the assumption that each proton has an equal contribution to the total charge density). Subsequently, the following relations can be derived

Because of the linearity, the corresponding Coulomb direct terms can be evaluated as

Based on this assumption, the self-interaction correction discussed above can be done simply by replacing $V_{dir}[\rho_{ch}^{Z}]$ with $V_{dir}[\rho_{ch}^{Z}]\times(Z-1)/Z$. Note that the Coulomb potential in Eq. (1) fully satisfies these properties with $\rho_{ch}^{Z}$ being the uniform distribution. It is interesting to remark that the factor $(Z-1)/Z$ goes to 1 as $Z$ goes to infinity, indicating that this correction has a greatest effect in light nuclei. We found that this replacement leads to lower proton energy levels, namely by about 500 keV in ¹⁶O and 200 keV in ²⁰⁸Pb, relative to the levels obtained with the pure Coulomb direct term $V_{dir}[\rho_{ch}^{Z}]$.

Throughout this subsection, the Coulomb functional obtained in this way will be referred to as “traditional functional”. In order to make a numerical test this traditional approach, it is useful to decompose the proton energy difference for a given orbital as

where $\Delta E^{T}_{F}$ is the single-particle energy difference between the traditional functional and the exact Fock term. On the right-hand-side, $\Delta E^{T}_{G}$, $\Delta E^{G}_{S}$ and $\Delta E^{S}_{F}$ are, respectively, the single-particle energy difference of the traditional treatment relative to GGA, GGA relative to Slater approximation, and Slater relative to the exact Fock term. Explicitly, $\Delta E^{T}_{F}=E^{T}-E^{F}$, $\Delta E^{T}_{G}=E^{T}-E^{G}$, $\Delta E^{G}_{S}=E^{G}-E^{S}$ and $\Delta E^{S}_{F}=E^{S}-E^{F}$. As the uniform distribution is inappropriate for the evaluation of Coulomb exchange term, the 2pF function will be selected instead for this test.

We notice that no exact treatment of Coulomb exchange term is performed in the present work, the data of $\Delta E^{S}_{F}$ for ¹⁶O, ⁴⁰Ca, ⁴⁸Ni and ²⁰⁸Pb are taken from Ref. [25]. Since the mass dependence of $\Delta E^{S}_{F}$ is not so strong, it is reasonable to use the values obtained for ⁴⁰Ca for ²⁸Si, ³²S and ⁴⁸Ca which were not considered in Ref. [25] because they reside in the neighbor in the nuclear chart. For the same reason, we use the $\Delta E^{S}_{F}$ values obtained for ²⁰⁸Pb for the remaining nuclei. One may argue here that the proton energy differences may depend on the nuclear component as well as the method for solving the Schrödinger equation, or they may vary significantly when transferring from a self-consistent to a phenomenological mean field. We have checked such dependence by looking at the term $\Delta E^{G}_{S}$, we found that its values obtained for ²⁰⁸Pb using WS potential are in the range between -1 and 19 keV, which are in remarkably good agreement with those calculated within the self-consistent Skyrme HF method.

Our numerical results for $\Delta E^{T}_{F}$ as well as for $\Delta E^{T}_{G}$ and $\Delta E^{G}_{S}$ are given in Table 3. It is seen that in light nuclei (around $Z=8$) negative $\Delta E^{T}_{G}$ values are obtained because the traditional self-interaction correction is stronger than the GGA exchange term. After $Z=8$, $\Delta E^{T}_{G}$ increases gradually with the atomic number and finally reaches a saturation at around $Z=80$. The saturated value of $\Delta E^{T}_{G}$ is about 300 keV. We notice also that the functional-driven energy differences are insensitive to the weakly bound effect, because the values obtained for ³²S do not differ significantly from those of the neighboring nuclei as observed in the two previous subsections. Combining these values with those of $\Delta E^{S}_{F}$ from the above-cited self-consistent calculations, we obtained that the total energy differences between the traditional correction and the exact treatment are in the range between $-130$ and 620 keV for light nuclei, and 136 and 800 keV for heavier nuclei. A larger value is expected if those induced by the charge density are added together. Therefore, the expression Eq. (1) for Coulomb potential should not be applied for a high precision calculation such as the shell-model description of the isospin-symmetry breaking.