11institutetext:  Department of Physics, Graduate School of Science, The
University of Tokyo, Tokyo 113-0033, Japan 22institutetext: RIKEN Nishina
Center, Wako 351-0198, Japan 33institutetext: 33email<EMAIL_ADDRESS>
# Relativistic correction of the Coulomb interaction for energies and radii in
doubly-magic nuclei

Tomoya Naito Present Address: RIKEN Interdisciplinary Theoretical and
Mathematical Sciences (iTHEMS), Wako 351-0198, Japan11 2 2 3 3

(Received: )

###### Abstract

Effects of the relativistic correction of the Coulomb interaction on doubly-
magic nuclei are discussed with Skyrme Hartree-Fock calculations. It is found
that the correction to the total energy is about $2.4\,\mathrm{MeV}$ for
${}^{208}\mathrm{Pb}$, while 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. Moreover, since the covariant (relativistic) density functional
theory often includes effects of magnetic interaction and retardation for
electromagnetic interaction, the method presented in this work enables one to
make a fair comparison between the non-relativistic and the relativistic
density functional calculations for electromagnetic energies.

###### pacs:

21.60.Jz Nuclear Density Functional Theory and extensions and 21.10.Sf Coulomb
energies, analogue states

##  1 Introduction

Atomic nuclei are composed of protons and neutrons, which interact via the
nuclear and electromagnetic (EM) interactions. Since the nuclear interaction
is much stronger than the EM one, the former determines most properties of the
atomic nuclei. Nevertheless, the EM interaction plays an important role when
one focuses on nuclear properties related to the isospin symmetry breaking
Okamoto1964Phys.Lett.11_150 ; Friedman1971Phys.Lett.B35_543 ;
Shlomo1982Phys.Scr.26_280 ; Suzuki1993Phys.Rev.C47_R1360 ;
Kaneko2012Phys.Rev.Lett.109_092504 ; Baczyk2018Phys.Lett.B778_178 ; Roca-
Maza2018Phys.Rev.Lett.120_202501 ; Baczyk2019J.Phys.G46_03LT01 ;
Loc2019Phys.Rev.C99_014311 ; Dong2019Phys.Rev.C99_014319 . For more details,
see reviews, e.g., Ref. Auerbach1983Phys.Rep.98_273 .

The Hartree-Fock method or the density functional theory (DFT)
Hohenberg1964Phys.Rev.136_B864 ; Kohn1965Phys.Rev.140_A1133 ;
Kohn1999Rev.Mod.Phys.71_1253  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 in the local density approximation (LDA)
Vosko1980Can.J.Phys.58_1200 ; Perdew1981Phys.Rev.B23_5048 ;
Yokota2021Phys.Rev.Research3_L012015 , and effects of the density gradient
have been included empirically or phenomenologically
Perdew1986Phys.Rev.B33_8800 ; Perdew1992Phys.Rev.B45_13244 ;
Perdew1996Phys.Rev.Lett.77_3865 . 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 EDF is one of the important current
issues Drut2010Prog.Part.Nucl.Phys.64_120 ; Dobaczewski2016J.Phys.G43_04LT01 ;
Yokota2019Phys.Rev.C99_024302 ; Salvioni2020J.Phys.G47_085107 ;
Accorto2021Phys.Rev.C103_044304 ; Marino2021Phys.Rev.C104_024315 ;
Naito:2021ozn . Therefore, most nuclear EDFs have been constructed
phenomenologically with fitting to the experimental data and different EDFs
may give different results Li2008Phys.Rep.464_113 ;
Danielewicz2009Nucl.Phys.A818_36 ; Danielewicz2014Nucl.Phys.A922_1 . The
accuracy of total energies obtained by widely used EDFs ranges approximately
$1$–$5\,\mathrm{MeV}$ Stoitsov2003Phys.Rev.C68_054312 ;
Kortelainen2010Phys.Rev.C82_024313 , while those of empirical mass formulae
and machine learning ranges $100$–$500\,\mathrm{keV}$
Wang2014Phys.Lett.B734_215 ; Niu2018Phys.Lett.B778_48 ; Niu2018Sci.Bull.63_759
; Yamaguchi2021Prog.Part.Nucl.Phys.120_103882 , which almost matches the
desired accuracy for the study of $r$-process nucleosynthesis
Lunney2003Rev.Mod.Phys.75_1021 .

To fit the parameters of nuclear EDF properly, the EM contribution to the
nuclear properties should be subtracted from experimental data properly.
Motivated by this issue, in a series of works Naito2018Phys.Rev.C97_044319 ;
Naito2019Phys.Rev.C99_024309 ; Naito2020Phys.Rev.C101_064311 , accurate
treatment of the Coulomb interaction in a Skyrme Hartree-Fock (SHF)
calculation Vautherin1972Phys.Rev.C5_626 , one of the branches of the nuclear
DFT Meng2006Prog.Part.Nucl.Phys.57_470 ; Niksic2011Prog.Part.Nucl.Phys.66_519
; Liang2015Phys.Rep.570_1 ; Nakatsukasa2016Rev.Mod.Phys.88_045004 ;
Colo2020Adv.Phys.X5_1740061 , 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
Naito2018Phys.Rev.C97_044319 ; Naito2019Phys.Rev.C99_024309 . On top of it,
more higher-order contributions are taken into account; the spatial charge
distributions of nucleons and the vacuum polarization between two protons
Auerbach1983Phys.Rep.98_273 ; Samaddar1986Nucl.Phys.A451_160 ;
Chamel2009Phys.Rev.C80_065804 ; Roca-Maza2018Phys.Rev.Lett.120_202501 ;
Naito2020Phys.Rev.C101_064311 . From the point of view of the quantum field
theory, the finite-size effects, i.e., effects originating from the spatial
charge distributions of nucleons, corresponds 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 $\alpha\approx 1/137$
(the fine-structure constant). Although these effects are just corrections for
the Coulomb interaction, which itself is already subdominant for the nuclear
total energies, it was 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}$ in $-8.2\,\mathrm{MeV}$, $+1.2\,\mathrm{MeV}$, and
$+3.7\,\mathrm{MeV}$, respectively Naito2020Phys.Rev.C101_064311 .

Then, it is natural to consider the corrections with the same order with
respect to $\alpha$; that is, 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 like the vacuum polarization,
while 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 inverse proportional to the nucleon mass Wiringa1995Phys.Rev.C51_38 . 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
Saue2011ChemPhysChem12_3077 ; Naito2020J.Phys.B53_215002 . Thus, in this work,
the finite-light-speed correction to 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 to the Coulomb interaction is often referred
to as the relativistic correction Saue2011ChemPhysChem12_3077 ;
Pershina2015Nucl.Phys.A944_578 ; Pershina2015Nucl.Phys.A944_578 ;
Eliav2015Nucl.Phys.A944_518 ; Schwerdtfeger2015Nucl.Phys.A944_551 . 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 Breit1929Phys.Rev.34_553 ;
Breit1930Phys.Rev.36_383 ; Breit1932Phys.Rev.39_616 . In atomic physics, the
Breit correction dominates more than the vacuum polarization of the Coulomb
field formed by the atomic nucleus Eides2001Phys.Rep.342_63 , and thus the
Breit correction should be considered once the vacuum polarization is
considered. In nuclear physics, the vacuum polarization discussed in Ref.
Naito2020Phys.Rev.C101_064311  originates from the interaction between two
protons, which is different from that in atomic physics. Hence, it is
interesting to compare which contribution dominates to 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 had been taken
into account for an exchange MacDonald1979J.Phys.C12_2977  or exchange-
correlation Kenny1996Phys.Rev.Lett.77_1099  EDF in the LDA level, and was
already taken into account fully for atomic DFT Naito2020J.Phys.B53_215002 .
It was found in Ref. Naito2020J.Phys.B53_215002  that the Breit correction
contributes to the total energies of atoms in almost the same magnitude but
the opposite to the density gradient effects. In the context of nuclear
physics, it was shown in Refs. LeBloas2011Phys.Rev.C84_014310 ; Roca-
Maza2016Phys.Rev.C94_044313 ; Naito2019Phys.Rev.C99_024309  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 (CDFT), 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
Meng2006Prog.Part.Nucl.Phys.57_470 ; Niksic2011Prog.Part.Nucl.Phys.66_519 ;
Niu2013Phys.Rev.C87_037301 ; Gu2013Phys.Rev.C87_041301 . Therefore, it is
interesting to take into account the relativistic correction to the Coulomb
interaction in SHF calculations to make a fair comparison with the CDFT
calculations for electromagnetic energies.

The relativistic correction to the Coulomb interaction was already discussed
in Refs. Gomez1973Phys.Lett.B44_231 ; Shlomo1976Phys.Lett.B60_244 . 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
Itoh1965Rev.Mod.Phys.37_159 ; Kenny1995Phys.Rev.A51_1898 . The aim of this
paper is modern calculation of Refs. Gomez1973Phys.Lett.B44_231 ;
Shlomo1976Phys.Lett.B60_244 : 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. Note
that it seems that “the Breit correction” refers the whole terms of the non-
relativistic reduced form of the correction [Eq. (5)] in the context of atomic
physics, while it refers only a term of the correction in the context of
nuclear physics. In order to avoid any confusion, hereinafter, I will use a
term “the relativistic correction”, instead of “the Breit correction”.

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

