Skyrme-type nuclear interaction as a tool for calculating the finite nuclear size correction to atomic energy levels and the bound-electron $g$ factor

A state-of-the-art approach for calculating the finite nuclear size correction to atomic energy levels and the bound-electron $g$ factor is introduced and demonstrated for a series of highly charged hydrogen-like ions. Firstly, self-consistent mean-field calculations based on the Skyrme-type nuclear interaction are employed in order to produce a realistic nuclear proton distribution. In the second step, the obtained nuclear charge density is used to construct the potential of an extended nucleus, and the Dirac equation is solved numerically. The ambiguity in the choice of a Skyrme parametrization is supressed by fine-tuning of only one parameter of the Skyrme force in order to accurately reproduce the experimental values of nuclear radii in each particular case. The homogeneously charged sphere approximation, the two-parameter Fermi distribution and experimental nuclear charge distributions are used for comparison with our approach, and the uncertainties of the presented calculations are estimated. In addition, suppression of the finite nuclear size effect for the specific differences of $g$ factors is demonstrated.


I. INTRODUCTION
Highly charged ions represent one of the simplest and most well-understood physical systems, and yet they still continue to provide an extremely rich scope of opportunities for fundamental research. They have been extensively used in past years for various high-precision tests of quantum electrodynamics, making it one of the most well-tested theories in physics [1][2][3][4][5]. Such a high accuracy has also been employed for a precise determination of the electron mass [6] and it has been proposed to be used for determination of the fine-structure constant [7][8][9] and even for search of its hypothetical variation [10][11][12][13][14]. Moreover, comparison between the experimental and theoretical results can be used to test theories beyond the Standard Model by setting bounds on parameters of new hypothetical forces [15]. Among various achievements in this field, the most prominent ones include measurements and calculations of the bound-electron g factor in highly charged ions to an extraordinary level of precision [16][17][18][19][20][21]. For all types of high-presicion spectroscopic measurements of highly charged ions [22][23][24] the essential and most fundamental quantities are atomic energy levels and corresponding transition energies which are needed to be known to a high level of accuracy from the theoretical side.
As the experimental precision is being improved, nuclear structure effects are also becoming observable and thus have to be calculated with an increasing accuracy. The largest correction of this kind is due to the finite nuclear size (FNS) effect. Analytical expressions for the FNS effect were presented in [25][26][27][28]. The FNS correction can also be calculated with a higher accuracy numerically by using the Fermi distribution as a model for nuclear charge density [29]. However, even this model does not describe any fine details of nuclear charge distributions that are unique for each nucleus. Hence, in order to perform more precise calculations of the FNS correction, it is necessary to use a more realistic nuclear-structure description and go beyond the simple Fermi model. As for other nuclear-structure corrections, we note that significant improvements in the evaluation of nuclear deformation and nuclear polarization effects have been made in recent years [30][31][32][33][34].
In this paper we present calculations of the FNS correction to atomic energy levels and the bound-electron g factor based on a more detailed description of nuclear charge distributions. The nuclear charge densities are calculated in the framework of the Hartree-Fock method based on the Skyrme-type nuclear interaction with adjustable parameters. We employ the skyrme_rpa program for this purpose [35]. The obtained data is then used to construct the potential of an extended nucleus and numerically solve the Dirac equation for an electron bound in this potential. The theoretically calculated nuclear charge densities are in a good agreement with the experimental ones. These results pave the way for a more accurate description of nuclear-structure effects in atomic systems.
The paper is organized as follows. After a brief description of the computational method we discuss the numerical results and their dependence on the Skyrme parameters. We compare our results to the FNS corrections obtained by using experimental nuclear charge distributions as well as simpler charge density models, such as the homogeneously charged sphere approximation and the Fermi distribution. Then we estimate the uncertainties of our calculations and also demonstrate suppression of the FNS effect for the specific differences [7,36] of g factors.
Relativistic system of units ( = c = 1) and Heaviside charge units (α = e 2 /4π, e < 0) are used throughout the paper. Three-vectors are denoted by bold letters.

