Relativistic correction of the Coulomb interaction in the local density approximation for energies and radii in doubly-magic nuclei

Abstract

Effects of the relativistic correction of the Coulomb interaction on doubly-magic nuclei are discussed with Skyrme Hartree–Fock calculations. The relativistic correction is treated by using the local density approximation. It is found that the correction to the total energy is about |$2.4 \, \mathrm{MeV}$| for |${}^{208}_{}\mathrm{Pb}$|⁠, while the proton and neutron radii do not change significantly. This difference is larger than the difference of the Coulomb exchange (Fock) energy calculated with the local density (Hartree–Fock–Slater) approximation and that with the exact treatment and the neutron finite-size effect. Effects of the correction are also compared to the correction due to the vacuum polarization. It is shown that the two contributions to the total energy are comparable in light nuclei, but the latter dominates in heavy nuclei, while the contribution of the relativistic correction to the total energy is non-negligible compared to the target accuracy of the DFT calculation.

1. Introduction

Atomic nuclei are composed of protons and neutrons, which interact via nuclear and electromagnetic (EM) interactions. Since the nuclear interaction is much stronger than the EM one, the former determines most of the properties of atomic nuclei. Nevertheless, the EM interaction plays an important role when one focuses on nuclear properties related to the isospin symmetry breaking [1–10]. For more details, see, e.g., Ref. [11].

The Hartree–Fock method or the density functional theory (DFT) [12–14] is a powerful method to calculate quantum many-fermion problems, including atomic nuclei. In DFT, an energy density functional (EDF), which contains information on interactions, governs the accuracy of calculation. Although the Coulomb interaction is rather simple, first-principle EDFs, even for electronic systems, have been derived approximately using the local density approximation (LDA) [15–17], and the effects of the density gradient have been included empirically or phenomenologically [18–20]. In nuclear physics, the situation is more complicated that the nuclear interaction in medium is still under discussion, and the derivation of ab initio nuclear EDFs is currently an important issue [21–27]. Therefore, most nuclear EDFs have been constructed phenomenologically by fitting to experimental data and different EDFs may give different results [28–30]. The accuracy of the total energies obtained by widely used EDFs ranges approximately from 1–|$5 \, \mathrm{MeV}$| [31,32], while those of empirical mass formulae and machine learning ranges from 100–|$500 \, \mathrm{keV}$| [33–36], which almost matches the desired accuracy for the study of r-process nucleosynthesis [37].

To fit the parameters of nuclear EDFs properly, the EM contribution to the nuclear properties should be subtracted from experimental data properly. Motivated by this issue, in a series of works [38–40], accurate treatment of the Coulomb interaction in a Skyrme Hartree–Fock (SHF) calculation [41], one of the branches of nuclear DFT [42–46], was proposed. In the series of works, the generalized gradient approximation (GGA) for the Coulomb exchange energy was introduced to reproduce the exact-Fock energy with low numerical costs [38,39]. In addition, more higher-order contributions were taken into account: the spatial charge distributions of nucleons and the vacuum polarization between two protons [7,11,40,47,48]. From the point of view of quantum field theory, the finite-size effects, i.e., effects originating from the spatial charge distributions of nucleons, correspond to the correction of vertices and the vacuum polarization corresponds to the higher-order correction of the Coulomb interaction in terms of the coupling constant α ≈ 1/137 (the fine-structure constant). Although these effects are just corrections of the Coulomb interaction, which itself is already subdominant for the nuclear total energies, it has been shown that they are non-negligible compared to the desired accuracy; for instance, the proton finite-size effect, the neutron one, and the vacuum polarization contribute to the total energy of |${}^{208}_{}\mathrm{Pb}$| by |$-8.2 \, \mathrm{MeV}$|⁠, |$+1.2 \, \mathrm{MeV}$|⁠, and |$+3.7 \, \mathrm{MeV}$|⁠, respectively [40].

Then, it is natural to consider corrections with the same order with respect to α, i.e., the finite-light-speed correction, the self-energy, and the two-photon exchange process. Although the self-energy contributes to the energy of electrons in atoms by the same magnitude as vacuum polarization, it may be weak for the energy of nucleons in atomic nuclei since the nucleon mass is heavy; the two-photon exchange may also be weak since it is also inversely proportional to the nucleon mass [49]. In contrast, the finite-light-speed effect does not depend on the nucleon mass, and it is known to be non-negligible in atomic structure [50,51]. Thus, in this work, the finite-light-speed correction of the Coulomb interaction is considered in the nuclear structure calculation in order to investigate whether its effect is small. Even if its effect is small, it is worthwhile to confirm it.

The finite-light-speed correction of the Coulomb interaction is often referred to as the relativistic correction [50,52,52–54]. The Coulomb interaction is instantaneous, while the velocity of photons, which mediate the EM interaction, is finite. Once the finite light speed is considered, the Coulomb interaction should be, therefore, corrected, and this correction is known as the Breit correction [55–57]. In atomic physics, the Breit correction dominates more than the vacuum polarization of the Coulomb field formed by the atomic nucleus [58], and thus the Breit correction should be considered once the vacuum polarization is considered. In nuclear physics, the vacuum polarization discussed in Ref. [40] originates from the interaction between two protons, which is different from that in atomic physics. Hence, it is interesting to compare which contribution is dominant in the nuclear total energy, the Breit correction or the vacuum polarization.

Another interesting point is a comparison between the Breit correction and the correction due to the density gradient, i.e., the difference between the LDA Coulomb energy and the GGA Coulomb energy. The Breit correction has been taken into account for an exchange [59] or exchange-correlation [60] EDF in the LDA level, and has already been taken into account fully for atomic DFT [51]. It was found in Ref. [51] that the Breit correction contributes to the total energies of atoms by almost the same magnitude as—but in an opposite direction to—the density gradient effect. In the context of nuclear physics, it was shown in Refs. [39,61,62] that difference between the Coulomb exchange energy calculated by the exact-Fock term or GGA and that in LDA is non-negligible in nuclear binding energies. Thus, it should be tested to what extent the Breit correction contributes to nuclear binding energies.

The covariant density functional theory, the nuclear DFT in the relativistic scheme, is based on the covariant Lagrangian with nucleon, meson, and photon fields. Hence, effects beyond the Coulomb interaction, such as magnetic interaction and retardation, have often been considered [42,43,63,64]. Therefore, it is interesting to take into account the relativistic correction of the Coulomb interaction in non-relativistic SHF calculations.

The relativistic correction of the Coulomb interaction has already been discussed in Refs. [65,66]. As discussed later, the correction used in these works is the same as the non-relativistic reduced form of the Breit correction for electronic systems [67,68], while anomalous magnetic moments of nucleons are, of course, considered. The aim of this paper is the modern calculation of Refs. [65,66]: deriving the EDF for such corrections and performing the self-consistent calculation. I shall discuss contributions to proton and neutron radii, total energies, and single-particle energies. The LDA is used to consider the relativistic correction; some terms of the correction vanish in the approximation, as will be discussed. Note that it seems that the “Breit correction” refers to the whole terms of the non-relativistic reduced form of the correction (5) in the context of atomic physics, while it refers to only a term of the correction in the context of nuclear physics. In order to avoid any confusion, hereinafter, I will use the term “relativistic correction”, instead of “Breit correction”.

This paper is organized as follows: In Sect. 2, the theoretical framework for the relativistic correction is shown. In Sect. 3, the relativistic correction to the total energy is estimated analytically, and in Sect. 4, numerical results are shown. Finally, in Sect. 5, the conclusion and perspectives are given.

2. Theoretical framework

In this section, I show the theoretical framework for the relativistic correction and the DFT. Throughout this paper, a unit with 4πε₀ = 1 is used and c denotes the speed of light.

Quantum many-body systems can be described by using the Hamiltonian

$$\begin{eqnarray*} H = T + V_{{{{{\rm{ext}}}}}} + V_{{{{{\rm{int}}}}}}, \end{eqnarray*} $$

(1)

where T, |$V_{{{{{\rm{ext}}}}}}$|⁠, and |$V_{{{{{\rm{int}}}}}}$| denote the kinetic energy operator, an external potential, and an interaction, respectively. In the case of electronic systems, |$V_{{{{{\rm{int}}}}}}$| is identical to the EM interaction |$V_{{{{{\rm{EM}}}}}}$|⁠, while in the case of atomic nuclei, |$V_{{{{{\rm{int}}}}}}$| includes a nuclear interaction |$V_{{{{{\rm{nucl}}}}}}$| and the EM interaction |$V_{{{{{\rm{EM}}}}}}$| with |$V_{{{{{\rm{ext}}}}}} \equiv 0$|⁠. In most works, the Coulomb interaction between protons

$$\begin{eqnarray*} V_{{{{{\rm{Coul}}}}}} \left( {\boldsymbol {r}}_j, {\boldsymbol {r}}_k \right) = \frac{e^2}{r_{jk}} \end{eqnarray*} $$

(2)