##  2 Theoretical Framework

In this section, I show the theoretical framework for the relativistic
correction and the DFT. Throughout this paper, an unit with
$4\pi\varepsilon_{0}=1$ is used and $c$ denotes the speed of light.

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

$H=T+V_{\scriptstyle{\text{{{{ext}}}}}}+V_{\scriptstyle{\text{{{{int}}}}}},$
(1)

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

$V_{\scriptstyle{\text{{{{Coul}}}}}}\left({\bm{r}}_{j},{\bm{r}}_{k}\right)=\frac{e^{2}}{r_{jk}}$
(2)

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

###  2.1 Breit correction for electronic systems

In electronic systems, the next-leading-order interaction
($O\left(1/c^{2}\right)$) in the Coulomb gauge is called the Breit correction,
whose form is Breit1929Phys.Rev.34_553 ; Breit1930Phys.Rev.36_383 ;
Breit1932Phys.Rev.39_616 ; Gorceix1988Phys.Rev.A37_1087 ;
Saue2011ChemPhysChem12_3077

$V_{\scriptstyle{\text{{{{Breit}}}}}}\left({\bm{r}}_{j},{\bm{r}}_{k}\right)=-\left[\frac{c{\bm{\alpha}}_{j}\cdot
c{\bm{\alpha}}_{k}}{2c^{2}r_{jk}}+\frac{\left(c{\bm{\alpha}}_{j}\cdot{\bm{r}}_{jk}\right)\left(c{\bm{\alpha}}_{k}\cdot{\bm{r}}_{jk}\right)}{2c^{2}r_{jk}^{3}}\right],$
(3)

where $\alpha$ is $4\times 4$ matrix defined by

${\bm{\alpha}}=\left(\begin{pmatrix}O_{2}&\sigma_{x}\\\
\sigma_{x}&O_{2}\end{pmatrix},\begin{pmatrix}O_{2}&\sigma_{y}\\\
\sigma_{y}&O_{2}\end{pmatrix},\begin{pmatrix}O_{2}&\sigma_{z}\\\
\sigma_{z}&O_{2}\end{pmatrix}\right)$ (4)

with the Pauli matrices $\sigma_{x}$, $\sigma_{y}$, and $\sigma_{z}$ and the
$2\times 2$ zero matrix $O_{2}$. In isolated atoms, the vacuum polarization of
the Coulomb field appears in higher order than the Breit correction
Eides2001Phys.Rep.342_63 . The Breit correction includes the magnetic
(current-current) interaction and the retardation in the Coulomb gauge
Gorceix1988Phys.Rev.A37_1087 . In the Lorenz (covariant) gauge, the next-
leading-order interaction corresponds to the magnetic interaction and the
retardation Sakurai1967AdvancedQuantumMechanics_AddisonWesley ;
Gorceix1988Phys.Rev.A37_1087 ; Lindroth1989Phys.Rev.A39_3794 .

After the Foldy-Wouthuysen-Tani transformation
Pryce1948Proc.RoyalSoc.Lond.A195_62 ; Foldy1950Phys.Rev.78_29 ;
Tani1951Prog.Theor.Phys.6_267 ; Foldy1952Phys.Rev.87_688 , the Breit
correction [Eq. (3)] in the non-relativistic scheme is obtained as
Itoh1965Rev.Mod.Phys.37_159 ; Kenny1995Phys.Rev.A51_1898

$\displaystyle\tilde{V}_{\scriptstyle{\text{{{{Breit}}}}}}\left({\bm{r}}_{j},{\bm{s}}_{j},{\bm{r}}_{k},{\bm{s}}_{j}\right)$
$\displaystyle=-\frac{\pi\hbar^{2}e^{2}}{m_{e}^{2}c^{2}}\delta\left({\bm{r}}_{jk}\right)$
(5a)
$\displaystyle\quad-\frac{e^{2}}{2m_{e}^{2}c^{2}}{\bm{p}}_{j}\cdot\left(\frac{1}{r_{jk}}+\frac{{\bm{r}}_{jk}{\bm{r}}_{jk}}{r_{jk}^{3}}\right)\cdot{\bm{p}}_{k}$
(5b)
$\displaystyle\quad-\frac{8\pi\hbar^{2}e^{2}}{3m_{e}^{2}c^{2}}\delta\left({\bm{r}}_{jk}\right){\bm{s}}_{j}\cdot{\bm{s}}_{k}$
(5c)
$\displaystyle\quad-\frac{\hbar^{2}e^{2}}{m_{e}^{2}c^{2}}{\bm{s}}_{j}\cdot\left(\frac{3{\bm{r}}_{jk}{\bm{r}}_{jk}}{r_{jk}^{5}}-\frac{1}{r_{jk}^{3}}\right)\cdot{\bm{s}}_{k}$
(5d)
$\displaystyle\quad+\frac{\hbar^{2}e^{2}}{m_{e}^{2}c^{2}}\frac{1}{r_{jk}^{3}}{\bm{s}}_{j}\cdot\left[{\bm{r}}_{jk}\times\left(2{\bm{p}}_{k}-{\bm{p}}_{j}\right)\right],$
(5e)

where $m_{e}\approx 0.511\,\mathrm{MeV}/c^{2}$ is the electron mass
Zyla2020Prog.Theor.Exp.Phys.2020_083C01  and ${\bm{p}}_{j}$ is the momentum of
the particle $j$. The first [Eq. (5a)] and second [Eq. (5b)] terms correspond
to the Darwin term for the electron-electron Coulomb interaction, which is
related to Zitterbewegung Sakurai1967AdvancedQuantumMechanics_AddisonWesley ,
and the retardation of the Coulomb interaction, respectively; the third and
fourth terms [Eqs. (5c) and (5d)] correspond to the spin-magnetic
interactions; and the fifth term [Eq. (5e)] corresponds to the electron-
electron spin-orbit interaction Itoh1965Rev.Mod.Phys.37_159 ;
Kenny1995Phys.Rev.A51_1898 ; Kenny1996Phys.Rev.Lett.77_1099 . Thus, in the
Schrödinger scheme, if one is eager to consider the Breit correction, Eq. (5)
should be considered on top of the Coulomb interaction [Eq. (2)], i.e.,