A. Skyrme interaction and nuclear charge density
In its standard form, the Skyrme interaction between two nucleons with spatial coordinates r 1 and r 2 can be expressed as [37]: are the Pauli spin matrices, and ρ is the total nucleon density. Here we note that t j , χ j (j ∈ {0, 1, 2, 3}), W 0 and λ are adjustable parameters of the Skyrme force [35]. Next, in order to derive the Hartree-Fock (HF) equations, single-particle wave functions {φ q i (x)} are introduced, where x denotes the set of spatial and spin coordinates, and the superscript q is used to distinguish between the neutron (q = "n") and proton (q = "p") orbitals. The many-body ground-state wave function is built out of these functions as a Slater determinant, and then by means of the variational principle one can obtain the HF equations of the general form: where the Hamiltonian H itself depends on the single particle wave functions. The explicit form of the equations as well as their detailed derivation can be found in various articles, for example in [38]. The HF equations are solved iteratively until self-consistency to a predefined accuracy is achieved. The obtained orbitals can then be used to construct point nucleon densities, in particular, the proton density: Finally, in order to obtain the nuclear charge distribution, the proton density is convoluted with the Gaussian form factor f p (r) to allow for the finite size of the proton [38]: We note here that in the following we assume spherical symmetry of nuclear charge distributions.
Other expressions for the nuclear charge density ρ c (r) are used in this paper for comparison purposes, and they include: 1) the homogeniously charged sphere approximation ("Sphere"): 2) Fermi distribution ("Fermi"): with the radius parameter c and the diffuseness parameter a = (2.3/4 ln3) fm [29]; 3) model-independent analyses of experimental scattering data [39]: a) expansion into a sum of spherical Bessel functions j 0 of order zero ("Bessel"): where R is the cutoff radius; b) expansion into a sum of Gaussians ("Gauss"): where R i and Q i are the positions and the amplitudes of the Gaussians, respectively, and the parameter γ is related to the root-mean-square radius R G of the Gaussians as follows: In this paper, the constants ρ sphere 0 and ρ Fermi 0 as well as the coefficients a ν and Q i in Eqs. (6) -(9) are chosen to fulfil the following normalization condition: 4π ∞ 0 ρ c (r)r 2 dr = 1.

B. Dirac equation
Once the nuclear charge density ρ c (r) is known, one can construct the potential describing the interaction between an electron and the nucleus as follows [40]: where Z is the nuclear charge, and α is the fine structure constant. This potential then enters the Dirac equation which determines the energy levels E and the fourcomponent wave functions ψ(r) of a bound electron [41]: where α and β are the usual Dirac matrices, and m e is the electron mass.
For an arbitraty central potential the electron wave function splits into radial and angular parts as: where n is the principal quantum number, κ is the relativistic angular momentum quantum number, and m is the total magnetic quantum number. The spherical spinors Ω ±κm (θ, ϕ) are the same for any central potential and are well known [42]. Hence, the problem can be reduced to the following set of radial Dirac equations: where the radial functions G(r) and F (r) satisfy the normalization condition: The radial wave functions G(r) and F (r) can then be found analytically for the Coulomb potential [41] or in general case numerically, for example, by expanding them in terms of B-splines and solving the resulting generalized matrix eigenvalue equations [42]. In our basis-set numerical calculations of the radial wave functions we used the dual-kinetic-balance approach [43].
In order to obtain the FNS correction to atomic energy levels, the numerically calculated values E ext [nκ] (in the case of an extended nucleus) are compared to the exact analytical solution E point [nκ] for the Coulomb potential V (r) = −Zα/r (i.e. point-like nucleus): Most generally, a g factor relates the electron's magnetic moment µ (in units of Bohr magneton µ B = |e|/2m e ) to its angular momentum M : where l is the orbital angular momentum, and s is the spin angular momentum. In the Dirac theory, i.e. without taking into account the radiative corrections, g s = 2 for a free electron, and g l is known to be exactly 1 [44]. Thus, the interaction Hamiltonian H int for an electron in an external homogeneous magnetic field B = (0, 0, B z ) can be expressed as: The corresponding first-order Zeeman splitting ∆E can then be written by introducing a new g factor, also called Landé g factor: We note that Eq. (17) is written in such a way to have the same form as for the simpler case where l = 0, and it can be considered as a definition of the Landé g factor of a bound electron.
On the other hand, the electromagnetic four-potential A µ can be chosen in the form (0, A(r) = [B × r]/2), and an application of the minimal coupling principle to the Dirac equation (11) implies: In this way, first-order perturbation theory gives: A calculation of the matrix element in Eq. (19) using the wave functions of the form (12) [45] and then taking into account Eq. (17) yields the following general formula for the g factor: where j = |κ|−1/2 is the total angular momentum quantum number.
In the case of the Coulomb potential V (r) = −Zα/r an exact analytical calculation can be performed [46], and the result reads: Finally, the FNS correction to the g factor for a state nκ is obtained by taking the difference between (20) and (21): Other contributions to the g factor are summarized e.g. in Ref. [29,47,48].