is used for |$V_{{{{{\rm{EM}}}}}}$|⁠, where |${\boldsymbol {r}}_j$| and |${\boldsymbol {r}}_k$| are the spatial coordinates of the protons j and k, |${\boldsymbol {r}}_{jk} = {\boldsymbol {r}}_j - {\boldsymbol {r}}_k$|⁠, and |$r_{jk} = \left| {\boldsymbol {r}}_{jk} \right|$|⁠. According to quantum electrodynamics, the leading-order instantaneous interaction, in which a photon mediates between two particles with c = ∞, is the Coulomb interaction.

2.1. Breit correction for electronic systems

In electronic systems, the next-leading-order interaction (O(1/c²)) in the Coulomb gauge is called the Breit correction, whose form is [50,55–57,69]

$$\begin{eqnarray*} V_{{{{{\rm{Breit}}}}}} \left( {\boldsymbol {r}}_j, {\boldsymbol {r}}_k \right) = - \left[ \frac{c {\boldsymbol {\alpha }}_j \cdot c {\boldsymbol {\alpha }}_k}{2c^2 r_{jk}} + \frac{\left( c {\boldsymbol {\alpha }}_j \cdot {\boldsymbol {r}}_{jk} \right) \left( c {\boldsymbol {\alpha }}_k \cdot {\boldsymbol {r}}_{jk} \right)}{2c^2 r_{jk}^3} \right], \end{eqnarray*} $$

(3)

where α is 4 × 4 matrix defined by

$$\begin{eqnarray*} {\boldsymbol {\alpha }} = \left( {\begin{pmatrix}O_2 & \quad \sigma _x \\ \sigma _x & \quad O_2 \end{pmatrix}}, {\begin{pmatrix}O_2 & \quad \sigma _y \\ \sigma _y & \quad O_2 \end{pmatrix}}, {\begin{pmatrix}O_2 & \quad \sigma _z \\ \sigma _z & \quad O_2 \end{pmatrix}} \right) \end{eqnarray*} $$

(4)

with the Pauli matrices σ_x, σ_y, and σ_z and the 2 × 2 zero matrix O₂. In isolated atoms, the vacuum polarization of the Coulomb field appears at a higher order than the Breit correction [58]. The Breit correction includes the magnetic (current–current) interaction and the retardation in the Coulomb gauge [69]. In the Lorenz (covariant) gauge, the next-leading-order interaction corresponds to the magnetic interaction and the retardation [69–71].

After the Foldy–Wouthuysen–Tani transformation [72–75], the Breit correction (3) in the non-relativistic scheme is obtained as [67,68]

$$\begin{eqnarray*} \tilde{V}_{{{{{\rm{Breit}}}}}} \left( {\boldsymbol {r}}_j, {\boldsymbol {s}}_j, {\boldsymbol {r}}_k, {\boldsymbol {s}}_k \right) &=& - \frac{\pi \hbar ^2 e^2}{m_e^2 c^2} \delta \left( {\boldsymbol {r}}_{jk} \right) \\ &&-\,\, \frac{e^2}{2 m_e^2 c^2} {\boldsymbol {p}}_j \cdot \left( \frac{1}{r_{jk}} + \frac{{\boldsymbol {r}}_{jk} {\boldsymbol {r}}_{jk}}{r_{jk}^3} \right) \cdot {\boldsymbol {p}}_k \\ &&-\,\, \frac{8 \pi \hbar ^2 e^2}{3 m_e^2 c^2} \delta \left( {\boldsymbol {r}}_{jk} \right) {\boldsymbol {s}}_j \cdot {\boldsymbol {s}}_k \\ &&-\,\, \frac{\hbar ^2 e^2}{m_e^2 c^2} {\boldsymbol {s}}_j \cdot \left( \frac{3 {\boldsymbol {r}}_{jk} {\boldsymbol {r}}_{jk}}{r_{jk}^5} - \frac{1}{r_{jk}^3} \right) \cdot {\boldsymbol {s}}_k \\ &&+\,\, \frac{\hbar ^2 e^2}{m_e^2 c^2} \frac{1}{r_{jk}^3} {\boldsymbol {s}}_j \cdot \left[ {\boldsymbol {r}}_{jk} \times \left( 2 {\boldsymbol {p}}_k - {\boldsymbol {p}}_j \right) \right], \end{eqnarray*} $$

(5a),(5b),(5c),(5d),(5e)

where |$m_e \approx 0.511 \, \mathrm{MeV} / c^2$| is the electron mass [76] and |${\boldsymbol {p}}_j$| is the momentum of the particle j. The first (5a) and second (5b) terms correspond to the Darwin term for the electron–electron Coulomb interaction, which is related to Zitterbewegung [70], and the retardation of the Coulomb interaction, respectively; the third and fourth terms (5c), (5d) correspond to the spin-magnetic interactions; and the fifth term (5e) corresponds to the electron–electron spin-orbit interaction [60,67,68]. Thus, in the Schrödinger scheme, if one is eager to consider the Breit correction, Eq. (5) should be considered above the Coulomb interaction (2), i.e.,

$$\begin{eqnarray*} V_{{{{{\rm{EM}}}}}} = \frac{1}{2} \sum _{j, \, k} \left[ V_{{{{{\rm{Coul}}}}}} \left( {\boldsymbol {r}}_j, {\boldsymbol {r}}_k \right) + \tilde{V}_{{{{{\rm{Breit}}}}}} \left( {\boldsymbol {r}}_j, {\boldsymbol {s}}_j, {\boldsymbol {r}}_k, {\boldsymbol {s}}_k \right) \right], \end{eqnarray*} $$

(6)

which will be referred to as the Coulomb–Breit interaction.

2.2. Beyond Coulomb interaction for nuclear systems

Corrections beyond the Coulomb interaction for nuclear systems were discussed several decades ago [65,66]. The correction was

$$\begin{eqnarray*} \tilde{V}_{{{{{\rm{Rel}}}}}} \left( {\boldsymbol {r}}_j, {\boldsymbol {s}}_j, {\boldsymbol {r}}_k, {\boldsymbol {s}}_k \right) &=& - \frac{e_{\tau _j} e_{\tau _k}}{2 M^2 c^2} {\boldsymbol {p}}_j \cdot \left( \frac{1}{r_{jk}} + \frac{{\boldsymbol {r}}_{jk} {\boldsymbol {r}}_{jk}}{r_{jk}^3} \right) \cdot {\boldsymbol {p}}_k \\ &&-\,\, \frac{2 \hbar \mu _{\tau _j} e_{\tau _k}}{M c^2} \frac{1}{r_{jk}^3} {\boldsymbol {s}}_j \cdot \left( {\boldsymbol {r}}_{jk} \times {\boldsymbol {p}}_k \right) + \frac{2 \hbar \mu _{\tau _k} e_{\tau _j}}{M c^2} \frac{1}{r_{jk}^3} {\boldsymbol {s}}_j \cdot \left( {\boldsymbol {r}}_{jk} \times {\boldsymbol {p}}_j \right) \\ &&-\,\, \frac{4 \mu _{\tau _j} \mu _{\tau _k}}{c^2} {\boldsymbol {s}}_j \cdot \left( \frac{3 {\boldsymbol {r}}_{jk} {\boldsymbol {r}}_{jk}}{r_{jk}^5} - \frac{1}{r_{jk}^3} \right) \cdot {\boldsymbol {s}}_k \\ &&-\,\, \frac{32 \pi }{3} \frac{\mu _{\tau _j} \mu _{\tau _k}}{c^2} \delta \left( {\boldsymbol {r}}_{jk} \right) {\boldsymbol {s}}_j \cdot {\boldsymbol {s}}_k \\ &&-\,\, \frac{e_{\tau _k} \hbar }{Mc^2} \left( 2 \mu _{\tau _j} - \frac{e_{\tau _j} \hbar }{2 M} \right) \frac{1}{r_{jk}^3} {\boldsymbol {s}}_j \cdot \left( {\boldsymbol {r}}_{jk} \times {\boldsymbol {p}}_j \right) \\ &&+\,\, \frac{e_{\tau _j} \hbar }{Mc^2} \left( 2 \mu _{\tau _k} - \frac{e_{\tau _k} \hbar }{2 M} \right) \frac{1}{r_{jk}^3} {\boldsymbol {s}}_k \cdot \left( {\boldsymbol {r}}_{jk} \times {\boldsymbol {p}}_k \right) \\ &&-\,\, \frac{\hbar \pi }{Mc^2} \left[ e_{\tau _k} \left( 2 \mu _{\tau _j} - \frac{e_{\tau _j} \hbar }{2M} \right) + e_{\tau _j} \left( 2 \mu _{\tau _k} - \frac{e_{\tau _k} \hbar }{2M} \right) \right] \delta \left( {\boldsymbol {r}}_{jk} \right) \\ &=& - \delta _{p \tau _j} \delta _{p \tau _k} \frac{e^2}{2 M^2 c^2} {\boldsymbol {p}}_j \cdot \left( \frac{1}{r_{jk}} + \frac{{\boldsymbol {r}}_{jk} {\boldsymbol {r}}_{jk}}{r_{jk}^3} \right) \cdot {\boldsymbol {p}}_k \nonumber \\ &&-\,\, \frac{\hbar ^2 e^2}{M^2 c^2} \frac{1}{r_{jk}^3} {\boldsymbol {s}}_j \cdot \left[ {\boldsymbol {r}}_{jk} \times \left( \tilde{\mu }_{\tau _j} \delta _{p \tau _k} {\boldsymbol {p}}_k - \tilde{\mu }_{\tau _k} \delta _{p \tau _j} {\boldsymbol {p}}_j \right) \right] \nonumber \\ &&-\,\, \tilde{\mu }_{\tau _j} \tilde{\mu }_{\tau _k} \frac{\hbar ^2 e^2}{M^2 c^2} {\boldsymbol {s}}_j \cdot \left( \frac{3 {\boldsymbol {r}}_{jk} {\boldsymbol {r}}_{jk}}{r_{jk}^5} - \frac{1}{r_{jk}^3} \right) \cdot {\boldsymbol {s}}_k \nonumber \\ &&-\,\, \tilde{\mu }_{\tau _j} \tilde{\mu }_{\tau _k} \frac{8 \pi \hbar ^2 e^2}{3 M^2 c^2} \delta \left( {\boldsymbol {r}}_{jk} \right) {\boldsymbol {s}}_j \cdot {\boldsymbol {s}}_k \nonumber \\ &&-\,\, \delta _{p \tau _k} \frac{\hbar ^2 e^2}{2 M^2 c^2} \left( 2 \tilde{\mu }_{\tau _j} - \delta _{p \tau _j} \right) \frac{1}{r_{jk}^3} {\boldsymbol {s}}_j \cdot \left( {\boldsymbol {r}}_{jk} \times {\boldsymbol {p}}_j \right) \nonumber \\ &&+\,\, \delta _{p \tau _j} \frac{\hbar ^2 e^2}{2 M^2 c^2} \left( 2 \tilde{\mu }_{\tau _k} - \delta _{p \tau _k} \right) \frac{1}{r_{jk}^3} {\boldsymbol {s}}_k \cdot \left( {\boldsymbol {r}}_{jk} \times {\boldsymbol {p}}_k \right) \nonumber \\ &&-\,\, \frac{\pi \hbar ^2 e^2}{2 M^2 c^2} \left[ \delta _{p \tau _k} \left( 2 \tilde{\mu }_{\tau _j} - \delta _{p \tau _j} \right) + \delta _{p \tau _j} \left( 2 \tilde{\mu }_{\tau _k} - \delta _{p \tau _k} \right) \right] \delta \left( {\boldsymbol {r}}_{jk} \right), \nonumber \end{eqnarray*} $$