$V_{\scriptstyle{\text{{{{EM}}}}}}=\frac{1}{2}\sum_{j,\,k}\left[V_{\scriptstyle{\text{{{{Coul}}}}}}\left({\bm{r}}_{j},{\bm{r}}_{k}\right)+\tilde{V}_{\scriptstyle{\text{{{{Breit}}}}}}\left({\bm{r}}_{j},{\bm{s}}_{j},{\bm{r}}_{k},{\bm{s}}_{k}\right)\right],$
(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 Gomez1973Phys.Lett.B44_231 ; Shlomo1976Phys.Lett.B60_244 .
The correction was

$\displaystyle\tilde{V}_{\scriptstyle{\text{{{{Rel}}}}}}\left({\bm{r}}_{j},{\bm{s}}_{j},{\bm{r}}_{k},{\bm{s}}_{j}\right)$
$\displaystyle=-\frac{e_{\tau_{j}}e_{\tau_{k}}}{2M^{2}c^{2}}{\bm{p}}_{j}\cdot\left(\frac{1}{r_{jk}}+\frac{{\bm{r}}_{jk}{\bm{r}}_{jk}}{r_{jk}^{3}}\right)\cdot{\bm{p}}_{k}$
(7a)
$\displaystyle\quad-\frac{2\hbar\mu_{\tau_{j}}e_{\tau_{k}}}{Mc^{2}}\frac{1}{r_{jk}^{3}}{\bm{s}}_{j}\cdot\left({\bm{r}}_{jk}\times{\bm{p}}_{k}\right)+\frac{2\hbar\mu_{\tau_{k}}e_{\tau_{j}}}{Mc^{2}}\frac{1}{r_{jk}^{3}}{\bm{s}}_{j}\cdot\left({\bm{r}}_{jk}\times{\bm{p}}_{j}\right)$
(7b)
$\displaystyle\quad-\frac{4\mu_{\tau_{j}}\mu_{\tau_{k}}}{c^{2}}{\bm{s}}_{j}\cdot\left(\frac{3{\bm{r}}_{jk}{\bm{r}}_{jk}}{r_{jk}^{5}}-\frac{1}{r_{jk}^{3}}\right)\cdot{\bm{s}}_{k}$
(7c)
$\displaystyle\quad-\frac{32\pi}{3}\frac{\mu_{\tau_{j}}\mu_{\tau_{k}}}{c^{2}}\delta\left({\bm{r}}_{jk}\right){\bm{s}}_{j}\cdot{\bm{s}}_{k}$
(7d)
$\displaystyle\quad-\frac{e_{\tau_{k}}\hbar}{Mc^{2}}\left(2\mu_{\tau_{j}}-\frac{e_{\tau_{j}}\hbar}{2M}\right)\frac{1}{r_{jk}^{3}}{\bm{s}}_{j}\cdot\left({\bm{r}}_{jk}\times{\bm{p}}_{j}\right)$
(7e)
$\displaystyle\quad+\frac{e_{\tau_{j}}\hbar}{Mc^{2}}\left(2\mu_{\tau_{k}}-\frac{e_{\tau_{k}}\hbar}{2M}\right)\frac{1}{r_{jk}^{3}}{\bm{s}}_{k}\cdot\left({\bm{r}}_{jk}\times{\bm{p}}_{k}\right)$
(7f)
$\displaystyle\quad-\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({\bm{r}}_{jk}\right)$
(7g)
$\displaystyle=-\delta_{p\tau_{j}}\delta_{p\tau_{k}}\frac{e^{2}}{2M^{2}c^{2}}{\bm{p}}_{j}\cdot\left(\frac{1}{r_{jk}}+\frac{{\bm{r}}_{jk}{\bm{r}}_{jk}}{r_{jk}^{3}}\right)\cdot{\bm{p}}_{k}$
$\displaystyle\quad-\frac{\hbar^{2}e^{2}}{M^{2}c^{2}}\frac{1}{r_{jk}^{3}}{\bm{s}}_{j}\cdot\left[{\bm{r}}_{jk}\times\left(\tilde{\mu}_{\tau_{j}}\delta_{p\tau_{k}}{\bm{p}}_{k}-\tilde{\mu}_{\tau_{k}}\delta_{p\tau_{j}}{\bm{p}}_{j}\right)\right]$
$\displaystyle\quad-\tilde{\mu}_{\tau_{j}}\tilde{\mu}_{\tau_{k}}\frac{\hbar^{2}e^{2}}{M^{2}c^{2}}{\bm{s}}_{j}\cdot\left(\frac{3{\bm{r}}_{jk}{\bm{r}}_{jk}}{r_{jk}^{5}}-\frac{1}{r_{jk}^{3}}\right)\cdot{\bm{s}}_{k}$
$\displaystyle\quad-\tilde{\mu}_{\tau_{j}}\tilde{\mu}_{\tau_{k}}\frac{8\pi\hbar^{2}e^{2}}{3M^{2}c^{2}}\delta\left({\bm{r}}_{jk}\right){\bm{s}}_{j}\cdot{\bm{s}}_{k}$
$\displaystyle\quad-\delta_{p\tau_{k}}\frac{\hbar^{2}e^{2}}{2M^{2}c^{2}}\left(2\tilde{\mu}_{\tau_{j}}-\delta_{p\tau_{j}}\right)\frac{1}{r_{jk}^{3}}{\bm{s}}_{j}\cdot\left({\bm{r}}_{jk}\times{\bm{p}}_{j}\right)$
$\displaystyle\quad+\delta_{p\tau_{j}}\frac{\hbar^{2}e^{2}}{2M^{2}c^{2}}\left(2\tilde{\mu}_{\tau_{k}}-\delta_{p\tau_{k}}\right)\frac{1}{r_{jk}^{3}}{\bm{s}}_{k}\cdot\left({\bm{r}}_{jk}\times{\bm{p}}_{k}\right)$
$\displaystyle\quad-\frac{\pi\hbar^{2}e^{2}}{2M^{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({\bm{r}}_{jk}\right),$

where $M\approx 938.919\,\mathrm{MeV}/c^{2}$ is the nucleon mass
Zyla2020Prog.Theor.Exp.Phys.2020_083C01 , $\delta_{p\tau}$ is the Kronecker
delta, which is equal to $1$ for $\tau=p$ and $0$ for $\tau=n$, $e_{\tau}$ and
$\mu_{\tau}$ are the charge and magnetic moment of a nucleon $\tau$, and
$\tilde{\mu}_{\tau}$ is the magnetic moment of the unit of the nuclear
magneton, i.e., Zyla2020Prog.Theor.Exp.Phys.2020_083C01

$\displaystyle\mu_{\tau}$ $\displaystyle=\tilde{\mu}_{\tau}\frac{e\hbar}{2M},$
(8a) $\displaystyle\tilde{\mu}_{p}$ $\displaystyle=2.79284734463,$ (8b)
$\displaystyle\tilde{\mu}_{n}$ $\displaystyle=-1.91304273.$ (8c)

If one assumes that these two particles are electrons instead of nucleons, and
accordingly, $m_{e}$, $-e$, and $\mu_{e}$, are used instead of $M$,
$e_{\tau}$, and $\mu_{\tau}$, one can find that Eq. (7) is identical to Eq.
(5). Here, the relation between the $g$-factor and the magnetic moment
$\mu_{e}$

$\mu_{e}=\frac{e\hbar}{4m_{e}}g_{e}$ (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) are cancelled out of each
other.

Table 1:  Correspondence between Eqs. (5) and (7). For convenience, correspondence to Eq. (2) of Ref. Itoh1965Rev.Mod.Phys.37_159  is also shown. Eq. (7) | Eq. (5) | Eq. (2) of Ref. Itoh1965Rev.Mod.Phys.37_159  | 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

$E_{\scriptstyle{\text{{{{gs}}}}}}=T_{0}+E_{\scriptstyle{\text{{{{nucl}}}}}}\left[\rho_{p}^{\scriptstyle{\text{{{{gs}}}}}},\rho_{n}^{\scriptstyle{\text{{{{gs}}}}}}\right]+E_{\scriptstyle{\text{{{{EM}}}}}}\left[\rho_{p}^{\scriptstyle{\text{{{{gs}}}}}}\right],$
(10)

where $T_{0}$ is the Kohn-Sham kinetic energy,
$E_{\scriptstyle{\text{{{{nucl}}}}}}$ and $E_{\scriptstyle{\text{{{{EM}}}}}}$
are EDFs of nuclear and EM interactions, respectively, and
$\rho^{p}_{\scriptstyle{\text{{{{gs}}}}}}$ and
$\rho^{n}_{\scriptstyle{\text{{{{gs}}}}}}$ are the ground-state densities of
protons and neutrons. If one considers the finite charge distributions of
nucleons are considered, $E_{\scriptstyle{\text{{{{EM}}}}}}$ should depend on
the ground-state charge distribution, instead of
$\rho_{p}^{\scriptstyle{\text{{{{gs}}}}}}$ Naito2020Phys.Rev.C101_064311 .

If one only considers the Coulomb interaction to
$V_{\scriptstyle{\text{{{{EM}}}}}}$, $E_{\scriptstyle{\text{{{{EM}}}}}}$ in
LDA reads Slater1951Phys.Rev.81_385

$\displaystyle E_{\scriptstyle{\text{{{{EM}}}}}}\left[\rho\right]$
$\displaystyle=E_{\scriptstyle{\text{{{{Coul}}}}}}^{\scriptstyle{\text{{{{H}}}}}}\left[\rho\right]+E_{\scriptstyle{\text{{{{Coul}}}}}}^{\scriptstyle{\text{{{{x}}}}}}\left[\rho\right],$
(11a) $\displaystyle
E_{\scriptstyle{\text{{{{Coul}}}}}}^{\scriptstyle{\text{{{{H}}}}}}\left[\rho\right]$
$\displaystyle=\frac{e^{2}}{2}\iint\frac{\rho\left({\bm{r}}\right)\rho\left({\bm{r}}^{\prime}\right)}{\left|{\bm{r}}-{\bm{r}}^{\prime}\right|}\,d{\bm{r}}\,d{\bm{r}}^{\prime},$
(11b) $\displaystyle
E_{\scriptstyle{\text{{{{Coul}}}}}}^{\scriptstyle{\text{{{{x}}}}}}\left[\rho\right]$
$\displaystyle=-\frac{3e^{2}}{4}\left(\frac{3}{\pi}\right)^{1/3}\int\left[\rho\left({\bm{r}}\right)\right]^{4/3}\,d{\bm{r}},$
(11c)

where the density $\rho$ is the electron density $\rho_{e}$ for electronic
systems and the proton density $\rho_{p}$ for nuclear systems.

Recently, an EDF for the Breit correction [Eq. (5)] for electronic systems in
the LDA was developed in Refs. Kenny1996Phys.Rev.Lett.77_1099 ;
Naito2020J.Phys.B53_215002 . In this paper, I consider the relativistic
correction [Eq. (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 111 Effects
of nucleon finite size were discussed in Refs. Auerbach1983Phys.Rep.98_273 ;
Naito2020Phys.Rev.C101_064311 ..

The Hartree and LDA exchange functionals for the Breit correction in the
electronic systems [Eq. (5)] were proposed in Ref. Naito2020J.Phys.B53_215002
and Ref. MacDonald1979J.Phys.C12_2977 , respectively. Their forms are

$\displaystyle
E_{\scriptstyle{\text{{{{Breit}}}}}}^{\scriptstyle{\text{{{{H}}}}}}\left[\rho\right]$
$\displaystyle=-\frac{\pi\hbar^{2}e^{2}}{2m_{e}^{2}c^{2}}\int\left[\rho\left({\bm{r}}\right)\right]^{2}\,d{\bm{r}},$
(12a) $\displaystyle
E_{\scriptstyle{\text{{{{Breit}}}}}}^{\scriptstyle{\text{{{{x}}}}}}\left[\rho\right]$
$\displaystyle=\frac{3\pi\hbar^{2}e^{2}}{2m_{e}^{2}c^{2}}\int\left[\rho\left({\bm{r}}\right)\right]^{2}\,d{\bm{r}},$
(12b) respectively, and accordingly, the total correction is $\displaystyle
E_{\scriptstyle{\text{{{{Breit}}}}}}^{\scriptstyle{\text{{{{tot}}}}}}\left[\rho\right]$
$\displaystyle=E_{\scriptstyle{\text{{{{Breit}}}}}}^{\scriptstyle{\text{{{{H}}}}}}\left[\rho\right]+E_{\scriptstyle{\text{{{{Breit}}}}}}^{\scriptstyle{\text{{{{x}}}}}}\left[\rho\right]$
$\displaystyle=\frac{\pi\hbar^{2}e^{2}}{m_{e}^{2}c^{2}}\int\left[\rho\left({\bm{r}}\right)\right]^{2}\,d{\bm{r}}.$
(12c)

The idea of the LDA is that the energy density is approximated to that of
homogeneous systems; hence, Eqs. (12a) and (12b) contain only effects
originating from Eqs. (5a) and (5b) Kenny1996Phys.Rev.Lett.77_1099 . In
addition, due to the symmetry, Eq. (12a) does not contain the effect of the
retardation [Eq. (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). Note that the spin-orbit interaction originating
from the magnetic form factors of nucleons [Eqs. (7b), (7e), and (7f)] have
been discussed in Refs. Roca-Maza2018Phys.Rev.Lett.120_202501 ; Roca-
Maza2018EPJWebConf.194_01002 ; Naito2020Phys.Rev.C101_064311 ; the spin-
magnetic interaction [Eqs. (7c) and (7d)] has not been considered yet for the
self-consistent calculation and it cannot be considered in the LDA, which is
left for the future study.

Hereinafter, I will construct the LDA EDF for Eq. (7). Referring the result on
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

$\displaystyle
E_{\scriptstyle{\text{{{{Darwin}}}}}}^{\scriptstyle{\text{{{{H}}}}}}\left[\rho_{p},\rho_{n}\right]$
$\displaystyle=\int\left\\{\frac{C_{pp}}{2}\left[\rho_{p}\left({\bm{r}}\right)\right]^{2}+C_{pn}\rho_{p}\left({\bm{r}}\right)\rho_{n}\left({\bm{r}}\right)\right.$
$\displaystyle\qquad\left.+\frac{C_{nn}}{2}\left[\rho_{n}\left({\bm{r}}\right)\right]^{2}\right\\}\,d{\bm{r}}$
$\displaystyle=\int\left\\{\frac{C_{pp}}{2}\left[\rho_{p}\left({\bm{r}}\right)\right]^{2}+C_{pn}\rho_{p}\left({\bm{r}}\right)\rho_{n}\left({\bm{r}}\right)\right\\}\,d{\bm{r}},$
(13a) $\displaystyle
E_{\scriptstyle{\text{{{{Darwin}}}}}}^{\scriptstyle{\text{{{{x}}}}}}\left[\rho_{p},\rho_{n}\right]$
$\displaystyle=-\int\left\\{\frac{C_{pp}}{4}\left[\rho_{p}\left({\bm{r}}\right)\right]^{2}+\frac{C_{nn}}{4}\left[\rho_{p}\left({\bm{r}}\right)\right]^{2}\right\\}\,d{\bm{r}}$
$\displaystyle=-\int\frac{C_{pp}}{4}\left[\rho_{p}\left({\bm{r}}\right)\right]^{2}\,d{\bm{r}},$
(13b)

where

$C_{\tau_{j}\tau_{k}}=-\frac{\pi\hbar^{2}e^{2}}{2M^{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]$
(14)

and accordingly

$\displaystyle C_{pp}$
$\displaystyle=-\frac{\pi\hbar^{2}e^{2}}{M^{2}c^{2}}\left(2\tilde{\mu}_{p}-1\right),$
(15a) $\displaystyle C_{pn}$
$\displaystyle=-\frac{\pi\hbar^{2}e^{2}}{M^{2}c^{2}}\tilde{\mu}_{n},$ (15b)
$\displaystyle C_{nn}$ $\displaystyle=0.$ (15c)

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

$\displaystyle\frac{3\pi\hbar^{2}e^{2}}{2m_{e}^{2}c^{2}}\int\left[\rho\left({\bm{r}}\right)\right]^{2}\,d{\bm{r}}-\frac{1}{4}\frac{\pi\hbar^{2}e^{2}}{m_{e}^{2}c^{2}}\int\left[\rho\left({\bm{r}}\right)\right]^{2}\,d{\bm{r}}$
$\displaystyle=\frac{5\pi\hbar^{2}e^{2}}{4m_{e}^{2}c^{2}}\int\left[\rho\left({\bm{r}}\right)\right]^{2}\,d{\bm{r}},$
(16)

since the Darwin EDF for electronic systems is

$-\frac{1}{4}\frac{\pi\hbar^{2}e^{2}}{m_{e}c^{2}}\int\left[\rho\left({\bm{r}}\right)\right]^{2}\,d{\bm{r}}.$
(17)

To apply the nuclear systems, with replacing the electron mass $m_{e}$ to the
nucleon mass $M$, one obtains

$E_{\scriptstyle{\text{{{{retar}}}}}}^{\scriptstyle{\text{{{{x}}}}}}\left[\rho_{p},\rho_{n}\right]=\frac{5\pi\hbar^{2}e^{2}}{4M^{2}c^{2}}\int\left[\rho_{p}\left({\bm{r}}\right)\right]^{2}\,d{\bm{r}}.$
(18)

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

$\displaystyle
E_{\scriptstyle{\text{{{{Rel}}}}}}^{\scriptstyle{\text{{{{tot}}}}}}\left[\rho_{p},\rho_{n}\right]$
$\displaystyle=\int\left\\{\frac{C_{pp}}{2}\left[\rho_{p}\left({\bm{r}}\right)\right]^{2}+C_{pn}\rho_{p}\left({\bm{r}}\right)\rho_{n}\left({\bm{r}}\right)\right\\}\,d{\bm{r}}$
$\displaystyle\quad-\int\frac{C_{pp}}{4}\left[\rho_{p}\left({\bm{r}}\right)\right]^{2}\,d{\bm{r}}+\frac{5\pi\hbar^{2}e^{2}}{4M^{2}c^{2}}\int\left[\rho_{p}\left({\bm{r}}\right)\right]^{2}\,d{\bm{r}}$
$\displaystyle=\left(\frac{C_{pp}}{2}-\frac{C_{pp}}{4}+\frac{5\pi\hbar^{2}e^{2}}{4M^{2}c^{2}}\right)\int\left[\rho_{p}\left({\bm{r}}\right)\right]^{2}\,d{\bm{r}}$
$\displaystyle\quad+C_{pn}\int\rho_{p}\left({\bm{r}}\right)\rho_{n}\left({\bm{r}}\right)\,d{\bm{r}}$
$\displaystyle=\left(3-\tilde{\mu}_{p}\right)\frac{\pi\hbar^{2}e^{2}}{2M^{2}c^{2}}\int\left[\rho_{p}\left({\bm{r}}\right)\right]^{2}\,d{\bm{r}}$
$\displaystyle\quad-\tilde{\mu}_{n}\frac{\pi\hbar^{2}e^{2}}{M^{2}c^{2}}\int\rho_{p}\left({\bm{r}}\right)\rho_{n}\left({\bm{r}}\right)\,d{\bm{r}}$
(19)
$\displaystyle\approx\frac{\pi\hbar^{2}e^{2}}{M^{2}c^{2}}\int\left[0.1035\rho_{p}\left({\bm{r}}\right)+1.913\rho_{n}\left({\bm{r}}\right)\right]\rho_{p}\left({\bm{r}}\right)\,d{\bm{r}}.$
(20)

##  3 Analytical Estimation

Before performing numerical many-body calculations with the relativistic
correction, I estimate contribution of the correction to the total energy with
a simple model and discuss the systematic behavior. As in Ref.
Naito2020Phys.Rev.C101_064311 , I assume hard-sphere distributions for both
proton and neutron distribution $\rho_{p}$ and $\rho_{n}$:

$\rho_{\tau}\left(r\right)=\begin{cases}\rho_{0}^{\tau}&r\leq R_{\tau},\\\
0&r>R_{\tau},\end{cases}$ (21)

where $R_{\tau}$ is the proton ($\tau=p$) or neutron ($\tau=n$) radius defined
by

$R_{\tau}=\left(\frac{3N_{\tau}}{4\pi\rho_{0}^{\tau}}\right)^{1/3}$ (22)

with $N_{p}=Z$ and $N_{n}=N$. For simplicity, the 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 Naito2020Phys.Rev.C101_064311

$\displaystyle
E_{\scriptstyle{\text{{{{Coul}}}}}}^{\scriptstyle{\text{{{{H}}}}}}$
$\displaystyle=\frac{3e^{2}}{5}\left(\frac{4\pi\rho_{0}^{p}}{3}\right)^{1/3}Z^{5/3}$
$\displaystyle\simeq 0.60Z^{5/3}\,\mathrm{MeV},$ (23a) $\displaystyle
E_{\scriptstyle{\text{{{{Coul}}}}}}^{\scriptstyle{\text{{{{x}}}}}}$
$\displaystyle=-\frac{3e^{2}}{4}\left(\frac{3\rho_{0}^{p}}{\pi}\right)^{1/3}Z$
$\displaystyle\simeq-0.46Z\,\mathrm{MeV},$ (23b)

respectively.

Since relations

$\displaystyle\int\left[\rho_{p}\left({\bm{r}}\right)\right]^{2}\,d{\bm{r}}$
$\displaystyle=\frac{4\pi R_{p}^{3}}{3}\left(\rho_{0}^{p}\right)^{2}$
$\displaystyle=\rho_{0}^{p}Z,$ (24a)
$\displaystyle\int\rho_{p}\left({\bm{r}}\right)\rho_{n}\left({\bm{r}}\right)\,d{\bm{r}}$
$\displaystyle=\frac{4\pi R_{<}^{3}}{3}\rho_{p}^{0}\rho_{n}^{0}$
$\displaystyle=\begin{cases}\rho_{n}^{0}Z&\text{for $R_{p}<R_{n}$},\\\
\rho_{p}^{0}N&\text{for $R_{p}>R_{n}$},\end{cases}$ (24b)

hold with $R_{<}=R_{p}$ for $R_{p}<R_{n}$ and $R_{<}=R_{n}$ for $R_{p}>R_{n}$,
one, accordingly, obtains

$\displaystyle
E_{\scriptstyle{\text{{{{Darwin}}}}}}^{\scriptstyle{\text{{{{H}}}}}}$
$\displaystyle=\frac{C_{pp}}{2}\frac{4\pi
R_{p}^{3}}{3}\left(\rho_{0}^{p}\right)^{2}+C_{pn}\frac{4\pi
R_{<}^{3}}{3}\rho_{p}^{0}\rho_{n}^{0}$
$\displaystyle=-\frac{2\pi^{2}\hbar^{2}e^{2}}{3M^{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_{<}^{3}\rho_{p}^{0}\rho_{n}^{0}\right],$
(25a) $\displaystyle
E_{\scriptstyle{\text{{{{Darwin}}}}}}^{\scriptstyle{\text{{{{x}}}}}}$
$\displaystyle=-\frac{C_{pp}}{4}\frac{4\pi
R_{p}^{3}}{3}\left(\rho_{0}^{p}\right)^{2}$
$\displaystyle=\frac{\pi^{2}\hbar^{2}e^{2}}{3M^{2}c^{2}}\left(2\tilde{\mu}_{p}-1\right)R_{p}^{3}\left(\rho_{0}^{p}\right)^{2}$
$\displaystyle=\frac{\pi\hbar^{2}e^{2}}{4M^{2}c^{2}}\left(2\tilde{\mu}_{p}-1\right)\rho_{0}^{p}Z$
$\displaystyle\simeq 0.018Z\,\mathrm{MeV},$ (25b) $\displaystyle
E_{\scriptstyle{\text{{{{retar}}}}}}^{\scriptstyle{\text{{{{x}}}}}}$
$\displaystyle=\frac{5\pi^{2}\hbar^{2}e^{2}}{3M^{2}c^{2}}R_{p}^{3}\left(\rho_{0}^{p}\right)^{2}$
$\displaystyle=\frac{5\pi\hbar^{2}e^{2}}{4M^{2}c^{2}}\rho_{0}^{p}Z$
$\displaystyle\simeq 0.020Z\,\mathrm{MeV},$ (25c) $\displaystyle
E_{\scriptstyle{\text{{{{Rel}}}}}}^{\scriptstyle{\text{{{{tot}}}}}}$
$\displaystyle=\frac{2\pi^{2}\hbar^{2}e^{2}}{3M^{2}c^{2}}\left[\left(3-\tilde{\mu}_{p}\right)R_{p}^{3}\left(\rho_{0}^{p}\right)^{2}-2\tilde{\mu}_{n}R_{<}^{3}\rho_{p}^{0}\rho_{n}^{0}\right].$
(25d)

In the case of $R_{p}<R_{n}$, Eqs. (25a) and (25d) read

$\displaystyle
E_{\scriptstyle{\text{{{{Darwin}}}}}}^{\scriptstyle{\text{{{{H}}}}}}$
$\displaystyle=-\frac{2\pi^{2}\hbar^{2}e^{2}}{3M^{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}$
$\displaystyle=-\frac{\pi\hbar^{2}e^{2}}{2M^{2}c^{2}}\left[\left(2\tilde{\mu}_{p}-1\right)\rho_{0}^{p}+2\tilde{\mu}_{n}\rho_{n}^{0}\right]Z$
$\displaystyle\simeq-0.006Z\,\mathrm{MeV},$ (26a) $\displaystyle
E_{\scriptstyle{\text{{{{Rel}}}}}}^{\scriptstyle{\text{{{{tot}}}}}}$
$\displaystyle=\frac{2\pi^{2}\hbar^{2}e^{2}}{3M^{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]$
$\displaystyle=\frac{\pi\hbar^{2}e^{2}}{2M^{2}c^{2}}\left[\left(3-\tilde{\mu}_{p}\right)\rho_{0}^{p}-2\tilde{\mu}_{n}\rho_{n}^{0}\right]Z$
$\displaystyle\simeq 0.032Z\,\mathrm{MeV};$ (26b)

in the case of $R_{p}>R_{n}$, they read

$\displaystyle
E_{\scriptstyle{\text{{{{Darwin}}}}}}^{\scriptstyle{\text{{{{H}}}}}}$
$\displaystyle=-\frac{2\pi^{2}\hbar^{2}e^{2}}{3M^{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]$
$\displaystyle=-\frac{\pi\hbar^{2}e^{2}}{2M^{2}c^{2}}\rho_{0}^{p}\left[\left(2\tilde{\mu}_{p}-1\right)Z+2\tilde{\mu}_{n}N\right]$
$\displaystyle\simeq-0.037Z+0.031N\,\mathrm{MeV},$ (27a) $\displaystyle
E_{\scriptstyle{\text{{{{Rel}}}}}}^{\scriptstyle{\text{{{{tot}}}}}}$
$\displaystyle=\frac{2\pi^{2}\hbar^{2}e^{2}}{3M^{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]$
$\displaystyle=\frac{\pi\hbar^{2}e^{2}}{2M^{2}c^{2}}\rho_{0}^{p}\left[\left(3-\tilde{\mu}_{p}\right)Z-2\tilde{\mu}_{n}N\right]$
$\displaystyle\simeq 0.002Z+0.031N\,\mathrm{MeV}.$ (27b)

Comparing Eqs. (23b) and (26b), one can find that the relativistic correction
to the total energy
$\left|E_{\scriptstyle{\text{{{{Rel}}}}}}^{\scriptstyle{\text{{{{tot}}}}}}\right|$
is $15$ times smaller than the Coulomb exchange energy
$\left|E_{\scriptstyle{\text{{{{Coul}}}}}}^{\scriptstyle{\text{{{{x}}}}}}\right|$.
This means that although
$E_{\scriptstyle{\text{{{{Rel}}}}}}^{\scriptstyle{\text{{{{tot}}}}}}$ is small
compared to the other contributions, one can estimate that
$E_{\scriptstyle{\text{{{{Rel}}}}}}^{\scriptstyle{\text{{{{tot}}}}}}$ for
${}^{208}\mathrm{Pb}$ may be around $2\,\mathrm{MeV}$, which is a non-
negligible contribution compared to the desired accuracy of
$E_{\scriptstyle{\text{{{{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
Eides2001Phys.Rep.342_63 . On the other hand, in nuclear structure, both the
vacuum polarization and the relativistic correction are corrections to the
proton-proton Coulomb interaction. Thus, it is worthwhile to compare them. In
Ref. Naito2020Phys.Rev.C101_064311 , the vacuum polarization was taken into
account in a nuclear structure calculation using the Uehling potential
Uehling1935Phys.Rev.48_55 ; WayneFullerton1976Phys.Rev.A13_1283 . The
contribution of the vacuum polarization was estimated as

$E_{\scriptstyle{\text{{{{VP}}}}}}\simeq 0.0047Z^{5/3}\,\mathrm{MeV}.$ (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_{\scriptstyle{\text{{{{Rel}}}}}}^{\scriptstyle{\text{{{{tot}}}}}}\approx
0.64\,\mathrm{MeV}$ and $E_{\scriptstyle{\text{{{{VP}}}}}}\approx
0.69\,\mathrm{MeV}$) but the latter dominates in heavy nuclei (e.g., for
$Z=100$,
$E_{\scriptstyle{\text{{{{Rel}}}}}}^{\scriptstyle{\text{{{{tot}}}}}}\approx
3.2\,\mathrm{MeV}$ and $E_{\scriptstyle{\text{{{{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 [Eq. (19)] is implemented to the self-consistent
Skyrme Hartree-Fock calculation code skyrme_rpa
Colo2013Comput.Phys.Commun.184_142 . In this work, doubly-magic nuclei are
focused on; thus, the pairing interaction is not considered, and the spherical
symmetry is assumed. The calculations are performed with a box of
$15\,\mathrm{fm}$ with $0.1\,\mathrm{fm}$ mesh. The SAMi functional Roca-
Maza2012Phys.Rev.C86_031306  is adopted for the nuclear EDF.

First of all, I shall mention effects to nuclear radii. Neither the proton nor
the neutron radii are found to be changed 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 labelled “Coulomb (NR-LDA)” and “Coulomb (R-LDA)” show results
with the non-relativistic Coulomb LDA [Eq. (11)] and the relativistic Coulomb
LDA [Eqs. (11) and (19)], respectively. For comparison, total energies
calculated with the vacuum polarization on top of the relativistic Coulomb LDA
and with the non-relativistic Coulomb GGA Naito2019Phys.Rev.C99_024309  are
shown in a column labelled “C. (R-LDA) & V.P.” and “Coulomb (NR-GGA)”,
respectively. Both the relativistic correction and the vacuum polarization
make the atomic nuclei less 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_{\scriptstyle{\text{{{{tot}}}}}}$ of ${}^{48}\mathrm{Ca}$ and
${}^{48}\mathrm{Ni}$ by $0.5\,\mathrm{MeV}$, the mass difference is not
changed.

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

I also compare the relativistic correction and the vacuum polarization. The
vacuum polarization also makes the atomic nuclei less bound, and 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 Sec. 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 difference of the
second and the third columns and the third and the fourth columns in Table 2,
respectively. The proton number $Z$ dependence of these contributions are
shown in Fig. 1. These energies are fitted with $E=aZ^{b}$, where $a$ and $b$
are adjustable parameters. The results are

$\displaystyle
E_{\scriptstyle{\text{{{{Rel}}}}}}^{\scriptstyle{\text{{{{tot}}}}}}$
$\displaystyle\simeq 0.01660Z^{1.12}\,\mathrm{MeV},$ (29a) $\displaystyle
E_{\scriptstyle{\text{{{{VP}}}}}}$ $\displaystyle\simeq
0.00398Z^{1.55}\,\mathrm{MeV},$ (29b)

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

Table 2: Ground-state energies calculated with the SAMi functional Roca-Maza2012Phys.Rev.C86_031306  shown in $\mathrm{MeV}$. The columns labelled “Coulomb (NR-LDA)” and “Coulomb (R-LDA)” show results calculated with the non-relativistic Coulomb LDA [Eq. (11)] and the relativistic Coulomb LDA [Eqs. (11) and (19)], respectively. For comparison, total energies calculated with the vacuum polarization on top of the relativistic Coulomb LDA and with the non-relativistic Coulomb GGA are shown in a column labelled “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$
Figure 1: Contribution of relativistic correction and vacuum polarization to
the total energies,
$E_{\scriptstyle{\text{{{{Rel}}}}}}^{\scriptstyle{\text{{{{tot}}}}}}$ and
$E_{\scriptstyle{\text{{{{VP}}}}}}$ shown as functions of proton number $Z$ in
circle and square points, respectively. Empirical fits with Eqs. (29) are also
shown by the solid and the dashed lines, respectively.

The breakdown of the relativistic correction [Eq. (19)] is shown in Table 3.
The columns labelled
$E_{\scriptstyle{\text{{{{Darwin}}}}}}^{\scriptstyle{\text{{{{H}}}}}}$,
$E_{\scriptstyle{\text{{{{Darwin}}}}}}^{\scriptstyle{\text{{{{x}}}}}}$, and
$E_{\scriptstyle{\text{{{{retar}}}}}}^{\scriptstyle{\text{{{{x}}}}}}$ are
Darwin Hartree energy [Eq. (13a)], Darwin exchange one [Eq. (13b)], and
retardation exchange one [Eq. (18)], respectively. As shown, in $N\leq Z$
nuclei,
$E_{\scriptstyle{\text{{{{Darwin}}}}}}^{\scriptstyle{\text{{{{H}}}}}}<0$
holds, while in $N>Z$ nuclei,
$E_{\scriptstyle{\text{{{{Darwin}}}}}}^{\scriptstyle{\text{{{{H}}}}}}>0$
holds. The Darwin Hartree energy reads

$\displaystyle
E_{\scriptstyle{\text{{{{Darwin}}}}}}^{\scriptstyle{\text{{{{H}}}}}}$
$\displaystyle=-\frac{\pi\hbar^{2}e^{2}}{M^{2}c^{2}}\int\left\\{\frac{2\tilde{\mu}_{p}-1}{2}\rho_{p}\left({\bm{r}}\right)+\tilde{\mu}_{n}\rho_{n}\left({\bm{r}}\right)\right\\}\rho_{p}\left({\bm{r}}\right)\,d{\bm{r}}$
$\displaystyle\approx-\frac{\pi\hbar^{2}e^{2}}{M^{2}c^{2}}\int\left[2.292\rho_{p}\left({\bm{r}}\right)-1.913\rho_{n}\left({\bm{r}}\right)\right]\rho_{p}\left({\bm{r}}\right)\,d{\bm{r}};$
(30)

hence,
$\left[2.292\rho_{p}\left({\bm{r}}\right)-1.913\rho_{n}\left({\bm{r}}\right)\right]$
determines the sign of
$E_{\scriptstyle{\text{{{{Darwin}}}}}}^{\scriptstyle{\text{{{{H}}}}}}$. In
$N\leq Z$ nuclei,
$\rho_{p}\left({\bm{r}}\right)\gtrsim\rho_{n}\left({\bm{r}}\right)$ holds and,
accordingly, it is negative (i.e.,
$E_{\scriptstyle{\text{{{{Darwin}}}}}}^{\scriptstyle{\text{{{{H}}}}}}>0$),
whereas, in $N<Z$ nuclei,
$\rho_{p}\left({\bm{r}}\right)\lesssim\rho_{n}\left({\bm{r}}\right)$ holds and
accordingly, it can be positive (i.e.,
$E_{\scriptstyle{\text{{{{Darwin}}}}}}^{\scriptstyle{\text{{{{H}}}}}}<0$) in
contrast to Eq. (25a). The Darwin and retardation exchange energies are also
fitted with $E=aZ^{b}$, where the results are

$\displaystyle
E_{\scriptstyle{\text{{{{Darwin}}}}}}^{\scriptstyle{\text{{{{x}}}}}}$
$\displaystyle\simeq 0.01312Z^{0.97}\,\mathrm{MeV},$ (31a) $\displaystyle
E_{\scriptstyle{\text{{{{retar}}}}}}^{\scriptstyle{\text{{{{x}}}}}}$
$\displaystyle\simeq 0.01431Z^{0.97}\,\mathrm{MeV},$ (31b)

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

Table 3: Breakdown of the relativistic correction [Eq. (19)] shown in $\mathrm{MeV}$. The columns labelled $E_{\scriptstyle{\text{{{{Darwin}}}}}}^{\scriptstyle{\text{{{{H}}}}}}$, $E_{\scriptstyle{\text{{{{Darwin}}}}}}^{\scriptstyle{\text{{{{x}}}}}}$, and $E_{\scriptstyle{\text{{{{retar}}}}}}^{\scriptstyle{\text{{{{x}}}}}}$ are Darwin Hartree energy [Eq. (13a)], Darwin exchange energy [Eq. (13b)], and retardation exchange energy [Eq. (18)], respectively. Nuclei | $E_{\scriptstyle{\text{{{{Darwin}}}}}}^{\scriptstyle{\text{{{{H}}}}}}$ | $E_{\scriptstyle{\text{{{{Darwin}}}}}}^{\scriptstyle{\text{{{{x}}}}}}$ | $E_{\scriptstyle{\text{{{{retar}}}}}}^{\scriptstyle{\text{{{{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 exact Coulomb
exchange energy

$\Delta
E_{\scriptstyle{\text{{{{x}}}}}}=\frac{E_{\scriptstyle{\text{{{{x}}}}}}^{\scriptstyle{\text{{{{LDA}}}}}}-E_{\scriptstyle{\text{{{{x}}}}}}^{\scriptstyle{\text{{{{exact}}}}}}}{E_{\scriptstyle{\text{{{{x}}}}}}^{\scriptstyle{\text{{{{exact}}}}}}}$
(32)

is shown in Fig. 2. Here, the exact Coulomb exchange energy is calculated by
using the modified GGA-PBE functional with $\lambda=1.25$, which reproduces
the exact Coulomb exchange energy less than $500\,\mathrm{keV}$ error
Naito2019Phys.Rev.C99_024309 . The non-relativistic and relativistic Coulomb
LDA exchange energies are calculated by using
$E_{\scriptstyle{\text{{{{Coul}}}}}}^{\scriptstyle{\text{{{{x}}}}}}$ [Eq.
(11c)] and
$E_{\scriptstyle{\text{{{{Coul}}}}}}^{\scriptstyle{\text{{{{x}}}}}}+E_{\scriptstyle{\text{{{{Darwin}}}}}}^{\scriptstyle{\text{{{{x}}}}}}+E_{\scriptstyle{\text{{{{retar}}}}}}^{\scriptstyle{\text{{{{x}}}}}}$
[Eqs. (13b) and (18)], respectively. It is seen that the absolute value of the
relative change $\left|\Delta E_{\scriptstyle{\text{{{{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_{\scriptstyle{\text{{{{x}}}}}}\approx-3\,\%$ in heavy nuclei
($A\gtrsim 100$). The absolute value of the relative change $\left|\Delta
E_{\scriptstyle{\text{{{{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_{\scriptstyle{\text{{{{x}}}}}}\approx-9\,\%$ in heavy nuclei. The Darwin
contribution and the retardation one are almost the same size.

This difference has been studied in the relativistic (covariant) density
functional theory Gu2013Phys.Rev.C87_041301 ; Niu2013Phys.Rev.C87_037301 . The
exact Coulomb exchange energy in this paper corresponds to Eq. (4) in Ref.
Gu2013Phys.Rev.C87_041301 , which does not seem to include the relativistic
correction. They also show that difference between the non-relativistic and
relativistic Coulomb LDA exchange energies is approximately $1\,\mathrm{MeV}$.
With considering such difference,
$E_{\scriptstyle{\text{{{{Coul}}}}}}^{\scriptstyle{\text{{{{x}}}}}}+E_{\scriptstyle{\text{{{{Rel}}}}}}^{\scriptstyle{\text{{{{x}}}}}}$
calculated in this work is almost consistent qualitatively with, but almost
twice of Figs. 1 and 4 of Ref. Gu2013Phys.Rev.C87_041301 , because Ref.
Gu2013Phys.Rev.C87_041301  does not include
$E_{\scriptstyle{\text{{{{Darwin}}}}}}^{\scriptstyle{\text{{{{x}}}}}}$. If one
considers
$E_{\scriptstyle{\text{{{{Darwin}}}}}}^{\scriptstyle{\text{{{{x}}}}}}$ only,
this work is quantitatively consistent with Figs. 1 and 4 of Ref.
Gu2013Phys.Rev.C87_041301 .

Figure 2: Ratio of the non-relativistic (square) or the relativistic (circle)
Coulomb LDA exchange energies to the exact Coulomb exchange energy, $\Delta
E_{\scriptstyle{\text{{{{x}}}}}}$ [Eq. (32)], shown as functions of mass
number $A$. The Darwin contribution only (up-triangle) and the retardation
contribution only (down-triangle) are also shown.

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 on top
of the relativistic LDA are shown. It is seen that the proton single-particle
energies change by about $0.04\,\mathrm{MeV}$ in inner shells while
$0.01\,\mathrm{MeV}$ in outer shells, due to the relativistic correction. In
contrast, the change due to the vacuum polarization does not depend on
orbitals strongly, and its value is about $0.1\,\mathrm{MeV}$. This difference
is due to the nature of interaction; the relativistic correction is the point-
coupling interaction, while the vacuum polarization is 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 while $0.01\,\mathrm{MeV}$ in outer
shells, due to the relativistic correction. In contrast, the vacuum
polarization changes neutron single-particle energies only by less than
$0.01\,\mathrm{MeV}$.

Table 4: Single-particle energies of ${}^{208}\mathrm{Pb}$ shown in $\mathrm{MeV}$ calculated with the SAMi functional. States | Coulomb | C. (R-LDA) & V.P.
---|---|---
NR-LDA | R-LDA |
$\pi 1s_{1/2}$ | $-44.980$ | $-44.936$ | $-44.829$
$\pi 1p_{3/2}$ | $-39.387$ | $-39.346$ | $-39.246$
$\pi 1p_{1/2}$ | $-39.107$ | $-39.067$ | $-38.967$
$\pi 1d_{5/2}$ | $-32.482$ | $-32.446$ | $-32.352$
$\pi 1d_{3/2}$ | $-31.815$ | $-31.779$ | $-31.686$
$\pi 2s_{1/2}$ | $-28.509$ | $-28.475$ | $-28.379$
$\pi 1f_{7/2}$ | $-24.692$ | $-24.661$ | $-24.572$
$\pi 1f_{5/2}$ | $-23.353$ | $-23.323$ | $-23.235$
$\pi 2p_{3/2}$ | $-19.411$ | $-19.383$ | $-19.292$
$\pi 2p_{1/2}$ | $-18.626$ | $-18.598$ | $-18.508$
$\pi 1g_{9/2}$ | $-16.338$ | $-16.312$ | $-16.229$
$\pi 1g_{7/2}$ | $-14.019$ | $-13.995$ | $-13.913$
$\pi 2d_{5/2}$ | $-10.255$ | $-10.233$ | $-10.148$
$\pi 2d_{3/2}$ | $-8.846$ | $-8.826$ | $-8.742$
$\pi 3s_{1/2}$ | $-7.673$ | $-7.653$ | $-7.568$
$\pi 1h_{11/2}$ | $-7.663$ | $-7.642$ | $-7.563$
$\nu 1s_{1/2}$ | $-59.291$ | $-59.253$ | $-59.238$
$\nu 1p_{3/2}$ | $-52.953$ | $-52.919$ | $-52.908$
$\nu 1p_{1/2}$ | $-52.656$ | $-52.622$ | $-52.612$
$\nu 1d_{5/2}$ | $-45.375$ | $-45.345$ | $-45.338$
$\nu 1d_{3/2}$ | $-44.744$ | $-44.715$ | $-44.709$
$\nu 2s_{1/2}$ | $-41.962$ | $-41.935$ | $-41.928$
$\nu 1f_{7/2}$ | $-36.904$ | $-36.880$ | $-36.876$
$\nu 1f_{5/2}$ | $-35.702$ | $-35.679$ | $-35.677$
$\nu 2p_{3/2}$ | $-32.094$ | $-32.072$ | $-32.069$
$\nu 2p_{1/2}$ | $-31.344$ | $-31.323$ | $-31.321$
$\nu 1g_{9/2}$ | $-27.863$ | $-27.843$ | $-27.843$
$\nu 1g_{7/2}$ | $-25.791$ | $-25.773$ | $-25.775$
$\nu 2d_{5/2}$ | $-22.277$ | $-22.261$ | $-22.262$
$\nu 2d_{3/2}$ | $-20.843$ | $-20.829$ | $-20.830$
$\nu 3s_{1/2}$ | $-19.959$ | $-19.945$ | $-19.946$
$\nu 1h_{11/2}$ | $-18.533$ | $-18.519$ | $-18.521$
$\nu 1h_{9/2}$ | $-15.306$ | $-15.295$ | $-15.300$
$\nu 2f_{7/2}$ | $-12.672$ | $-12.661$ | $-12.664$
$\nu 2f_{5/2}$ | $-10.622$ | $-10.614$ | $-10.619$
$\nu 3p_{3/2}$ | $-9.859$ | $-9.851$ | $-9.854$
$\nu 3p_{1/2}$ | $-9.063$ | $-9.055$ | $-9.059$
$\nu 1i_{13/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, while only the
first and second terms were considered in this paper because of the LDA.

It was found that the relativistic correction makes ${}^{40}\mathrm{Ca}$ less
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 LeBloas2011Phys.Rev.C84_014310 ;
Roca-Maza2016Phys.Rev.C94_044313 ; Naito2018Phys.Rev.C97_044319 ;
Naito2019Phys.Rev.C99_024309 , 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 was pointed out
also in the context of atomic physics Naito2020J.Phys.B53_215002 . In
addition, since the desired accuracy of the nuclear EDF is $100\,\mathrm{keV}$
order, 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 sub-
dominant to the Breit correction Eides2001Phys.Rep.342_63 . Although the
contribution of vacuum polarization to the total energy is usually tiny, it is
non-negligible to discuss mirror nuclei mass difference or the isobaric analog
energy Auerbach1983Phys.Rep.98_273 ; Roca-Maza2018Phys.Rev.Lett.120_202501 ;
Naito2020Phys.Rev.C101_064311 , whereas the relativistic correction does not
contribute to such properties.

In relativistic Hartree-Fock calculation, the current-current Hartree
interaction is often considered, and the relativistic correction is also
considered to the Fock term exactly or by using LDA. In this work, the
relativistic correction is consistently considered in the Hartree-Fock-Slater
approximation. The magnetic (current-current) interaction and the retardation
appear in the next-leading order of the electromagnetic interaction in the
Lorenz gauge, which corresponds to the Breit correction in the Coulomb gauge
Gorceix1988Phys.Rev.A37_1087 . Hence, the difference between them is which
gauge is introduced in quantum electrodynamics calculation
Sakurai1967AdvancedQuantumMechanics_AddisonWesley ;
Gorceix1988Phys.Rev.A37_1087 ; Lindroth1989Phys.Rev.A39_3794 . Thus, it can be
concluded that by introducing the relativistic correction to non-relativistic
Hartree-Fock calculations as proposed in this work, one can directly compare
electromagnetic energies obtained by non-relativistic and relativistic
calculations.

In this calculation, the spin-orbit and spin-magnetic interactions were not
considered. The former has been discussed in Refs. Roca-
Maza2018Phys.Rev.Lett.120_202501 ; Roca-Maza2018EPJWebConf.194_01002 ;
Naito2020Phys.Rev.C101_064311 , 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 calculation in both atomic and nuclear structure
calculations, which remains as a future perspective.

###### 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 the RIKEN iTHEMS program, the JSPS Grant-in-Aid for JSPS Fellows
under Grant No. 19J20543, and the RIKEN Special Postdoctoral Researchers
Program. The numerical calculations were performed on cluster computers at the
RIKEN iTHEMS program.

## References

  * (1) K. Okamoto, Phys. Lett. 11, 150 (1964)
  * (2) E. Friedman, Phys. Lett. B 35, 543 (1971)
  * (3) S. Shlomo, W.G. Love, Phys. Scr. 26, 280 (1982)
  * (4) T. Suzuki, H. Sagawa, N. Van Giai, Phys. Rev. C 47, R1360 (1993)
  * (5) K. Kaneko, T. Mizusaki, Y. Sun, S. Tazaki, G. de Angelis, Phys. Rev. Lett. 109, 092504 (2012)
  * (6) P. Bączyk, J. Dobaczewski, M. Konieczka, W. Satuła, T. Nakatsukasa, K. Sato, Phys. Lett. B 778, 178 (2018)
  * (7) X. Roca-Maza, G. Colò, H. Sagawa, Phys. Rev. Lett. 120, 202501 (2018)
  * (8) P. Bączyk, W. Satuła, J. Dobaczewski, M. Konieczka, J. Phys. G 46, 03LT01 (2019)
  * (9) B.M. Loc, N. Auerbach, G. Colò, Phys. Rev. C 99, 014311 (2019)
  * (10) J.M. Dong, J.Z. Gu, Y.H. Zhang, W. Zuo, L.J. Wang, Y.A. Litvinov, Y. Sun, Phys. Rev. C 99, 014319 (2019)
  * (11) N. Auerbach, Phys. Rep. 98, 273 (1983)
  * (12) P. Hohenberg, W. Kohn, Phys. Rev. 136, B864 (1964)
  * (13) W. Kohn, L.J. Sham, Phys. Rev. 140, A1133 (1965)
  * (14) W. Kohn, Rev. Mod. Phys. 71, 1253 (1999)
  * (15) S.H. Vosko, L. Wilk, M. Nusair, Can. J. Phys. 58, 1200 (1980)
  * (16) J.P. Perdew, A. Zunger, Phys. Rev. B 23, 5048 (1981)
  * (17) T. Yokota, T. Naito, Phys. Rev. Research 3, L012015 (2021)
  * (18) J.P. Perdew, Y. Wang, Phys. Rev. B 33, 8800 (1986)
  * (19) J.P. Perdew, Y. Wang, Phys. Rev. B 45, 13244 (1992)
  * (20) J.P. Perdew, K. Burke, M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996)
  * (21) J. Drut, R. Furnstahl, L. Platter, Prog. Part. Nucl. Phys. 64, 120 (2010)
  * (22) J. Dobaczewski, J. Phys. G 43, 04LT01 (2016)
  * (23) T. Yokota, K. Yoshida, T. Kunihiro, Phys. Rev. C 99, 024302 (2019)
  * (24) G. Salvioni, J. Dobaczewski, C. Barbieri, G. Carlsson, A. Idini, A. Pastore, J. Phys. G 47, 085107 (2020)
  * (25) G. Accorto, T. Naito, H. Liang, T. Nikšić, D. Vretenar, Phys. Rev. C 103, 044304 (2021)
  * (26) F. Marino, C. Barbieri, A. Carbone, G. Colò, A. Lovato, F. Pederiva, X. Roca-Maza, E. Vigezzi, Phys. Rev. C 104, 024315 (2021)
  * (27) T. Naito, G. Colò, H. Liang, X. Roca-Maza, H. Sagawa (2021), 2107.14436
  * (28) B.A. Li, L.W. Chen, C.M. Ko, Phys. Rep. 464, 113 (2008)
  * (29) P. Danielewicz, J. Lee, Nucl. Phys. A 818, 36 (2009)
  * (30) P. Danielewicz, J. Lee, Nucl. Phys. A 922, 1 (2014)
  * (31) M.V. Stoitsov, J. Dobaczewski, W. Nazarewicz, S. Pittel, D.J. Dean, Phys. Rev. C 68, 054312 (2003)
  * (32) M. Kortelainen, T. Lesinski, J. Moré, W. Nazarewicz, J. Sarich, N. Schunck, M.V. Stoitsov, S. Wild, Phys. Rev. C 82, 024313 (2010)
  * (33) N. Wang, M. Liu, X. Wu, J. Meng, Phys. Lett. B 734, 215 (2014)
  * (34) Z.M. Niu, H.Z. Liang, Phys. Lett. B 778, 48 (2018)
  * (35) Z. Niu, H. Liang, B. Sun, Y. Niu, J. Guo, J. Meng, Sci. Bull. 63, 759 (2018)
  * (36) T. Yamaguchi, H. Koura, Y.A. Litvinov, M. Wang, Prog. Part. Nucl. Phys. 120, 103882 (2021)
  * (37) D. Lunney, J.M. Pearson, C. Thibault, Rev. Mod. Phys. 75, 1021 (2003)
  * (38) T. Naito, R. Akashi, H. Liang, Phys. Rev. C 97, 044319 (2018)
  * (39) T. Naito, X. Roca-Maza, G. Colò, H. Liang, Phys. Rev. C 99, 024309 (2019)
  * (40) T. Naito, X. Roca-Maza, G. Colò, H. Liang, Phys. Rev. C 101, 064311 (2020)
  * (41) D. Vautherin, D.M. Brink, Phys. Rev. C 5, 626 (1972)
  * (42) J. Meng, H. Toki, S. Zhou, S. Zhang, W. Long, L. Geng, Prog. Part. Nucl. Phys. 57, 470 (2006)
  * (43) T. Nikšić, D. Vretenar, P. Ring, Prog. Part. Nucl. Phys. 66, 519 (2011)
  * (44) H. Liang, J. Meng, S.G. Zhou, Phys. Rep. 570, 1 (2015)
  * (45) T. Nakatsukasa, K. Matsuyanagi, M. Matsuo, K. Yabana, Rev. Mod. Phys. 88, 045004 (2016)
  * (46) G. Colò, Adv. Phys.:X 5, 1740061 (2020)
  * (47) S.K. Samaddar, M.M. Majumdar, B.C. Samanta, J.N. De, Nucl. Phys. A 451, 160 (1986)
  * (48) N. Chamel, S. Goriely, J.M. Pearson, Phys. Rev. C 80, 065804 (2009)
  * (49) R.B. Wiringa, V.G.J. Stoks, R. Schiavilla, Phys. Rev. C 51, 38 (1995)
  * (50) T. Saue, ChemPhysChem 12, 3077 (2011)
  * (51) T. Naito, R. Akashi, H. Liang, S. Tsuneyuki, J. Phys. B 53, 215002 (2020)
  * (52) V. Pershina, Nucl. Phys. A 944, 578 (2015)
  * (53) E. Eliav, S. Fritzsche, U. Kaldor, Nucl. Phys. A 944, 518 (2015)
  * (54) P. Schwerdtfeger, L.F. Pašteka, A. Punnett, P.O. Bowman, Nucl. Phys. A 944, 551 (2015)
  * (55) G. Breit, Phys. Rev. 34, 553 (1929)
  * (56) G. Breit, Phys. Rev. 36, 383 (1930)
  * (57) G. Breit, Phys. Rev. 39, 616 (1932)
  * (58) M.I. Eides, H. Grotch, V.A. Shelyuto, Phys. Rep. 342, 63 (2001)
  * (59) A.H. MacDonald, S.H. Vosko, J. Phys. C 12, 2977 (1979)
  * (60) S.D. Kenny, G. Rajagopal, R.J. Needs, W.K. Leung, M.J. Godfrey, A.J. Williamson, W.M.C. Foulkes, Phys. Rev. Lett. 77, 1099 (1996)
  * (61) J. Le Bloas, M.H. Koh, P. Quentin, L. Bonneau, J.I.A. Ithnin, Phys. Rev. C 84, 014310 (2011)
  * (62) X. Roca-Maza, L.G. Cao, G. Colò, H. Sagawa, Phys. Rev. C 94, 044313 (2016)
  * (63) Z.M. Niu, Q. Liu, Y.F. Niu, W.H. Long, J.Y. Guo, Phys. Rev. C 87, 037301 (2013)
  * (64) H.Q. Gu, H. Liang, W.H. Long, N. Van Giai, J. Meng, Phys. Rev. C 87, 041301 (2013)
  * (65) J.M.G. Gomez, J. Sesma, Phys. Lett. B 44, 231 (1973)
  * (66) S. Shlomo, Phys. Lett. B 60, 244 (1976)
  * (67) T. Itoh, Rev. Mod. Phys. 37, 159 (1965)
  * (68) S.D. Kenny, G. Rajagopal, R.J. Needs, Phys. Rev. A 51, 1898 (1995)
  * (69) O. Gorceix, P. Indelicato, Phys. Rev. A 37, 1087 (1988)
  * (70) J.J. Sakurai, _Advanced Quantum Mechanics_ (Addison Wesley, Boston, 1967)
  * (71) E. Lindroth, A.M. Mårtensson-Pendrill, Phys. Rev. A 39, 3794 (1989)
  * (72) M.H.L. Pryce, Proc. Royal Soc. Lond. A 195, 62 (1948)
  * (73) L.L. Foldy, S.A. Wouthuysen, Phys. Rev. 78, 29 (1950)
  * (74) S. Tani, Prog. Theor. Phys. 6, 267 (1951)
  * (75) L.L. Foldy, Phys. Rev. 87, 688 (1952)
  * (76) P.A. Zyla, R.M. Barnett, J. Beringer, O. Dahl, D.A. Dwyer, D.E. Groom, C.J. Lin, K.S. Lugovsky, E. Pianori, D.J. Robinson et al. (Particle Data Group), Prog. Theor. Exp. Phys. 2020, 083C01 (2020)
  * (77) J.C. Slater, Phys. Rev. 81, 385 (1951)
  * (78) X. Roca-Maza, G. Colò, H. Sagawa, EPJ Web Conf. 194, 01002 (2018)
  * (79) E.A. Uehling, Phys. Rev. 48, 55 (1935)
  * (80) L. Wayne Fullerton, G.A. Rinker, Phys. Rev. A 13, 1283 (1976)
  * (81) G. Colò, L. Cao, N. Van Giai, L. Capelli, Comput. Phys. Commun. 184, 142 (2013)
  * (82) X. Roca-Maza, G. Colò, H. Sagawa, Phys. Rev. C 86, 031306 (2012)