A. Choice of Skyrme parametrization
First, let us discuss the influence of the choice of a Skyrme parameter set on the computational results. For this purpose, we consider the FNS correction to the ground-state (1s 1/2 ) energy and g factor for hydrogenlike lead 208 82 Pb 81+ . In Table I three different widely used parametrizations (known as LNS, SLy5 and SKP [35]) are compared. In the first three columns of Table I some selected Skyrme parameters are presented in order to illustrate large differences between these parameter sets. These differences can be seen even more clearly by comparing the values of root-mean-square (RMS) nuclear radius obtained by using each of the parameter sets. As a result, the FNS corrections presented in the last two columns also vary significantly in such a way that the results can turn out to be larger or smaller than the FNS corrections obtained in the homogeneously charged sphere approximation (using the RMS radius value of 5.5012 fm for 208 82 Pb [49]). However, it is well known that the magnitude of the FNS correction is highly influenced by the value of RMS nuclear radius [25,26]. Hence, it is natural to adjust Skyrme parameters to reproduce the experimental RMS radius beforehand and only then calculate the FNS correction. We found that the RMS radius is most sensitive with respect to varying the Skyrme parameter t 0 . In Table II the results of such adjustments in t 0 (to obtain r 2 = 5.5012 fm) are shown. It can be seen that once the value of RMS radius is fixed, the calculated magnitudes of the FNS corrections become stable, despite the significant differences between the parameter sets. We tested this observation on a wide range of ions and parametrizations, and we found that the ambiguity in the choice of a Skyrme parameter set was largely suppressed in all cases simply by adjusting the RMS nuclear radius.
All the FNS corrections, presented in the following discussion, were obtained by using the SLy5 parameter set, one of the most widely used parametrizations of the Skyrme force, and the parameter t 0 was adjusted to reproduce the experimental values of RMS nuclear radii in each particular case. In Table III (7) 208 82 Pb 81+ . The FNS corrections in the "Bessel" and "Gauss" rows correspond to experimental charge densities. Such densities are obtained by expanding ρ c (r) into sums of spherical zero-order Bessel functions or Gaussians according to Eqs. (8) - (9) and fitting the expansion coefficients (as well as any other parameters) to the measured cross sections. All the values of the fitting parameters used in this paper were taken from Ref. [39]. We note that for 208 82 Pb nucleus two sets of the "Bessel" coefficients are known [50,51], and for simplicity we present here the results only for the parameters from the Ref. [51]. The parameters of all the theoretical charge distributions were adjusted to yield the following experimental values of RMS nuclear radii:  [49].
One peculiar feature can be immediately seen from the results presented in Table III: the values obtained in the "Fermi" and "Skyrme" models agree with each other much better than with the "experimental values". The explanation for this observation comes from the fact that the value of RMS nuclear radius turns out to be a crucial input parameter, and the experimental charge distributions do not reproduce RMS radii to the current level of precision (as used in the "Sphere", "Fermi" and "Skyrme" calculations). This interesting effect can be seen more clearly from the Figure 1, where we compare different nuclear charge distributions employed in the calculations for 40 20 Ca 19+ ion. It is instructive to note that despite the fact that the Skyrme and experimental charge distributions are in an excellent agreement with each other, the difference in the corresponding FNS corrections is larger than even between the "Skyrme" and "Sphere" values. This surprising result simply comes from the fact that the experimental "Gauss" distribution yields r 2 = 3.4797 fm instead of 3.4776 fm, and it emphasizes the great influence of the RMS nuclear radius on the magnitude of the FNS effect. The observation described above suggests a straightforward way to estimate the calculation uncertainties for the FNS corrections. Since the RMS nuclear radius turns out to be the main source of uncertainty, one can simply vary the value of the RMS radius within its experimental error bars in the Skyrme model (by varying the t 0 parameter) and calculate the corresponding variation in ∆E FNS or ∆g FNS . The calculation uncertainties obtained in such a manner are presented in Tables III and IV. C. g factor and cancellation of the FNS effect in specific differences In general, the same trends as described above for the energy levels hold true also in the case of the FNS corrections to the bound-electron g factor. In this last section we additionally consider the specific differences of the g factors in 1s 1/2 and 2s 1/2 states, as well as in 1s 1/2 and 2p 1/2 states, for all the charge distributions mentioned above. These quantities were introduced [7,36] with the aim of supressing the FNS effect. Thus, one can expect the specific differences to have more stable values with respect to the choice of nuclear charge distribution. The specific differences are defined as follows: By expanding the analytical (second-order perturbation theory) expression for ∆g FNS from Ref. [26] in powers of (Zα), we obtain: The calculated values of ∆g FNS = g ext −g point , together with the FNS corrections to the g factors in 1s 1/2 , 2s 1/2 and 2p 1/2 states for 40 20 Ca 19+ , 116 50 Sn 49+ and 208 82 Pb 81+ , are shown in Table IV. It can be seen that for the specific differences g s and g p the FNS effect is indeed suppressed by several orders of magnitude.
However, we also note here that instead of using the analytical expression for ∆g FNS , one could alternatively evaluate ξ s and ξ p numerically, for example, in the framework of the homogeneously charged sphere approximation. In this approach, using the new values of ξ s and ξ p for other nuclear models leads to an even stronger suppression of the FNS effect for the specific differences. For example, in the case of 208 82 Pb 81+ ion, the corrections ∆g FNS [1s 1/2 , 2s 1/2 ] and ∆g FNS [1s 1/2 , 2p 1/2 ] within the Skyrme model become only −1.1 × 10 −9 and 5.4 × 10 −10 , respectively, which is 2 − 3 orders of magnitude smaller than the corresponding values given in Table IV. This shows that in the case of heavy ions a direct numerical calculation of ξ s and ξ p (from the best available nuclear model) should be preferred over using analytical formulas.