(7a),(7b),(7c),(7d),(7e),(7f),(7g)

where |$M \approx 938.919 \, \mathrm{MeV} / c^2$| is the nucleon mass [76], δ_pτ is the Kronecker delta, which is equal to 1 for τ = p and 0 for τ = n, e_τ and μ_τ are the charge and magnetic moment of a nucleon τ, and |$\tilde{\mu }_{\tau }$| is the magnetic moment of the unit of the nuclear magneton, i.e. [76],

$$\begin{eqnarray*} \mu _{\tau } &=& \tilde{\mu }_{\tau } \frac{e \hbar }{2M}, \end{eqnarray*} $$

(8a)

$$\begin{eqnarray*} \tilde{\mu }_p &=& 2.792\,847\,344\,63, \end{eqnarray*} $$

(8b)

$$\begin{eqnarray*} \tilde{\mu }_n &=& -1.913\,042\,73. \end{eqnarray*} $$

(8c)

If one assumes that these two particles are electrons instead of nucleons and, accordingly, m_e, −e, and μ_e are used instead of M, e_τ, and μ_τ, one can find that Eq. (7) is identical to Eq. (5). Here, the relation between the g-factor and the magnetic moment μ_e

$$\begin{eqnarray*} \mu _e = \frac{e \hbar }{4 m_e} g_e \end{eqnarray*} $$

(9)

and g_e = 2 is used, apart from the second term of Eq. (7b), in which the g-factor for the orbital motion of electrons g_e = 1 is used. In more detail, the correspondence between Eqs. (5) and (7) is shown in Table 1. In the case of electronic systems, Eqs. (7e) and (7f) cancel each other out.

Table 1.

Correspondence between Eqs. (5) and (7). For convenience, correspondence to Eq. (2) of Ref. [67] is also shown.

Eq. (7)	Eq. (5)	Eq. (2) of Ref. [67]	Origin
Eq. (7a)	Eq. (5b)	Eq. (2b)	Retardation
Eq. (7b)	Eq. (5e)	Eqs. (2c) and (2d)	Spin-orbit interaction
Eq. (7c)	Eq. (5d)	Eq. (2e)	Spin-magnetic interaction (non-mutually penetrating)
Eq. (7d)	Eq. (5c)	Eq. (2f)	Spin-magnetic interaction (mutually penetrating)
Eq. (7g)	Eq. (5a)	Eq. (2g)	Darwin term

2.3. Density functional theory for relativistic correction

In nuclear DFT, the ground-state energy is obtained by

$$\begin{eqnarray*} E_{{{{{\rm{gs}}}}}} = T_0 + E_{{{{{\rm{nucl}}}}}} \left[ \rho _p^{{{{{\rm{gs}}}}}}, \rho _n^{{{{{\rm{gs}}}}}} \right] + E_{{{{{\rm{EM}}}}}} \left[ \rho _p^{{{{{\rm{gs}}}}}} \right], \end{eqnarray*} $$

(10)

where T₀ is the Kohn–Sham kinetic energy, |$E_{{{{{\rm{nucl}}}}}}$| and |$E_{{{{{\rm{EM}}}}}}$| are EDFs of nuclear and EM interactions, respectively, and |$\rho ^p_{{{{{\rm{gs}}}}}}$| and |$\rho ^n_{{{{{\rm{gs}}}}}}$| are the ground-state densities of protons and neutrons. If the finite charge distributions of nucleons are considered, |$E_{{{{{\rm{EM}}}}}}$| should depend on the ground-state charge distribution, instead of |$\rho _p^{{{{{\rm{gs}}}}}}$| [40].

If one only considers the Coulomb interaction to |$V_{{{{{\rm{EM}}}}}}$|⁠, |$E_{{{{{\rm{EM}}}}}}$| in LDA reads [77]

$$\begin{eqnarray*} E_{{{{{\rm{EM}}}}}} \left[ \rho \right] &=& E_{{{{{\rm{Coul}}}}}}^{{{{{\rm{H}}}}}} \left[ \rho \right] + E_{{{{{\rm{Coul}}}}}}^{{{{{\rm{x}}}}}} \left[ \rho \right], \end{eqnarray*} $$

(11a)

$$\begin{eqnarray*} E_{{{{{\rm{Coul}}}}}}^{{{{{\rm{H}}}}}} \left[ \rho \right] &=& \frac{e^2}{2} \iint \frac{\rho \left( {\boldsymbol {r}} \right) \rho \left( {\boldsymbol {r}}^{\prime } \right)}{\left| {\boldsymbol {r}} - {\boldsymbol {r}}^{\prime } \right|} \, d {\boldsymbol {r}} \, d {\boldsymbol {r}}^{\prime }, \end{eqnarray*} $$

(11b)

$$\begin{eqnarray*} E_{{{{{\rm{Coul}}}}}}^{{{{{\rm{x}}}}}} \left[ \rho \right] &=& - \frac{3 e^2}{4} \left( \frac{3}{\pi } \right)^{1/3} \int \left[ \rho \left( {\boldsymbol {r}} \right) \right]^{4/3} \, d {\boldsymbol {r}}, \end{eqnarray*} $$

(11c)

where the density ρ is the electron density ρ_e for electronic systems and the proton density ρ_p for nuclear systems.

Recently, an EDF for the Breit correction (5) for electronic systems in the LDA was developed in Refs. [51,60]. In this paper, I consider the relativistic correction (7) in the nuclear DFT using the idea of such an EDF and the correspondence shown in Table 1. This newly developed EDF enables one to perform the self-consistent calculation. To isolate the relativistic correction, I shall not consider the nucleon finite-size effects¹.

The Hartree and LDA exchange EDFs for the Breit correction in the electronic systems (5) were proposed in Refs. [51] and [59], respectively. Their forms are

$$\begin{eqnarray*} E_{{{{{\rm{Breit}}}}}}^{{{{{\rm{H}}}}}} \left[ \rho \right] &=& - \frac{\pi \hbar ^2 e^2}{2 m_e^2 c^2} \int \left[ \rho \left( {\boldsymbol {r}} \right) \right]^2 \, d {\boldsymbol {r}}, \end{eqnarray*} $$

(12a)

