# Hydrogen and hydrogen-like-ion bound states and hyperfine splittings: finite
nuclear size effects

Igor Kuzmenko1,2, Tetyana Kuzmenko1, Y. Avishai1,3, Y. B. Band1,2,4
1Department of Physics, Ben-Gurion University of the Negev, Beer-Sheva 84105,
Israel
2Department of Chemistry, Ben-Gurion University of the Negev, Beer-Sheva
84105, Israel
3 Yukawa Institute for Theoretical Physics, Kyoto, Japan
4The Ilse Katz Center for Nano-Science, Ben-Gurion University of the Negev,
Beer-Sheva 84105, Israel

###### Abstract

Using the Dirac equation, we study corrections to electron binding energies
and hyperfine splittings of atomic hydrogen and hydrogen-like ions due to
finite nuclear size (FNS) effects, relativistic QED radiative corrections and
nuclear recoil corrections. Three models for the charge distribution and the
magnetic moment distribution within the nucleus are considered. Calculations
are carried for light atoms (H, He and K) and heavy atoms (Rb, Cs, Pb, Bi, U).
The FNS corrections to the ground-state energy are shown to be smaller than
the electron-nucleus reduced mass corrections, and comparable to the
relativistic QED radiative corrections for the light nuclei, but much larger
than both these corrections for heavy nuclei. Comparison is made with an
experiment on the $1s$-$2s$ transition frequency for hydrogen. FNS corrections
to the ground state hyperfine splitting are comparable in size to the
relativistic QED radiative corrections for light nuclei, but are larger for
heavy nuclei.

##  I Introduction

Hydrogen is the most abundant element in the universe, and atomic hydrogen has
the simplest electronic structure. If the proton in atomic hydrogen is taken
as a point charge, analytic expressions for wave functions and bound-state
energies exist within both the non-relativistic framework (using the
Schrödinger equation) and the relativistic one (using the Dirac equation) [1,
2, 3, 4, 5]. This is also true for hydrogen-like ions (HLI) if the nucleus is
taken as a point charge. Moreover, the ground state hyperfine splitting of
these systems can also be obtained analytically (if the nuclear $g$-factor is
known) [6, 7, 8]. However, nuclei are not point-like; even the proton is a
composite particle with finite radius, $r_{p}=0.8414(19)\times 10^{-15}\mbox{
m}=1.5900(36)\times 10^{-5}$ Bohr [9]. Here we consider the effects of a
finite nuclear size (FNS) $r_{N}$ (and the nuclear charge distribution),
nuclear recoil corrections and QED radiative corrections on the electronic
properties and hyperfine splittings of hydrogen and hydrogen-like ions
(H&HLI).

Finite nuclear size (FNS) corrections (i.e., finite nuclear radius
corrections) to the binding energies of H&HLI have been extensively studied,
see e.g., Refs. [10, 11, 12, 13, 14, 15, 16, 17, 6, 7, 8, 20, 21, 22, 23] and
references quoted therein. Theoretical investigations often use the
Schrödinger equation with relativistic corrections, e.g., spin-orbit coupling,
obtained using $1/c$ expansions [4]. One can apply Rayleigh–Schrödinger
perturbation theory to calculate nuclear structure corrections for the ground
state hyperfine splitting as well as the electronic transition energies [15].
However, the ground state binding energies, $\epsilon$, of H&HLI with large
nuclear charge, $Ze$, are not small compared to $m_{e}c^{2}$, where $m_{e}$ is
the electron mass, hence fully relativistic calculations are necessary. For
example, for the uranium hydrogen-like ion, $\epsilon\approx
0.2588\,m_{e}c^{2}$. Hence, using the Schrödinger equation (even with spin-
orbit interaction) is not satisfactory for heavy nuclei; one must use the
Dirac equation. Moreover, even for light nuclei, one should use relativistic
calculations to obtain high accuracy. In Ref. [15], the authors solve the
Schrödinger equation and add relativistic corrections such as spin-orbit
interaction. They assume that the nuclear charge is distributed homogeneously
inside the nucleus and find the correction to the ground-state energy
$(2/3)Z\epsilon_{H}r_{N}^{2}/a_{0}^{2}\ll Z\epsilon_{H}$, where in Gaussian
units $\epsilon_{H}=e^{2}/a_{0}$ is the Hartree energy, $r_{N}$ is the nuclear
charge radius, and $a_{0}$ is the Bohr radius. However, for high $Z$ nuclei,
the relativistic Dirac equation should be applied. In Refs. [16, 17], the
authors solve the Dirac equation for hydrogen having a FNS with the charge
distributed on the surface of the nucleus, approximate the position dependence
of the wave function inside the nucleus by a polynomial, and find the binding
energies for the ground and excited states. Note that perturbation theory for
the FNS corrections using the Dirac equation is a good approximation for
hydrogen and light  nuclei. Such an approximation leads to a closed form
expression for the FNS correction in the non-relativistic limit [see Eq.
(71)]. However, perturbation theory is not valid for high $Z$ nuclei.

H&HLI systems are traditionally used in ultra-high precision experimental
tests of quantum mechanics, as well as of nuclear structure, since
complications due to many-electron effects are absent [18, 19]. It will be
crucial to compare the results of our calculations with these experiments.
Furthermore, accurate experiments of muonic hydrogen atoms (muonic hydrogen
contains a negatively charged muon instead of an electron) have been reported,
see e.g., [12, 8].

Here, within a relativistic Dirac equation formulation, we calculate the FNS
corrections to the ground ($1S_{1/2}$) and excited $S$ states energies, and to
the ground state hyperfine energy splittings, for a number of elements and
their isotopes for H&HLI. We use three different models for the spherically
symmetric charge and magnetic moment distributions within the nucleus: In
model (a) the nuclear charge and magnetic moment are distributed homogeneously
on the surface of a sphere with radius $r_{a}$. In model (b) the charge and
magnetic moment are distributed homogeneously inside the nucleus with radius
$r_{b}$. In model (c) we use the so called two-parameter Fermi model (TPFM)
for the nuclear charge density and magnetic moment density distributions [21,
22, 24], which fall off smoothly near the surface of nucleus. We use the same
distribution for both nuclear charge and nuclear magnetic moment because of
lack of parameters for the magnetic moment density.

We also calculate FNS corrections to the ground state hyperfine splitting of
H&HLI due to (1) the correction for the charge distribution of the nucleus and
(2) the correction for the magnetic moment distribution of the nucleus, i.e.,
the Bohr-Weisskopf correction [13, 14]. Measurements of the ground state
hyperfine splitting of the H&HLI are among the most accurately measured
quantities [6, 7, 8]. These measurements have been used to determine both
nuclear radii and nuclear $g$-factors. To compare theoretical calculations for
FNS corrections with experimental measurements, one needs to incorporate
relativistic quantum electrodynamics (QED) radiative corrections [20]. FNS
corrections to the ground-state hyperfine splitting have been analyzed in
Refs. [6, 7, 8, 13, 14, 20, 21, 22, 26]. In Ref. [8], the authors report
measurements of the hyperfine splitting of the ground state of hydrogen used
to estimate the proton radius. In Refs. [7, 13, 14, 20, 21], the authors apply
QED radiative corrections and calculate the FNS corrections to the ground-
state hyperfine splitting of H&HLI. References [22, 26] studied the hyperfine
splitting for high-$Z$ HLI.

The paper is organized as follows. In Sec. II, we formulate the Dirac equation
for the H&HLI with potential energy of the electron determined using models
(a) and (b). Electron bound states for H&HLI are studied in Sec. III. The
solution of the Dirac equation for a point-like nucleus is presented in Sec.
III.1. The analytic solution of the Dirac equation for model (a) is given in
Sec III.2 and the numerical solution for model (b) is given in Sec III.3. FNS
effects on the ground-state energy for both models (a) and (b) are calculated
in Sec. IV within perturbation theory. Section V discusses the TPFM for the
nuclear charge distribution. Ground state hyperfine splittings are calculated
in Sec. VI. Section VII estimates the uncertainty of the hyperfine splitting
due to uncertainty in the nuclear $g$-factors, $\delta g_{I}$, and nuclear
radii $\delta r_{N}$ and Sec. VIII compares the results of models (a), (b) and
(c) with experiments. In Sec. IX we summarize our main results. Several
technical details are presented in the Appendices. Specifically, in Appendix A
we discuss the root mean squared nuclear charge radius, $r_{N}$, in Appendix B
we calculate the matrix elements of $\mathbf{e}_{r}\times\bm{\alpha}$, where
$\mathbf{e}_{r}$ is the basis vector in the direction of the position vector
$\mathbf{r}$ from the nucleus to the electron and $\bm{\alpha}$ is the vector
of the Dirac matrices. Appendix D discusses the magnetic moment distribution
for nuclei with one nucleon outside a closed nuclear shell.

##  II Models for FNS potentials

The stationary Dirac equation for H&HLI is given by,

$\Big{[}-i\hbar c\bm{\alpha}\cdot\nabla+V(r)+\beta
m_{e}c^{2}\Big{]}\psi(\mathbf{r})=\varepsilon\psi(\mathbf{r}),$ (1)

where $m_{e}$ is the electron mass, $\bm{\alpha}$ and $\beta$ are the
4$\times$4 Dirac matrices, and $V(r)=V_{C}(r)=-Ze^{2}/r$ is the Coulomb
potential wherein $Z$ is the nuclear charge. In model (a), the nuclear charge
is distributed on the surface of the sphere of radius $r_{a}$ and therefore
the potential is constant within the nucleus, and in model (b) the nuclear
charge density $\varrho_{b}(r)$ is uniformly distributed within a sphere of
radius $r_{b}$,

$\varrho_{b}(r)=\frac{3Ze}{4\pi r_{b}^{3}}\,\Theta(r_{b}-r),$ (2)

where $\Theta(\bullet)$ is the Heaviside theta function. Therefore the
potential for model (b) is quadratic inside the nucleus. Explicit expressions
for the potentials in models (a) and (b) are:

$\displaystyle V_{a}(r)=\begin{cases}\displaystyle-\frac{Ze^{2}}{r_{a}},\ \
r\leq r_{a},\\\ \displaystyle-\frac{Ze^{2}}{r},\ \ r>r_{a},\end{cases}$ (3)
$\displaystyle
V_{b}(r)=\left\\{\begin{array}[]{cc}\displaystyle-\frac{Ze^{2}}{2r_{b}}\Big{(}3-\frac{r^{2}}{r_{b}^{2}}\Big{)},&r\leq
r_{b},\\\ \displaystyle-\frac{Ze^{2}}{r},&r>r_{b},\end{array}\right.$ (6)

The nuclear radii $r_{a}$ and $r_{b}$ are expressed in terms of root mean
square (RMS) nuclear charge radius $r_{N}$ as [10, 27]:

$r_{a}=r_{N},\quad r_{b}=\sqrt{\frac{5}{3}}\,r_{N}.$ (7)

See Appendix A for details.

The difference between the FNS potential and the Coulomb potential is the
repulsive potential

$W_{\nu}(r)=V_{\nu}(r)-V_{C}(r)\quad\nu=a,b.$ (8)

Note that $V_{\nu}(r)=V_{C}(r)+W_{\nu}(r)$ is less attractive than $V_{C}(r)$
in both models. We also define the difference potential
$W_{ba}(r)=V_{b}(r)-V_{a}(r)$,

$W_{ba}(r)=\left\\{\begin{array}[]{cc}\displaystyle\frac{Ze^{2}}{r_{a}}-\frac{Ze^{2}}{2r_{b}}\Big{(}3-\frac{r^{2}}{r_{b}^{2}}\Big{)},&r\leq
r_{a},\\\
\displaystyle\frac{Ze^{2}}{r}-\frac{Ze^{2}}{2r_{b}}\Big{(}3-\frac{r^{2}}{r_{b}^{2}}\Big{)},&r_{a}<r\leq
r_{b},\\\ 0&r>r_{b},\end{array}\right\\}$ (9)

which will be useful in our perturbation analysis below. $W_{ba}(r)$ needs to
be added to $V_{a}(r)$ to get $V_{b}(r)$. Figure 1 plots the potentials
$V_{a}(r)$ and $V_{b}(r)$ given in Eqs. (3) and (6), the repulsive potentials
$W_{a,b}(r)$, and the difference potential $W_{ba}(r)$ versus $r$.

Figure 1: The potentials $V_{\nu=a}(r)$ and $V_{\nu=b}(r)$ versus $r$ for a
hydrogen atom defined in Eqs. (3) and (6) (solid blue and red curves
respectively). The repulsive potentials $W_{\nu}(r)=V_{\nu}(r)-V_{C}(r)$,
(blue and red dashed curves respectively), and the difference potentials
$W_{ab}(r)$ defined in Eq. (9) (dashed-dotted purple curve). Note that for
$r>r_{\nu}$, $W_{\nu}(r)=0$. Here $r_{a}$ and $r_{b}$ are given by Eq. (7).

The nuclear radii $r_{a},r_{b}$ and their uncertainties $\delta r_{a},\delta
r_{b}$, are expressed in terms of the RMS nuclear charge radius $r_{N}$ and
its uncertainty $\delta r_{N}$ as

$\displaystyle r_{a}\pm\delta r_{a}=r_{N}\pm\delta r_{N},$ $\displaystyle
r_{b}\pm\delta r_{b}=\sqrt{\frac{5}{3}}\,\big{(}r_{N}\pm\delta r_{N}\big{)}.$
(10)

The uncertainties $\delta r_{\nu}$ give rise to the uncertainties in
$V_{\nu}(r)$:

$\displaystyle V_{\nu}(r)\pm\delta
V_{\nu}(r)=V_{\nu}(r)\pm\bigg{|}\frac{\partial V_{\nu}(r)}{\partial
r_{\nu}}\bigg{|}\,\delta r_{\nu}.$ (11)

Explicitly,

$\delta V_{a}(r)=\frac{Ze^{2}\delta
r_{a}}{r_{a}^{2}}\Theta\big{(}r_{a}-r\big{)},$ (12)

and

$\delta V_{b}(r)=\frac{3Ze^{2}\delta
r_{b}}{2r_{b}^{2}}\,\bigg{(}1-\frac{r^{2}}{r_{b}^{2}}\bigg{)}\,\Theta\big{(}r_{b}-r\big{)}.$
(13)

In Sec. V we will introduce the TPFM for the nuclear charge distribution [24,
21, 22], and we will refer to it as model (c) [see Eq. (82)]. In this model
the nuclear charge distribution near $r_{N}$ is smooth. Figure 2 plots the
charge distribution for three different smoothing parameters and the charge
distribution of model (b) (see the dashed curve). The shape of the TPFM charge
distribution is sometimes called the Woods-Saxon shape. The potential of the
TPFM charge distribution is also referred to as the Woods-Saxon potential. It
can be determined analytically and then numerical solution of the Dirac
equation can be obtained.

The potential for this model can be determined analytically and then numerical
solution of the Dirac equation eigenstates with this potential by numerically
solving the Dirac equation for heavy nuclei.

Figure 2: The charge density $\varrho$ versus $r/r_{N}$ for the TPFM [see Eq.
(82)] for three different smoothing parameters (solid colored curves) and the
model (b) charge density [see Eq. (2)] (dashed black curve). The integral of
$\varrho$ over all space equals $Ze$.

##  III Dirac Equation Bound states

The Dirac $4$-component spinor wave function $\psi(\mathbf{r})$ can be written
in the form [29],

$\psi(\mathbf{r})=\left(\begin{array}[]{c}g(r)\,\Omega_{jlm}\\\
if(r)\,\Omega_{jl^{\prime}m}\end{array}\right),$ (14)

where $\Omega_{jlm}$ is a normalized spherical harmonic spinor which is an
eigenfunction of the total electronic angular momentum operator squared,
$\mathbf{J}^{2}$, the electronic orbital angular momentum operator squared,
$\mathbf{L}^{2}$, the electronic spin angular momentum operator squared,
$\mathbf{S}^{2}$ and the $z$-projection of the total electronic angular
momentum operator, $J_{z}$. These have eigenvalues $j(j+1)$, $l(l+1)$,
$\frac{3}{4}$ and $m$, where $l=j\pm\frac{1}{2}$ and $l^{\prime}=2j-l$. The
radial wave functions $g(r)$ and $f(r)$ satisfy the equations,

$\displaystyle\hbar
c\,\Big{[}\big{(}rg(r)\big{)}^{\prime}+\frac{\kappa}{r}\,\big{(}rg(r)\big{)}\Big{]}=$
$\displaystyle\big{[}m_{e}c^{2}-V(r)+\varepsilon\big{]}\big{(}rf(r)\big{)},$
(15a) $\displaystyle\hbar
c\,\Big{[}\big{(}rf(r)\big{)}^{\prime}-\frac{\kappa}{r}\,\big{(}rf(r)\big{)}\Big{]}=$
$\displaystyle\big{[}m_{e}c^{2}+V(r)-\varepsilon\big{]}\big{(}rg(r)\big{)},$
(15b)

where $\kappa=-1$ for the $s$-state.  $S$ or $s$ ?

###  III.1 Solution of the Dirac equation (15) for the point-like nucleus

The solution of Eqs. (15) for the point-like nucleus (i.e., for $r_{N}=0$) is
well known [29]; for an $s$-state,

$g(r)=\mathcal{N}\sqrt{m_{e}c^{2}+\varepsilon}\,\Big{[}Q_{1}(\rho)+Q_{2}(\rho)\Big{]}\,\rho^{\gamma-1}e^{-\rho/2},$
(16a)
$f(r)=-\mathcal{N}\sqrt{m_{e}c^{2}-\varepsilon}\,\Big{[}Q_{1}(\rho)-Q_{2}(\rho)\Big{]}\,\rho^{\gamma-1}e^{-\rho/2},$
(16b)

where $\mathcal{N}$ is a normalization constant found from the equation,
$\int\limits_{0}^{\infty}\Big{(}\big{|}g(r)\big{|}^{2}+\big{|}f(r)\big{|}^{2}\Big{)}r^{2}dr=1$.
Here we define the quantities

$\displaystyle\rho=\frac{2\lambda r}{\hbar
c},\quad\lambda=\sqrt{m_{e}^{2}c^{4}-\varepsilon^{2}},\quad\gamma=\sqrt{\kappa^{2}-Z^{2}\alpha^{2}},$
(17)

and $\alpha=\frac{e^{2}}{\hbar c}$ is the fine structure constant (where the
right hand side of the equation is given in Gaussian units). The functions
$Q_{1}(\rho)$ and $Q_{2}(\rho)$ are [29]

$\displaystyle Q_{1}(\rho)$ $\displaystyle=$ $\displaystyle
A_{1}\,F\bigg{(}\gamma-\frac{Z\alpha\varepsilon}{\lambda},\,2\gamma+1,\,\rho\bigg{)},$
(18a) $\displaystyle Q_{2}(\rho)$ $\displaystyle=$ $\displaystyle
A_{2}\,F\bigg{(}\gamma+1-\frac{Z\alpha\varepsilon}{\lambda},\,2\gamma+1,\,\rho\bigg{)},$
(18b)

where $F(a,b,\rho)$ is the Kummer confluent hypergeometric function. The
normalization constants $A_{1}$ and $A_{2}$ are related as follows:

$A_{2}=-\frac{\gamma\lambda-Z\alpha\varepsilon}{\kappa\lambda-Z\alpha
m_{e}c^{2}}\,A_{1}.$ (19)

In the special case,

$\gamma-\frac{Z\alpha\varepsilon}{\lambda}=-n_{r},$ (20)

where $n_{r}$ is a positive integer, the Kummer function reduces to a Laguerre
polynomial,

$F\big{(}-n_{r},2\gamma+1,\rho\big{)}=\frac{\Gamma(2\gamma+1)\Gamma(n_{r}+1)}{\Gamma(2\gamma+n_{r}+1)}\,L_{n_{r}}^{(2\gamma)}(\rho),$

why there is no number?
and thus, the functions $g(r)$ and $f(r)$ in Eq. (16) decay exponentially for
$r\to\infty$. In the special case,

$\gamma-\frac{Z\alpha\varepsilon}{\lambda}=0,\ \ \mbox{and}\quad\kappa<0,$

$A_{2}=0$ and $Q_{1}(\rho)$ decay exponentially for $r\to\infty$. Otherwise,
when

$\gamma-\frac{Z\alpha\varepsilon}{\lambda}\neq-n_{r},$

$Q_{1}(\rho)$ and $Q_{2}(\rho)$ diverge as $e^{\rho}$ when $\rho\to\infty$,
and thus $g(r)$ and $f(r)$ in Eq. (16) diverge as $e^{\rho/2}$.

Let us recall the expressions for the relativistic energies
$\varepsilon_{n}^{(0)}$ of a H&HLI system calculated for a point-like nucleus
[29],

$\displaystyle\varepsilon_{n}^{(0)}=m_{e}c^{2}+hcR_{\infty}\frac{2\mathcal{G}(n)}{1+\mu_{ep}},$
(21)

where, for arbitrary $Z$, $\mu_{ep}\equiv m_{e}/M(Z,A)$ is the electron to
nuclear mass ratio, $M(Z,A)$ is the mass of the nucleus comprising $Z$ protons
and $A-Z$ neutrons, and $cR_{\infty}$ is the Rydberg frequency, which is
related to the Hartree energy $\epsilon_{H}$,
$cR_{\infty}=\frac{\epsilon_{H}}{2h}[=\frac{e^{4}m_{e}}{4\pi\hbar^{3}}$ in
Gaussian units]. In Eq. (21) $n=1,2,3,\ldots$ is the principal quantum number,
and

$\displaystyle\mathcal{G}(n)$ $\displaystyle=$
$\displaystyle\frac{1}{Z^{2}\alpha^{2}}\,\bigg{[}\frac{1}{\sqrt{1+X_{n}^{2}}}-1\bigg{]}$
(22) $\displaystyle=$
$\displaystyle-\frac{1}{Z^{2}\alpha^{2}}\,\frac{X_{n}^{2}}{\sqrt{1+X_{n}^{2}}\,\big{[}\sqrt{1+X_{n}^{2}}+1\big{]}},$

where $X_{n}=\frac{Z\alpha}{\sqrt{1-Z^{2}\alpha^{2}}+n-1}$. We found that the
first form on the right-hand side of Eq. (22) can lead to significant
numerical errors, so for numerical calculations, the second form is
preferable. As noted in Refs. [30, 31, 32, 33], the concept of a reduced mass
in the Dirac equation is ambiguous to some extent. Nevertheless, the necessity
of introducing the reduced mass as the substitution for electron mass in the
Dirac energy and the Dirac equation in order to obtain theoretical results
that are close to experimental results has been stressed in the literature
(see Refs. [30, 31, 32, 33]). The ground state energy in Eq. (21) ($n=1$) is

$\varepsilon_{1}^{(0)}(Z,A)=m_{e}c^{2}+mc^{2}(\sqrt{1-Z^{2}\alpha^{2}}-1),$
(23)

where $m\equiv m(Z,A)$ is the electron-nucleus reduced mass,

$m\equiv m(Z,A)=\frac{m_{e}}{1+\mu_{ep}},$ (24)

and the ground-state wave functions in Eq. (16) are,

$\displaystyle g(r)$ $\displaystyle=$
$\displaystyle\frac{2Z^{3/2}\sqrt{1+\gamma}}{a_{B}^{3/2}\sqrt{\Gamma(2\gamma+1)}}\,\rho^{\gamma-1}e^{-\rho/2},$
(25a) $\displaystyle f(r)$ $\displaystyle=$
$\displaystyle-\frac{2Z^{3/2}\sqrt{1-\gamma}}{a_{B}^{3/2}\sqrt{\Gamma(2\gamma+1)}}\,\rho^{\gamma-1}e^{-\rho/2},$
(25b)

where $a_{B}$ is taken to be the effective Bohr radius, including the reduced
mass effect, and is given (in Gaussian units) by

$a_{B}=\frac{\hbar^{2}}{me^{2}}.$ (26)

The nuclear radii $r_{N}$ are tabulated in Table 1 for a number of nuclei.

Table 1: RMS nuclear charge radii $r_{N}$ and their uncertainties $\delta r_{N}$ (in Bohr radius units $a_{B}$). Data for the proton and deuteron are taken from Ref. [10], and data for the other isotopes are taken from Ref. [27]. Here (and below) $E^{-n}$ is shorthand notation for $\times 10^{-n}$. Isotope ($Z,A$) | $r_{N}$ | $\delta r_{N}$
---|---|---
H(1, 1) | $1.5900E^{-5}$ | $3.6E^{-8}$
H(1, 2) | $4.02132E^{-5}$ | $1.40E^{-8}$
H(1, 3) | $3.3242E^{-5}$ | $6.86E^{-7}$
He(2, 3) | $3.7154E^{-5}$ | $5.7E^{-8}$
He(2, 4) | $3.1662E^{-5}$ | $5.3E^{-8}$
K(19, 39) | $6.4910E^{-5}$ | $3.6E^{-8}$
K(19, 40) | $6.4971E^{-5}$ | $5.3E^{-8}$
K(19, 41) | $6.5230E^{-5}$ | $1.04E^{-7}$
Rb(37, 85) | $7.9437E^{-5}$ | $4.5E^{-8}$
Rb(37, 87) | $7.9348E^{-5}$ | $3.4E^{-8}$
Cs(55, 133) | $9.0784E^{-5}$ | $8.7E^{-8}$
Cs(55, 135) | $9.0833E^{-5}$ | $8.9E^{-8}$
Pb(82, 204) | $1.03563E^{-4}$ | $2.6E^{-8}$
Pb(82, 206) | $1.03750E^{-4}$ | $2.6E^{-8}$
Pb(82, 207) | $1.03827E^{-4}$ | $2.6E^{-8}$
Pb(82, 208) | $1.03958E^{-4}$ | $2.5E^{-8}$
Bi(83, 209) | $1.04334E^{-4}$ | $4.9E^{-8}$
U(92, 235) | $1.10241E^{-4}$ | $7.7E^{-8}$
U(92, 236) | $1.10419E^{-4}$ | $7.2E^{-8}$
U(92, 238) | $1.10683E^{-4}$ | $6.2E^{-8}$

The energies $\varepsilon_{1}^{(0)}$ and $\varepsilon_{2}^{(0)}$ in Eq. (21)
are plotted in Fig. 3 as functions of $Z$. For $Z\alpha<1$,
$\varepsilon_{1}^{(0)}$ and $\varepsilon_{2}^{(0)}$ satisfy the inequalities,

$0<\varepsilon_{1}^{(0)}<\varepsilon_{2}^{(0)}<m_{e}c^{2}.$

For $Z=Z_{\max}=137$, $\varepsilon_{1}^{(0)}=0.02292\,m_{e}c^{2}$ and
$\varepsilon_{2}^{(0)}=0.71516\,m_{e}c^{2}$. For $Z>Z_{\max}$, the parameter
$\gamma$ in Eq. (17) becomes imaginary, and the wave function in Eq. (25) is
singular as $r\to 0$ [29]. Indeed, the Dirac equation (15) with the Coulomb
potential of a point-like nucleus is meaningful only for $Z<137$.

Figure 3: Energies of the $1s$ bound state $\varepsilon_{1}^{(0)}$ (solid red
curve) and the excited $2s$ state $\varepsilon_{2}^{(0)}$ (dashed blue curve)
for the H&HLI isotopes defined in Eq. (21) as functions of the proton number
$Z$. The gold and turquoise regions show the continuous energy spectrum with
$|\varepsilon|>m_{e}c^{2}$. The horizontal line is $\varepsilon=0$, and the
vertical line indicates $Z_{\max}=137$.

###  III.2 Dirac equation (15) solution for Model (a)

Equations (15) for model (a) potential are solved separately for $r>r_{N}$ and
for $r<r_{N}$, to obtain the corresponding wave functions $f_{\pm}(r)$ and
$g_{\pm}(r)$ up to as yet unknown multiplicative constants $A_{\pm}$. The
matching conditions at $r=r_{N}$ are,

$A_{-}g_{-}(r_{N})=A_{+}g_{+}(r_{N}),\quad
A_{-}f_{-}(r_{N})=A_{+}f_{+}(r_{N}).$ (27)

The set of equations (27) has nontrivial solutions when the determinant

$\mathfrak{D}(\varepsilon)\equiv
g_{-}(r_{N})f_{+}(r_{N})-g_{+}(r_{N})f_{-}(r_{N})=0.$ (28)

Note that the condition above is necessary but not sufficient to assure that
$g_{-}=g_{+}$ and $f_{-}=f_{+}$ at $r_{N}$.

####  III.2.1 Solution of Dirac equation (15) for $r>r_{N}$

The solution of Eqs. (15) can be written in the form [29]

$f_{+}(r)=\sqrt{\frac{m_{e}c^{2}-\varepsilon}{m_{e}c^{2}}}\,\Big{[}Q_{1}(\rho)+Q_{2}(\rho)\Big{]}\,\rho^{\gamma-1}e^{-\rho/2},$
(29a)
$g_{+}(r)=\sqrt{\frac{m_{e}c^{2}+\varepsilon}{m_{e}c^{2}}}\,\Big{[}Q_{2}(\rho)-Q_{1}(\rho)\Big{]}\,\rho^{\gamma-1}e^{-\rho/2},$
(29b)

where $\rho$, $\lambda$ and $\gamma$ are given by Eq. (17). The functions
$Q_{1}(\rho)$ and $Q_{2}(\rho)$ satisfy the equations,

$\displaystyle\rho
Q^{\prime}_{1}(\rho)+\bigg{(}\gamma-\frac{Z\alpha\varepsilon}{\lambda}\bigg{)}Q_{1}(\rho)$
$\displaystyle-\bigg{(}\kappa-\frac{Z\alpha
m_{e}c^{2}}{\lambda}\bigg{)}Q_{2}(\rho)$ $\displaystyle=$ $\displaystyle 0,$
(30a) $\displaystyle\rho
Q^{\prime}_{2}(\rho)+\bigg{(}\gamma+\frac{Z\alpha\varepsilon}{\lambda}-\rho\bigg{)}Q_{2}(\rho)$
$\displaystyle-\bigg{(}\kappa+\frac{Z\alpha
m_{e}c^{2}}{\lambda}\bigg{)}Q_{1}(\rho)$ $\displaystyle=$ $\displaystyle 0,$
(30b)

where the prime denotes differentiation with respect to $\rho$, and
$\kappa=-1$ for the $s$ states.  $S$ or $s$ ? Substituting $\kappa=-1$ into
Eq. (30a), we can express $Q_{2}(\rho)$ in terms of $Q_{1}(\rho)$ and
$Q^{\prime}_{1}(\rho)$ as follows,

$\displaystyle Q_{2}(\rho)$ $\displaystyle=$
$\displaystyle-\frac{\lambda}{\lambda+Z\alpha m_{e}c^{2}}\,\Bigg{[}\rho
Q^{\prime}_{1}(\rho)$ (31)
$\displaystyle+\bigg{(}\gamma-\frac{Z\alpha\varepsilon}{\lambda}\bigg{)}\,Q_{1}(\rho)\Bigg{]}.$

Substituting this equation into Eq. (30b), we find

$\rho
Q^{\prime\prime}_{1}(\rho)+\big{(}2\gamma+1-\rho\big{)}Q^{\prime}_{1}(\rho)-\bigg{(}\gamma-\frac{Z\alpha\varepsilon}{\lambda}\bigg{)}Q_{1}(\rho)=0.$
(32)

The solution of Eq. (32) that vanishes at $\rho\to\infty$ is [29]

$Q_{1}(\rho)=AU\bigg{(}\gamma-\frac{Z\alpha\varepsilon}{\lambda},\,2\gamma+1,\,\rho\bigg{)},$
(33)

where $U(a,b,\rho)$ is the confluent hypergeometric function.

Similarly, substituting in Eq. (30b) $\kappa=-1$, we can express $Q_{1}(\rho)$
in terms of $Q_{2}(\rho)$ and $Q^{\prime}_{2}(\rho)$,

$\displaystyle Q_{1}(\rho)$ $\displaystyle=$
$\displaystyle-\frac{\lambda}{\lambda-Z\alpha m_{e}c^{2}}\,\Bigg{[}\rho
Q^{\prime}_{2}(\rho)$ (34)
$\displaystyle+\bigg{(}\gamma+\frac{Z\alpha\varepsilon}{\lambda}-\rho\bigg{)}\,Q_{2}(\rho)\Bigg{]}.$

Substituting Eq. (34) into Eq. (30a), we get

$\displaystyle\rho
Q^{\prime\prime}_{2}(\rho)+\big{(}2\gamma+1-\rho\big{)}Q^{\prime}_{2}(\rho)$
$\displaystyle-$
$\displaystyle\bigg{(}\gamma+1-\frac{Z\alpha\varepsilon}{\lambda}\bigg{)}Q_{2}(\rho)$
$\displaystyle=$ $\displaystyle 0.$ (35)

Solution of Eq. (35) vanishing at $\rho\to\infty$ is [29]

$Q_{2}(\rho)=BU\bigg{(}\gamma+1-\frac{Z\alpha\varepsilon}{\lambda},\,2\gamma+1,\,\rho\bigg{)}.$
(36)

Substituting Eqs. (33) and (36) into Eq. (31), we get

$A=\frac{\lambda}{m_{e}c^{2}},\quad B=\frac{\lambda-Z\alpha
m_{e}c^{2}}{m_{e}c^{2}}.$ (37)

Alternatively, we can substitute Eqs. (33) and (36) into Eq. (34), and get Eq.
(37). Finally, the solution of the Dirac equation (15) for $r>r_{N}$ is,

$\displaystyle f_{+}(r)$ $\displaystyle=$
$\displaystyle\sqrt{\frac{m_{e}c^{2}-\varepsilon}{m_{e}c^{2}}}\,\rho^{\gamma-1}e^{-\rho/2}\times$
$\displaystyle\Bigg{[}\frac{\lambda}{m_{e}c^{2}}\,U\bigg{(}\gamma-\frac{Z\alpha\varepsilon}{\lambda},2\gamma+1,\rho\bigg{)}+$
$\displaystyle\frac{\lambda-Z\alpha
m_{e}c^{2}}{m_{e}c^{2}}\,U\bigg{(}\gamma+1-\frac{Z\alpha\varepsilon}{\lambda},2\gamma+1,\rho\bigg{)}\Bigg{]},$
$\displaystyle g_{+}(r)$ $\displaystyle=$
$\displaystyle-\sqrt{\frac{m_{e}c^{2}+\varepsilon}{m_{e}c^{2}}}\,\rho^{\gamma-1}e^{-\rho/2}\times$
$\displaystyle\Bigg{[}\frac{\lambda}{m_{e}c^{2}}\,U\bigg{(}\gamma-\frac{Z\alpha\varepsilon}{\lambda},2\gamma+1,\rho\bigg{)}-$
$\displaystyle\frac{\lambda-Z\alpha
m_{e}c^{2}}{m_{e}c^{2}}\,U\bigg{(}\gamma+1-\frac{Z\alpha\varepsilon}{\lambda},2\gamma+1,\rho\bigg{)}\Bigg{]}.$

####  III.2.2 Dirac equation (15) solution for $r<r_{N}$

An analytic solution for the case of constant potential inside the nucleus is
well-known [2]. Equation (15a) allows us to express $(rf_{-}(r))$ in terms of
$(rg_{-}(r))$ and $(rg_{-}(r))^{\prime}$,

$\displaystyle\big{(}rf_{-}(r)\big{)}$ $\displaystyle=$
$\displaystyle\frac{\hbar c}{m_{e}c^{2}-V(r_{N})+\varepsilon}\times$ (39)
$\displaystyle\bigg{[}\big{(}rg_{-}(r)\big{)}^{\prime}-\frac{1}{r}\,\big{(}rg_{-}(r)\big{)}\bigg{]}.$

Substituting Eq. (39) into Eq. (15b), we get

$\big{(}rg_{-}(r)\big{)}^{\prime\prime}+k^{2}(\varepsilon)\big{(}rg_{-}(r)\big{)}=0,$
(40)

where

$k(\varepsilon)=\frac{1}{\hbar
c}\sqrt{(\varepsilon-V(r_{N}))^{2}-m_{e}^{2}c^{4}}.$ (41)

The nonsingular solution of Eq. (40) at $r=0$ is,

$g_{-}(r)=C\,\frac{\sin(k(\varepsilon)r)}{k(\varepsilon)r},$ (42)

where $C$ is a normalization constant.

Equation (15b) allows us to express $(rg_{-}(r))$ in terms of $(rf_{-}(r))$
and $(rf_{-}(r))^{\prime}$,

$\displaystyle\big{(}rg_{-}(r)\big{)}$ $\displaystyle=$
$\displaystyle\frac{\hbar c}{m_{e}c^{2}+V(r_{N})-\varepsilon}\times$ (43)
$\displaystyle\bigg{[}\big{(}rf_{-}(r)\big{)}^{\prime}+\frac{1}{r}\,\big{(}rf_{-}(r)\big{)}\bigg{]}.$

Substituting Eq. (43) into Eq. (15a), we get

$\big{(}rf_{-}(r)\big{)}^{\prime\prime}+\bigg{(}k^{2}(\varepsilon)-\frac{2}{r^{2}}\bigg{)}\big{(}rf_{-}(r)\big{)}=0.$
(44)

The nonsingular solution of Eq. (44) at $r=0$ is,

$f_{-}(r)=D\,\frac{k(\varepsilon)r\cos(k(\varepsilon)r)-\sin(k(\varepsilon)r)}{k^{2}(\varepsilon)r^{2}},$
(45)

where $D$ is a normalization constant.

Substituting Eqs. (42) and (45) into Eq. (39), we get

$\frac{D}{C}=\sqrt{\frac{\varepsilon-V(r_{N})-m_{e}c^{2}}{\varepsilon-V(r_{N})+m_{e}c^{2}}}.$
(46)

Alternatively, we can substitute Eqs. (42) and (45) into Eq. (43), and get Eq.
(46). Finally, the solution of the Dirac equation (15) for $r<r_{N}$ is,

$\displaystyle f_{-}(r)$ $\displaystyle=$
$\displaystyle\sqrt{\frac{\varepsilon-V(r_{N})-m_{e}c^{2}}{m_{e}c^{2}}}\times$
(47a)
$\displaystyle\frac{k(\varepsilon)r\cos(k(\varepsilon)r)-\sin(k(\varepsilon)r)}{k^{2}(\varepsilon)r^{2}},$
$\displaystyle g_{-}(r)$ $\displaystyle=$
$\displaystyle\sqrt{\frac{\varepsilon-V(r_{N})+m_{e}c^{2}}{m_{e}c^{2}}}\,\frac{\sin(k(\varepsilon)r)}{k(\varepsilon)r}.$
(47b)

####  III.2.3 Boundary conditions at $r=r_{N}$

The functions $f_{\pm}(r)$ and $g_{\pm}(r)$ satisfy the boundary conditions
specified in Eq. (28) at $r=r_{N}$. Clearly, $f_{\pm}(r_{N})$ and
$g_{\pm}(r_{N})$ in Eqs. (38) and (47), and hence $\mathfrak{D}$ in Eq. (28)
depend on $\varepsilon$. Solving the equation $\mathfrak{D}(\varepsilon)=0$
for $\varepsilon$ yields the energies of the hydrogen atom. These are
tabulated in Table 2. Figure 4(a) plots the function

$\tilde{\mathfrak{D}}(\varepsilon)=\mathfrak{L}(\varepsilon)\mathfrak{D}(\varepsilon),$
(48)

versus $\varepsilon$ for $Z=1$, $r_{N}=0.8414\times 10^{-13}$ cm (the proton
RMS radius [10]) and $|\varepsilon|<0.9\,m_{e}c^{2}$, where
$\mathfrak{D}(\varepsilon)$ is given in Eq. (28), and

$\mathfrak{L}(\varepsilon)=\frac{m_{e}^{2}c^{4}-\varepsilon^{2}}{m_{e}^{2}c^{4}}.$

$\mathfrak{L}(\varepsilon)$ is positive for $|\varepsilon|<m_{e}c^{2}$. Both
$\tilde{\mathfrak{D}}(\varepsilon)$ and $\mathfrak{L}(\varepsilon)$ are
dimensionless. The wave functions in Eqs. (38) and (47) are dimensionless and
will be normalized in what follows [by multiplying by the factor
$\mathcal{N}_{n}$ given by Eq. (49), where the dimensions of $\mathcal{N}_{n}$
is $[\mathcal{N}_{n}]=\mathrm{Length}^{-3/2}$].
$\tilde{\mathfrak{D}}(\varepsilon)$ is positive, and there are no bound states
in the interval $-m_{e}c^{2}<\varepsilon<0.99\,m_{e}c^{2}$. Figure 4 shows
$\tilde{\mathfrak{D}}(\varepsilon)$ for $\varepsilon$ close to and below
$m_{e}c^{2}$. $\tilde{\mathfrak{D}}(\varepsilon)=0$ when $\varepsilon$ is very
close to the Dirac energies $\varepsilon_{n}^{(0)}$ of the hydrogen atom given
by Eq. (21) with $Z=A=1$, where the subscript $(0)$ indicates the point charge
nucleus limit.

Figure 4: (a) The function $\tilde{\mathfrak{D}}(\varepsilon)$ defined in Eq. (48) versus the energy $\varepsilon$ for $|\varepsilon|<0.999\,m_{e}c^{2}$. (b) Log-linear plot of $\tilde{\mathfrak{D}}(\varepsilon)$ in Eq. (48) as a function of the binding energy $\epsilon=mc^{2}-\varepsilon$ (blue curve) for $\varepsilon$ close to $mc^{2}$ (the range includes the 11 lowest $s$ 1H states). The red dots show the Dirac binding energies $\epsilon_{n}^{(0)}$ of the hydrogen atom with zero proton radius (i.e., the Dirac energies). Here $\tilde{\mathfrak{D}}(\varepsilon)$ is dimensionless, and the binding energy in (b) is in Hartree.  $n$ | 1 | 2 | 3
---|---|---|---
$\epsilon_{n}$ | $0.499734640577$ | $0.1249340759790$ | $0.0555260711867$
$\epsilon_{n}^{(0)}$ | $0.499734640745$ | $0.1249340760015$ | $0.0555260711927$
Table 2: Hydrogen binding energies $\epsilon_{n}(r_{N})$ in Hartree calculated
for model (a) using Eq. (28), and the hydrogen atom Dirac energies
$\epsilon_{n}^{(0)}$ for $r_{N}=0$ (i.e., for point-like proton) in Eq. (21).
$r_{N}$ is given in Table 1.

The binding energies $\epsilon_{n}(r_{N})$ (with $n=1,2,3$ and $r_{N}$ being
the nuclear radius of 1H) calculated numerically from Eq. (28) and
$\epsilon_{n}^{(0)}=m_{e}c^{2}-\varepsilon_{n}^{(0)}$, with
$\varepsilon_{n}^{(0)}$ given in Eq. (21), are listed in Table 2. It is seen
that
$\big{|}\epsilon_{1}(r_{N})-\epsilon_{1}^{(0)}\big{|}~{}\approx~{}1.680\times
10^{-10}$ Hartree, which is clearly very much smaller than
$\epsilon_{1}^{(0)}$. The Dirac energies $\varepsilon_{n}(r_{N})$ are pushed
up by the repulsive potential
$W_{a}(r)=[V_{C}(r_{N})-V_{C}(r)]\Theta(r_{N}-r)$ [see Eq. (8)]. Note the
following inequalities:

$\big{|}\epsilon_{n+1}(r_{N})-\epsilon_{n+1}^{(0)}\big{|}~{}\ll~{}\big{|}\epsilon_{n}(r_{N})-\epsilon_{n}^{(0)}\big{|},$

where $n=1,2,3,\ldots$

####  III.2.4 Wave functions

We now consider normalization of the wave functions $g_{n}$ and $f_{n}$. Here
$g_{n}(r)$ is given by Eqs. (LABEL:eq:g-res-plus) for $r>r_{N}$ and by Eq.
(47b) for $r<r_{N}$, and $f_{n}(r)$ is given by Eq. (LABEL:eq:f-res-plus) for
$r>r_{N}$ and by Eq. (47a) for $r<r_{N}$. The normalization constants
$\mathcal{N}_{n}$ are chosen such that

$\mathcal{N}_{n}^{2}\int\limits_{0}^{\infty}[g_{n}^{2}(r)+f_{n}^{2}(r)]r^{2}dr=1.$
(49)

Figure 5 plots $\mathcal{N}_{1}(g_{1}(r),f_{1}(r))$ for the $1s$ state. The
functions $(g_{1}(r),f_{1}(r))$ and their derivatives
$(g^{\prime}_{1}(r),f^{\prime}_{1}(r))$ are continuous at $r=r_{N}$, where
$r_{N}=1.658\times 10^{-5}a_{B}$. $g_{1}(r)$ has no nodes, has a maximum at
$r=0$, and vanishes as $r\to\infty$. $f_{1}(r)$ vanishes at $r=0$ and is
negative for finite $r$; it reaches its minimum and tends to zero as
$r\to\infty$.

Figure 5: Log-linear plot of (a) $\mathcal{N}g_{1}(r)$ and (b)
$\mathcal{N}f_{1}(r)$ for the $1s$ state of hydrogen atom with the energy
$\varepsilon_{1}$ (in Hartree). Here $g_{1}(r)$ is given by Eq.
(LABEL:eq:g-res-plus) for $r>r_{N}$ (purple) and by Eq. (47b) for $r<r_{N}$
(blue). $f_{1}(r)$ is given by Eq. (LABEL:eq:f-res-plus) for $r>r_{N}$
(purple) and by Eq. (47a) for $r<r_{N}$ (blue). The normalization constant
$\mathcal{N}$ is found from Eq. (49). The red dots show $g_{1}(r_{N})$ and
$f_{1}(r_{N})$.

Figure 6 is a log-linear plot for the $2s$ normalized wave function
$\mathcal{N}_{2}(g_{2}(r),f_{2}(r))$. The functions $(g_{2}(r),f_{2}(r))$ and
their derivatives $(g^{\prime}_{2}(r),f^{\prime}_{2}(r))$ are continuous at
$r=r_{N}$. $g_{2}(r)$ has a maximum at $r=0$, and one node. It vanishes at
$r\gg a_{B}$. $f_{2}(r)$ has two nodes: one node at $r=0$, and another node at
finite $r$; it vanishes as $r\to\infty$.

Figure 6: The log-linear plot of (a) $\mathcal{N}_{2}g_{2}(r)$ and (b)
$\mathcal{N}_{2}f_{2}(r)$ for the $2s$ state of hydrogen atom with the energy
$\varepsilon_{2}$ (in Hartree). Here $g_{2}(r)$ is given by Eq.
(LABEL:eq:g-res-plus) for $r>r_{N}$ (purple) and by Eq. (47b) for $r<r_{N}$
(blue). $f_{2}(r)$ is given by Eq. (LABEL:eq:f-res-plus) for $r>r_{N}$
(purple) and by Eq. (47a) for $r<r_{N}$ (blue). The normalization constant
$\mathcal{N}_{2}$ is found from Eq. (49). The red dots show the wave functions
at $r=r_{N}$.

####  III.2.5 $1s$ and $2s$ state energies for various isotopes

We now solve the determinantal equation (28) for model (a) (constant potential
inside the nucleus) and find the FNS corrections to the energies of the ground
$1s$ state and the first excited $2s$ state for various isotopes.

Table 3: FNS corrections $\varepsilon_{1,a}$ and $\varepsilon_{1,b}$ (in Hartree) to the energy of the ground $1s$ state of H&HLI. The numbers in the brackets indicate the uncertainties $\delta\varepsilon_{1,a}$ and $\delta\varepsilon_{1,b}$ in the calculated ground-state energy due to the uncertainty in the nuclear radius $r_{N}$ in Table 1 [see Eqs. (51) and (61)]. The concise notation used for errors is explained in Ref. [34]. Isotope ($Z,A$) | $\varepsilon_{1,a}(\delta\varepsilon_{1,a})$ | $\varepsilon_{1,b}(\delta\varepsilon_{1,b})$
---|---|---
H(1, 1) | $1.6837(76)E^{-10}$ | $1.6837(76)E^{-10}$
H(1, 2) | $1.07773(75)E^{-9}$ | $1.07773(75)E^{-9}$
H(1, 3) | $7.37(30)E^{-10}$ | $7.37(30)E^{-10}$
He(2, 3) | $1.4746(45)E^{-8}$ | $1.4746(45)E^{-8}$
He(2, 4) | $1.0711(36)E^{-8}$ | $1.0711(36)E^{-8}$
K(19, 39) | $4.1987(46)E^{-4}$ | $4.1974(46)E^{-4}$
K(19, 40) | $4.2064(68)E^{-4}$ | $4.2051(68)E^{-4}$
K(19, 41) | $4.240(13)E^{-4}$ | $4.238(13)E^{-4}$
Rb(37, 85) | $1.2547(14)E^{-2}$ | $1.2534(14)E^{-2}$
Rb(37, 87) | $1.2520(12)E^{-2}$ | $12507(12)E^{-2}$
Cs(55, 133) | $0.13135(23)$ | $0.13107(23)$
Cs(55, 135) | $0.13148(23)$ | $0.13120(23)$
Pb(82, 204) | $2.4687(10)$ | $2.4586(10)$
Pb(82, 206) | $2.4758(10)$ | $2.4657(10)$
Pb(82, 207) | $2.4787(10)$ | $2.4686(10)$
Pb(82, 208) | $2.48367(93)$ | $2.47351(93)$
Bi(83, 209) | $2.7559(20)$ | $2.7444(20)$
U(92, 235) | $7.3062(75)$ | $7.2718(75)$
U(92, 236) | $7.3234(69)$ | $7.2889(69)$
U(92, 238) | $7.3490(60)$ | $7.3143(60)$

Tables 3 and 4 show the FNS corrections

$\varepsilon_{n,a}=\varepsilon_{n}-\varepsilon_{n}^{(0)},\quad n=1,2.$ (50)

Here $\varepsilon_{1}$ and $\varepsilon_{2}$ are calculated numerically using
Eq. (28), and $\varepsilon_{1}^{(0)}$ and $\varepsilon_{2}^{(0)}$ are given in
analytic form by Eq. (21). $\varepsilon_{1,a}$ and $\varepsilon_{2,a}$ are
functions of the nuclear radius $r_{a}=r_{N}$, where $r_{N}$ is tabulated in
Table 1. The uncertainty in $r_{N}$ gives rise to uncertainty
$\delta\varepsilon_{n,a}$ in $\varepsilon_{n,a}$, which can be estimated as,

$\delta\varepsilon_{n,a}=\bigg{|}\frac{\partial\varepsilon_{n,a}}{\partial
r_{N}}\bigg{|}\,\delta r_{N}.$ (51)

The FNS correction $\varepsilon_{1,a}$ to the ground-state energy, and the
uncertainty $\delta\varepsilon_{1,a}$ are tabulated in Table 3, and the FNS
correction $\varepsilon_{2,a}$ to the excited-state energy, and the
uncertainty $\delta\varepsilon_{2,a}$ are tabulated in Table 4. The
corrections $\varepsilon_{1,a}$ and $\varepsilon_{2,a}$ satisfy the
inequalities $\varepsilon_{1,a}\gg\varepsilon_{2,a}>0$. Note that, for
isotopes of atoms with a given $Z$, $\varepsilon_{1,a}$ and
$\varepsilon_{2,a}$ do not always increase with $A$. This is because the
nuclear radii $r_{N}$ in Table 1 for some isotopes are not monotonic with $A$
(see $r_{N}$ for He and Rb). The repulsive potential $W_{a}(r)$ in Eq. (8)
increases with the nuclear radius $r_{N}$, and $r_{N}$ increases with $Z$, see
Table 1, hence the energy correction $\varepsilon_{a}$ in Table 3 increases
with $Z$.

Table 4: FNS corrections $\varepsilon_{2,a}$ and $\varepsilon_{2,b}$ (in Hartree) to the energy of the excited $2s$ state of H&HLI. The numbers in brackets indicate the uncertainties $\delta\varepsilon_{2,a}$ and $\delta\varepsilon_{2,b}$ in the calculated ground-state energy due to the uncertainty in the nuclear radius $r_{N}$ given in Table 1 [see Eqs. (51) and (61)]. Isotope ($Z,A$) | $\varepsilon_{2,a}(\delta\varepsilon_{2,a})$ | $\varepsilon_{2,b}(\delta\varepsilon_{2,b})$
---|---|---
H(1, 1) | $2.1047(95)E^{-11}$ | $2.1047(95)E^{-11}$
H(1, 2) | $1.34722(94)E^{-10}$ | $1.34722(94)E^{-10}$
H(1, 3) | $9.21(38)E^{-11}$ | $9.21(38)E^{-11}$
He(2, 3) | $1.8436(56)E^{-9}$ | $1.8436(56)E^{-9}$
He(2, 4) | $1.3391(45)E^{-9}$ | $1.3391(45)E^{-9}$
K(19, 39) | $5.3381(58)E^{-5}$ | $5.3365(58)E^{-5}$
K(19, 40) | $5.3480(86)E^{-5}$ | $5.3463(86)E^{-5}$
K(19, 41) | $5.390(17)E^{-5}$ | $5389(17)E^{-5}$
Rb(37, 85) | $1.6732(18)E^{-3}$ | $1.6714(18)E^{-3}$
Rb(37, 87) | $1.6696(16)E^{-3}$ | $1.6678(16)E^{-3}$
Cs(55, 133) | $1.8969(33)E^{-2}$ | $1.8928(33)E^{-2}$
Cs(55, 135) | $1.8988(34)E^{-2}$ | $1.8947(34)E^{-2}$
Pb(82, 204) | $0.42868(17)$ | $0.42692(17)$
Pb(82, 206) | $0.42991(17)$ | $0.42814(17)$
Pb(82, 207) | $0.43041(17)$ | $0.42864(17)$
Pb(82, 208) | $0.43127(16)$ | $0.42950(16)$
Bi(83, 209) | $0.48259(36)$ | $0.48057(36)$
U(92, 235) | $1.3880(14)$ | $1.3814(14)$
U(92, 236) | $1.3912(13)$ | $1.3846(13)$
U(92, 238) | $1.3961(11)$ | $1.3895(11)$

Table 4 shows the calculated FNS corrections $\varepsilon_{2,a}$ and
$\varepsilon_{2,b}$. The FNS corrections to the $1s$-$2s$ hydrogen transition
frequency $\Delta\omega_{1s,2s}^{(a)}$ and $\Delta\omega_{1s,2s}^{(b)}$ are
defined as,

$\Delta\omega_{1s,2s}^{(a)}=\frac{1}{\hbar}\,\big{(}\varepsilon_{2,a}-\varepsilon_{1,a}\big{)},\quad\Delta\omega_{1s,2s}^{(b)}=\frac{1}{\hbar}\big{(}\varepsilon_{2,b}-\varepsilon_{1,b}\big{)}.$
(52)

Note that the method to accurately measure the isotope shifts of atomic
transitions reported in Ref. [38] can be used to experimentally probe the
isotope shifts calculated here.

Equation (21) for a point-like nucleus and the numerical calculations with the
nuclear charge distribution for models (a) and (b) take into account recoil
effects by replacing the electron mass $m_{e}$ by the reduced mass $m\equiv
m(Z,A)$ defined in Eq. (24). This naive expression can be used as a starting
point for calculating the recoil corrections to the electron energy levels. An
explicit expression for the $s$-state energies with the reduced mass
dependence is derived from quantum electrodynamics in Ref. [33],

$\displaystyle\varepsilon_{n}^{\mathrm{(Dirac)}}$ $\displaystyle=$
$\displaystyle m_{e}c^{2}+Z^{2}hcR_{\infty}\times$ (53)
$\displaystyle\Big{[}\frac{2\mathcal{G}(n)}{1+\mu_{ep}}-\frac{\mu_{ep}\,Z^{2}\alpha^{2}\,\mathcal{G}^{2}(n)}{\big{(}1+\mu_{ep}\big{)}^{3}}\Big{]}.$

The dimensionless function $\mathcal{G}(n)$ is given in Eq. (22). The first
two terms on the right hand side of Eq. (53) give the energy of the $s$ state
in Eq. (21), but the last term in this equation, which is proportional to
$\mathcal{G}^{2}(n)$ is proportional to $\mu_{ep}(Z\alpha)^{4}$. This term
vanishes in the limit of infinitely large nuclear mass, and hence is absent in
Eq. (21). For hydrogen, this correction term to the energy is

$hcR_{\infty}\frac{\mu_{e}\,\alpha^{2}\,\mathcal{G}^{2}(n)}{\big{(}1+\mu_{ep}\big{)}^{3}}\approx\frac{7.2\times
10^{-9}}{n^{2}}~{}\epsilon_{H}.$

The Dirac energy in Eq. (53) can be expanded in a power series in
$Z^{2}\alpha^{2}$, as in Eq. (3.4) of Ref. [33]:

$\displaystyle\tilde{\varepsilon}_{n}^{\mathrm{(Dirac)}}$ $\displaystyle=$
$\displaystyle
m_{e}c^{2}-\frac{hcR_{\infty}}{1+\mu_{ep}}\,\bigg{\\{}\frac{Z^{2}}{n^{2}}+\frac{Z^{4}\alpha^{2}}{n^{3}}\,\Big{(}1-\frac{3}{4n}$
(54)
$\displaystyle+\frac{\mu_{ep}}{(1+\mu_{ep})^{2}}\,\frac{1}{4n}\Big{)}\bigg{\\}}.$

Equation (54) is a power series expansion of Eq. (53) in $\alpha^{2}$, up to
$\alpha^{4}$ (recall that a factor of $\alpha^{2}$ is contained in
$R_{\infty}$). The $1s$ energy obtained using (54) differs from that obtained
using (53) by about a MHz.

###  III.3 Dirac equation (15) solutions for model (b)

Here we find the wave functions and binding energies of H&HLI with nuclear
charge distributed homogeneously inside a sphere of radius $r_{b}$ defined in
Eq. (7). The potential $V_{b}(r)$ in Eq. (6) is harmonic inside the nucleus.
In this case an analytic solution is unknown. We rewrite Eq. (15) with
$\kappa=-1$ as,

$\displaystyle\frac{\hbar
c\,r^{2}g^{\prime}_{-}(r)}{mc^{2}-V_{b}(r)+\varepsilon}=r^{2}f_{b,-}(r),$
(55a)
$\displaystyle\big{(}r^{2}f_{b,-}(r)\big{)}^{\prime}=\frac{mc^{2}+V_{b}(r)-\varepsilon}{\hbar
c}\,r^{2}g_{b,-}(r).$ (55b)

Differentiating Eq. (55a) and using Eq. (55b) we get the following second-
order differential equation for $g_{b,-}(r)$:

$\displaystyle\frac{1}{r^{2}}\,\frac{d}{dr}\bigg{(}\frac{r^{2}}{mc^{2}-V_{b}(r)+\varepsilon}\,\frac{dg_{b,-}(r)}{dr}\bigg{)}$
$\displaystyle=\frac{1}{\hbar^{2}c^{2}}\,\big{[}mc^{2}+V_{b}(r)-\varepsilon\big{]}g_{b,-}(r).$
(56)

We solve Eq. (56) numerically and find the radial wave function $g_{b,-}(r)$
for $r<r_{b}$. For this purpose, we need the binding energy of the atom (which
we do not know). In order to find the energy and the wave function, we apply
the iteration method described below. The radial wave function $f_{b,-}(r)$
can be found from the equation,

$f_{b,-}(r)=\frac{\hbar c\,g^{\prime}_{b,-}(r)}{mc^{2}-V_{b}(r)+\varepsilon}.$
(57)

The ground-state wave functions satisfy the boundary conditions at $r=0$,

$f_{b,-}(0)=0,\quad g_{b,-}(0)=1.$ (58)

The radial wave functions $g_{+}(r)$ and $f_{+}(r)$ with $r>r_{b}$ are given
in Eq. (38), and the radial wave functions $g_{b}(r),f_{b}(r)$ satisfy the
boundary conditions at $r=r_{b}$,

$A_{-}g_{b,-}(r_{b})=A_{+}g_{+}(r_{b}),\quad
A_{-}f_{b,-}(r_{b})=A_{+}f_{+}(r_{b}).$ (59)

Here $A_{\pm}$ are (unknown) normalization constants. The set of equations
(59) has nontrivial solutions when the determinant $\mathfrak{D}(\varepsilon)$
vanishes, where

$\mathfrak{D}(\varepsilon)=g_{b,-}(r_{b})f_{+}(r_{b})-g_{+}(r_{b})f_{-}(r_{b}).$
(60)

Appendix C explains the iteration method for finding the ground state and
excited $s$ state energies using Eq. (60).

The uncertainty $\delta\varepsilon_{n,b}$ in the model (b) energy
$\varepsilon_{n,b}$ of the $ns$ state is estimated as,

$\delta\varepsilon_{n,b}=\int\limits_{0}^{r_{b}}\delta
V_{b}(r)\Big{[}\big{|}g_{n,-}(r)\big{|}^{2}+\big{|}f_{n,-}(r)\big{|}^{2}\Big{]}r^{2}dr,$
(61)

where $\delta V_{b}(r)$ is given in Eq. (13), and $n=1$ for the ground state.

The FNS correction $\varepsilon_{1,b}$ to the ground-state energy, and the
uncertainty $\delta\varepsilon_{1,b}$ are tabulated in Table 3, and the FNS
correction $\varepsilon_{2,b}$ to the excited-state energy, and the
uncertainty $\delta\varepsilon_{2,b}$ are tabulated in Table 4. Note that for
light isotopes, i.e., for H, He and K, $\varepsilon_{n,a}$ and
$\varepsilon_{n,b}$ [$n=1,2$] and their uncertainties
$\delta\varepsilon_{n,a}$ and $\delta\varepsilon_{n,b}$ satisfy the relations

$\delta\varepsilon_{n,a}\approx\delta\varepsilon_{n,b}>\varepsilon_{n,a}-\varepsilon_{n,b},$

whereas for the heavy isotopes,

$\delta\varepsilon_{n,a}\approx\delta\varepsilon_{n,b}<\varepsilon_{n,a}-\varepsilon_{n,b}.$

Therefore, the comparison of the calculated and measured ground-state binding
energies of heavy HLI (such as Rb, Cs, Pb, Bi and U) can be used to determine
the nuclear charge distribution. The ground-state energy of the hydrogen atom,
and of light HLI (such as He and K) are almost the same for models (a) and
(b), hence, measuring the ground-state binding energies would not be an
accurate method for determining the nuclear charge distributions.

The numerical solution technique developed in this subsection can be
generalized to the case of other charge distributions inside the nucleus,
e.g., for the TPFM charge distribution [21, 22] introduced in Eq. (82) below.

We calculated the FNS corrections analytically for model (a) and numerically
for model (b). The results were tabulated in Table 3 for the ground $1s$ state
of H&HLI, and in Table 4 for the first excited $2s$ state. For comparison, we
will calculate in Sec. IV below the FNS correction using perturbation theory.
We shall show that the perturbation theory is a good approximation for the FNS
corrections for H and light HLI, compare Tables 7 and 3. For heavy HLI, the
perturbation theory results in Table 7 are dissimilar to the results in Table
3, i.e., the perturbation theory is poor. References [10] and [11] (2018
CODATA and 2002 CODATA) and many other papers use perturbation theory to
calculate FNS corrections for hydrogen and light HLIs. These references expand
the FNS corrections in a power series in $r_{N}$. In contrast, we numerically
calculate the perturbation theory.

###  III.4 Finite-nuclear size corrections in 2018 CODATA

Finite-nuclear-size and nuclear-polarizability corrections are ordered by
powers in $\alpha$, following Ref. [10]:

$\varepsilon_{n}^{(2018)}=\sum_{i=4}^{\infty}\varepsilon_{n}^{(i)},$ (62)

where $n$ is the principal quantum number, and index $i$ indicates the power
of $Z\alpha$ in the expansion. The first and lowest-order contributions to
$s$-state energies are

$\varepsilon_{n}^{(4)}=\frac{2}{3}\,m_{e}c^{2}\,\frac{(Z\alpha)^{4}}{n^{3}}\,\bigg{(}\frac{m}{m_{e}}\bigg{)}^{3}\,\bigg{(}\frac{r_{N}}{\lambdabar_{C}}\bigg{)}^{2},$
(63)

where $\lambdabar_{C}=\hbar/(m_{e}c)$ is the reduced Compton wavelength. The
numerical value of $\lambdabar_{C}$ is given in Ref. [10] as,
$\lambdabar_{C}=3.8615926796(12)\times 10^{-12}$ m. Equation (63) is a non-
relativistic approximation to the FNS correction to the $s$-state energy of
H&HLI. Reference [10] gives relativistic corrections $\varepsilon_{n}^{(i)}$
with $i\geq 5$ just for hydrogen and deuterium atoms, but not for tritium
atoms or for hydrogen-like ions.

Table 5: FNS corrections $\varepsilon_{1}^{(4)}$ (in Hartree) to the energy of the ground $1s$ state of H&HLI. The numbers in the brackets indicate the uncertainties $\delta\varepsilon_{1}^{(4)}$ in the calculated ground-state energy due to the uncertainty in the nuclear radius $r_{N}$ in Table 1. $\Delta\varepsilon_{1,\nu}$ with $\nu=a,b$ is given in Eq. (64).  Isotope ($Z,A$) | $\varepsilon_{1}^{(4)}(\delta\varepsilon_{1}^{(4)})$ | $\Delta\varepsilon_{1,a}^{(4)}$ | $\Delta\varepsilon_{1,b}^{(4)}$
---|---|---|---
H(1, 1) | $1.6827(76)E^{-10}$ | $10E^{-13}$ | $10E^{-13}$
H(1, 2) | $1.07719(75)E^{-9}$ | $5.4E^{-13}$ | $5.4E^{-13}$
H(1, 3) | $7.36(30)E^{-10}$ | $4E^{-13}$ | $4E^{-13}$
He(2, 3) | $1.4716(45)E^{-8}$ | $3.0E^{-11}$ | $3.0E^{-11}$
He(2, 4) | $1.0689(36)E^{-8}$ | $2.2E^{-11}$ | $2.2E^{-11}$
K(19, 39) | $3.6604(40)E^{-4}$ | $5.383E^{-5}$ | $5.370E^{-5}$
K(19, 40) | $36673(60)E^{-4}$ | $5.392E^{-5}$ | $5.379E^{-5}$
K(19, 41) | $3.697(12)E^{-4}$ | $5.43E^{-5}$ | $5.42E^{-5}$
Rb(37, 85) | $7.884(9)E^{-3}$ | $4.663E^{-3}$ | $4.650E^{-3}$
Rb(37, 87) | $7.866(8)E^{-3}$ | $4.653E^{-3}$ | $4.640E^{-3}$
Cs(55, 133) | $5.028(10)E^{-2}$ | $8.108E^{-2}$ | $8.079E^{-2}$
Cs(55, 135) | $5.033(10)E^{-2}$ | $8.115E^{-2}$ | $8.087E^{-2}$
Pb(82, 204) | $0.32327(17)$ | $2.1455$ | $2.1354$
Pb(82, 206) | $0.32444(17)$ | $2.1514$ | $2.1412$
Pb(82, 207) | $0.32493(17)$ | $2.1538$ | $2.1437$
Pb(82, 208) | $0.32574(16)$ | $2.15793$ | $2.14777$
Bi(83, 209) | $0.34440(32)$ | $2.4115$ | $2.4000$
U(92, 235) | $0.58042(82)$ | $6.7258$ | $6.6913$
U(92, 236) | $0.58229(76)$ | $6.7411$ | $6.7066$
U(92, 238) | $0.58509(66)$ | $6.7639$ | $6.7293$

For H&HLIs the FNS correction $\varepsilon_{n}^{(4)}$ to the $ns$-state energy
and its uncertainty $\delta\varepsilon_{n}^{(4)}$ due to the uncertainty of
the nuclear RMS charge radius are tabulated in Table 5 for the ground $n=1$
state, and in Table 6 for the first excited $n=2$ state. These tables also
show $\Delta\varepsilon^{(4)}_{n,a}$ and $\Delta\varepsilon^{(4)}_{n,b}$ where

$\Delta\varepsilon_{n,\nu}^{(4)}=\varepsilon_{n}^{(4)}-\varepsilon_{n,\nu},$
(64)

and $\nu=a,b$ for the models (a) and (b) of the nuclear charge distribution.
We will now compare $\varepsilon_{1}^{(4)}\approx\varepsilon_{1}^{(2018)}$
with $\varepsilon_{1,a}$ and $\varepsilon_{1,b}$ in Table 3, and
$\varepsilon_{2}^{(4)}\approx\varepsilon_{2}^{(2018)}$ with
$\varepsilon_{2,a}$ and $\varepsilon_{2,b}$ in Table 4 for H&HLIs. Moreover,
for hydrogen and light HLIs, we compare $\Delta\varepsilon_{n,\nu}^{(4)}$ with
the uncertainty $\delta\varepsilon_{n,\nu}$ of the FNS correction in Tablea 3
and 4. For AH atoms and He ions,
$\Delta\varepsilon_{n,\nu}^{(4)}<\delta\varepsilon_{n,\nu}$, and therefore the
non-relativistic approximation in Eq. (63) is a good approximation to the FNS
correction to the $1s$ state energy ($n=1$) and to the $2s$ state energy
($n=2$). For heavier HLI,
$\Delta\varepsilon_{n,\nu}^{(4)}>\delta\varepsilon_{n,\nu}$, and the non-
relativistic approximation fails.

Table 6: FNS corrections $\varepsilon_{2}^{(4)}$ (in Hartree) to the energy of the first excited $2s$ state of H&HLI. The numbers in the brackets indicate the uncertainties $\delta\varepsilon_{2}^{(4)}$ in the calculated ground-state energy due to the uncertainty in the nuclear radius $r_{N}$ in Table 1. $\Delta\varepsilon_{2,\nu}$ with $\nu=a,b$ is given in Eq. (64).  Isotope ($Z,A$) | $\varepsilon_{2}^{(4)}(\delta\varepsilon_{2}^{(4)})$ | $\Delta\varepsilon_{2,a}^{(4)}$ | $\Delta\varepsilon_{2,b}^{(4)}$
---|---|---|---
H(1, 1) | $2.1034(75)E^{-11}$ | $1.3E^{-14}$ | $1.3E^{-14}$
H(1, 2) | $1.34648(94)E^{-10}$ | $7.4E^{-14}$ | $7.4E^{-14}$
H(1, 3) | $9.20(38)E^{-11}$ | $5E^{-14}$ | $5E^{-14}$
He(2, 3) | $1.8396(56)E^{-9}$ | $4.1E^{-12}$ | $4.1E^{-12}$
He(2, 4) | $1.3361(45)E^{-9}$ | $3.0E^{-12}$ | $3.0E^{-12}$
K(19, 39) | $4.5755(51)E^{-5}$ | $7.626E^{-6}$ | $7.610E^{-6}$
K(19, 40) | $4.5841(75)E^{-5}$ | $7.639E^{-6}$ | $7.623E^{-6}$
K(19, 41) | $4.621(15)E^{-5}$ | $7.70E^{-6}$ | $7.68E^{-6}$
Rb(37, 85) | $9.855(11)E^{-4}$ | $6.877E^{-4}$ | $6.859E^{-4}$
Rb(37, 87) | $9.833(10)E^{-4}$ | $6.863E^{-4}$ | $6.846E^{-4}$
Cs(55, 133) | $6.285(12)E^{-3}$ | $1.2684E^{-2}$ | $1.2643E^{-2}$
Cs(55, 135) | $6292(12)E^{-3}$ | $1.2696E^{-2}$ | $1.2655E^{-2}$
Pb(82, 204) | $4.0409(21)E^{-2}$ | $0.38827$ | $0.38651$
Pb(82, 206) | $4.0555(21)E^{-2}$ | $0.38935$ | $0.38758$
Pb(82, 207) | $4.0616(21)E^{-2}$ | $0.38980$ | $0.38803$
Pb(82, 208) | $4.0718(19)E^{-2}$ | $0.39055$ | $0.38878$
Bi(83, 209) | $4.3050(41)E^{-2}$ | $0.43954$ | $0.43752$
U(92, 235) | $7.255(10)E^{-2}$ | $1.3154$ | $1.3088$
U(92, 236) | $7.2787(95)E^{-2}$ | $1.3184$ | $1.3118$
U(92, 238) | $7.3136(82)E^{-2}$ | $1.3229$ | $1.3163$

##  IV Perturbation Theory Corrections due to FNS

In this section we apply perturbation theory for models (a) and (b) and find
the corrections to the ground-state energies of H&HLI (assuming the potentials
$W_{\nu}(r)$ and $W_{ba}(r)$ are small).

To be specific, we will concentrate on the ground electronic state, $n=1$, and
therefore the subscript $n$ on energies and wave functions will often be
omitted below. Since an analytic solution to the Dirac equation for a harmonic
potential is not known, we use perturbation theory in order to obtain FNS
corrections for model (b), but we also numerically solve the Dirac equation
(see previous section).

Perturbation theory for models (a) and (b) with weak $W_{\nu}(r)$ is detailed
in Sec. IV.1. It starts from the Hamiltonian, $H_{\nu}=-i\hbar
c\bm{\alpha}\cdot\nabla+V_{\nu}(r)+\beta mc^{2}$, rewriten as,

$H_{\nu}=H_{0}+V_{\nu}(r)-V_{C}(r)\equiv H_{0}+W_{\nu}(r),$ (65)

where $H_{0}=-i\hbar c\bm{\alpha}\cdot\nabla+V_{C}(r)+\beta mc^{2}$ is the
Dirac Hamiltonian for a point-like nucleus with $V_{C}(r)=-Ze^{2}/r$.

Perturbation theory for model (b) alone with weak $W_{ba}(r)$ is detailed in
Sec. IV.2. It starts from the Hamiltonoian $H_{b}$ in the form,

$H_{b}=H_{a}+W_{ba}(r),$ (66)

where $H_{a}$ is given in Eq. (65), and $W_{ba}(r)=V_{b}(r)-V_{a}(r)$ is given
in Eq. (9).

###  IV.1 Perturbation theory with $W_{\nu}(r)$

The first-order correction to the energy is,

$\varepsilon^{(1)}_{n,\nu}=\mathcal{N}^{2}\int\limits_{0}^{r_{\nu}}W_{\nu}(r)\,\big{[}g^{2}(r)+f^{2}(r)\big{]}r^{2}dr,$
(67)

where $f(r)$ and $g(r)$ are defined in Eq. (25). Note that the singularities
at $r=0$ of $f(r)$ and $g(r)$ are integrable [29] and $W_{\nu}(r)$ is given by
Eq. (8). Explicitly,

$\displaystyle\varepsilon^{(1)}_{n,\nu}$ $\displaystyle=$
$\displaystyle\frac{\hbar^{3}c^{3}\mathcal{N}^{2}}{8\lambda^{3}}\int\limits_{0}^{\rho_{\nu}}W_{\nu}\bigg{(}\frac{\hbar
c\rho}{2\lambda}\bigg{)}\Big{\\{}m_{e}c^{2}\,\big{[}Q_{1}^{2}(\rho)+Q_{2}^{2}(\rho)\big{]}$
(68)
$\displaystyle+2\varepsilon\,Q_{1}(\rho)Q_{2}(\rho)\Big{\\}}\rho^{2\gamma}e^{-\rho}d\rho,$

where

$\displaystyle\mathcal{N}^{-2}$ $\displaystyle=$
$\displaystyle\frac{\hbar^{3}c^{3}}{8\lambda^{3}}\int\limits_{0}^{\infty}\Big{\\{}m_{e}c^{2}\,\big{[}Q_{1}^{2}(\rho)+Q_{2}^{2}(\rho)\big{]}$
(69)
$\displaystyle+2\varepsilon\,Q_{1}(\rho)Q_{2}(\rho)\Big{\\}}\rho^{2\gamma}e^{-\rho}d\rho,$

and the dimensionless nuclear radius is

$\rho_{\nu}=\frac{2\lambda r_{\nu}}{\hbar c}=\frac{2Zr_{\nu}}{a_{B}}.$ (70)

Equation (17) was used to obtain the second equality in Eq. (70). The nuclear
radii $r_{N}$ are tabulated in Table 1 for a number of elements and isotopes.
Let us restore the $n$ dependence in $\varepsilon^{(1)}_{n,\nu}$. The FNS
corrections $\varepsilon^{(1)}_{1,\nu}$ to the ground-state energy are given
in Table 7. The dependence of $\varepsilon^{(1)}_{1,\nu}$ on $A$ and $Z$
follows the behavior of the nuclear radius $r_{N}(A,Z)$; because $r_{N}$
increases with $Z$, so does $\varepsilon^{(1)}_{1,\nu}$. Dependence of
$\varepsilon^{(1)}_{1,\nu}$ on $A$ is less simple. For most of the elements in
Table 1, $r_{N}$ increases with $A$, and so does $\varepsilon^{(1)}_{1,\nu}$,
[see, e.g., $\varepsilon^{(1)}_{1,a}$ and $\varepsilon^{(1)}_{1,b}$ for
isotopes of K, Cs, Pb and U]. However, for some elements, $r_{N}$ decreases
with $A$, and so does $\varepsilon^{(1)}_{1,\nu}$, [see, e.g.,
$\varepsilon^{(1)}_{1,a}$ and $\varepsilon^{(1)}_{1,b}$ for isotopes of He and
Rb]. Moreover, for some elements, $r_{N}$ is not a monotonic function of $A$,
hence $\varepsilon^{(1)}_{1,\nu}$ is also not a monotonic function of $A$. For
example, the nuclear radii $r_{N}({}^{A}\text{H})$ for the isotopes
${}^{A}\text{H}$ satisfy the inequalities,

$r_{N}\big{(}{}^{1}\text{H}\big{)}<r_{N}\big{(}{}^{3}\text{H}\big{)}<r_{N}\big{(}{}^{2}\text{H}\big{)}.$

Thus $\varepsilon^{(1)}_{1,\nu}({}^{A}\text{H})$ satisfy

$\varepsilon^{(1)}_{1,\nu}\big{(}{}^{1}\text{H}\big{)}<\varepsilon^{(1)}_{1,\nu}\big{(}{}^{3}\text{H}\big{)}<\varepsilon^{(1)}_{1,\nu}\big{(}{}^{2}\text{H}\big{)}.$

Comparing $\varepsilon^{(1)}_{1,a}$ in Table 7 with $\varepsilon_{a}$ in Table
3 shows that for light isotopes, perturbation theory predictions are very
close to the non-perturbative numerical calculations, but for the heavy
isotopes, the differences between $\varepsilon_{1,a}$ and
$\varepsilon^{(1)}_{1,a}$ are large (see the difference between the values of
$\varepsilon_{1,a}$ in Table 3 and $\varepsilon^{(1)}_{1,a}$ in Table 7).

Table 7: First-order perturbation theory correction due to FNS in Eq. (68), for $\varepsilon^{(1)}_{1,\nu}$ [with $\nu=a,b$]. The nuclear radii $r_{a}$ and $r_{b}$ are given in Eq. (7), and the RMS nuclear charge radii $r_{N}$ are listed in Table 1. Isotope ($Z,A$) |  $\varepsilon^{(1)}_{1,a}$ (Hartree) |  $\varepsilon^{(1)}_{1,b}$ (Hartree)
---|---|---
H (1, 1) | $1.683743E^{-10}$ | $1.683736E^{-10}$
H (1, 2) | $1.077785E^{-9}$ | $1.077778E^{-9}$
H (1, 3) | $7.36714E^{-10}$ | $7.36710E^{-10}$
He (2, 3) | $1.47484E^{-8}$ | $1.47482E^{-8}$
He (2, 4) | $1.07128E^{-8}$ | $1.07126E^{-8}$
K (19, 39) | $4.24722E^{-4}$ | $4.24232E^{-4}$
K (19, 40) | $4.25506E^{-4}$ | $4.25015E^{-4}$
K (19, 41) | $4.28869E^{-4}$ | $4.28374E^{-4}$
Rb (37, 85) | $1.31164E^{-2}$ | $1.30612E^{-2}$
Rb (37, 87) | $1.30882E^{-2}$ | $1.30331E^{-2}$
Cs (55, 133) | $0.145819$ | $0.144503$
Cs (55, 135) | $0.145964$ | $0.144646$
Pb (82, 204) | $3.25838$ | $3.19659$
Pb (82, 206) | $3.26777$ | $3.20580$
Pb (82, 207) | $3.27167$ | $3.20962$
Pb (82, 208) | $3.27822$ | $3.21605$
Bi(83, 209) | $3.67177$ | $3.60062$
U (92, 235) | $10.7543$ | $10.5056$
U (92, 236) | $10.7799$ | $10.5305$
U (92, 238) | $10.8180$ | $10.5677$

In the non-relativistic approximation, the perturbation theory of the FNS
correction, $\varepsilon^{(1)}_{n,\nu}$, takes the form of the simple equation
[11, 10]

$\varepsilon^{(1)}_{n,\nu}\approx\frac{2m_{e}c^{2}}{3}\,\frac{(Z\alpha)^{4}}{n^{3}}\,\frac{m^{3}}{m_{e}^{3}}\,\frac{r_{N}^{2}}{\lambdabar_{C}^{2}}.$
(71)

We verified that this is an excellent approximation for hydrogen. The
perturbation theory result in Eq. (71) is obtained from Eq. (68) by using the
approximations: $Q_{1}(\rho)\approx Q_{1}(0)$, $Q_{2}(\rho)\approx Q_{2}(0)$
and $\gamma\approx 1$. The first and second approximations are based on the
inequalities $0<\rho<\rho_{\nu}\ll 1$ [see Eq. (68)] which are valid for H and
all HLIs. The third approximation is based on the inequality $Z\alpha\ll 1$
which is satisfied for H and light HLIs. For the heavy HLIs, $Z\alpha$ is not
a small parameter, and Eq. (71) should be modified. The FNS correction to the
ground state energy of all HLIs is

$\displaystyle\varepsilon^{(1)}_{1,a}$ $\displaystyle\approx$
$\displaystyle\frac{2m_{e}c^{2}}{\gamma(2\gamma+1)}\,\frac{(Z\alpha)^{4}}{n^{3}}\,\frac{m^{3}}{m_{e}^{3}}\,\frac{r_{N}^{2\gamma}}{\lambdabar_{C}^{2\gamma}},$
(72a) $\displaystyle\varepsilon^{(1)}_{1,b}$ $\displaystyle\approx$
$\displaystyle\frac{5^{\gamma}\,3^{1-\gamma}\,2\,m_{e}c^{2}}{\gamma(2\gamma+1)(2\gamma+3)}\,\frac{(Z\alpha)^{4}}{n^{3}}\,\frac{m^{3}}{m_{e}^{3}}\,\frac{r_{N}^{2\gamma}}{\lambdabar_{C}^{2\gamma}}.$
(72b)

Note that the FNS correction to the $1s$ state energy of H and light HLIs in
Eq. (71) depends just on the nuclear RMS charge radius $r_{N}$, and not on the
details of the nuclear charge distribution. However, for heavy HLIs with high
$Z$, the FNS correction to the $s$ state energy is sensitive to the details of
the nuclear charge distribution, and Eqs. (72a) and (72b) give different FNS
corrections for the models (a) and (b). Therefore, it is important to develop
and use more accurate models for the nuclear charge distribution in order to
obtain better comparisons with the experimental energies of the $s$ state and
the frequencies of the $1s$-$2s$ and $1s$-$3s$ transitions.

###  IV.2 Perturbation theory with $W_{ba}(r)$

The first order correction $\varepsilon_{ba}^{(1)}$ is given by,

$\displaystyle\varepsilon_{ba}^{(1)}$ $\displaystyle=$
$\displaystyle\frac{4\pi
A_{-}^{2}}{\mathcal{N}_{a}}\int\limits_{0}^{r_{a}}W_{ba}(r)\Big{[}g_{-}^{2}(r)+f_{-}^{2}(r)\Big{]}r^{2}dr+$
(73) $\displaystyle\frac{4\pi
A_{+}^{2}}{\mathcal{N}_{a}}\int\limits_{r_{a}}^{r_{b}}W_{ba}(r)\Big{[}g_{+}^{2}(r)+f_{+}^{2}(r)\Big{]}r^{2}dr,$

where $g_{\pm}(r)$ and $f_{\pm}(r)$ are the radial wave functions for the
model (a). The functions $g_{+}(r)$ and $f_{+}(r)$ are given in Eq. (38), and
$g_{-}(r)$ and $f_{-}(r)$ are given in Eq. (47). The constants $A_{\pm}$ are
found from Eq. (27), and the normalization constant $\mathcal{N}_{a}$ is,

$\displaystyle\mathcal{N}_{a}$ $\displaystyle=$ $\displaystyle 4\pi
A_{-}^{2}\int\limits_{0}^{r_{a}}\Big{[}g_{-}^{2}(r)+f_{-}^{2}(r)\Big{]}r^{2}dr+$
(74) $\displaystyle 4\pi
A_{+}^{2}\int\limits_{r_{a}}^{\infty}\Big{[}g_{+}^{2}(r)+f_{+}^{2}(r)\Big{]}r^{2}dr.$

In order to calculate $\varepsilon_{ba}^{(1)}$, it is useful to rewrite Eq.
(73) employing the electron density $\varrho_{a}(r)$ given by,

$\varrho_{a}(r)=\left\\{\begin{array}[]{cc}\displaystyle\frac{4\pi
A_{-}^{2}}{\mathcal{N}}\,\Big{[}g_{-}^{2}(r)+f_{-}^{2}(r)\Big{]},&r<r_{a},\\\
\displaystyle\frac{4\pi
A_{+}^{2}}{\mathcal{N}}\,\Big{[}g_{+}^{2}(r)+f_{+}^{2}(r)\Big{]},&r>r_{a},\end{array}\right.$
(75)

hence Eq. (73) can be written as

$\varepsilon_{ba}^{(1)}=\int\limits_{0}^{r_{b}}W_{ba}(r)\varrho_{a}(r)r^{2}dr.$
(76)

The two-component wave function in Eq. (14) is continuous at $r=r_{a}$, i.e.,

$A_{-}\,\big{(}g_{-}(r_{a}),\,f_{-}(r_{a})\big{)}=A_{+}\,\big{(}g_{+}(r_{a}),\,f_{+}(r_{a})\big{)},$

and thus the electron density $\varrho_{a}(r)$ in Eq. (75) is also continuous.
Moreover, $\varrho_{a}(r)$ is smooth within the interval $0<r<r_{b}$, i.e.,
for any $r\in(0,r_{b})$, the following inequality is satisfied:
$\big{|}\varrho(r)-\varrho(r_{a})\big{|}\ll\varrho(r_{a})$. Note that,
although one may expect that Eq. (76) can be approximated as

$\varepsilon_{ba}^{(1)}\approx\varrho(r_{a})\int\limits_{0}^{r_{b}}W_{ba}(r)r^{2}dr,$
(77)

it can be easily verified that $\int\limits_{0}^{r_{b}}W_{ba}(r)r^{2}dr=0$,
hence the approximation in Eq. (77) is not valid.

Table 8: Perturbation theory corrections $\varepsilon_{ba}^{(1)}$ [in Hartree] to the ground state of H&HLI systems [see Eq. (78)]. The dimensionless parameter $\delta_{ba}$ is defined in Eq. (79). The energy correction $\varepsilon_{1,ba}$ and its uncertainty $\delta\varepsilon_{1,ba}$ due to the uncertainty in the nuclear RMS nuclear charge radius [in units of Hartree] are calculated non-perturbatively from Eqs. (80) and (81). Isotope ($Z,A$) | $\varepsilon_{ba}^{(1)}$ | $\delta_{ba}$ | $\varepsilon_{1,ba}(\delta\varepsilon_{1,ba})$
---|---|---|---
H (1, 1) | $-2.39058E^{-16}$ | $-4.78371E^{-16}$ | $-3.26(15)E^{-16}$
H (1, 2) | $-2.37631E^{-15}$ | $-4.75385E^{-15}$ | $-2.48(3)E^{-15}$
H (1, 3) | $-1.91704E^{-15}$ | $-3.83472E^{-15}$ | $-1.57(11)E^{-15}$
He (2, 3) | $-8.14903E^{-14}$ | $-4.07504E^{-14}$ | $-8.59(4)E^{-14}$
He (2, 4) | $-5.4718E^{-14}$ | $-2.73613E^{-14}$ | $-5.80(3)E^{-14}$
K (19, 39) | $-1.20571E^{-7}$ | $-6.64768E^{-10}$ | $-1.2837(15)E^{-7}$
K (19, 40) | $-1.20812E^{-7}$ | $-6.66093E^{-10}$ | $-1.286(2)E^{-7}$
K (19, 41) | $-1.21844E^{-7}$ | $-6.71783E^{-10}$ | $-1.297(4)E^{-7}$
Rb (37, 85) | $-1.24451E^{-5}$ | $-1.78437E^{-8}$ | $-1.3283(15)E^{-5}$
Rb (37, 87) | $-1.24169E^{-5}$ | $-1.78034E^{-8}$ | $-1.3254(13)E^{-5}$
Cs (55, 133) | $-2.65761E^{-4}$ | $-1.68324E^{-7}$ | $-2.842(5)E^{-4}$
Cs (55, 135) | $-2.66035E^{-4}$ | $-1.68497E^{-7}$ | $-2.845(5)E^{-4}$
Pb (82, 204) | $-9.41203E^{-3}$ | $-2.52128E^{-6}$ | $-1.0096(4)E^{-2}$
Pb (82, 206) | $-9.43977E^{-3}$ | $-2.52871E^{-6}$ | $-1.0125(4)E^{-2}$
Pb (82, 207) | $-9.45126E^{-3}$ | $-2.53179E^{-6}$ | $-1.0138(4)E^{-2}$
Pb (82, 208) | $-947062E^{-3}$ | $-2.53698E^{-6}$ | $-1.0159(4)E^{-2}$
Bi (83, 209) | $-1.06803E^{-2}$ | $-2.78397E^{-6}$ | $-1.1458(9)E^{-2}$
U (92, 235) | $-3.20879E^{-2}$ | $-6.60084E^{-6}$ | $-3.447(4)E^{-2}$
U (92, 236) | $-3.21652E^{-2}$ | $-6.61674E^{-6}$ | $-3.455(3)E^{-2}$
U (92, 238) | $-3.22805E^{-2}$ | $-6.64045E^{-6}$ | $-3.467(3)E^{-2}$

The correction $\varepsilon_{ba}^{(1)}$ in Eq. (76) can be written as

$\varepsilon_{ba}^{(1)}=\int\limits_{0}^{r_{b}}W_{ba}(r)\,\big{[}\varrho(r)-\varrho(0)\big{]}\,r^{2}dr.$
(78)

The correction to the ground state energy of model (b) is given in Eq. (78)
and tabulated in Table 8, which also tabulates the dimensionless corrections
to the ground-state binding energy,

$\delta_{ba}=\frac{\varepsilon_{ba}^{(1)}}{mc^{2}-\varepsilon_{1}^{(0)}},$
(79)

where $\varepsilon_{1}^{(0)}$ is given in Eq. (23). For comparison, Table 8
shows $\varepsilon_{1,ba}$ and its uncertainty $\delta\varepsilon_{1,ba}$
calculated non-perturbatively from the equations,

$\displaystyle\varepsilon_{1,ba}$ $\displaystyle=$
$\displaystyle\varepsilon_{1,b}-\varepsilon_{1,a},$ (80)
$\displaystyle\delta\varepsilon_{1,ba}$ $\displaystyle=$
$\displaystyle\bigg{|}\frac{\partial\varepsilon_{1,ba}}{\partial
r_{N}}\bigg{|}\,\delta r_{N}$ (81) $\displaystyle=$
$\displaystyle\Big{|}\delta\varepsilon_{1,b}-\delta\varepsilon_{1,a}\Big{|},$

where $\delta\varepsilon_{1,a}$ and $\delta\varepsilon_{1,b}$ are given by

$\delta\varepsilon_{1,a}=\bigg{|}\frac{\partial\varepsilon_{1,a}}{\partial
r_{N}}\bigg{|}\,\delta
r_{N},\quad\delta\varepsilon_{1,b}=\bigg{|}\frac{\partial\varepsilon_{1,b}}{\partial
r_{N}}\bigg{|}\,\delta r_{N}.$

$\varepsilon_{1,a}$, $\delta\varepsilon_{1,a}$ and $\varepsilon_{1,b}$,
$\delta\varepsilon_{1,b}$ are tabulated in Table 3. Note that the first-order
perturbation correction $\varepsilon_{ba}^{(1)}$ is close to the non-
perturbative correction $\varepsilon_{1,ba}$ for all the isotopes. In other
words, we conclude that the first-order perturbation correction
$\varepsilon_{ba}^{(1)}$ gives accurate results for both light and heavy
isotopes.

In summary, for light nuclei, perturbation theory gives results very close to
the non-perturbative numerical calculations for the difference between
$\varepsilon_{1,a}$ and $\varepsilon^{(1)}_{1,a}$, but for heavy isotopes the
difference is large (see the large differences for the heavy isotopes between
the values of $\varepsilon_{1,a}$ in Table 3 analytically calculated using
model (a) and $\varepsilon^{(1)}_{1,a}$ calculated by perturbation theory
using model (a) in Table 7).

##  V Two-parameter Fermi model for the nuclear charge distribution

In model (c) the distribution of charge within the nucleus is taken to be the
TPFM [24]. References [21, 22, 25] studied the FNS corrections to the
hyperfine splitting of the ground state of H&HLI using a TPFM to parameterize
the distribution of the nuclear charge. The TPFM is more realistic than either
model (a) or (b), at least for heavier nuclei; but it is unclear how realistic
it is for light nuclei. The TPFM charge distribution is,

$\varrho_{F}(r)=\frac{Ze}{8\pi
a_{F}^{3}\big{|}\mathrm{Li}_{3}(-e^{c_{F}/a_{F}})\big{|}}\,\bigg{[}\exp\bigg{(}\frac{r-c_{F}}{a_{F}}\bigg{)}+1\Bigg{]}^{-1},$
(82)

where $\mathrm{Li}_{n}(\bullet)$ is the polylogarithm function [47] of order
$n$. Here $c_{F}$ is the half-density nuclear radius, and $a_{F}$ is the
distribution thickness parameter. The charge density is normalized as follows:

$4\pi\int\limits_{0}^{\infty}\varrho_{F}(r)r^{2}dr=Ze.$

The nuclear charge radius is $r_{N}=\sqrt{\langle r^{2}\rangle_{N}}$ is
obtained by the expression,

$\displaystyle r_{N}^{2}\equiv\langle r^{2}\rangle_{N}$ $\displaystyle=$
$\displaystyle\frac{4\pi}{Ze}\int\limits_{0}^{\infty}\varrho_{F}(r)r^{4}dr$
(83) $\displaystyle=$ $\displaystyle
12a_{F}^{2}\,\frac{\mathrm{Li}_{5}(-e^{c_{F}/a_{F}})}{\mathrm{Li}_{3}(-e^{c_{F}/a_{F}})}.$

Note that the arguments of the polylogarithm functions contain an exponential
function with argument $c_{F}/a_{F}$. We find $c_{F}$ as a function of $r_{N}$
and $a_{F}$ using Eq. (83). The uncertainty $\delta c_{F}$ in $c_{F}$ due to
the uncertainties $\delta r_{N}$ and $\delta a_{F}$ in $r_{N}$ and $a_{F}$ is,

$\delta c_{F}=\sqrt{\bigg{(}\frac{\partial c_{F}}{\partial r_{N}}\,\delta
r_{N}\bigg{)}^{2}+\bigg{(}\frac{\partial c_{F}}{\partial a_{F}}\,\delta
a_{F}\bigg{)}^{2}}.$ (84)

The nuclear charge radii $r_{N}$ and their uncertainties $\delta r_{N}$ for
many nuclei are published in Ref. [24], and the parameters $a_{F}$ and the
uncertainties $\delta a_{F}$ in $a_{F}$ are published in Ref. [27]. As
$a_{F}\to 0$, the charge distribution inside the nucleus becomes,

$\lim_{a_{F}\to 0}\varrho_{F}(r)=\frac{3Ze}{4\pi
r_{b}^{3}}\,\Theta(r_{b}-r),\quad\mathrm{and\,\,\,}r_{b}=\lim_{a_{F}\to
0}c_{F}.$

Figure 2 shows the charge density distribution versus $r/r_{N}$ for three
values of the width parameter, $a_{F}=0.02\,r_{N},0.1\,r_{N},0.15\,r_{N}$
(red, green and blue, respectively). Clearly, as $a_{F}\to 0$, the charge
distribution approaches the limit given above. Moreover, in the limit as
$a_{F}\to 0$, $c_{F}\to\sqrt{5/3}\,r_{N}$. The dashed black curve in Fig. 2
shows the charge density in this limit.

The TPFM charge distribution yields the following potential (for simplicity we
have used $c$ and $a$ for $a_{F}$ and $c_{F}$ in this equation):

$\displaystyle
V_{F}(r)=-\frac{Ze^{2}}{r}+\frac{4\pi}{r}\,\int\limits_{r}^{\infty}\varrho_{F}(r^{\prime})\,r^{\prime}(r^{\prime}-r)\,dr^{\prime}=$
$\displaystyle
Ze^{2}\left[\frac{\text{Li}_{3}\left(-e^{\frac{c}{a}-\frac{r}{a}}\right)}{r\,\text{Li}_{3}\left(-e^{\frac{c}{a}}\right)}+\frac{\text{Li}_{2}\left(-e^{\frac{c}{a}-\frac{r}{a}}\right)}{2a\,\text{Li}_{3}\left(-e^{\frac{c}{a}}\right)}-\frac{1}{r}\right].$
(85)

Table 9: The half-density nuclear radii $c_{F}$ (in Bohr), the distribution thickness parameter $a_{F}$ (in Bohr), and the uncertainties $\delta c_{F}$ and $\delta a_{F}$ in $c_{F}$ and $a_{F}$ for a number of heavy isotopes. $a_{F}$ and $\delta a_{F}$ are taken from Ref. [24], $c_{F}$ is calculated from Eq. (83), $\delta c_{F}$ is calculated from Eq. (84), and $r_{N}$ and $\delta r_{N}$ are tabulated in Table 1. Isotope ($Z,A$) |  $c_{F}$ ($\delta c_{F}$) |  $a_{F}$ ($\delta a_{F}$)
---|---|---
Pb (82, 206) | $1.2449(29)E^{-4}$ | $1.030(15)E^{-5}$
Pb (82, 207) | $1.2456(36)E^{-4}$ | $1.032(19)E^{-5}$
Bi (83, 209) | $1.278(12)E^{-4}$ | $8.84(74)E^{-6}$
U (92, 238) | $1.3194(61)E^{-4}$ | $1.143(30)E^{-5}$

We numerically solve Dirac equation (15) for the heavy nuclei 206Pb, 207Pb,
209Bi, and 238U to obtain radial wave functions $g_{n}(r)$ and $f_{n}(r)$ and
eigen-energies $\varepsilon_{n}$ with potential (85) with $c_{F}$ and $a_{F}$
as tabulated in Table 9. The uncertainty $\delta\varepsilon_{n}$ in the energy
$\varepsilon_{n}$ of the $ns$ state can be estimated as follows,

$\delta\varepsilon_{n}=\sqrt{\Big{(}\delta\varepsilon_{n}^{(r)}\Big{)}^{2}+\Big{(}\delta\varepsilon_{n}^{(a)}\Big{)}^{2}},$
(86)

where $\delta\varepsilon_{n}^{(r)}$ is the uncertainty in $\varepsilon_{n}$
due to the uncertainty in $r_{N}$, and $\delta\varepsilon_{n}^{(a)}$ is the
uncertainty in $\varepsilon_{n}$ due to the uncertainty in $a_{F}$.
$\delta\varepsilon_{n}^{(r)}$ and $\delta\varepsilon_{n}^{(a)}$ can be
approximated as,

$\displaystyle\delta\varepsilon_{n}^{(r)}$ $\displaystyle\approx$
$\displaystyle\frac{1}{2}\bigg{|}\varepsilon_{n}\Big{(}c_{F}(r_{N}+\delta
r_{N},a_{F}),a_{F}\Big{)}$
$\displaystyle-\varepsilon_{n}\Big{(}c_{F}(r_{N}-\delta
r_{N},a_{F}),a_{F}\Big{)}\bigg{|},$ $\displaystyle\delta\varepsilon_{n}^{(a)}$
$\displaystyle\approx$
$\displaystyle\frac{1}{2}\bigg{|}\varepsilon_{n}\Big{(}c_{F}(r_{N},a_{F}+\delta
a_{F}),a_{F}+\delta a_{F}\Big{)}$
$\displaystyle-\varepsilon_{n}\Big{(}c_{F}(r_{N},a_{F}-\delta
a_{F}),a_{F}-\delta a_{F}\Big{)}\bigg{|},$

where $\varepsilon_{n}(c_{F}(r_{N}\pm\delta r_{N},a_{F}),a_{F})$ is the $ns$
state energy for the nuclear RMS charge radius $r_{N}\pm\delta r_{N}$, and
$\varepsilon_{n}(c_{F}(r_{N},a_{F}\pm\delta a_{F}),a_{F}\pm\delta a_{F})$ is
the $ns$ state energy for the distribution thickness parameter $a_{F}\pm\delta
a_{F}$. The FNS energy correction $\varepsilon_{n,c}$ is

$\varepsilon_{n,c}=\varepsilon_{n}-\varepsilon_{n}^{(0)},$ (87)

where $\varepsilon_{n}^{(0)}$ is the energy of the $ns$ state for a point-like
nucleus. $\varepsilon_{n,c}$ and $\delta\varepsilon_{n}$ are tabulated in
Table 10 for heavy isotopes. Comparing $\varepsilon_{n,c}$ with
$\varepsilon_{n,a}$ and $\varepsilon_{n,b}$ in Tables 3 and 4, one sees that
$\varepsilon_{n,a}>\varepsilon_{n,b}>\varepsilon_{n,c}>0$ and
$\varepsilon_{n,a}-\varepsilon_{n,b}>\varepsilon_{n,b}-\varepsilon_{n,c}$. For
206Pb, 207Pb and 238U, the uncertainty $\delta\varepsilon_{n}$ in
$\varepsilon_{n}$ is the same for the models (a), (b) and (c). This is because
$\delta\varepsilon_{n}^{(a)}\ll\varepsilon_{n}^{(r)}$, see Eq. (86). For
209Bi, $\delta\varepsilon_{n}^{(a)}$ and $\varepsilon_{n}^{(r)}$ are of the
same order of magnitude, and the uncertainty $\delta\varepsilon_{n}$ for the
model (c) is larger than the uncertainty for models (a) and (b).

Table 10: The FNS corrections $\varepsilon_{n,c}$ to the energies $\varepsilon_{n}$ of the $ns$ states in Eq. (87), and the uncertainties $\delta\varepsilon_{n}$ in $\varepsilon_{n}$ in Eq. (86). Isotope ($Z,A$) |  $\varepsilon_{1,c}$ ($\delta\varepsilon_{1}$) |  $\varepsilon_{2,c}$ ($\delta\varepsilon_{2}$)
---|---|---
Pb (82, 206) | $2.4605(10)$ | $0.42724(17)$
Pb (82, 207) | $2.4634(10)$ | $0.42774(17)$
Bi (83, 209) | $2.740(29)$ | $0.4798(50)$
U (92, 238) | $7.2952(60)$ | $1.3858(11)$

##  VI Ground state hyperfine splitting

Within the Dirac equation formalism, the hyperfine interaction between the
magnetic moments of the electron (located at space point ${\bf r}$) and the
nucleus (located at the origin) can be expressed in terms of the vector
potential $\mathbf{A}(\mathbf{r})$. The interaction Hamiltonian is

$H_{\mathrm{mag}}=[c]\,e\,\bm{\alpha}\cdot\mathbf{A}(\mathbf{r}),$ (88)

where the factor $[c]$ is present in SI units but is unity in Gaussian units,
and $\mathbf{A}(\mathbf{r})$ is the vector potential generated by the magnetic
moment of the nucleus with magnetic moment density
$\bm{\mu}(\mathbf{r}^{\prime})$,

$\mathbf{A}_{\nu}(\mathbf{r})=\left[{\frac{{\mu_{0}}}{{4\pi}}}\right]\int\frac{{\bm{\mu}_{\nu}(\mathbf{r}^{\prime})}\times(\mathbf{r}-\mathbf{r}^{\prime})}{|\mathbf{r}-\mathbf{r}^{\prime}|^{3}}\,d^{3}\mathbf{r}^{\prime},$
(89)

where integration is over all space. The subscript $\nu$ specifies the models
for the magnetic moment distribution: $\nu=0,a,b,c$ are the point-like,
uniform surface, uniform volume, and TPFM. Here the factor
$\left[{\frac{{\mu_{0}}}{{4\pi}}}\right]$ is present in SI units but is unity
in Gaussian units. We restrict ourselves by considering isotropic models where
the magnetic moment density $\bm{\mu}(r^{\prime})$ depends only on the
distance $r^{\prime}=|\mathbf{r}^{\prime}|$ from the center of mass of the
nucleus. For the point-like distribution, $\nu=0$,

$\bm{\mu}_{0}(r)=g_{I}\mu_{N}\mathbf{I}\,\delta^{3}(\mathbf{r}),$ (90)

where $\mu_{N}$ is the nuclear magneton, $\mu_{N}=e\hbar/(2m_{p}c)$ in
Gaussian units, $m_{p}$ is the proton mass, $g_{I}$ is the nuclear $g$-factor,
$\mathbf{I}$ is the nuclear spin operator, and $\delta^{3}(\mathbf{r})$ is the
three dimensional Dirac delta function. For $\nu=a$,

$\bm{\mu}_{a}(r)=\frac{g_{I}\mu_{N}}{4\pi
r_{a}^{2}}\,\mathbf{I}\,\delta(r-r_{a}),$ (91)

where $\delta(\bullet)$ is the Dirac $\delta$ function. For $\nu=b$,

$\bm{\mu}_{b}(r)=\frac{3g_{I}\mu_{N}}{4\pi
r_{b}^{3}}\,\mathbf{I}\,\Theta(r_{b}-r).$ (92)

For the Fermi model, $\nu=c$,

$\bm{\mu}_{c}(r)=\frac{\vartheta_{0}}{e^{(r-c_{F})/a_{F}}+1}\,g_{I}\mu_{N}\mathbf{I}.$
(93)

The factor $\vartheta_{0}$ is determined from the normalization condition,

$4\pi\vartheta_{0}\int\limits_{0}^{\infty}\frac{r^{2}dr}{e^{(r-c_{F})/a_{F}}+1}=1.$

This yields,

$\vartheta_{0}=-\frac{1}{8\pi a_{F}^{3}\mathrm{Li}_{3}(-e^{c_{F}/a_{F}})}.$

Substituting $\bm{\mu}_{\nu}(r^{\prime})$ into Eq. (89), we get,

$\mathbf{A}_{\nu}(\mathbf{r})=\left[{\frac{{\mu_{0}}}{{4\pi}}}\right]\frac{g_{I}\mu_{N}}{r^{2}}\,\big{[}\mathbf{I}\times\mathbf{e}_{r}\big{]}\,\mathcal{M}_{\nu}(r),$
(94)

where

$\displaystyle\mathcal{M}_{0}(r)$ $\displaystyle=$ $\displaystyle 1,$ (95)
$\displaystyle\mathcal{M}_{a}(r)$ $\displaystyle=$
$\displaystyle\Theta(r-r_{a}),$ (96) $\displaystyle\mathcal{M}_{b}(r)$
$\displaystyle=$
$\displaystyle\frac{r^{3}}{r_{b}^{3}}\,\Theta(r_{b}-r)+\Theta(r-r_{b}).$ (97)
$\displaystyle\mathcal{M}_{c}(r)$ $\displaystyle=$ $\displaystyle
4\pi\vartheta_{0}\int\limits_{0}^{r}\frac{{r^{\prime}}^{2}dr^{\prime}}{e^{(r^{\prime}-c_{F})/a_{F}}+1}$
(98) $\displaystyle=$
$\displaystyle\frac{1}{\mathrm{Li}_{3}\big{(}-e^{c_{F}/a_{F}}\big{)}}\,\Bigg{[}\frac{r^{2}}{2a^{2}}\ln\big{(}1+e^{(c_{F}-r)/a_{F}}\big{)}$
$\displaystyle+\mathrm{Li}_{3}\big{(}-e^{c_{F}/a_{F}}\big{)}-\mathrm{Li}_{3}\big{(}-e^{(c_{F}-r)/a_{F}}\big{)}$
$\displaystyle-\frac{r}{a_{F}}\mathrm{Li}_{2}\big{(}-e^{(c_{F}-r)/a_{F}}\big{)}\Bigg{]}.$

The function $\mathcal{M}_{\nu}(r)$ is plotted versus $r$ in Fig. 7 for the
four magnetic moment distribution models. Substituting Eq. (94) into Eq. (88),
we obtain

$H_{\mathrm{hf}}^{(\nu)}=\left[\frac{\mu_{0}c}{4\pi}\right]\frac{eg_{I}\mu_{N}}{r^{2}}\,\mathcal{M}_{\nu}(r)\,\Big{(}\mathbf{I}\cdot\big{[}\mathbf{e}_{r}\times\bm{\alpha}\big{]}\Big{)}.$
(99)

Figure 7: $\mathcal{M}_{\nu}(r)$ versus $r$ for 207Pb nucleus for different
magnetic moment distribution models: The point-like nucleus, $\nu=0$ (red
dotted curve), distributed homogeneously on a sphere of radius $r_{a}$,
$\nu=a$ (green dot-dashed curve), distributed homogeneously inside a sphere of
radius $r_{b}$, $\nu=b$ (dashed blue curve), and the TPFM in Eq. (98), $\nu=c$
(solid purple curve).

Assuming that the hyperfine splitting of the ground state is much smaller than
the excitation energy of the $2p$ state, we can restrict ourselves to first-
order perturbation theory,

$\epsilon_{\mathrm{hf}}^{(\nu)}(F)=\big{\langle}\Psi_{F,m_{F}}^{(\nu)}\big{|}H_{\mathrm{hf}}^{(\nu)}\big{|}\Psi_{F,m_{F}}^{(\nu)}\big{\rangle},$
(100)

where the ground-state wave function $|\Psi_{F,m_{F}}\rangle$ is

$\displaystyle\langle{\bf r}\big{|}\Psi_{F,m_{F}}^{(\nu)}\big{\rangle}$
$\displaystyle=$
$\displaystyle\sum_{m_{J},m_{F}}\big{\langle}\frac{1}{2},m_{J};\,I,m_{I}\big{|}F,m_{F}\big{\rangle}$
(101) $\displaystyle\times\big{\langle}{\bf
r}\big{|}\nu;\,\frac{1}{2},m_{J};\,I,m_{I}\big{\rangle},$

$|\nu;\frac{1}{2},m_{J};I,m_{I}\rangle$ is the atomic angular ket with total
electronic angular momentum vector ${\bf J}$ ($J=1/2$) and projection
$m_{J}=\pm\tfrac{1}{2}$, nuclear spin ${\bf I}$, and projection $m_{I}$, and

$\big{\langle}\mathbf{r}\big{|}\nu;\,\frac{1}{2},m_{J};\,I,m_{I}\big{\rangle}=\chi_{I,m_{I}}\,\psi_{m_{J}}^{(\nu)}(\mathbf{r}),$

$\chi_{I,m_{I}}$ is a nuclear spin wave function, and
$\psi_{m_{J}}^{(\nu)}(\mathbf{r})$ is a spatial wave function,

$\psi_{m_{J}}^{(\nu)}(\mathbf{r})=\left(\begin{array}[]{c}g_{\nu}(r)\Omega_{\frac{1}{2}0m_{J}}\\\
if_{\nu}(r)\Omega_{\frac{1}{2}1m_{J}}\end{array}\right).$ (102)

For $\nu=0$, the radial wave functions $g_{0}(r)$ and $f_{0}(r)$ are given in
Eq. (25). For $\nu=a$, the radial wave functions $g_{a}(r)$ and $f_{a}(r)$ are

$\displaystyle g_{a}(r)$ $\displaystyle=$ $\displaystyle
g_{+}(r)\Theta(r-r_{a})+g_{-}(r)\Theta(r_{a}-r),$ (103a) $\displaystyle
f_{a}(r)$ $\displaystyle=$ $\displaystyle
f_{+}(r)\Theta(r-r_{a})+f_{-}(r)\Theta(r_{a}-r),$ (103b)

where $g_{+}(r)$ and $f_{+}(r)$ are given in Eq. (29), and $g_{-}(r)$ and
$f_{-}(r)$ are given in Eq. (47). For $\nu=b$, the wave functions
$(g_{b}(r),f_{b}(r))$ are numerically calculated as described in Sec. III.3,
and for $\nu=c$ the wave functions must be numerically calculated using the
Dirac equation for a potential which includes a two-parameter Fermi nuclear
model distribution.

Substituting Eq. (101) into Eq. (100), we find
$\epsilon_{\mathrm{hf}}^{(\nu)}(F)$,

$\epsilon_{\mathrm{hf}}^{(\nu)}(F)=\frac{\Delta_{\nu}}{2}\,\big{[}F(F+1)-I(I+1)-S(S+1)\big{]},$
(104)

where $S=1/2$. The hyperfine splitting
$\Delta_{\nu}=\epsilon_{\mathrm{hf}}^{(\nu)}(I+1/2)-\epsilon_{\mathrm{hf}}^{(\nu)}(I-1/2)$
is

$\displaystyle\Delta_{\nu}$ $\displaystyle=$
$\displaystyle-\left[\frac{\mu_{0}c}{4\pi}\right]\Big{(}I+\frac{1}{2}\Big{)}\,\frac{8eg_{I}\mu_{N}}{3}$
(105)
$\displaystyle\times\int\limits_{0}^{\infty}\mathcal{M}_{\nu}(r)g_{\nu}(r)f_{\nu}(r)dr.$

The expression for the hyperfine splitting of the ground state of H&HLI given
in Ref. [21] is equivalent to

$\displaystyle\Delta$ $\displaystyle=$
$\displaystyle\frac{4mc^{2}}{3}\,\frac{\alpha(Z\alpha)^{3}}{\gamma(2\gamma-1)}\,g_{I}\,\frac{m_{e}}{m_{p}}\,\Big{(}I+\frac{1}{2}\Big{)}$
(106)
$\displaystyle\times\bigg{[}(1-\xi)(1-\eta)+\zeta_{\mathrm{QED}}\bigg{]}.$

Note the use of the electron reduced mass $m$ in Eq. (106), as discussed in
Refs. [30, 31, 32, 33]. In the following we shall write Eq. (106) in the form

$\Delta=\Delta_{0}\bigg{[}(1-\xi)(1-\eta)+\zeta_{\mathrm{QED}}\bigg{]},$ (107)

where

$\Delta_{0}=\frac{4mc^{2}}{3}\,\alpha\big{(}Z\alpha\big{)}^{3}\,g_{I}\,\frac{m_{e}}{m_{p}}\,\Big{(}I+\frac{1}{2}\Big{)}\,\frac{1}{\gamma(2\gamma-1)}.$
(108)

$\Delta_{0}$ is the well-known uncorrected relativistic hyperfine splitting
[21]. Here $\xi$ and $\eta$ are the FNS corrections to the hyperfine splitting
due to the charge distribution and the magnetic moment distribution (Bohr-
Weisskopf correction) respectively [see Eqs. (111), (115) and (119), and the
text below for models (a), (b) and (c)], and the quantity
$\zeta_{\mathrm{QED}}$ is the QED radiative correction to the hyperfine
splitting. Relativistic QED radiative corrections are addressed in Ref. [21],
where a quantity $\chi_{\mathrm{QED}}$ is defined. We prefer to present the
QED radiative correction in terms of the quantity
$\zeta_{\mathrm{QED}}\equiv\gamma(2\gamma-1)\chi_{\mathrm{QED}}$ because
$\zeta_{\mathrm{QED}}$ can be directly compared with $\xi$ and $\eta$, as is
clear from Eq. (106). Note that we do not calculate the QED radiative
corrections here but rather use values reported in the literature.

Below, we calculate hyperfine splitting for the four models $\nu=0,a,b,c$ for
the magnetic moment distribution. We should point out that the nuclear
magnetic moment distribution can be very different from the nuclear charge
distribution. This is particularly true for nuclei with one nucleon outside a
closed nuclear shell. Appendix D discusses this case. In our calculations we
used the same distribution for both nuclear charge and nuclear magnetic
moment.

###  VI.1 Point-like charge and magnetic moment

Numerical values of $\Delta_{0}$ are tabulated in Table 11 for a number of
isotopes. For some ions, $g_{I}$ is negative, and the hyperfine interaction is
ferromagnetic, hence $\Delta_{0}<0$. Figure 8 plots the magnitude of the
uncorrected hyperfine splitting versus $A$.

Table 11: Nuclear $g$-factors $g_{I}$, nuclear spin $I$, and ground-state hyperfine splitting $\Delta_{0}$ calculated using Eq. (108) [in Hartree], for point-like nuclei. The nuclear $g$ factors for 1H, 2H, 3H and 3He nuclei are taken from 2018 CODATA, i.e., Ref. [10], and those for 39K, 40K, 41K, 85Rb, 87Rb, 133Cs are taken from Ref. [56], those for 135Cs and 235U are taken from Ref. [58], those for 207Pb and 209Bi are taken from Ref. [57]. Isotope ($Z,A$) | $g_{I}$ | $I$ | $\Delta_{0}$
---|---|---|---
H (1, 1) | $5.5856946893(16)$ | $1/2$ | $2.1589180648(6)E^{-7}$
H (1, 2) | $0.8574382338(22)$ | $1$ | $4.972458067(13)E^{-8}$
H (1, 3) | $5.957924931(12)$ | $1/2$ | $2.3036226460(46)E^{-7}$
He (2, 3) | $-4.255250615(50)$ | $1/2$ | $1.316544463(15)E^{-8}$
K (19, 39) | $0.26061413(22)$ | $3/2$ | $1.4234088(12)E^{-4}$
K (19, 40) | $-0.324063(62)$ | $4$ | $-3.98239(77)E^{-4}$
K (19, 41) | $0.14304731(15)$ | $3/2$ | $7.812889(8)E^{-5}$
Rb (37, 85) | $0.5391679(11)$ | $5/2$ | $3.554396(7)E^{-3}$
Rb (37, 87) | $1.8272315(18)$ | $3/2$ | $8.030529(8)E^{-3}$
Cs (55, 133) | $0.732356746(95)$ | $7/2$ | $2.4735911(32)E^{-2}$
Cs (55, 135) | $0.78069(6)$ | $7/2$ | $2.63683(19)E^{-2}$
Pb (82, 207) | $1.16438(4)$ | $1/2$ | $5.1434(2)E^{-2}$
Bi (83, 209) | $0.91347(4)$ | $9/2$ | $2.14590(10)E^{-1}$
U (92, 235) | $-0.109(9)$ | $7/2$ | $-3.66(29)E^{-2}$
Figure 8: Absolute value of the hyperfine splitting $|\Delta_{0}|$ (in
Hartree) in Eq. (108) for a point-like nucleus versus $A$. Positive values of
the hyperfine splitting $\Delta_{0}$ (positive $g_{I}$) are marked by blue
symbols, and negative values of $\Delta_{0}$ (negative $g_{I}$) are shown in
red.

###  VI.2 Charge and magnetic moment distributions for model (a)

Substituting Eqs. (96) and (103) into Eq. (105) with $\nu=a$, we find

$\displaystyle\Delta_{a}$ $\displaystyle=$
$\displaystyle-\frac{4mc^{2}}{3}\,\alpha\big{(}Z\alpha\big{)}^{3}\,g_{I}\,\frac{m_{e}}{m_{p}}\Big{(}I+\frac{1}{2}\Big{)}$
(109) $\displaystyle\times
a_{B}^{2}\int\limits_{r_{a}}^{\infty}g_{+}(r)f_{+}(r)dr.$

As in Eq. (107), it is convenient to write $\Delta_{a}$ in the form

$\Delta_{a}=\Delta_{0}\big{(}1-\xi_{a}\big{)}\big{(}1-\eta_{a}\big{)},$ (110)

where $\Delta_{0}$ is given in Eq. (108). The corrections $\xi_{a}$ and
$\eta_{a}$ are calculate within model (a). Explicitly,

$\xi_{a}=\frac{\mathcal{G}_{0}-\mathcal{G}_{a}}{\mathcal{G}_{0}},\quad\eta_{a}=\frac{\mathcal{G}_{-}^{(a)}}{\mathcal{G}_{a}},$
(111)

where

$\displaystyle\mathcal{G}_{0}$ $\displaystyle=$
$\displaystyle\frac{a_{B}^{2}}{Z^{3}\alpha}\int\limits_{0}^{\infty}g_{0}(r)f_{0}(r)dr$
(112a) $\displaystyle=$ $\displaystyle-\frac{1}{\gamma(2\gamma-1)},$
$\displaystyle\mathcal{G}_{-}^{(a)}$ $\displaystyle=$
$\displaystyle\frac{a_{B}^{2}}{Z^{3}\alpha}\int\limits_{0}^{r_{a}}g_{a,-}(r)f_{a,-}(r)dr,$
(112b) $\displaystyle\mathcal{G}_{a}$ $\displaystyle=$
$\displaystyle\mathcal{G}_{-}^{(a)}+\frac{a_{B}^{2}}{Z^{3}\alpha}\int\limits_{r_{a}}^{\infty}g_{+}(r)f_{+}(r)dr.$
(112c)

Table 12: The corrections $\xi_{a}$ and $\xi_{b}$ to the hyperfine splitting for models (a) [Eq. (111)], and (b) [Eq. (115)]. Isotope ($Z,A$) | $\xi_{a}$ | $\xi_{b}$
---|---|---
H (1, 1) | $1.1092(2)E^{-3}$ | $1.1106(2)E^{-3}$
H (1, 2) | $5.9912(2)E^{-4}$ | $6.0261(2)E^{-4}$
H (1, 3) | $4.087(9)E^{-4}$ | $4.116(10)E^{-4}$
He (2, 3) | $5.126(2)E^{-4}$ | $5.0786(10)E^{-4}$
K (19, 39) | $2.8513(15)E^{-3}$ | $2.7643(7)E^{-3}$
K (19, 40) | $2.853(2)E^{-3}$ | $2.7661(10)E^{-3}$
K (19, 41) | $2.864(4)E^{-3}$ | $2.776(2)E^{-3}$
Rb (37, 85) | $9.173(5)E^{-3}$ | $8.911(2)E^{-3}$
Rb (37, 87) | $9.163(4)E^{-3}$ | $8.902(2)E^{-3}$
Cs (55, 133) | $2.4526(19)E^{-2}$ | $2.3913(9)E^{-2}$
Cs (55, 135) | $2.454(2)E^{-2}$ | $2.3924(9)E^{-2}$
Pb (82, 207) | $0.107053(17)$ | $0.105226(7)$
Bi (83, 209) | $0.11338(3)$ | $0.111483(13)$
U (92, 235) | $0.19317(6)$ | $0.19062(3)$
Table 13: The corrections $\eta_{a}$ and $\eta_{b}$ to the hyperfine splitting for models (a) [Eq. (111)], and (b) [Eq. (115)]. Isotope ($Z,A$) | $\eta_{a}$ | $\eta_{b}$
---|---|---
H (1, 1) | $1.060(11)E^{-5}$ | $1.847(19)E^{-5}$
H (1, 2) | $2.6812(9)E^{-5}$ | $4.6729(16)E^{-5}$
H (1, 3) | $2.22(5)E^{-5}$ | $3.86(8)E^{-5}$
He (2, 3) | $4.961(8)E^{-5}$ | $8.646(13)E^{-5}$
K (19, 39) | $9.083(5)E^{-4}$ | $1.5807(9)E^{-3}$
K (19, 40) | $9.092(7)E^{-4}$ | $1.5822(13)E^{-3}$
K (19, 41) | $9.127(14)E^{-4}$ | $1.588(2)E^{-3}$
Rb (37, 85) | $2.6911(14)E^{-3}$ | $4.665(3)E^{-3}$
Rb (37, 87) | $2.6883(12)E^{-3}$ | $4.660(2)E^{-3}$
Cs (55, 133) | $6.179(5)E^{-3}$ | $1.0645(9)E^{-2}$
Cs (55, 135) | $6.182(5)E^{-3}$ | $1.0650(9)E^{-2}$
Pb (82, 207) | $1.8471(3)E^{-2}$ | $3.1386(6)E^{-2}$
Bi (83, 209) | $1.9209(6)E^{-2}$ | $3.2619(11)E^{-2}$
U (92, 235) | $2.7368(11)E^{-2}$ | $4.620(2)E^{-2}$

The corrections $\xi_{a}$ and $\xi_{b}$ are tabulated in Table 12, and
$\eta_{a}$ and $\eta_{b}$ are tabulated in Table 13. Here the numbers in the
parenthesis are the uncertainties. For example, $\xi_{a}$ for the 1H isotope
is

$\xi_{a}=1.1092(2)E^{-3}\equiv(1.1092\pm 0.0002)\times 10^{-3}.$

For all the isotopes, $|\xi_{a}-\xi_{b}|$ is larger than the uncertainties in
$\xi_{a}$ and $\xi_{b}$, and $|\eta_{a}-\eta_{b}|$ is much larger than the
uncertainties in $\eta_{a}$ and $\eta_{b}$. Thus the ground-state hyperfine
splitting is very sensitive to the nuclear charge and magnetic moment
distribution.

###  VI.3 Charge and magnetic moment distribution inside the nucleus for model
(b)

Substituting Eqs. (97) and (103) into Eq. (105) with $\nu=b$, we get

$\displaystyle\Delta_{b}$ $\displaystyle=$
$\displaystyle-\frac{4mc^{2}}{3}\,\alpha\big{(}Z\alpha\big{)}^{3}\,g_{I}\,\frac{m_{e}}{m_{p}}\Big{(}I+\frac{1}{2}\Big{)}$
(113)
$\displaystyle\times\bigg{\\{}a_{B}^{2}\int\limits_{0}^{r_{b}}\frac{r^{3}}{r_{b}^{3}}\,g_{b,-}(r)f_{b,-}(r)dr$
$\displaystyle+a_{B}^{2}\int\limits_{r_{b}}^{\infty}g_{+}(r)f_{+}(r)dr\bigg{\\}}.$

We write $\Delta_{b}$ in the form

$\Delta_{b}=\Delta_{0}\big{(}1-\xi_{b}\big{)}\big{(}1-\eta_{b}\big{)},$ (114)

where $\Delta_{0}$ is defined in Eq. (108),

$\xi_{b}=\frac{\mathcal{G}_{0}-\mathcal{G}_{b}}{\mathcal{G}_{0}},\quad\eta_{b}=\frac{\mathcal{G}_{-}^{(b)}-\mathcal{F}_{-}^{(b)}}{\mathcal{G}_{b}},$
(115)

$\mathcal{G}_{0}$ is given in Eq. (112a), and

$\displaystyle\mathcal{F}_{-}^{(b)}$ $\displaystyle=$
$\displaystyle\frac{a_{B}^{2}}{Z^{3}\alpha}\int\limits_{0}^{r_{b}}\frac{r^{3}}{r_{b}^{3}}\,g_{b,-}(r)f_{b,-}(r)dr,$
(116a) $\displaystyle\mathcal{G}_{-}^{(b)}$ $\displaystyle=$
$\displaystyle\frac{a_{B}^{2}}{Z^{3}\alpha}\int\limits_{0}^{r_{b}}g_{b,-}(r)f_{b,-}(r)dr,$
(116b) $\displaystyle\mathcal{G}_{b}$ $\displaystyle=$
$\displaystyle\mathcal{G}_{-}^{(b)}+\frac{a_{B}^{2}}{Z^{3}\alpha}\int\limits_{r_{b}}^{\infty}g_{+}(r)f_{+}(r)dr.$
(116c)

Figure 9: Corrections to the ground-state hyperfine splitting due to FNS
distributions of the nuclear charge and magnetic moment as a function of the
atomic mass number $A$. The corrections $\xi_{a}$ and $\eta_{a}$ [red and
orange] are defined in Eq. (111) and the correction $\xi_{b}$ and $\eta_{b}$
[blue and green] is defined in Eq. (115) [blue]. The solid curves that fit the
data for $\xi_{a}$, $\eta_{a}$ and $\xi_{b}$, $\eta_{b}$ are meant to guide
the eye (see the text).
Redo figure.

###  VI.4 FNS corrections to hyperfine splitting in models (a) and (b)

The FNS corrections to the hyperfine splitting $\xi_{a}$, $\xi_{b}$,
$\eta_{a}$ and $\eta_{b}$ are tabulated in Tables 12 and 13 and are plotted in
Fig. 9 as functions of $A$. Numerical analysis shows that $\ln\xi_{a}$,
$\ln\eta_{a}$, $\ln\xi_{b}$ and $\ln\eta_{b}$ can be well-fit by

$\displaystyle\ln\xi_{a}(A)$ $\displaystyle=$
$\displaystyle-6.79743-1.20214\,\ln A$ (117a)
$\displaystyle+0.384078\,\ln^{2}A,$ $\displaystyle\ln\eta_{a}(A)$
$\displaystyle=$ $\displaystyle-11.2423+0.661711\,\ln A$ (117b)
$\displaystyle+0.128055\,\ln^{2}A,$ $\displaystyle\ln\xi_{b}(A)$
$\displaystyle=$ $\displaystyle-6.78386-1.22495\,\ln A$ (117c)
$\displaystyle+0.38719\,\ln^{2}A,$ $\displaystyle\ln\eta_{b}(A)$
$\displaystyle=$ $\displaystyle-10.6924+0.673483\,\ln A$ (117d)
$\displaystyle+0.125229\,\ln^{2}A.$

FNS corrections to the hyperfine splitting of the ground state of heavy HLI
were calculated in Ref. [22]. The corrections $\xi_{a}$ for 207Pb and 209Bi in
Table 12 are not far from the values calculated in Ref. [22],
$\xi(\mathrm{Pb})=0.1049$, and $\xi(\mathrm{Bi})=0.1111$, the values of
$\eta_{a}$ for 209Bi in Table 13 not far from the value reported in [22],
$\epsilon=0.0118$, and $\eta_{b}$ for 207Pb not too far away from
$\epsilon=0.0419$ in Ref. [22].

Within the assumption that $W_{ba}(r)$ in Eq. (9) is small [i.e., $W_{ba}(r)$
is much weaker than $W_{a}(r)$ or $W_{b}(r)$ in Eq. (8), see Fig. 1], we
showed that the correction $\xi$ is almost the same for both models (a) and
(b), whereas models (a) and (b) predict very different values for the
corrections $\eta_{\nu}$ due to the magnetic moment distribution.

###  VI.5 Charge and magnetic moment distribution inside the nucleus for model
(c)

For the TPFM, model (c),

$\Delta_{c}=\Delta_{0}\big{(}1-\xi_{c}\big{)}\big{(}1-\eta_{c}\big{)},$ (118)

where $\Delta_{0}$ is given in (108),

$\xi_{c}=\frac{\mathcal{G}_{0}-\mathcal{G}_{c}}{\mathcal{G}_{0}},\quad\eta_{c}=\frac{\mathcal{F}_{c}}{\mathcal{G}_{c}},$
(119)

and

$\displaystyle\mathcal{G}_{c}$ $\displaystyle=$
$\displaystyle\frac{a_{B}^{2}}{Z^{3}\alpha}\int\limits_{0}^{\infty}g_{c}(r)f_{c}(r)dr,$
(120a) $\displaystyle\mathcal{F}_{c}$ $\displaystyle=$
$\displaystyle\frac{a_{B}^{2}}{Z^{3}\alpha}\int\limits_{0}^{\infty}\Big{[}1-\mathcal{M}_{c}(r)\Big{]}g_{c}(r)f_{c}(r)dr.$
(120b)

###  VI.6 Relativistic QED radiative corrections for hyperfine splitting

We now include the relativistic QED radiative corrections
$\zeta_{\mathrm{QED}}=\gamma(2\gamma-1)\chi_{\mathrm{QED}}$ given in Eq.
(107). The quantity $\chi_{\mathrm{QED}}$ is calculated in Refs. [22, 25] for
some elements. The corrections $\zeta_{\mathrm{QED}}$ are tabulated in Table
14.

Table 14: Radiative QED correction $\zeta_{\mathrm{QED}}=\gamma(2\gamma-1)\chi_{\mathrm{QED}}$ to the hyperfine splitting in Eq. (107), where $\gamma$ is given in Eq. (17) and $\chi_{\mathrm{QED}}$ is calculated in Refs. [22, 25] for some elements. The fit $\zeta_{\mathrm{QED}}^{\mathrm{(fit)}}$ is calculated using Eq. (121). $Z$ | $\zeta_{\mathrm{QED}}$ | $\zeta_{\mathrm{QED}}^{\mathrm{(fit)}}$
---|---|---
$1$ | $0.00105756$ | $0.00106387$
$2$ |  | $0.000946973$
$3$ | $0.000835844$ | $0.000831409$
$5$ | $0.000609149$ | $0.000603679$
$7$ | $0.000382069$ | $0.00038002$
$10$ | $0.0000466108$ | $0.0000516092$
$19$ |  | $-0.000884945$
$37$ |  | $-0.00253464$
$49$ | $-0.00348472$ | $-0.00345706$
$53$ | $-0.00373754$ | $-0.00373092$
$55$ |  | $-0.00386134$
$57$ | $-0.0039463$ | $-0.00398737$
$63$ | $-0.00427325$ | $-0.00433865$
$67$ | $-0.00467521$ | $-0.00455007$
$71$ | $-0.00461933$ | $-0.00474287$
$75$ | $-0.00501339$ | $-0.00491674$
$82$ | $-0.0051733$ | $-0.0051745$
$83$ | $-0.0051822$ | $-0.00520641$
92 |  | $-0.00543714$

We fit $\zeta_{\mathrm{QED}}$ by the formula [25],

$\displaystyle\zeta_{\mathrm{QED}}^{\mathrm{(fit)}}(Z)=\frac{\alpha}{\pi}\Big{[}A_{0}-A_{1}\,\alpha
Z+A_{2}\,\big{(}\alpha Z)^{2}$ $\displaystyle~{}~{}~{}~{}~{}+A_{3}\,(\alpha
Z)^{2}\ln(\alpha Z)+A_{4}\big{(}\alpha Z\ln(\alpha Z)\big{)}^{2}\Big{]}.$
(121)

The fit gives,

$\displaystyle A_{0}$ $\displaystyle=$ $\displaystyle 0.509092$ $\displaystyle
A_{1}$ $\displaystyle=$ $\displaystyle 7.08715,$ $\displaystyle A_{2}$
$\displaystyle=$ $\displaystyle 4.43975,$ $\displaystyle A_{3}$
$\displaystyle=$ $\displaystyle 0.696079,$ $\displaystyle A_{4}$
$\displaystyle=$ $\displaystyle 0.448402,$

and the mean standard deviation $\sigma_{\mathrm{QED}}$ is

$\displaystyle\sigma_{\mathrm{QED}}$ $\displaystyle=$
$\displaystyle(\frac{1}{N_{Z}}\sum_{n_{Z}}^{N_{Z}}\big{(}\zeta_{\mathrm{QED}}(Z)-\zeta_{\mathrm{QED}}^{\mathrm{(fit)}}(Z)\big{)}^{2})^{1/2}$
$\displaystyle=$ $\displaystyle 0.000058397,$

where $\zeta_{\mathrm{QED}}$ was taken from Table 14, $N_{Z}=14$ is the number
of elements for which $\zeta_{\mathrm{QED}}$ is tabulated in the table, and
$n_{Z}$ is a running index over the elements used. The QED radiative
correction $\zeta_{\mathrm{QED}}$ and its fit using Eq. (121) are shown in
Fig. 10. The fit (blue curve) is close to the red dots, hence we conclude that
the approximation in Eq. (121) is satisfactory. The QED radiative corrections
$\zeta_{\mathrm{QED}}^{\mathrm{(fit)}}(Z)$ in Eq. (121) are tabulated in Table
14 for H&HLI. For heavy isotopes, $\zeta_{\mathrm{QED}}(Z)$ satisfies the
inequalities: $\zeta_{\mathrm{QED}}<0$, and
$\eta_{a}<|\zeta_{\mathrm{QED}}|<\eta_{b}\ll\xi$, where $\xi_{a}$ and
$\xi_{b}$ are tabulated in Table 12, and $\eta_{a}$ and $\eta_{b}$ are
tabulated in Table 13. Therefore, we conclude that the main contribution to
the hyperfine splitting is the correction $\xi$ due to the finite nuclear
charge distribution, which is much larger than the relativistic QED radiative
corrections and the correction of the magnetic momentum distribution due to
FNS. For light elements (e.g., H and He), $\zeta_{\mathrm{QED}}>0$ and
satisfies the inequalities,
$\zeta_{\mathrm{QED}}>\xi_{a},\xi_{b}>\eta_{a},\eta_{b}$, i.e., the QED
radiative corrections to the ground-state hyperfine splitting are larger than
the FNS effects. For heavy elements, $\zeta_{\mathrm{QED}}<0$ and satisfies
the inequalities, $|\zeta_{\mathrm{QED}}|<\eta_{a},\eta_{b}<\xi_{a},\xi_{b}$;
The FNS corrections to the hyperfine splitting due to the charge distribution
are the largest correction using either model (a) or (b).

Figure 10: Radiative QED correction (red dots)
$\zeta_{\mathrm{QED}}=\gamma(2\gamma-1)\chi_{\mathrm{QED}}$ to the hyperfine
splitting in Eq. (107), versus nuclear charge $Z$ where $\chi_{\mathrm{QED}}$
is calculated in Refs. [22, 25]. The blue curve is the fit of
$\zeta_{\mathrm{QED}}$ in Eq. (121).

##  VII Estimates of the uncertainty of the hyperfine splitting in models (a),
(b) and (c)

The ground-state hyperfine splittings $\Delta_{a}$, Eq. (110), and
$\Delta_{b}$, Eq. (114), depend on the nuclear g factor $g_{I}$ and the RMS
nuclear charge radius $r_{N}$. Both $g_{I}$ and $r_{N}$ are measured in
experiment with their uncertainties $\delta g_{I}$ and $\delta r_{N}$, see
Tables 1 and 11. The uncertainties $\delta g_{I}$ and $\delta r_{N}$ give rise
to uncertainties $\delta\Delta_{a}$ and $\delta\Delta_{b}$ that are given by

$\delta\Delta_{\nu}=\sqrt{\bigg{(}\frac{\partial\Delta_{\nu}}{\partial
g_{I}}\,\delta g_{I}\bigg{)}^{2}+\bigg{(}\frac{\partial\Delta_{\nu}}{\partial
r_{N}}\,\delta r_{N}\bigg{)}^{2}}\,.$ (122)

Using Eq. (110) and (114), we can write Eq. (122) as,

$\delta\Delta_{\nu}=\Delta_{\nu}\,\sqrt{\bigg{(}\frac{\delta
g_{I}}{g_{I}}\bigg{)}^{2}+\bigg{(}\frac{\delta\xi_{\nu}}{1-\xi_{\nu}}+\frac{\delta\eta_{\nu}}{1-\eta_{\nu}}\bigg{)}^{2}}\,.$
(123)

The ground-state hyperfine splitting $\Delta_{c}$ in Eq. (118) depends on
$g_{I}$, $r_{N}$ and the distribution thickness parameter $a_{F}$. The
parameters $g_{I}$, $r_{N}$ and $a_{F}$, and their uncertainties $\delta
g_{I}$, $\delta r_{N}$ and $\delta a_{F}$, are measured in experiment, see
Tables 1, 9 and 11. The uncertainty $\delta\Delta_{c}$ in $\Delta_{c}$ duet to
the uncertainties $\delta g_{I}$, $\delta r_{N}$ and $\delta a_{F}$ is
estimated as,

$\delta\Delta_{c}=\sqrt{\bigg{(}\frac{\partial\Delta_{c}}{\partial
g_{I}}\,\delta g_{I}\bigg{)}^{2}+\bigg{(}\frac{\partial\Delta_{c}}{\partial
r_{N}}\,\delta r_{N}\bigg{)}^{2}+\bigg{(}\frac{\partial\Delta_{c}}{\partial
a_{F}}\,\delta a_{F}\bigg{)}^{2}}\,.$ (124)

Using Eq. (118), we can rewrite Eq. (124) in the form,

$\displaystyle\delta\Delta_{c}$ $\displaystyle=$
$\displaystyle\Bigg{[}\bigg{(}\frac{\delta
g_{I}}{g_{I}}\bigg{)}^{2}+\bigg{(}\frac{\delta\xi_{c}^{(r)}}{1-\xi_{c}}+\frac{\delta\eta_{c}^{(r)}}{1-\eta_{c}}\bigg{)}^{2}$
(125)
$\displaystyle+\bigg{(}\frac{\delta\xi_{c}^{(a)}}{1-\xi_{c}}+\frac{\delta\eta_{c}^{(a)}}{1-\eta_{c}}\bigg{)}^{2}\Bigg{]}^{1/2},$

where $\delta\xi_{c}^{(r)}$ and $\delta\xi_{c}^{(a)}$ are the uncertainties in
$\xi_{c}$ due to the uncertainties $\delta r_{N}$ and $\delta a_{F}$ in
$r_{N}$ and $a_{F}$,

$\delta\xi_{c}^{(r)}=\bigg{|}\frac{\partial\xi_{c}}{\partial
r_{N}}\bigg{|}\,\delta
r_{N},\quad\delta\xi_{c}^{(a)}=\bigg{|}\frac{\partial\xi_{c}}{\partial
a_{F}}\bigg{|}\,\delta a_{F},$

and $\delta\eta_{c}^{(r)}$ and $\delta\eta_{c}^{(a)}$ are the uncertainties in
$\xi_{c}$ due to the uncertainties $\delta r_{N}$ and $\delta a_{F}$ in
$r_{N}$ and $a_{F}$,

$\delta\eta_{c}^{(r)}=\bigg{|}\frac{\partial\eta_{c}}{\partial
r_{N}}\bigg{|}\,\delta
r_{N},\quad\delta\eta_{c}^{(a)}=\bigg{|}\frac{\partial\eta_{c}}{\partial
a_{F}}\bigg{|}\,\delta a_{F}.$

##  VIII Hyperfine splitting: Comparison with experiment

We now compare the results of experiments on the ground state hyperfine
splitting of H&HLI with the results of our calculations using model (a) and
(b). We present results of the calculations without and with relativistic QED
radiative corrections.

The percent uncertainties of the theoretical results compared with the
experimental values are defined as,

$\displaystyle\tau_{\nu}=\frac{\Delta_{\mathrm{exp}}-\tilde{\Delta}_{\nu}}{\Delta_{\mathrm{exp}}},$
(126)

where $\nu=a,b,c$, and $\tilde{\Delta}_{\nu}$ is a hyperfine splitting with
the QED radiative correction given by

$\tilde{\Delta}_{\nu}=\Delta_{0}\,\Big{[}\big{(}1-\xi_{\nu}\big{)}\big{(}1-\eta_{\nu}\big{)}+\zeta_{\mathrm{QED}}\Big{]}=\Delta_{\nu}+\Delta_{0}\,\zeta_{\mathrm{QED}}.$
(127)

The experimental hyperfine splitting $\Delta_{\mathrm{exp}}$ for hydrogen
atoms and helium ions are taken from Refs. [48, 49, 50, 51, 52], values for
hydrogenic lead ions are taken from Ref. [53], and values for hydrogenic
bismuth ions are taken from Ref. [54]. The quantities $\tilde{\Delta}_{a}$,
$\tilde{\Delta}_{b}$, $\Delta_{\mathrm{exp}}$, $\tau_{a}$ and $\tau_{b}$ are
tabulated in Table 15. For 1H atoms and 3He ions, $\tau_{a}<0$ and
$\tau_{b}>0$, and
$\tilde{\Delta}_{a}<\Delta_{\mathrm{exp}}<\tilde{\Delta}_{b}$. Moreover,
$|\tau_{a}|$ and $\tau_{b}$ are smaller than the FNS corrections $\eta_{\nu}$
and $\xi_{\nu}$ and the QED radiative correction $\zeta_{\mathrm{QED}}$. For
2H and 3H atoms, both $\tau_{a}$ and $\tau_{b}$ are positive, and
$\tilde{\Delta}_{a}$ and $\tilde{\Delta}_{b}$ are below
$\Delta_{\mathrm{exp}}$. Moreover, $\tau_{a}$ and $\tau_{b}$ exceed
$\eta_{\nu}$, but are below $\xi_{\nu}$ and $\zeta_{\mathrm{QED}}$. For 207Pb
ions, $\tau_{a}<0$ and $\tau_{b}>0$, and for 209Bi ions, $\tau_{a}>0$ and
$\tau_{b}>0$. For both 207Pb and 209Bi, $|\tau_{a}|$ and $\tau_{b}$ are below
$\eta_{\nu}$ and $\xi_{\nu}$, but above $\zeta_{\mathrm{QED}}$.

Table 15: Experimental values of the hyperfine splitting $\Delta_{\mathrm{exp}}$ [in Hartree], as compared with the calculated results $\tilde{\Delta}_{a}$ and $\tilde{\Delta}_{b}$ in Eq. (127) using models (a) and (b) with relativistic QED radiative corrections, and the percent uncertainties ($\tilde{\tau}_{a},\tilde{\tau}_{b}$) of theoretical results compared with $\Delta_{\mathrm{exp}}$. The numbers in brackets show uncertainty $\delta\Delta_{\mathrm{exp}}$ in $\Delta_{\mathrm{exp}}$, and $\delta\Delta_{\nu}$ in $\Delta_{\nu}$ with $\nu=a,b$, see Eq. (123). The experimental values for hydrogen atoms and helium ions are taken from Refs. [48, 49, 50, 51, 52], values for hydrogenic lead ions are taken from Ref. [53] and values for hydrogenic bismuth ions are taken from Ref. [54]. The percent uncertainties $\tau_{\nu}$ of the theoretical results compared with the experimental values [see Eq. (126)]. The model (c) results for heavy ions are not presented in this table; they are presented below. Isotope ($Z,A$) | $\tilde{\Delta}_{a}$ | $\delta\Delta_{a}/\Delta_{a}$ | $\tilde{\Delta}_{b}$ | $\delta\Delta_{b}/\Delta_{b}$ | $\Delta_{\mathrm{exp}}$ | $\tau_{a}$ | $\tau_{b}$
---|---|---|---|---|---|---|---
H (1, 1) | $2.1587837(7)E^{-7}$ | $3E^{-7}$ | $2.1587637(8)E^{-7}$ | $4E^{-7}$ | $2.158775054569(2)E^{-7}$ | $-4.0E^{-6}$ | $5.2E^{-6}$
H (1, 2) | $4.97460440(14)E^{-8}$ | $3E^{-8}$ | $4.97448807(18)E^{-8}$ | $4E^{-8}$ | $4.9756850998(3)E^{-8}$ | $2.2E^{-4}$ | $2.4E^{-4}$
H (1, 3) | $2.305066(3)E^{-7}$ | $1.4E^{-6}$ | $2.305022(4)E^{-7}$ | $1.8E^{-6}$ | $2.30512816232(1)E^{-7}$ | $2.7E^{-5}$ | $4.6E^{-5}$
He (2, 3) | $-1.3170511(4)E^{-6}$ | $3E^{-7}$ | $-1.3170088(3)E^{-6}$ | $2E^{-7}$ | $-1.317031330040(2)E^{-6}$ | $-1.5E^{-5}$ | $1.7E^{-5}$
K(19, 39) | $1.416801(3)E^{-4}$ | $2E^{-6}$ | $1.415971(3)E^{-4}$ | $1.8E^{-6}$ |  |  |
K(19, 40) | $-3.9639(8)E^{-4}$ | $1.9E^{-4}$ | $-3.9616(8)E^{-4}$ | $1.9E^{-4}$ |  |  |
K(19, 41) | $7.77649(4)E^{-5}$ | $5E^{-6}$ | $7.77191(3)E^{-5}$ | $4E^{-6}$ |  |  |
Rb(37, 85) | $3.50330(2)E^{-3}$ | $7E^{-6}$ | $3.497280(19)E^{-3}$ | $5E^{-6}$ |  |  |
Rb(37, 87) | $7.91520(4)E^{-3}$ | $4E^{-6}$ | $7.90160(3)E^{-3}$ | $4E^{-6}$ |  |  |
Cs(55, 133) | $2.38846(6)E^{-2}$ | $2E^{-5}$ | $2.37919(4)E^{-2}$ | $1.8E^{-5}$ |  |  |
Cs(55, 135) | $2.54604(19)E^{-2}$ | $8E^{-5}$ | $2.53615(19)E^{-2}$ | $8E^{-5}$ |  |  |
Pb (82, 207) | $4.5117(19)E^{-2}$ | $4E^{-5}$ | $4.46840(17)E^{-2}$ | $4E^{-5}$ | $4.4684(7)E^{-2}$ | $-9.7E^{-3}$ | $8.3E^{-3}$
Bi (83, 209) | $1.85493(12)E^{-1}$ | $1.2E^{-5}$ | $1.83335E^{-1}$ | $1.0E^{-5}$ | $1.86871(3)E^{-1}$ | $1.4E^{-3}$ | $3.5E^{-3}$
U(92, 235) | $-2.85(22)E^{-2}$ | $8E^{-2}$ | $-3.46(27)$ | $8E^{-2}$ |  |  |

We now compare our results for model (c) with experiment. Calculations for the
207Pb ion yield:

$\xi_{c}=0.11105(8),\quad\eta_{c}=0.02086(3).$

Taking into account that $\Delta_{0}=5.1434(2)\times 10^{-2}$, we find

$\Delta_{c}=4.4858(2)\times 10^{-2},$

hence

$\tau_{c}=\frac{\Delta_{\mathrm{exp}}-\Delta_{c}}{\Delta_{\mathrm{exp}}}=-0.0039.$

Comparing $\tau_{c}$ with $\tau_{a}=-0.0097$ and $\tau_{b}=0.0083$ in Table
15, one can see that model (c) results are better for the ground-state
hyperfine splitting than for models (a) and (b). However,
$|\Delta_{\mathrm{exp}}-\Delta_{c}|>\delta\Delta_{c}$.

Calculations for 209Bi ion give:

$\xi_{c}=0.104649(27),\quad\eta_{c}=0.020142(10).$

Taking into account that $\Delta_{0}=2.14590(10)E^{-2}$, we get

$\Delta_{c}=0.185665(15)E^{-2},$

such that

$\tau_{c}=\frac{\Delta_{\mathrm{exp}}-\Delta_{c}}{\Delta_{\mathrm{exp}}}=0.0064526.$

Comparing $\tau_{c}$ with $\tau_{a}=0.0014$ and $\tau_{b}=0.0035$ in Table 15,
one sees that model (a) gives somewhat better results for the ground-state
hyperfine splitting than models (b) and (c) for 209Bi.

##  IX Summary and Conclusions

We calculated FNS corrections for the ground state ($1s$) and the first
excited state ($2s$) of H&HLI using the Dirac equation with three different
models for the atomic nucleus charge distribution. For model (a) [nuclear
charge on the surface of the nucleus], we present the analytic solution of the
Dirac equation, whereas for model (b) [nuclear charge uniformly distributed
inside the nucleus] and model (c) [the two-parameter Fermi model – TPFM], we
use both numerical calculations and first-order perturbation theory (the
latter is not valid for large $Z$ nuclei). The FNS corrections to the ground
state energy and the first excited state energy are smaller than the electron-
nucleus reduced mass correction and the relativistic QED radiative corrections
for light nuclei (i.e., low $Z$ nuclei) [33]. Nevertheless, it is important to
obtain accurate FNS corrections so that a good comparison can be made with the
extremely high-accuracy experimental results for hydrogen [18, 19] and
deuterium [69, 70]. A detailed study of the hydrogen $1s$-$2s$ transition has
been reported in Ref. [28]. For heavy nuclei, the electron-nucleus reduced
mass correction is small because $m_{e}/M(Z,A)$ is small, hence the FNS
correction is larger than the electron-nucleus reduced mass correction.
Furthermore, for heavy nuclei, $Z\alpha$ is not a small parameter, hence
higher order relativistic QED radiative corrections, which can be neglected
for hydrogen and light nuclei, are required, and therefore it is not clear
whether FNS corrections are more important than relativistic QED radiative
corrections since the latter have not been calculated.

We find that for H atoms and light HLI, the difference of the FNS corrections
obtained using model (a) and model (b), $\varepsilon_{n,a}-\varepsilon_{n,b}$,
are small; the difference between model (a) and model (b) energies
$\varepsilon_{n}$ for $n=1$ (the $1s$ state) and for $n=2$ (the $2s$ state) is
smaller than the uncertainties $\delta\varepsilon_{n,a}$ and
$\delta\varepsilon_{n,b}$ in $\varepsilon_{n,a}$ and $\varepsilon_{n,b}$, see
Tables 3 and 4. Hence, measurement of the $1s$-$2s$ transition frequency
cannot be used to determine details of the nuclear charge distribution for low
$Z$ nuclei [the perturbation theory expression in Eq. (71) is a good
approximation for low $Z$ nuclei, and perturbation theory expression does not
depend on the details of the charge distribution, but only on $r_{N}$].
However, for heavy HLIs, the difference,
$\varepsilon_{n,a}-\varepsilon_{n,b}$, is larger than the uncertainties
$\delta\varepsilon_{n,a}$ and $\delta\varepsilon_{n,b}$ in $\varepsilon_{n,a}$
and $\varepsilon_{n,b}$. Thus measuring the $1s$-$2s$ transition frequency for
heavy nuclei HLIs can be used to determine not only the nuclear charge radii
but the details of the nuclear charge distribution. Note that QED radiative
corrections are available only for $Z\alpha\ll 1$. Hence, there is a problem
applying the QED corrections for heavy nuclei.

Experimental results for the $1s$-$3s$ transition frequency for hydrogen and
deuterium have also been reported, see e.g., Ref. [65]. We intend to calculate
this transition frequency to compare with these experiments in future work.

We have also calculated the FNS corrections for the ground-state hyperfine
splitting of H&HLI using models (a), (b) and (c). After calculating the
hyperfine splitting for point-like nuclei (model $\nu=0$), we calculate the
corrections $\xi_{\nu}$ (where $\nu=a,b,c$) due to the nuclear charge
distribution, and the corrections $\eta_{\nu}$ due to the magnetic moment
distribution. The percent uncertainties $\tau_{\nu}$ of the theoretical
results are compared with the experimental values [see Eq. (126)]. Both
$|\tau_{a}|$ and $|\tau_{b}|$ are smaller than the correction $\xi_{\nu}$ to
the hyperfine splitting. For 1H atoms and 3He ions, $|\tau_{\nu}|<\eta_{\nu}$.
For 2H, 3H atoms, and 207Pb and 209Bi HLIs, $|\tau_{\nu}|>\eta_{\nu}$. We
calculate also the ground-state hyperfine splitting for 207Pb and 209Bi HLIs
with model (c) for the nuclear charge and magnetic moment distribution. For
the 207Pb HLI, model (c) gives a somewhat better result for the ground-state
hyperfine splitting than models (a) and (b). For the 209Bi HLI, model (a)
gives a somewhat better result for the ground-state hyperfine splitting than
models (b) and (c).

We show that $\xi_{a}>\eta_{a}$ for all the isotopes (see tables 12 and 13),
i.e., the nuclear charge distribution correction to the ground-state hyperfine
splitting is larger than the magnetic moment distribution correction. On the
other hand, the difference $|\xi_{a}-\xi_{b}|$ is much larger than the
uncertainties $\delta\xi_{a}$ and $\delta\xi_{b}$ in $\xi_{a}$ and $\xi_{b}$,
and the difference $|\eta_{a}-\eta_{b}|$ is much larger than the uncertainties
$\delta\eta_{a}$ and $\delta\eta_{b}$ in $\eta_{a}$ and $\eta_{b}$. Hence, it
is crucial to test more accurate models of the nuclear charge and magnetic
moment distributions in order to obtain results closer to experimental
hyperfine splittings. The work that has been invested into calculating QED
radiative corrections, can only bear fruit if accurate FNS corrections are
developed. Given the fact that the FNS effects for hyperfine splittings are as
big or bigger than the QED radiative corrections for all the nuclei studied,
how can one obtain more accurate nuclear charge and magnetic moment
distributions in order to improve the FNS corrections? Several approaches are
possible. (1) One can use more elaborate models for the charge and magnetic
moment distributions. For example, one can use the three-parameter Fermi
model, which has an additional term, $(1+w_{F}\,r^{2}/c_{F}^{2})$, multiplying
the TPFM [24] or the double three-parameter Fermi model [63]. (2) Direct
numerical fits to nuclear charge and magnetic moment distributions can be
employed using data from electron scattering from nuclei [64]. (3) The charge
and magnetic moment distributions are in general not equal, see e.g., Appendix
D, hence distributions should be determined for each from experiment.

Acknowledgement: We would like to thank Shalom Shlomo for useful
correspondence.

##  Appendix A Nuclear charge radius

The nuclear charge radius $r_{N}=\sqrt{\langle r^{2}\rangle_{N}}$ is defined
as a RMS charge radius which is a measure of the normalized charge
distribution $\rho_{N}(\mathbf{r})$,

$\langle r^{2}\rangle_{N}=\int r^{2}\rho_{N}(r)d^{3}\mathbf{r},$ (128)

where $\rho_{N}(\mathbf{r})$ is normalized by the condition

$\int\rho_{N}(\mathbf{r})d^{3}\mathbf{r}=1.$

We now apply Eq. (128) and calculate the nuclear radius for models (a) and
(b).

Model (a): When the charge is distributed uniformly on a sphere of radius
$r_{a}$, the normalized proton distribution
$\rho_{N}(\mathbf{r})\equiv\rho_{a}(r)$ is,

$\rho_{a}(r)=\frac{Ze}{4\pi r_{a}^{2}}\,\delta(r-r_{a}),$ (129)

Substituting Eq. (129) into Eq. (128), we get

$r_{a}=r_{N}.$ (130)

Model (b): When the charge is distributed uniformly inside a sphere of radius
$r_{b}$, the normalized proton distribution
$\rho_{N}(\mathbf{r})\equiv\rho_{b}(r)$ is,

$\rho_{b}(r)=\frac{3Ze}{4\pi r_{b}^{3}}\,\Theta(r_{b}-r),$ (131)

where $\Theta(r_{b}-r)$ is equal to $1$ for $r<r_{b}$, $\frac{1}{2}$ for
$r=r_{b}$ and $0$ for $r>r_{b}$. Substituting Eq. (131) into Eq. (128), we get
$r_{N}=r_{b}\sqrt{3/5}$, or

$r_{b}=r_{N}\,\sqrt{\frac{5}{3}}.$ (132)

##  Appendix B Matrix elements of $\mathbf{e}_{r}\times\bm{\alpha}$

The vector operator $\mathbf{e}_{r}\times\bm{\alpha}$ can be written as,

$\displaystyle\mathbf{e}_{r}\times\bm{\alpha}$ $\displaystyle=$
$\displaystyle\left(\begin{array}[]{cc}0&\mathbf{e}_{r}\times\bm{\sigma}\\\
\mathbf{e}_{r}\times\bm{\sigma}&0\end{array}\right),$

where the 2$\times$2 dimensional matrices $\mathbf{e}_{r}\times\bm{\sigma}$
are given by

$\displaystyle\mathbf{e}_{r}\times\bm{\sigma}$ $\displaystyle=$
$\displaystyle\left(\begin{array}[]{cc}\sin\theta\sin\phi&i\cos\theta\\\
-i\cos\theta&-\sin\theta\sin\phi\end{array}\right)\mathbf{e}_{x}+$
$\displaystyle\left(\begin{array}[]{cc}-\sin\theta\cos\phi&\cos\theta\\\
\cos\theta&\sin\theta\cos\phi\end{array}\right)\mathbf{e}_{y}+$
$\displaystyle\left(\begin{array}[]{cc}0&-ie^{-i\phi}\\\
ie^{i\phi}&0\end{array}\right)\mathbf{e}_{z}.$

The matrix elements of $\mathbf{e}_{r}\times\bm{\alpha}$ with the wave
functions $\psi_{m_{J}}(\mathbf{r})$ in Eq. (102) take the form,

$\displaystyle\int\psi_{m_{J}}^{{\dagger}}(\mathbf{r})\,\big{[}\mathbf{e}_{r}\times\bm{\alpha}\big{]}\,\psi_{m^{\prime}_{J}}(\mathbf{r})\sin\theta
d\theta d\phi\times$
$\displaystyle~{}~{}~{}=g(r)f(r)\int\limits_{0}^{\pi}\sin\theta
d\theta\int\limits_{0}^{2\pi}d\phi\times$
$\displaystyle~{}~{}~{}~{}~{}~{}\Big{\\{}i\,\Omega_{\frac{1}{2}0m_{J}}^{{\dagger}}\big{[}\mathbf{e}_{r}\times\bm{\sigma}\big{]}\Omega_{\frac{1}{2}1m^{\prime}_{J}}+\mathrm{c.c.}\Big{\\}}$
$\displaystyle~{}~{}~{}=-\frac{8}{3}\,g(r)f(r)\,\mathbf{J}_{m_{J},m^{\prime}_{J}},$
(137)

where $\mathbf{J}_{m_{J},m^{\prime}_{J}}$ are matrix elements of the orbital
angular momentum operator $\mathbf{J}$.

##  Appendix C Iteration method for determining $s$ state energies using Eq.
(60)

We shall now explain the iteration method, using the case of the 3He ion for
specificity. The ground-state energy of the 3He ion for model (a) is
$\varepsilon_{1}^{(0)}=1.47460\times 10^{-8}\epsilon_{H}$ (see Table 3), where
$\varepsilon_{1}^{(0)}$ is the ground-state energy of the He ion with the
point-like nucleus given in Eq. (23). The energy using model (b) is close to
the energy for model (a), see Table 8 below. Thus we take
$\mathcal{E}_{1}=\varepsilon_{1}^{(0)}-4\times 10^{-8}\epsilon_{H}$ as a first
iteration for the energy. Solving Eq. (56) numerically for $\mathcal{E}$, we
find $g_{b,-}(r),f_{b,-}(r)$ for $r<r_{b}$. Substituting
$g_{b,-}(r_{b}),f_{b,-}(r_{b})$, and $g_{+}(r_{b}),f_{+}(r_{b})$ with
$\mathcal{E}_{1}$ [see in Eq. (38)] into Eq. (60), we get
$\mathfrak{D}_{1}\equiv\mathfrak{D}(\mathcal{E}_{1})=-5.16387\times 10^{-5}$.
The second iteration is $\mathcal{E}_{2}=\varepsilon_{1}^{(0)}-1.4\times
10^{-8}\epsilon_{H}$. Solving Eq. (56) numerically for $\mathcal{E}_{2}$, we
find $g_{b,-}(r),f_{b,-}(r)$ for $r<r_{b}$. Substituting
$g_{b,-}(r_{b}),f_{b,-}(r_{b})$, and $g_{+}(r_{b}),f_{+}(r_{b})$ with
$\mathcal{E}_{2}$ [see in Eq. (38)] into Eq. (60), we get
$\mathfrak{D}_{2}\equiv\mathfrak{D}(\mathcal{E}_{2})=1.52535\times 10^{-6}$.
In order to find the third iteration for the ground-state energy
$\varepsilon$, we plot a graph of $\mathfrak{D}(\varepsilon)$ versus
$\varepsilon$, with the two points, $(\mathcal{E}_{1},\mathfrak{D}_{1})$, and
$(\mathcal{E}_{2},\mathfrak{D}_{2})$. We then connect these points by a line,
and find the crossing point $\mathcal{E}_{3}$ of this line with the $x$ axis,

$\mathcal{E}_{3}=\frac{\mathcal{E}_{1}\mathfrak{D}_{2}-\mathcal{E}_{2}\mathfrak{D}_{1}}{\mathfrak{D}_{2}-\mathfrak{D}_{1}}=\varepsilon_{1}^{(0)}-1.474597673\times
10^{-8}\,\epsilon_{H}.$

Solving Eq. (56) numerically for $\mathcal{E}_{3}$, we find
$g_{b,-}(r),f_{b,-}(r)$ for $r<r_{b}$. Substituting
$g_{b,-}(r_{b}),f_{b,-}(r_{b})$, and $g_{+}(r_{b}),f_{+}(r_{b})$ with
$\mathcal{E}_{3}$ [see in Eq. (38)] into Eq. (60), we get
$\mathfrak{D}_{3}\equiv\mathfrak{D}(\mathcal{E}_{3})=4.59235\times 10^{-13}$.
The fourth iteration for the energy is found from the equation,

$\mathcal{E}_{4}=\frac{\mathcal{E}_{2}\mathfrak{D}_{3}-\mathcal{E}_{3}\mathfrak{D}_{2}}{\mathfrak{D}_{3}-\mathfrak{D}_{2}}=\varepsilon_{1}^{(0)}-1.474597699\times
10^{-8}\,\epsilon_{H},$

such that $\mathcal{E}_{3}-\mathcal{E}_{4}=2.6\times 10^{-16}\epsilon_{H}$.
Solving Eq. (56) numerically for $\mathcal{E}_{4}$, we find
$g_{b,-}(r),f_{b,-}(r)$ for $r<r_{b}$. Substituting
$g_{b,-}(r_{b}),f_{b,-}(r_{b})$, and $g_{+}(r_{b}),f_{+}(r_{b})$ with
$\mathcal{E}_{4}$ [see in Eq. (38)] into Eq. (60), we get
$\mathfrak{D}_{4}\equiv\mathfrak{D}(\mathcal{E}_{4})=-9.41473\times 10^{-13}$.
Since $\mathcal{E}_{3}-\mathcal{E}_{4}$ and $\mathfrak{D}_{4}$ are very small,
$\mathcal{E}_{4}$ is practically equal to the 3He ground state energy. —————–

##  Appendix D Nuclear magnetic moment distribution for special nuclei

This appendix discusses the magnetic moment distribution for nuclei with one
nucleon outside a closed nuclear shell.

###  D.1 Nuclear shell model

According to the nuclear shell model [59, 60], each nucleon in a nucleus moves
in a self-consistent field due to the other nucleons. The nuclear self-
consistent potential decreases rapidly outside the volume bounded by the
surface of the nucleus. The quantum state of the nucleus is described by
specifying the states of the individual nucleons. The self-consistent field is
spherically symmetric, and the center of symmetry is the center of mass of the
nucleus. Hence, the quantum states of the individual nucleons are parametrized
by a radial quantum number $n_{r}$, angular momentum quantum number $L$,
orbital angular momentum $J=L\pm 1/2$, and projection $M_{J}$ of the vector
$\mathbf{J}=\mathbf{L}+\mathbf{S}$. Nucleons fill shells according to the
Pauli principle such that the nucleon states are distributed among the groups
shown in Table 16 [59, 60]. For each group, the total number of proton or
neutron occupation is shown. Hence, the occupation of a group is completed
when the number $Z$ of protons or the number $N$ of neutrons is equal to one
of the numbers: 2, 8, 28, 50, 82, 126, $\ldots$. These numbers are called
magic numbers [60]. Nuclei with both $Z$ and $N$ being magic numbers, are
called double magic nuclei and have nuclear spin $I=0$.

shells | nucleons
---|---
$1s_{1/2}$ | 2
$1p_{3/2}$, $1p_{1/2}$ | 6
$1d_{5/2}$, $2s_{1/2}$, $1d_{3/2}$, $1f_{7/2}$ | 20
$2p_{3/2}$, $1f_{5/2}$, $2p_{1/2}$, $1g_{9/2}$ | 22
$2d_{5/2}$, $1g_{7/2}$, $1h_{11/2}$, $2d_{3/2}$, $3s_{1/2}$ | 32
$1h_{9/2}$, $2f_{7/2}$, $2i_{13/2}$, $2f_{5/2}$, $3p_{3/2}$, $3p_{1/2}$ | 44
Table 16: Groups of shells and the number of nucleon completing each group.
Here the quantum number $n_{r}$ before the letter is a radial quantum number,
and the letters $s$, $p$, $d$, $f$, $g$, $h$, $i$, $\dots$ refer to the
angular momentum quantum number, $L_{j}=0,1,2,3,4,5,6,\ldots$, where $j=p,n$
for protons and neutrons, and the subscript on the letters indicates the total
orbital angular momentum of the nucleon, $J_{j}=L_{j}\pm 1/2$.

A nucleon moves in an effective attractive potential
$V_{\mathrm{nucl}}^{(j)}(r)$ formed by all the other nucleons. The potential
$V_{\mathrm{nucl}}^{(j)}(r)$ is given by [59]

$V_{\mathrm{nucl}}^{(j)}(r)=-\frac{V_{j}}{e^{(r-R_{0})/a_{F}}+1}+U_{\mathrm{rep}}(r)\delta_{j,p},$
(138)

where the first term on the right hand side of Eq. (138) is the Woods-Saxon
attractive potential, and the second term, $U_{\mathrm{rep}}(r)$, is the
electrostatic repulsive interaction of a proton with the other $Z-1$ protons
in the nucleus. The Kronecker delta indicates that the Coulomb interaction is
present for protons and absent for neutrons. Here $R_{0}=1.25\,A^{1/3}$, and
$V_{j}$ is the depth of the potential well for protons ($j=p$) and neutrons
($j=n$). Note that $R_{0}>c_{F}$, where $c_{F}$ is the half-density nuclear
radius in the two-parameter Fermi nuclear charge distribution, Eq. (82), and
the difference $R_{0}-c_{F}$ is due to the finite range of the nuclear force.
The potential depth $V_{j}$ is

$\displaystyle V_{j}$ $\displaystyle=$ $\displaystyle
U_{0}+\eta_{j}U_{1}\,\frac{N-Z}{A},$ (139)

where $U_{0}=57~{}\text{MeV}$, $U_{1}=27~{}\text{MeV}$, $Z$ is the number of
protons in the nucleus, $N$ is the number of neutrons in the nucleus, and
$A=Z+N$ is the total number of nucleons in the nucleus. The second term in the
right hand side of Eq. (139) is the symmetry energy arising when $Z\neq N$,
where $\eta_{p}=1$ and $\eta_{n}=-1$ [59]. The second term in the right hand
side of Eq. (138), $U_{\mathrm{rep}}(r)$, is the electrostatic repulsion of a
proton with the other $Z-1$ protons in the nucleus. For the model (b) for the
charge distribution, $U_{\mathrm{rep}}(r)$ is

$U_{\mathrm{rep}}(r)=\left\\{\begin{array}[]{ccc}\displaystyle\frac{(Z-1)e^{2}}{2r_{b}}\bigg{(}3-\frac{r^{2}}{r_{b}^{2}}\bigg{)},&\text{for}&r\leq
r_{b},\\\
\displaystyle\frac{(Z-1)e^{2}}{r},&\text{for}&r>r_{b}.\end{array}\right.$

In addition, spin-orbit interaction of a nucleon in the self-consistent
potential is present, which is weaker than the strong interaction and can be
taken into account using perturbation theory.

A nucleon wave function $\Psi_{n_{r},L_{j},M_{L}}(\mathbf{r})$ in a shell with
the radial quantum number $n_{r}$, the angular momentum quantum number $L_{j}$
and the projection $M_{L}$ of the angular momentum $\mathbf{L}_{j}$ on the $z$
axis is found from the Schrödinger equation,

$\left[-\frac{\hbar^{2}}{2M_{j}}\,\nabla^{2}+V_{\mathrm{nucl}}^{(j)}(r)\right]\Psi_{n_{r},L_{j},M_{L}}(\mathbf{r})=\epsilon\Psi_{n_{r},L_{j},M_{L}}(\mathbf{r}),$
(140)

where $M_{j}$ is the nucleon mass. $\Psi_{n_{r},L_{j},M_{L}}(\mathbf{r})$ can
be written as,

$\Psi_{n_{r},L_{j},M_{L}}(\mathbf{r})=\psi_{n_{r},L_{j}}(r)Y_{L_{j},M_{L}}(\theta,\phi),$

where $\psi_{n_{r},L_{j}}(r)$ is a radial wave function and
$Y_{L_{j},M_{L}}(\theta,\phi)$ is a spherical harmonic.
$\psi_{n_{r},L_{j}}(r)$ is found from the equation,

$\displaystyle-\frac{\hbar^{2}}{2M_{j}r^{2}}\,\frac{d}{dr}\bigg{(}r^{2}\,\frac{d\psi_{n_{r},L}(r)}{dr}\bigg{)}+\frac{\hbar^{2}L_{j}(L_{j}+1)}{2M_{j}r^{2}}\,\psi_{n_{r},L}(r)$
$\displaystyle\;\;+V_{\mathrm{nucl}}^{(j)}(r)\psi_{n_{r},L_{j}}(r)=\epsilon\psi_{n_{r},L_{j}}(r),\quad$
(141)

where the potential $V_{\mathrm{nucl}}^{(j)}(r)$ is given in Eq. (138).

We shall now apply the shell model to describe the nuclear magnetic moment
distribution for 207Pb and 209Bi isotopes. For 207Pb, $R_{0}=7.394$ fm,
$c_{F}=6.591$ fm, and the difference is $R_{0}-c_{F}=0.803$ fm. For 209Bi,
$R_{0}=7.418$ fm, $c_{F}=6.765$ fm, and the difference is $R_{0}-c_{F}=0.653$
fm.

###  D.2 ${}^{207}\text{Pb}$ isotope nuclear density

The nucleus of 207Pb consists of $Z=82$ protons and $N=125$ neutrons. Filling
of the nuclear shells by the nucleons is illustrated in Fig. 11(a): The number
of protons is a magic number, i.e., the proton close the lowest shells and
form a singlet state. The number of neutrons is one less than the magic number
126\. The lowest 5 shells are closed by the neutrons, and the outer shell is
partially filled with 43 neutrons. We consider the outer shell as a closed
shell, partially filled by 1 hole, and derive the wave function for this hole.
The nuclear spin $I$ is equal to the orbital angular momentum of the hole. The
nuclear spin of 207Pb nucleus is $I=1/2$. The hole is on the $3p_{1/2}$ orbit
with the radial quantum number $n_{r}=3$, angular momentum quantum number
$L=1$, and the orbital angular momentum quantum number $I=1/2$.

Figure 11: (a) Protons (red) and neutrons (blue) in the 207Pb nucleus fill the
nuclear shells. The closed shells consist of 2, 6, 20, 22, 32 and 44 nucleons.
The number of protons $Z=82$ is a magic number, and the number of neutrons
$N=125$ is one less than the magic number 126. The empty circle shows the
vacancy (a hole) in the outer shell. The nuclear spin $I=1/2$ is the orbital
angular momentum of the hole on the $3p_{1/2}$ shell. (b) Nuclear magnetic
moment density for 207Pb isotope.

Wave function $\Psi_{3,1,M_{L}}(\mathbf{r})$ of a hole with the radial quantum
number $n_{r}=3$, angular momentum number $L=1$ and the projection $M_{L}$ of
$\mathbf{L}$ on the $z$ axis can be written as,

$\Psi_{n_{r},L,M_{L}}(\mathbf{r})=\psi_{n_{r},L}(r)Y_{L,M_{L}}(\theta,\phi),$

where $\psi_{n_{r},L}(r)$ is a radial wave function and
$Y_{L,M_{L}}(\theta,\phi)$ is a spherical harmonics. $\psi_{n_{r},L}(r)$ is
found from the equation,

$\displaystyle-\frac{\hbar^{2}}{2M_{n}r^{2}}\,\frac{d}{dr}\bigg{(}r^{2}\,\frac{d\psi_{n_{r},L}(r)}{dr}\bigg{)}+\frac{\hbar^{2}L(L+1)}{2M_{n}r^{2}}\,\psi_{n_{r},L}(r)$
$\displaystyle\;\;+V_{\mathrm{nucl}}^{(n)}(r)\psi_{n_{r},L}(r)=\epsilon\psi_{n_{r},L}(r),$
(142)

where $M_{n}$ is the neutron mass, and the potential
$V_{\mathrm{nucl}}^{(n)}(r)$ is given in Eq. (138) with $j=n$. We apply the
Mathematica command NDEigensystem to solve Eq. (141) and find radial wave
functions $\psi_{n_{r},1}(r)$ and energies of the $1p$, $2p$ and $3p$ states.

The magnetic moment density $\varrho_{m}(r)$ is

$\varrho_{m}(r)=\big{|}\psi_{3,1}(r)\big{|}^{2}.$ (143)

The magnetic moment density $\varrho_{m}(r)$ in Eq. (143) is plotted in Fig.
11. $\varrho_{m}(r)$ has three nodes, at $r=0$, $0.59821\,r_{N}$ and
$1.03036\,r_{N}$, and three maxima, $\varrho_{m}(r_{1})=11.7192\,r_{N}^{-3}$
at $r_{1}=0.277111\,r_{N}$, $\varrho_{m}(r_{2})=1.74229\,r_{N}^{-3}$ at
$r_{2}=0.79112\,r_{N}$, and $\varrho_{m}(r_{3})=0.824837\,r_{N}^{-3}$ at
$r_{3}=1.25533r_{N}$. Note that the highest peak position $r_{1}$ in
$\varrho(r)$ is far below $r_{N}$.