IV. CONCLUSIONS AND OUTLOOK
We have demonstrated the use of the Skyrme-Hartree-Fock nuclear-structure method as a tool for calculating the finite nuclear size effect in highly charged ions. We have shown that, despite the fact that various parametrizations of the Skyrme force differ from each other drastically, the ambiguity in the choice of a parameter set can be significantly suppressed by fixing the value of root-mean-square nuclear radius. For this purpose, we suggest adjusting a single Skyrme parameter that has the biggest influence on the value of RMS radius, namely, the parameter t 0 . In this way, the ambiguity associated with the choice of a Skyrme parametrization becomes smaller than the ambiguity stemming from uncertainties in values of nuclear radii.
Our results strongly emphasize the importance of the values of RMS nuclear radii in calculations of FNS corrections. We have demonstrated that in some cases the value of nuclear radius can be even more important than the shape of the nuclear charge distribution. In fact, the FNS corrections obtained by means of our approach and by simply using the Fermi distribution agree with each other within the uncertainties in values of nuclear radii. However, it is clear that the Skyrme model provides a more realistic and thus more reliable description of nuclear charge distributions, which will become crucial in the future when the values of nuclear radii are known to a higher level of accuracy. Table IV. Finite nuclear size (FNS) corrections ∆gFNS to the g factors in 1s 1/2 , 2s 1/2 and 2p 1/2 states for highly charged hydrogen-like ions 40 20 Ca 19+ , 116 50 Sn 49+ and 208 82 Pb 81+ . In the last two columns the magnitudes of the remaining FNS contribution to the specific differences g s and g p are presented. Different models of the nuclear charge distributions were used in the calculations. The presented calculation uncertainties are due to the experimental uncertainties in RMS nuclear radii [49]. The names of the distributions are explained in Section II A. 40 20