$$\begin{eqnarray*} E_{{{{{\rm{Breit}}}}}}^{{{{{\rm{x}}}}}} \left[ \rho \right] &=& \frac{3 \pi \hbar ^2 e^2}{2 m_e^2 c^2} \int \left[ \rho \left( {\boldsymbol {r}} \right) \right]^2 \, d {\boldsymbol {r}}, \end{eqnarray*} $$

(12b)

respectively, and, accordingly, the total correction is

$$\begin{eqnarray*} E_{{{{{\rm{Breit}}}}}}^{{{{{\rm{tot}}}}}} \left[ \rho \right] &=& E_{{{{{\rm{Breit}}}}}}^{{{{{\rm{H}}}}}} \left[ \rho \right] + E_{{{{{\rm{Breit}}}}}}^{{{{{\rm{x}}}}}} \left[ \rho \right] \nonumber \\ &=& \frac{\pi \hbar ^2 e^2}{m_e^2 c^2} \int \left[ \rho \left( {\boldsymbol {r}} \right) \right]^2 \, d {\boldsymbol {r}}. \end{eqnarray*} $$

(12c)

The idea of the LDA is that the energy density is approximated to that of homogeneous systems; hence, contributions of Eqs. (5c), (5d), and (5e) vanish. Consequently, Eqs. (12a) and (12b) contain only effects originating from Eqs. (5a) and (5b) [60]. In addition, due to the symmetry, Eq. (12a) does not contain the effect of the retardation (5b).

I shall construct an EDF for Eq. (7) in the LDA, using the knowledge of Eqs. (12a) and (12b). Because of the idea of the LDA, the EDF contains only the effect of Eqs. (7a) and (7g), as in the electronic systems. Note that the spin-orbit interaction originating from the magnetic form factors of nucleons (7b), (7e), (7f) have been discussed in Refs. [7,40,78]; the spin-magnetic interaction (7c), (7d) has not been considered yet for the self-consistent calculation and it cannot be considered in the LDA, which is left for future study.

Then, I will construct the LDA EDF for Eq. (7). Referring to the result for Eq. (12), one can find that Eq. (7g) contributes both the Hartree and the exchange terms, whereas Eq. (7a) contributes the exchange term only.

The Hartree and the exchange EDFs of Eq. (7g), respectively, read

$$\begin{eqnarray*} E_{{{{{\rm{Darwin}}}}}}^{{{{{\rm{H}}}}}} \left[ \rho _p, \rho _n \right] &=& \int \left\lbrace \frac{C_{pp}}{2} \left[ \rho _p \left( {\boldsymbol {r}} \right) \right]^2 + C_{pn} \rho _p \left( {\boldsymbol {r}} \right) \rho _n \left( {\boldsymbol {r}} \right) + \frac{C_{nn}}{2} \left[ \rho _n \left( {\boldsymbol {r}} \right) \right]^2 \right\rbrace \, d {\boldsymbol {r}} \nonumber \\ &=& \int \left\lbrace \frac{C_{pp}}{2} \left[ \rho _p \left( {\boldsymbol {r}} \right) \right]^2 + C_{pn} \rho _p \left( {\boldsymbol {r}} \right) \rho _n \left( {\boldsymbol {r}} \right) \right\rbrace \, d {\boldsymbol {r}}, \end{eqnarray*} $$

(13a)

$$\begin{eqnarray*} E_{{{{{\rm{Darwin}}}}}}^{{{{{\rm{x}}}}}} \left[ \rho _p, \rho _n \right] &=& - \int \left\lbrace \frac{C_{pp}}{4} \left[ \rho _p \left( {\boldsymbol {r}} \right) \right]^2 + \frac{C_{nn}}{4} \left[ \rho _p \left( {\boldsymbol {r}} \right) \right]^2 \right\rbrace \, d {\boldsymbol {r}} \nonumber \\ &=& - \int \frac{C_{pp}}{4} \left[ \rho _p \left( {\boldsymbol {r}} \right) \right]^2 \, d {\boldsymbol {r}} , \end{eqnarray*} $$

(13b)

where

$$\begin{eqnarray*} C_{\tau _j \tau _k} = - \frac{\pi \hbar ^2 e^2}{2 M^2 c^2} \left[ \delta _{p \tau _k} \left( 2 \tilde{\mu }_{\tau _j} - \delta _{p \tau _j} \right) + \delta _{p \tau _j} \left( 2 \tilde{\mu }_{\tau _k} - \delta _{p \tau _k} \right) \right] \end{eqnarray*} $$

(14)

and accordingly

$$\begin{eqnarray*} C_{pp} &=& - \frac{\pi \hbar ^2 e^2}{M^2 c^2} \left( 2 \tilde{\mu }_p - 1 \right), \end{eqnarray*} $$

(15a)

$$\begin{eqnarray*} C_{pn} &=& - \frac{\pi \hbar ^2 e^2}{M^2 c^2} \tilde{\mu }_n, \end{eqnarray*} $$

(15b)

$$\begin{eqnarray*} C_{nn} &=& 0. \end{eqnarray*} $$

(15c)

The exchange EDF for Eq. (7a) is constructed from Eq. (12b). Equation (12b) includes not only a contribution from Eq. (5b) but also that from Eq. (5a). Hence, the exchange EDF for Eq. (5b) reads

$$\begin{eqnarray*} \frac{3 \pi \hbar ^2 e^2}{2 m_e^2 c^2} \int \left[ \rho \left( {\boldsymbol {r}} \right) \right]^2 \, d {\boldsymbol {r}} - \frac{1}{4} \frac{\pi \hbar ^2 e^2}{m_e^2 c^2} \int \left[ \rho \left( {\boldsymbol {r}} \right) \right]^2 \, d {\boldsymbol {r}} = \frac{5 \pi \hbar ^2 e^2}{4 m_e^2 c^2} \int \left[ \rho \left( {\boldsymbol {r}} \right) \right]^2 \, d {\boldsymbol {r}}, \end{eqnarray*} $$

(16)

since the Darwin EDF for electronic systems is

$$\begin{eqnarray*} - \frac{1}{4} \frac{\pi \hbar ^2 e^2}{m_e c^2} \int \left[ \rho \left( {\boldsymbol {r}} \right) \right]^2 \, d {\boldsymbol {r}}. \end{eqnarray*} $$

(17)

To apply to nuclear systems, replacing the electron mass m_e with the nucleon mass M, one obtains

$$\begin{eqnarray*} E_{{{{{\rm{Retar}}}}}}^{{{{{\rm{x}}}}}} \left[ \rho _p, \rho _n \right] = \frac{5 \pi \hbar ^2 e^2}{4 M^2 c^2} \int \left[ \rho _p \left( {\boldsymbol {r}} \right) \right]^2 \, d {\boldsymbol {r}}. \end{eqnarray*} $$

(18)

Combining Eqs. (13) and (18), the total EDF for the correction reads

$$\begin{eqnarray*} E_{{{{{\rm{Rel}}}}}}^{{{{{\rm{tot}}}}}} \left[ \rho _p, \rho _n \right] &=& \int \left\lbrace \frac{C_{pp}}{2} \left[ \rho _p \left( {\boldsymbol {r}} \right) \right]^2 + C_{pn} \rho _p \left( {\boldsymbol {r}} \right) \rho _n \left( {\boldsymbol {r}} \right) \right\rbrace \, d {\boldsymbol {r}} - \int \frac{C_{pp}}{4} \left[ \rho _p \left( {\boldsymbol {r}} \right) \right]^2 \, d {\boldsymbol {r}}\nonumber \\ &&+\,\, \frac{5 \pi \hbar ^2 e^2}{4 M^2 c^2} \int \left[ \rho _p \left( {\boldsymbol {r}} \right) \right]^2 \, d {\boldsymbol {r}} \nonumber \\ &=& \left( \frac{C_{pp}}{2} - \frac{C_{pp}}{4} + \frac{5 \pi \hbar ^2 e^2}{4 M^2 c^2} \right) \int \left[ \rho _p \left( {\boldsymbol {r}} \right) \right]^2 \, d {\boldsymbol {r}} + C_{pn} \int \rho _p \left( {\boldsymbol {r}} \right) \rho _n \left( {\boldsymbol {r}} \right) \, d {\boldsymbol {r}} \nonumber \\ &=& \left( 3 - \tilde{\mu }_p \right) \frac{\pi \hbar ^2 e^2}{2 M^2 c^2} \int \left[ \rho _p \left( {\boldsymbol {r}} \right) \right]^2 \, d {\boldsymbol {r}} - \tilde{\mu }_n \frac{\pi \hbar ^2 e^2}{M^2 c^2} \int \rho _p \left( {\boldsymbol {r}} \right) \rho _n \left( {\boldsymbol {r}} \right) \, d {\boldsymbol {r}} \\ &\approx& \frac{\pi \hbar ^2 e^2}{M^2 c^2} \int \left[ 0.1035 \rho _p \left( {\boldsymbol {r}} \right) + 1.913 \rho _n \left( {\boldsymbol {r}} \right) \right] \rho _p \left( {\boldsymbol {r}} \right) \, d {\boldsymbol {r}}. \end{eqnarray*} $$

(19),(20)

3. Analytical estimation

Before performing numerical many-body calculations with the relativistic correction, I estimate the contribution of the correction to the total energy with a simple model and discuss the systematic behavior. As in Ref. [40], I assume hard-sphere distributions for both proton and neutron distributions ρ_p and ρ_n:

$$\begin{eqnarray*} \rho _{\tau } \left( r \right) = \left\lbrace \begin{array}{@{}l@{\quad }l@{}}\rho _0^{\tau } & r \le R_{\tau }, \\ 0 & r \gt R_{\tau }, \end{array}\right. \end{eqnarray*} $$

(21)

where R_τ is the proton (τ = p) or neutron (τ = n) radius defined by

$$\begin{eqnarray*} R_{\tau } = \left( \frac{3 N_{\tau }}{4 \pi \rho _0^{\tau }} \right)^{1/3} \end{eqnarray*} $$

(22)

with N_p = Z and N_n = N. For simplicity, half of the saturation density of atomic nuclei |$\rho _0^p = \rho _0^n \approx 0.08 \, \mathrm{fm}^{-3}$| is assumed. With this density, one can estimate the Coulomb Hartree and exchange energies as [40]

$$\begin{eqnarray*} E_{{{{{\rm{Coul}}}}}}^{{{{{\rm{H}}}}}} &=& \frac{3e^2}{5} \left( \frac{4 \pi \rho _0^p}{3} \right)^{1/3} Z^{5/3} \nonumber \\ &\simeq& 0.60 Z^{5/3} \, \mathrm{MeV}, \end{eqnarray*} $$

(23a)

$$\begin{eqnarray*} E_{{{{{\rm{Coul}}}}}}^{{{{{\rm{x}}}}}} &=& - \frac{3e^2}{4} \left( \frac{3 \rho _0^p}{\pi } \right)^{1/3} Z \nonumber \\ &\simeq& - 0.46 Z \, \mathrm{MeV}, \end{eqnarray*} $$

(23b)

respectively.

Since the relations

$$\begin{eqnarray*} \int \left[ \rho _p \left( {\boldsymbol {r}} \right) \right]^2 \, d {\boldsymbol {r}} &=& \frac{4 \pi R_p^3}{3} \left( \rho _0^p \right)^2 \nonumber \\ &=& \rho _0^p Z, \end{eqnarray*} $$

(24a)

$$\begin{eqnarray*} \int \rho _p \left( {\boldsymbol {r}} \right) \rho _n \left( {\boldsymbol {r}} \right) \, d {\boldsymbol {r}} &=& \frac{4 \pi R_{\lt }^3}{3} \rho _p^0 \rho _n^0 \nonumber \\ &=& \left\lbrace \begin{array}{@{}l@{\quad }l@{}}\rho _n^0 Z & \text{for $ R_p \lt R_n $}, \\ \rho _p^0 N & \text{for $ R_p \gt R_n $}, \end{array}\right. \end{eqnarray*} $$

(24b)

hold with R_< = R_p for R_p < R_n and R_< = R_n for R_p > R_n, one, accordingly, obtains

$$\begin{eqnarray*} E_{{{{{\rm{Darwin}}}}}}^{{{{{\rm{H}}}}}} &=& \frac{C_{pp}}{2} \frac{4 \pi R_p^3}{3} \left( \rho _0^p \right)^2 + C_{pn} \frac{4 \pi R_{\lt }^3}{3} \rho _p^0 \rho _n^0 \nonumber \\ &=& - \frac{2 \pi ^2 \hbar ^2 e^2}{3 M^2 c^2} \left[ \left( 2 \tilde{\mu }_p - 1 \right) R_p^3 \left( \rho _0^p \right)^2 + 2 \tilde{\mu }_n R_{\lt }^3 \rho _p^0 \rho _n^0 \right], \end{eqnarray*} $$

(25a)

$$\begin{eqnarray*} E_{{{{{\rm{Darwin}}}}}}^{{{{{\rm{x}}}}}} &=& - \frac{C_{pp}}{4} \frac{4 \pi R_p^3}{3} \left( \rho _0^p \right)^2 \nonumber \\ &=& \frac{\pi ^2 \hbar ^2 e^2}{3 M^2 c^2} \left( 2 \tilde{\mu }_p - 1 \right) R_p^3 \left( \rho _0^p \right)^2 \nonumber \\ &=& \frac{\pi \hbar ^2 e^2}{4 M^2 c^2} \left( 2 \tilde{\mu }_p - 1 \right) \rho _0^p Z \nonumber \\ &\simeq& 0.018 Z \, \mathrm{MeV}, \end{eqnarray*} $$

(25b)

$$\begin{eqnarray*} E_{{{{{\rm{Retar}}}}}}^{{{{{\rm{x}}}}}} &=& \frac{5 \pi ^2 \hbar ^2 e^2}{3 M^2 c^2} R_p^3 \left( \rho _0^p \right)^2 \nonumber \\ &=& \frac{5 \pi \hbar ^2 e^2}{4 M^2 c^2} \rho _0^p Z \nonumber \\ &\simeq& 0.020 Z \, \mathrm{MeV}, \end{eqnarray*} $$

(25c)

$$\begin{eqnarray*} E_{{{{{\rm{Rel}}}}}}^{{{{{\rm{tot}}}}}} &=& \frac{2 \pi ^2 \hbar ^2 e^2}{3 M^2 c^2} \left[ \left( 3 - \tilde{\mu }_p \right) R_p^3 \left( \rho _0^p \right)^2 - 2 \tilde{\mu }_n R_{\lt }^3 \rho _p^0 \rho _n^0 \right]. \end{eqnarray*} $$

(25d)

In the case of R_p < R_n, Eqs. (25a) and (25d) read

$$\begin{eqnarray*} E_{{{{{\rm{Darwin}}}}}}^{{{{{\rm{H}}}}}} &=& - \frac{2 \pi ^2 \hbar ^2 e^2}{3 M^2 c^2} \left[ \left( 2 \tilde{\mu }_p - 1 \right) \left( \rho _0^p \right)^2 + 2 \tilde{\mu }_n \rho _p^0 \rho _n^0 \right] R_p^3 \nonumber \\ &=& - \frac{\pi \hbar ^2 e^2}{2 M^2 c^2} \left[ \left( 2 \tilde{\mu }_p - 1 \right) \rho _0^p + 2 \tilde{\mu }_n \rho _n^0 \right] Z \nonumber \\ &\simeq& - 0.006 Z \, \mathrm{MeV}, \end{eqnarray*} $$

(26a)

$$\begin{eqnarray*} E_{{{{{\rm{Rel}}}}}}^{{{{{\rm{tot}}}}}} &=& \frac{2 \pi ^2 \hbar ^2 e^2}{3 M^2 c^2} \left[ \left( 3 - \tilde{\mu }_p \right) R_p^3 \left( \rho _0^p \right)^2 - 2 \tilde{\mu }_n R_p^3 \rho _p^0 \rho _n^0 \right] \nonumber \\ &=& \frac{\pi \hbar ^2 e^2}{2 M^2 c^2} \left[ \left( 3 - \tilde{\mu }_p \right) \rho _0^p - 2 \tilde{\mu }_n \rho _n^0 \right] Z \nonumber \\ &\simeq& 0.032 Z \, \mathrm{MeV} ; \end{eqnarray*} $$

(26b)

in the case of R_p > R_n, they read

$$\begin{eqnarray*} E_{{{{{\rm{Darwin}}}}}}^{{{{{\rm{H}}}}}} &=& - \frac{2 \pi ^2 \hbar ^2 e^2}{3 M^2 c^2} \left[ \left( 2 \tilde{\mu }_p - 1 \right) R_p^3 \left( \rho _0^p \right)^2 + 2 \tilde{\mu }_n R_n^3 \rho _p^0 \rho _n^0 \right] \nonumber \\ &=& - \frac{\pi \hbar ^2 e^2}{2 M^2 c^2} \rho _0^p \left[ \left( 2 \tilde{\mu }_p - 1 \right) Z + 2 \tilde{\mu }_n N \right] \nonumber \\ &\simeq& -0.037 Z + 0.031 N \, \mathrm{MeV}, \end{eqnarray*} $$

(27a)

$$\begin{eqnarray*} E_{{{{{\rm{Rel}}}}}}^{{{{{\rm{tot}}}}}} &=& \frac{2 \pi ^2 \hbar ^2 e^2}{3 M^2 c^2} \left[ \left( 3 - \tilde{\mu }_p \right) R_p^3 \left( \rho _0^p \right)^2 - 2 \tilde{\mu }_n R_n^3 \rho _p^0 \rho _n^0 \right] \nonumber \\ &=& \frac{\pi \hbar ^2 e^2}{2 M^2 c^2} \rho _0^p \left[ \left( 3 - \tilde{\mu }_p \right) Z - 2 \tilde{\mu }_n N \right] \nonumber \\ &\simeq& 0.002 Z + 0.031 N \, \mathrm{MeV}. \end{eqnarray*} $$

(27b)

Comparing Eqs. (23b) and (26b), one can find that the relativistic correction to the total energy |$\left| E_{{{{{\rm{Rel}}}}}}^{{{{{\rm{tot}}}}}} \right|$| is 15 times smaller than the Coulomb exchange energy |$\left| E_{{{{{\rm{Coul}}}}}}^{{{{{\rm{x}}}}}} \right|$|⁠. This means that although |$E_{{{{{\rm{Rel}}}}}}^{{{{{\rm{tot}}}}}}$| is small compared to the other contributions, one can estimate that |$E_{{{{{\rm{Rel}}}}}}^{{{{{\rm{tot}}}}}}$| for |${}^{208}_{}\mathrm{Pb}$| may be around |$2 \, \mathrm{MeV}$|⁠, which is a non-negligible contribution compared to the desired accuracy of |$E_{{{{{\rm{nucl}}}}}}$| (⁠|$O \left( 100 \right) \, \mathrm{keV}$|⁠).

In atomic physics, the vacuum polarization to the Coulomb potential formed by an atomic nucleus is the next-leading order to the Breit correction, i.e., the relativistic correction of electron–electron interaction [58]. On the other hand, in nuclear structure, both the vacuum polarization and the relativistic correction are corrections of the proton–proton Coulomb interaction. Thus, it is worthwhile to compare them. In Ref. [40], the vacuum polarization was taken into account in a nuclear structure calculation using the Uehling potential [79,80]. The contribution of the vacuum polarization was estimated as [40]

$$\begin{eqnarray*} E_{{{{{\rm{VP}}}}}} \simeq 0.0047 Z^{5/3} \, \mathrm{MeV}. \end{eqnarray*} $$

(28)

Hence, the relativistic correction and the vacuum polarization provide comparable contributions to the total energy in the light nuclei (e.g., for Z = 20, |$E_{{{{{\rm{Rel}}}}}}^{{{{{\rm{tot}}}}}} \approx 0.64 \, \mathrm{MeV}$| and |$E_{{{{{\rm{VP}}}}}} \approx 0.69 \, \mathrm{MeV}$|⁠) but the latter dominates in heavy nuclei (e.g., for Z = 100, |$E_{{{{{\rm{Rel}}}}}}^{{{{{\rm{tot}}}}}} \approx 3.2 \, \mathrm{MeV}$| and |$E_{{{{{\rm{VP}}}}}} \approx 10 \, \mathrm{MeV}$|⁠). In the next section, I shall confirm this by performing numerical calculations.

4. Skyrme Hartree–Fock calculation

To calculate the ground-state density distribution and energy numerically, one needs to use a many-body calculation technique. In this paper, the Skyrme Hartree–Fock method is used.

The EDF for the correction (19) is implemented with the self-consistent Skyrme Hartree–Fock calculation code skyrme_rpa [81]. In this work, doubly-magic nuclei are focused on; thus, the pairing interaction is not considered, and spherical symmetry is assumed. The calculations are performed with a box of |$15 \, \mathrm{fm}$| with |$0.1 \, \mathrm{fm}$| mesh. The SAMi EDF [82] is adopted for the nuclear EDF.

It is worthwhile to mention here that the non-relativistic Coulomb GGA exchange EDF, more precisely, the modified Coulomb GGA-Perdew–Burke–Ernzerhof (PBE) EDF with λ = 1.25 [39], reproduces the exact Coulomb Fock energy with a less than |$100 \, \mathrm{keV}$| error, which is much less than the relativistic correction, which I will discuss; thus, in this paper, the exact Coulomb Fock energy is approximated to the Coulomb GGA exchange energy, |$E_{{{{{\rm{Coul}}}}}}^{{{{{\rm{x-GGA}}}}}} \simeq E_{{{{{\rm{Coul}}}}}}^{{{{{\rm{x-exact}}}}}}$|⁠.

First of all, I shall mention the effects on nuclear radii. Neither the proton nor the neutron radii are found to be changed by more than |$0.002 \, \mathrm{fm}$| due to the relativistic correction and the vacuum polarization. Therefore, the effect of the correction on nuclear radii can be safely neglected.

Next, the ground-state energies for doubly-magic nuclei are shown in Table 2. The columns labeled “Coulomb (NR-LDA)” and “Coulomb (R-LDA)” show results with the non-relativistic Coulomb LDA (11) and the relativistic Coulomb LDA (11), (19), respectively. For comparison, total energies calculated with the vacuum polarization above the relativistic Coulomb LDA and with the non-relativistic Coulomb GGA [39] are shown in a column labeled “C. (R-LDA) & V.P.” and “Coulomb (NR-GGA)”, respectively. Both the relativistic correction and the vacuum polarization make the atomic nuclei less strongly bound. The contribution of the former is around |$0.42 \, \mathrm{MeV}$| in |${}^{40}_{}\mathrm{Ca}$| and |$2.37 \, \mathrm{MeV}$| in |${}^{208}_{}\mathrm{Pb}$|⁠. Although the relativistic correction changes |$\Delta E_{{{{{\rm{tot}}}}}}$| of |${}^{48}_{}\mathrm{Ca}$| and |${}^{48}_{}\mathrm{Ni}$| by |$0.5 \, \mathrm{MeV}$|⁠, the mass difference is not changed.

Table 2.

Ground-state energies calculated with the SAMi EDF [82] shown in MeV. The columns labeled “Coulomb (NR-LDA)” and “Coulomb (R-LDA)” show results calculated with the non-relativistic Coulomb LDA (11) and the relativistic Coulomb LDA (11), (19), respectively. For comparison, total energies calculated with the vacuum polarization above the relativistic Coulomb LDA and with the non-relativistic Coulomb GGA are shown in columns labeled “C. (R-LDA) & V.P.” and “Coulomb (NR-GGA)”, respectively.

Nuclei	Coulomb (NR-LDA)	Coulomb (R-LDA)	C. (R-LDA) & V.P.	Coulomb (NR-GGA)
|${}^{4}_{}\mathrm{He}$|	−27.5263	−27.5001	−27.4918	−27.6011
|${}^{16}_{}\mathrm{O}$|	−130.4800	−130.3340	−130.2444	−130.6630
|${}^{40}_{}\mathrm{Ca}$|	−347.0848	−346.6639	−346.2544	−347.4045
|${}^{48}_{}\mathrm{Ca}$|	−415.6148	−415.1104	−414.7093	−415.9293
|${}^{48}_{}\mathrm{Ni}$|	−352.6388	−352.1357	−351.4160	−353.0443
|${}^{100}_{}\mathrm{Sn}$|	−811.6641	−810.4970	−808.6876	−812.2390
|${}^{132}_{}\mathrm{Sn}$|	−1103.0881	−1101.6342	−1099.9484	−1103.6397
|${}^{208}_{}\mathrm{Pb}$|	−1636.6149	−1634.2456	−1630.5350	−1637.3582
|${}^{310}_{}\mathrm{126}$|	−2131.4146	−2127.8446	−2120.5577	−2132.3711

It was shown in Refs. [39,62] that the exact treatment, or its simplified treatment called GGA, of the Coulomb exchange energy makes |${}^{40}_{}\mathrm{Ca}$| more strongly bound by around |$0.32 \, \mathrm{MeV}$| and |${}^{208}_{}\mathrm{Pb}$| by around |$0.74\, \mathrm{MeV}$|⁠, which can be understood as the correction of the density gradient on the LDA. Hence, the relativistic correction of the Coulomb interaction is almost comparable with or even larger than the gradient correction, but in an opposite direction to the total energies, as in the case of atoms [51]. Thus, if both contributions are considered simultaneously, it is not even obvious whether the total energy decreases or increases. Indeed, this fact was pointed out in the context of atomic physics in Ref. [51], while no first-principles relativistic Coulomb GGA EDF is available, which remains an open question.

I also compare the relativistic correction with the vacuum polarization. The vacuum polarization also makes the atomic nuclei less strongly bound, and the contribution to the total energy is approximately |$0.4 \, \mathrm{MeV}$| in |${}^{40}_{}\mathrm{Ca}$| and |$3.7 \, \mathrm{MeV}$| in |${}^{208}_{}\mathrm{Pb}$|⁠. Thus, as discussed in Sect. 3, the vacuum polarization gives a comparable contribution to the total energy to the relativistic correction in light nuclei, but it contributes more in heavy nuclei.

The systematic behavior of the relativistic correction and the vacuum polarization shall be discussed next. The contributions of the relativistic correction and the vacuum polarization can be seen as the differences of the second and third columns and the third and fourth columns in Table 2, respectively. The proton-number Z dependence of these contributions is shown in Fig. 1. These energies are fitted with E = aZ^b, where a and b are adjustable parameters. The results are

$$\begin{eqnarray*} E_{{{{{\rm{Rel}}}}}}^{{{{{\rm{tot}}}}}} &\simeq& 0.016\,60 Z^{1.12} \, \mathrm{MeV}, \end{eqnarray*} $$

(29a)

$$\begin{eqnarray*} E_{{{{{\rm{VP}}}}}} &\simeq& 0.003\,98 Z^{1.55} \, \mathrm{MeV}, \end{eqnarray*} $$

(29b)

which are consistent with Eqs. (26b) and (28).

Fig. 1.

Contribution of relativistic correction and vacuum polarization to the total energies, |$E_{{{{{\rm{Rel}}}}}}^{{{{{\rm{tot}}}}}}$| and |$E_{{{{{\rm{VP}}}}}}$|⁠, shown as functions of proton number Z by circles and squares, respectively. Empirical fits with Eq. (29) are also shown by the solid and dashed lines, respectively.

The breakdown of the relativistic correction (19) is shown in Table 3. The columns labeled |$E_{{{{{\rm{Darwin}}}}}}^{{{{{\rm{H}}}}}}$|⁠, |$E_{{{{{\rm{Darwin}}}}}}^{{{{{\rm{x}}}}}}$|⁠, and |$E_{{{{{\rm{Retar}}}}}}^{{{{{\rm{x}}}}}}$| are the Darwin Hartree (13a), Darwin exchange (13b), and retardation exchange (18) energy, respectively. As shown, in N ≤ Z nuclei, |$E_{{{{{\rm{Darwin}}}}}}^{{{{{\rm{H}}}}}} \lt 0$| holds, while, in N > Z nuclei, |$E_{{{{{\rm{Darwin}}}}}}^{{{{{\rm{H}}}}}} \gt 0$| holds. The Darwin Hartree energy reads

$$\begin{eqnarray*} E_{{{{{\rm{Darwin}}}}}}^{{{{{\rm{H}}}}}} &=& - \frac{\pi \hbar ^2 e^2}{M^2 c^2} \int \left\lbrace \frac{2 \tilde{\mu }_p - 1}{2} \rho _p \left( {\boldsymbol {r}} \right) + \tilde{\mu }_n \rho _n \left( {\boldsymbol {r}} \right) \right\rbrace \rho _p \left( {\boldsymbol {r}} \right) \, d {\boldsymbol {r}} \nonumber \\ &\simeq& - \frac{\pi \hbar ^2 e^2}{M^2 c^2} \int \left[ 2.292 \rho _p \left( {\boldsymbol {r}} \right) - 1.913 \rho _n \left( {\boldsymbol {r}} \right) \right] \rho _p \left( {\boldsymbol {r}} \right) \, d {\boldsymbol {r}}; \end{eqnarray*} $$

(30)

hence, |$\left[ 2.292 \rho _p \left( {\boldsymbol {r}} \right) - 1.913 \rho _n \left( {\boldsymbol {r}} \right) \right]$| determines the sign of |$E_{{{{{\rm{Darwin}}}}}}^{{{{{\rm{H}}}}}}$|⁠. In N ≤ Z nuclei, |$\rho _p \left( {\boldsymbol {r}} \right) \gtrsim \rho _n \left( {\boldsymbol {r}} \right)$| holds and, accordingly, it is positive (i.e., |$E_{{{{{\rm{Darwin}}}}}}^{{{{{\rm{H}}}}}} \lt 0$|⁠), whereas, in N > Z nuclei, |$\rho _p \left( {\boldsymbol {r}} \right) \lesssim \rho _n \left( {\boldsymbol {r}} \right)$| holds and, accordingly, it can be negative (i.e., |$E_{{{{{\rm{Darwin}}}}}}^{{{{{\rm{H}}}}}} \gt 0$|⁠) in contrast to Eq. (25a). The Darwin and retardation exchange energies are also fitted with E = aZ^b, where the results are

$$\begin{eqnarray*} E_{{{{{\rm{Darwin}}}}}}^{{{{{\rm{x}}}}}} &\simeq& 0.013\,12 Z^{0.97} \, \mathrm{MeV}, \end{eqnarray*} $$

(31a)

$$\begin{eqnarray*} E_{{{{{\rm{Retar}}}}}}^{{{{{\rm{x}}}}}} &\simeq& 0.014\,31 Z^{0.97} \, \mathrm{MeV}, \end{eqnarray*} $$

(31b)

which are consistent with Eqs. (25b) and (25c).

Table 3.

Breakdown of the relativistic correction (19) shown in MeV. The columns labeled |$E_{{{{{\rm{Darwin}}}}}}^{{{{{\rm{H}}}}}}$|⁠, |$E_{{{{{\rm{Darwin}}}}}}^{{{{{\rm{x}}}}}}$|⁠, and |$E_{{{{{\rm{Retar}}}}}}^{{{{{\rm{x}}}}}}$| are Darwin Hartree energy (13a), Darwin exchange energy (13b), and retardation exchange energy (18), respectively.

Nuclei	|$E_{{{{{\rm{Darwin}}}}}}^{{{{{\rm{H}}}}}}$|	|$E_{{{{{\rm{Darwin}}}}}}^{{{{{\rm{x}}}}}}$|	|$E_{{{{{\rm{Retar}}}}}}^{{{{{\rm{x}}}}}}$|
|${}^{4}_{}\mathrm{He}$|	−0.0048	+0.0148	+0.0162
|${}^{16}_{}\mathrm{O}$|	−0.0256	+0.0821	+0.0895
|${}^{40}_{}\mathrm{Ca}$|	−0.0704	+0.2349	+0.2562
|${}^{48}_{}\mathrm{Ca}$|	+0.0379	+0.2230	+0.2432
|${}^{48}_{}\mathrm{Ni}$|	−0.2429	+0.3568	+0.3890
|${}^{100}_{}\mathrm{Sn}$|	−0.1869	+0.6475	+0.7060
|${}^{132}_{}\mathrm{Sn}$|	+0.3026	+0.5505	+0.6002
|${}^{208}_{}\mathrm{Pb}$|	+0.4221	+0.9311	+1.0153
|${}^{310}_{}\mathrm{126}$|	+0.5492	+1.4445	+1.5750

The exchange contribution to the EM energy shall be further discussed. The ratio of the Coulomb or relativistic exchange energy to the non-relativistic Coulomb GGA exchange energy

$$\begin{eqnarray*} \Delta E_{{{{{\rm{x}}}}}} = \frac{E^{{{{{\rm{x-LDA}}}}}} - E_{{{{{\rm{Coul}}}}}}^{{{{{\rm{x-GGA}}}}}}}{E_{{{{{\rm{Coul}}}}}}^{{{{{\rm{x-GGA}}}}}}} \simeq \frac{E^{{{{{\rm{x-LDA}}}}}} - E_{{{{{\rm{Coul}}}}}}^{{{{{\rm{x-exact}}}}}}}{E_{{{{{\rm{Coul}}}}}}^{{{{{\rm{x-exact}}}}}}} \end{eqnarray*} $$

(32)

is shown in Fig. 2. The non-relativistic and relativistic Coulomb LDA exchange energies, |$E^{{{{{\rm{x-LDA}}}}}}$|⁠, are calculated by using |$E_{{{{{\rm{Coul}}}}}}^{{{{{\rm{x}}}}}}$| (11c) and |$E_{{{{{\rm{Coul}}}}}}^{{{{{\rm{x}}}}}} + E_{{{{{\rm{Darwin}}}}}}^{{{{{\rm{x}}}}}} + E_{{{{{\rm{Retar}}}}}}^{{{{{\rm{x}}}}}}$| (13b), (18), respectively. It is seen that the absolute value of the relative change |$\left| \Delta E_{{{{{\rm{x}}}}}} \right|$| for the difference between the exact and LDA non-relativistic Coulomb exchange energy gets smaller as A increases in light nuclei and reaches almost constant |$\Delta E_{{{{{\rm{x}}}}}} \approx -3 \, \%$| in heavy nuclei (A ≳ 100). The absolute value of the relative change |$\left| \Delta E_{{{{{\rm{x}}}}}} \right|$| for the difference between the non-relativistic and relativistic LDA Coulomb exchange energy behaves similarly, while its value is larger and reaches almost constant |$\Delta E_{{{{{\rm{x}}}}}} \approx -9 \, \%$| in heavy nuclei. The Darwin contribution and the retardation one are almost the same size.

Fig. 2.

Ratio of the non-relativistic (square) or the relativistic (circle) Coulomb LDA exchange energies to the modified Coulomb GGA-PBE exchange energy with λ = 1.25, which reproduces the exact Coulomb Fock energy within sufficient accuracy, |$\Delta E_{{{{{\rm{x}}}}}}$| (32), shown as a function of the mass number A. The Darwin contribution only (up-triangle) and the retardation contribution only (down-triangle) are also shown.

This difference has been studied in the relativistic (covariant) DFT [63,64]. The exact Coulomb exchange energy in this paper corresponds to Eq. (4) in Ref. [64], which does not seem to include the relativistic correction. They also show that the difference between the non-relativistic and relativistic Coulomb LDA exchange energies is approximately |$1 \, \mathrm{MeV}$|⁠. Considering such difference, |$E_{{{{{\rm{Coul}}}}}}^{{{{{\rm{x}}}}}} + E_{{{{{\rm{Rel}}}}}}^{{{{{\rm{x}}}}}}$| calculated in this work is almost consistent qualitatively with, but almost twice that of, Figs. 1 and 4 of Ref. [64], because Ref. [64] does not include |$E_{{{{{\rm{Darwin}}}}}}^{{{{{\rm{x}}}}}}$|⁠. If one considers |$E_{{{{{\rm{Darwin}}}}}}^{{{{{\rm{x}}}}}}$| only, this work is quantitatively consistent with Figs. 1 and 4 of Ref. [64].

Lastly, a change of the single-particle energies is discussed, for which |${}^{208}_{}\mathrm{Pb}$| is taken as an example. The single-particle energies calculated with the non-relativistic and relativistic Coulomb LDA are shown in Table 4. For comparison, those calculated with the vacuum polarization above the relativistic LDA are shown. It is seen that the proton single-particle energies change by about |$0.04 \, \mathrm{MeV}$| in inner shells and |$0.01 \, \mathrm{MeV}$| in outer shells, due to the relativistic correction. In contrast, the change due to the vacuum polarization does not depend strongly on orbitals, and its value is about |$0.1 \, \mathrm{MeV}$|⁠. This difference is due to the nature of the interaction; the relativistic correction is a point-coupling interaction, while the vacuum polarization is a finite range. Since the Darwin Hartree term contributes to the neutron mean-field potential [see Eq. (20)], the neutron single-particle energies change by about |$0.04 \, \mathrm{MeV}$| in inner shells and |$0.01 \, \mathrm{MeV}$| in outer shells, due to the relativistic correction. In contrast, the vacuum polarization changes neutron single-particle energies by less than |$0.01 \, \mathrm{MeV}$| only.

Table 4.

Single-particle energies of |${}^{208}_{}\mathrm{Pb}$| shown in MeV calculated with the SAMi EDF.

States	Coulomb		C. (R-LDA) & V.P.
	NR-LDA	R-LDA
π1s1/2	−44.980	−44.936	−44.829
π1p3/2	−39.387	−39.346	−39.246
π1p1/2	−39.107	−39.067	−38.967
π1d5/2	−32.482	−32.446	−32.352
π1d3/2	−31.815	−31.779	−31.686
π2s1/2	−28.509	−28.475	−28.379
π1f7/2	−24.692	−24.661	−24.572
π1f5/2	−23.353	−23.323	−23.235
π2p3/2	−19.411	−19.383	−19.292
π2p1/2	−18.626	−18.598	−18.508
π1g9/2	−16.338	−16.312	−16.229
π1g7/2	−14.019	−13.995	−13.913
π2d5/2	−10.255	−10.233	−10.148
π2d3/2	−8.846	−8.826	−8.742
π3s1/2	−7.673	−7.653	−7.568
π1h11/2	−7.663	−7.642	−7.563
ν1s1/2	−59.291	−59.253	−59.238
ν1p3/2	−52.953	−52.919	−52.908
ν1p1/2	−52.656	−52.622	−52.612
ν1d5/2	−45.375	−45.345	−45.338
ν1d3/2	−44.744	−44.715	−44.709
ν2s1/2	−41.962	−41.935	−41.928
ν1f7/2	−36.904	−36.880	−36.876
ν1f5/2	−35.702	−35.679	−35.677
ν2p3/2	−32.094	−32.072	−32.069
ν2p1/2	−31.344	−31.323	−31.321
ν1g9/2	−27.863	−27.843	−27.843
ν1g7/2	−25.791	−25.773	−25.775
ν2d5/2	−22.277	−22.261	−22.262
ν2d3/2	−20.843	−20.829	−20.830
ν3s1/2	−19.959	−19.945	−19.946
ν1h11/2	−18.533	−18.519	−18.521
ν1h9/2	−15.306	−15.295	−15.300
ν2f7/2	−12.672	−12.661	−12.664
ν2f5/2	−10.622	−10.614	−10.619
ν3p3/2	−9.859	−9.851	−9.854
ν3p1/2	−9.063	−9.055	−9.059
ν1i13/2	−9.150	−9.141	−9.145

5. Conclusion and perspectives

In this paper, the relativistic correction of the Coulomb interaction was introduced in Skyrme Hartree–Fock calculations using the local density approximation (LDA). The correction contains the Darwin term, the retardation, the spin-orbit interaction, and the spin-magnetic interaction, but only the first and second terms were considered in this paper because the remaining terms vanish in the LDA.

It was found that the relativistic correction makes |${}^{40}_{}\mathrm{Ca}$| less strongly bound by |$0.4 \, \mathrm{MeV}$| and |${}^{208}_{}\mathrm{Pb}$| by |$2.4 \, \mathrm{MeV}$|⁠, which are comparable with or even larger than the difference between the exact Coulomb exchange (exact-Fock), or its simplified treatment called the generalized gradient approximation (GGA), and the LDA (Hartree–Fock–Slater approximation) Coulomb exchange energies [38,39,61,62], but with the opposite sign. Thus, once both the exact treatment and the relativistic correction are considered at the same time, even the sign of the total correction is not obvious, as has also been pointed out in the context of atomic physics [51]. In addition, since the desired accuracy of the nuclear EDF is of the order of |$100 \, \mathrm{keV}$|⁠, the relativistic correction is often non-negligible, as well as the correction due to the density gradient.

The correction was also compared with the vacuum polarization. Since the proton-number dependence of the vacuum polarization is stronger than that of the relativistic correction, these two contributions to the total energy are comparable in light nuclei, but the former dominates in heavy nuclei. This behavior is in contrast to atoms, in which the vacuum polarization is subdominant to the Breit correction [58]. Although the contribution of vacuum polarization to the total energy is usually tiny, it is non-negligible in discussions of mirror nuclei mass differences or the isobaric analog energy [7,11,40], whereas the relativistic correction does not contribute to such properties.

In this calculation, the spin-orbit and spin-magnetic interactions were not able to be considered since they vanish in the LDA. The former has been discussed in Refs. [7,40,78], while the latter has not been considered in the self-consistent calculation. Consideration of the latter is left for a future task. The relativistic Coulomb GGA functional is indispensable to achieve more accurate calculations in both atomic and nuclear structure calculations, which is also left for a future perspective.

The Coulomb correlation energy, which originates from many-body correlations [83–87], is also an important and interesting topic. In electronic systems, the Coulomb correlation energy can be defined without any ambiguity as the difference between the “exact” energy and the Hartree–Fock energy, which can be calculated by using the Monte Carlo calculation and some other many-body techniques [15–17,88–92]. Reference [60] derived the relativistic Coulomb exchange-correlation EDF in the LDA level, while Ref. [51] discussed that the relativistic correction of the Coulomb correlation energy is almost negligible. Meanwhile, in nuclear systems, there exist both the nuclear and Coulomb interactions; hence, there is still ambiguity in defining the Coulomb correlation energy, which has been discussed in several works [38,93–95]. Thus, it should be discussed whether the relativistic correction of the Coulomb correlation energy in nuclear systems is negligible; this is also left for future work.

Once all the relevant corrections of the Coulomb energy are taken into account for the nuclear DFT calculation above the non-relativistic Coulomb LDA exchange EDF, the correction may be more than |$1 \, \mathrm{MeV}$| for heavy nuclei, such as |${}^{208}_{}\mathrm{Pb}$|⁠. For instance, the nucleon finite-size effects and the vacuum polarization contribute to the binding energy of |${}^{208}_{}\mathrm{Pb}$| by about |$4 \, \mathrm{MeV}$| in total [40], the relativistic correction discussed in this paper is about |$2 \, \mathrm{MeV}$|⁠, and the Coulomb correlation energy may reach an order of |$1 \, \mathrm{MeV}$| [94]. The current accuracy of nuclear EDFs is about 1–|$5 \, \mathrm{MeV}$| in binding energies [31,96], which is the same order as or even smaller than the Coulomb corrections. Therefore, in order to achieve accurate calculation of binding energies using a nuclear EDF quantitatively, the nuclear EDF should, in principle, be refitted by considering all the relevant Coulomb corrections.

ACKNOWLEDGEMENTS

The author thanks Takumi Doi, Kenji Fukushima, Kouichi Hagino, Haozhao Liang, and Kenichi Yoshida for stimulating discussions. The author also thanks Kouichi Hagino for variable comments on the manuscript. The author acknowledges a JSPS Grant-in-Aid for JSPS Fellows under Grant No. 19J20543, a JSPS Grant-in-Aid for Research Activity Start-up under Grant No. 22K20372, the RIKEN iTHEMS program, the RIKEN Special Postdoctoral Researchers Program, and the Science and Technology Hub Collaborative Research Program from the RIKEN Cluster for Science, Technology and Innovation Hub (RCSTI). The numerical calculations were performed on cluster computers at the RIKEN iTHEMS program.

Footnotes

1

The effects of nucleon finite size are discussed in Refs. [11,40].