Full leading-order nuclear polarization in highly charged ions

The nuclear-polarization corrections to the energy levels of highly charged ions are systematically investigated to leading order in the fine-structure constant. To this end, the notion of effective photon propagators with nuclear-polarization insertions is employed, where the nuclear excitation spectrum is calculated by means of the Hartree-Fock-based random-phase approximation. The effective Skyrme force is used to describe the interaction between nucleons, and the model dependence is analyzed. To leading order, the formalism predicts two contributions given by the effective vacuum-polarization and self-energy diagrams. The existing ambiguity around the vacuum-polarization term is resolved by demonstrating that it is effectively absorbed in the standard finite-nuclear-size correction. The self-energy part is evaluated with the full electromagnetic electron-nucleus interaction taken into account, where the importance of the effects of the nuclear three-currents is emphasized.


I. INTRODUCTION
Ever since the first experimental detection of a nuclearstructure effect (i.e., the field shift) in the hyperfine spectra of the thallium isotopes 203 Tl and 205 Tl in 1931 [1], there have been extraordinary advances in high-precision atomic physics, both from the experimental and the theoretical sides.On the one hand, this has led to a multitude of remarkable applications such as stringent tests of quantum electrodynamics (QED) [2][3][4][5], determination of fundamental physical constants [6][7][8][9], including the search for their variation [10][11][12][13][14], as well as the quest for physics beyond the Standard Model [15][16][17].On the other hand, these applications are now becoming more and more limited by the ability to accurately take into account the increasingly more pronounced nuclear-structure effects.This is especially true for highly charged ions, which hold a special place in the field due to their simplicity and the extreme electromagnetic environments they provide.These advantages, however, come with the price of strong sensitivity to the size [18,19], the shape [20][21][22] and even the internal dynamic structure [23][24][25] of the nucleus.The latter is known as the nuclear-polarization (NP) effect, and it poses one of the biggest challenges in theoretical description of the spectra of highly charged ions.
On the fundamental level, the NP effect stems from the difference between a static nuclear charge distribution, which is a common approximation in atomic calculations, and the true dynamic nature of the nucleus with the intricate motion of individual protons and neutrons within it.Thus, the first calculations of the NP corrections to atomic energy levels of hydrogenlike ions were carried out by means of the ordinary second-order perturbation theory [26].This approach, however, is immediately challenged for relativistic systems by the presence of the negative continuum in the electronic Dirac spectrum.As with any other difficulties plaguing the relativistic single-particle quantum mechanics, this issue could only be resolved within a field-theoretical description [23][24][25], where the vacuum contribution arises naturally.Indeed, as it turned out, a naive extension of the summation incorporates the negative part of the Dirac spectrum with a wrong expression in the energy denominator [27], while omitting this contribution leads to a significant overestimation of the NP effect.
Despite the success and the power of the fieldtheoretical approach for NP, there are still aspects of it that have not received enough attention in the literature.Firstly, to leading order in the fine-structure constant, this formalism leads to two terms corresponding to the effective self-energy and vacuum-polarization diagrams; however, only the former is usually considered, while the latter is omitted entirely, with only a few exceptions in the literature [28][29][30].Additionally, the question of possible double counting in the effective vacuum polarization, which was raised in Ref. [28], needs further clarification.Secondly, in the calculations of the effective self-energy contribution, only the longitudinal, or Coulomb, part of the electron-nucleus interaction is often taken into account, whereas it has been shown in Ref. [31] that the transverse interaction is far from being negligible.Furthermore, from the nuclear side of the computations, the common practice of employing simple phenomenological sum rules for estimating the NP contributions from giant resonances [32], which are dominant in the case of spherical nuclei, does not allow to reliably assess the theoretical uncertainty of the NP corrections.To the best of our knowledge, there is only one fully microscopic study of NP for electronic systems, which was carried out for hydrogenlike 208 Pb 81+ in the framework of the randomphase approximation [31].However, these calculations were performed only for one set of the Migdal-force parameters, and the nuclear model dependence was not explored.
The aim of this paper is to present a complete overview and analysis of the NP effect in highly charged ions to leading order in the fine-structure constant α.That means considering both the self-energy and the vacuumpolarization effective diagrams, taking into account the full electromagnetic electron-nucleus interaction, incorporating a microscopic description of the nucleus as well as investigating the dominant theoretical uncertainty stemming from the nuclear model dependence.For this purpose, we first outline the general field-theoretical formalism and the computational techniques from both atomic and nuclear parts of the calculations.Then we demonstrate for both electronic and muonic systems that the effective vacuum polarization is already contained in the standard finite-nuclear-size effect and thus does not need to be included in the NP correction.Finally, we present our results for the effective self-energy diagram for a series of hydrogenlike ions and a wide range of parametrizations of the effective Skyrme force, which is used in this work for nuclear description.Relative contributions from different nuclear excitation modes, the importance of the transverse electron-nucleus interaction, and the model uncertainties are discussed.
Relativistic system of units (ℏ = c = 1) and Heaviside charge units (α = e 2 /4π, e < 0) are used throughout the paper.Four-vectors (x) are represented by regular typeface, while bold upright letters are used for threevectors (x), whose lengths are denoted by non-bold upright letters (|x| := x).

A. Modified photon propagator
The starting point of any perturbative calculation is the definition of a zeroth-order approximation.In the case of bound atomic systems, we begin with the following Lagrangian density L (0) (omitting the issue of renormalization in this brief overview for simplicity) [33,34]: where ψ(x) is the electron-positron field operator, m e is the electron mass, and γ µ are the Dirac matrices.A µ stat (x) denotes a classical time-averaged electromagnetic four-potential generated by the nucleus, while L free EM is the Lagrangian density for the free electromagnetic field operator Âµ free (x) [33]: where , the second gauge-fixing term corresponds to the Lorenz condition, and the value of the parameter ζ determines the gauge choice.In this subsection, we assume for simplicity ζ = 1.The introduction of the field A µ stat (x) in Eq. ( 1) corresponds to the external field approximation [34] and allows to describe bound states, which are of central interest in atomic physics.The usual interaction term is given by [33] and it is responsible for all kinds of QED corrections, which can be calculated in the so-called Furry picture [35] while treating the nucleus simply as a static charge distribution.
In order to treat the NP effect on the same fieldtheoretical footing, the total nuclear four-current density operator Ĵµ N is introduced as the following sum [23]: with the classical static part J µ N, stat corresponding to the average over the nuclear ground state and the fluctuating part Ĵµ N, fluc describing the intrinsic nuclear dynamics.In the same way as J µ N, stat is taken into account by introducing the corresponding classical field A µ stat in Eq. ( 1), a second-quantized photon field Âµ fluc can be associated with the fluctuating current Ĵµ N, fluc .In this view, a bound atomic electron interacts with the total electromagnetic field where the total quantum radiation field Âµ rad is defined as the sum of the free photon field Âµ free and the fluctuating part Âµ fluc generated by Ĵµ N, fluc .The latter is described by the following equation of motion: As a consequence, in order to keep utilizing the standard machinery of Wick's theorem and Feynman diagrams, one is led to the modified photon propagator instead of the usual free photon propagator given by where |0⟩ denotes the "vacuum" state in the presence of the external field A µ stat (x), which corresponds to the nucleus being in its ground state.Since the mixed terms arising from Eq. ( 7) vanish, only the ⟨0|T [ Âfluc µ (x) Âfluc ν (x ′ )]|0⟩ term needs to be evaluated.By using the fact that the free photon propagator is the Green's function of the free equation of motion: and performing integration by parts with vanishing boundary terms, one can show that It is important to note that the derivatives in the last line of Eq. ( 10) act not only on the fields Âµ fluc but also on the θ-functions from the time-ordered product: producing an additional term that we denote as iS ξζ N (x 1 , x 2 ): while the resulting two-point current correlation function defines the so-called NP tensor Going back to Eq. ( 7), the expression for the modified photon propagator can be written as defining the NP correction D NP µν (x, x ′ ) to the free photon propagator as follows: In the original paper [23] by Plunien et al., where the modified photon propagator (7) was first introduced, only the NP tensor appeared in the final expression.However, it was pointed out back in 1961 by Kenneth Johnson that a time-ordered product of two currents is in general not a covariant function [36].Thus, the role of the iS ξζ N term is to maintain the Lorentz covariance of the vacuum expectation value ⟨0|T [ Âfluc µ (x) Âfluc ν (x ′ )]|0⟩ and the modified photon propagator as a whole.In addition, Lowell Brown demonstrated a connection between the requirement of gauge invariance and the restored Lorenz covariance of a properly defined two-point current correlation operator [37].This requirement implies that which leads to where the second term is equal to zero due to the continuity equation of nuclear charge conservation.While the equal-time commutator in Eq. ( 17) vanishes for ζ = 0, it was shown by Julian Schwinger from fundamental principles of quantum field theory that charge and current (ζ = 1, 2, 3) densities cannot commute at a common time [38].It follows from Eq. ( 17) that these nonvanishing commutators, known as the Schwinger terms, must be cancelled by the divergence of iS ξζ N , if gauge invariance is to be satisfied.The contribution iS ξζ N is often called the "seagull" or "catastrophic" term, and this kind of cancellation is in fact a very general result in currentalgebra theories [39].
It is clear that the expression for iS ξζ N depends on a specific definition of Ĵµ N .For example, it can be shown that in the case of the non-relativistic nuclear chargecurrent density operators the seagull term takes on the following form [31]: where ρN (x 1 ) := Ĵ0 N (0, x 1 ), M p is the proton mass, and δ ξζ is the Kronecker delta extended to four dimensions with δ 00 = 0.

B. Nuclear-polarization insertion
As one can see from Eq. ( 14), every photon line in ordinary QED diagrams receives the NP correction in the form of Eq. ( 15).The part [Π ξζ N + S ξζ N ] in between the two free photon propagators is called the NP insertion, and it contains all the information about the intrinsic nuclear dynamics.In order to implement this correction in practical calculations, we employ the two-time Green's function method from Ref. [40].For this purpose, we need a Fourier-transformed version of D NP µν with respect to the time variables.First, by writing the time evolution of the Heisenberg operators Ĵµ N, fluc (t, x), inserting a complete set of nuclear excitations |λ⟩ in Eq. ( 13) (the ground state does not contribute due to the definition in Eq. ( 4)) and using the integral representation of the θ-function, it is easy to show that the NP tensor is homogeneous in time [24]: with where Ĵµ N (x) := Ĵµ N (0, x), and ω λ = E λ − E 0 are the nuclear excitation energies.Similarly, the Fouriertransformed version of the seagull term (18) reads Then, by using Eq. ( 19) together with one readily obtains the desired equivalent of Eq. ( 15): Thus, we can supplement the set of the Feynman rules in Ref. [40] with an additional one shown in Fig. 1.The NP insertion is represented by a shaded circle, and, by analogy with the rule for an internal photon line, the energy variable ω has to be integrated over.
Feynman rule for evaluating the NP correction to the photon propagator.
To leading order in the fine-structure constant, the NP insertion can be applied to both the self-energy (SE) and the vacuum-polarization (VP) corrections resulting in the two effective diagrams presented in Fig. 2. Analogously to the ordinary SE and VP effects, the following general expressions for the energy shifts of an electronic level |a⟩ due to the SE-NP and VP-NP diagrams are obtained as where the dressed electron propagator is constructed from a complete set of eigenstates and eigenvalues of the stationary single-particle Dirac equation: with V (x) = eA 0 (x) denoting the potential energy corresponding to the electrostatic potential of an atomic nucleus, β = γ 0 , and The SE-NP contribution has received much more attention in the literature because of its correspondence to the usual NP correction calculated earlier in secondorder perturbation theory.In the language of field theory, Eq. ( 24) can also be interpreted as a two-photon exchange between the bound electron and the nucleus, with the ladder and the cross diagrams in Fig. 2(a) corresponding to the first and the second terms in the expression (20) for the NP tensor, respectively.The seagull term (21) can be represented as a coupling of the electromagnetic currents to the nucleus at the same point, and its physical significance is to ensure gauge invariance of the calculated SE-NP corrections.In this work, we employ the formulas for the ladder, cross, and seagull terms in the momentum representation that can be found in Ref. [31].As for the VP-NP contribution, we present some more details on it in the following subsection.

C. Effective vacuum polarization
The trace part of Eq. ( 25) can be rewritten as [41] −ieTr where the current operator of the Dirac field is defined as ĵν (x ′ ) = e 2 [ ψ(x ′ ), γ ν ψ(x ′ )], and the limit x = x ′ is approached symmetrically from t = t ′ − 0 and t = t ′ + 0. To first order in Zα, this induced current is related to  to the photon propagator (wavy line).The former can be interpreted as a two-photon exchange between the bound atomic electron (double line) and the nucleus (single solid line), which includes the ladder, the cross and the seagull diagrams; while the latter can be understood as the interaction of the induced vacuum polarization with the bound electron via additional virtual nuclear excitations.
the static external current of the nucleus in momentum space as follows [41]: with the renormalized Uehling polarization function Since we consider electrostatic nuclear potentials, the only non-vanishing component in Eq. ( 29) is the nuclear ground-state charge density J0 N, stat (k) := ρN (k).As a result, only the component D NP 00 (0, x, x ′ ) contributes to the VP-NP energy shift: where ρ VP (x ′ ) := ⟨ ĵ0 (x ′ )⟩ (1) is the induced vacuumpolarization density.In this case, the NP correction to the photon propagator is obtained most conveniently in the Coulomb gauge: Moreover, in the case of spherically symmetric nuclear charge distributions ρN , the induced vacuum-polarization density ρ VP becomes spherical, too.This fact leads to a further simplification of Eqs.(31) and (32) when one takes into account the following expression for the nuclear matrix elements entering the NP tensor [42]: where ρ tr λ (x) is the radial charge transition density for a nuclear excited state |λ⟩ with the angular quantum numbers J and M .Thus, it is easy to see that only the monopole part of the Coulomb interaction in Eq. ( 32) survives, which also means that only the monopole nuclear excitations (the so-called breathing modes) contribute to the VP-NP correction.In this case, the NP tensor takes on an especially simple form: such that where Finally, we also recall that in central potentials V (x) the electron wave function factorizes into the radial and angular parts: where n is the principal quantum number, κ is the relativistic angular momentum number, m is the total magnetic number, and Ω ±κm (Ω x ) are the spherical spinors.
After collecting everything together, we obtain the following expression for the VP-NP energy shift: with We note that the VP-NP correction can also be expressed in terms of the reduced transition probabilities B(E0), which was done in Ref. [28].

D. Computational techniques
The calculations of the NP corrections require input from both atomic and nuclear physics.For the atomic part, we solve the stationary Dirac equation ( 27) numerically by confining the system to a spherical cavity of a large radius and expanding the radial part of the electron wavefunction in terms of B -splines within the dualkinetic-balance approach [43].The resulting generalized matrix eigenvalue equations are then readily solved by means of the LAPACK library.This approach is especially useful in the case of the SE-NP correction since it allows to reduce the infinite sum over the bound states and the integrals over the positive and negative continua to finite sums with no remainders to evaluate.At the same time, the convergence of the results is readily controlled by varying the size of the cavity and the number of B -splines used.As for the central nuclear potential V (x) corresponding to the zeroth approximation, it is sufficient to construct it from the simple two-parameter Fermi charge distribution where we use the standard value of the diffuseness parameter a = 2.3/[4 ln(3)] fm and adjust the half-density radius c such that the tabulated values [44] of the rootmean-square (RMS) radii are reproduced.
For the nuclear part of the computations, we employ the skyrme_rpa program [42] in order to obtain the entire nuclear excitation spectrum.In the first step, the selfconsistent Hartree-Fock mean field is built assuming the effective Skyrme-type interaction between the nucleons.Then, based on this mean field, the excited nuclear states are calculated in the framework of the random-phase approximation (RPA).Under the spherical symmetry assumption, the RPA equations are solved separately for given angular momentum and parity J π .In our calculations we take into account the 0 + , 1 − , 2 + , 3 − , 4 + , 5 − , and 1 + excitation modes.The completeness of the obtained spectra is controlled by the degree of fulfilment of the energy-weighted sum rule.In addition, we extended the code to include the transition densities for the nuclear three-currents, which are needed for the transverse part of the SE-NP corrections.Non-relativistic charge-current operators are used in the calculations of the nuclear matrix elements, and the corresponding expressions can be found in Ref. [45].

A. Effective vacuum polarization
In this section we present the numerical results of our NP calculations.First, we address the VP-NP diagram and its physical meaning, as the question was raised in Ref. [28] on whether this contribution might already be contained in the ordinary finite-nuclear-size (FNS) effect.Based on a scaling argument, the authors of Ref. [28] suggested that there is a part of the VP-NP correction that is qualitatively distinct from the standard FNS contribution; however, they did not provide a clear procedure of extracting it.At the same time, the diagrammatic representation in Fig. 2(b) suggests that the VP-NP effect corresponds to a correction to the Coulomb potential of the nucleus due to the interaction between the induced vacuum polarization and nuclear degrees of freedom.Therefore, in this work, we approach this issue by additionally considering the effect of the induced vacuum polarization on the nuclear ground state.
To this end, we modify the skyrme_rpa program by introducing the Uehling potential [46] V into the mean field generated by the protons with the density distribution ρ p (r ′ ) and estimate the resulting change in the RMS charge radius.The obtained enlargement of the nucleus can then be translated into the corresponding increase of the FNS effect (taken with respect to the tabulated RMS radius), which we denote as ∆FNS.
The VP-NP corrections calculated from the monopole excitation spectrum as described in Subsection II C also decrease the atomic binding energies and can be directly compared to the ∆FNS values.We present such a comparison in Table I for the example of the 208 Pb nucleus and three parametrizations of the Skyrme force, namely SKX, KDE0, and SkP [47][48][49].In order to examine the scaling behaviour, both the electron and the muon are considered as a bound atomic particle, and the results are given for the 1s 1/2 , 2s 1/2 , and 2p 1/2 states.Considering both leptons is of importance, given that the nuclear charge radii are determined experimentally from muonic spectroscopy [50] and subsequently used as input parameters in the description of electronic systems.Therefore, it is essential to treat these two cases consistently.
The extremely close agreement between the ∆FNS and VP-NP values, independent of the state and the bound particle considered, is clear evidence that the VP-NP effect simply corresponds to the change in the nuclear charge radius due to the induced vacuum polarization.Put differently, if one calculates the VP-NP correction and expresses the corresponding energy shift in terms of a shift of the RMS radius, the resulting radius shift will be the same for different states and also for electronic and muonic atoms.We note that a similar level of agreement was also found for other nuclei considered in this work.Since the experimentally obtained nuclear charge radii are inseparable from such QED corrections and thus already contain them, we conclude that the VP-NP contribution does not need to be taken into account as part of the NP effect.

B. Effective self-energy in Coulomb approximation
Before turning to the full SE-NP correction, let us first consider the approximation commonly used in the literature.Based on the argument that the velocities associated with nuclear dynamics are mainly non-relativistic, the contributions from the nuclear three-currents ĴN are often neglected, meaning that only the Π 00 N component of the NP tensor (the longitudinal, or Coulomb, part) is taken into account.We refer to this framework as the "Coulomb approximation".Moreover, the calculations can be simplified even further by expressing the SE-NP energy shift in terms of the experimentally measurable nuclear transition probabilities B(EL; L → 0) λ and excitation energies ω λ , resulting in the most widely used formula for the SE-NP correction [27]: where the radial functions F L are often assumed to be with R 0 being the radius of the nucleus as a homogeneously charged sphere.However, the experimental In this approach, the giant resonances are assumed to be concentrated in a single state for each multipolarity and isospin.Since the skyrme_rpa program also allows to compute the B(EL) λ values, it is of interest to compare the Coulomb SE-NP corrections calculated from such a microscopic nuclear theory with the results obtained from the experimental and EWSR nuclear parameters.Such a comparison for the ground state of the 208 Pb 81+ ion is presented in Table II.Eleven different Skyrme parametrizations [47][48][49][51][52][53][54][55][56][57][58] (nine of which have been shown to cover a wide range in the parameter space [59]) have been used in the microscopic calculations, while the last row has been calculated using the same B(EL) λ and ω λ values as in Ref. [27].The results turn out to be quite stable with respect to a nuclear model and are in remarkably good agreement with the extremely simple EWSR estimations.To appreciate this fact, we note that, unlike concentrating the giant resonances in just a few "states" in the EWSR approach, the microscopic RPA calculations produce numerous excitations, e.g., around 1500 for the 3 − mode alone.It is also immediately apparent that most of the nuclear model dependence comes from the 1 − contribution, which is the largest one due to the longer range of the dipole NP potential and which is almost entirely dominated by the giant dipole resonances.The use of Eq. ( 42) involves certain assumptions about the radial dependence of the Coulomb part of the NP Table III.SE-NP contributions (in meV) with different multipolarities for the 1s 1/2 state of 208 Pb 81+ obtained with the SLy5 parametrization of the Skyrme force.Two ways of implementing the Coulomb approximation are compared to the full calculation including the transverse electron-nucleus interaction.See Subsections III B and III C for more details.
[SLy5] tensor.It is not uncommon to employ the same radial functions (44), which correspond to harmonic surface vibrations, for all types of nuclear excitations, except for the special case of the breathing modes, where Eq. ( 43) is used instead.This simplification is justified by the small overlap between electronic wavefunctions and a nucleus such that the details of this radial dependence are of minor importance.Therefore, it is also of interest to test this approximation by using Π 00 N constructed directly from the RPA charge transition densities ρ tr λ (given by Eq. ( 33)) with a unique radial dependence for each nuclear excitation.We compare the {B(EL) λ +F L } and the ρ tr λ results for the SLy5 parametrization in the first two rows of Table III, again using the 1s 1/2 state of 208 Pb 81+ as an example.One can see that the more accurate and detailed treatment of the Coulomb SE-NP correction by means of ρ tr λ leads to slightly larger values, with most of the difference coming from the dominant 1 − contribution.Nevertheless, the overall agreement between the EWSR estimation from the last row of Table II and the microscopic calculation from the second row of Table III is noteworthy, given how much simpler the former is compared to the latter.

C. Full effective self-energy
The assumption of including only the Π 00 N component of the NP tensor turns out to be not as accurate as was initially thought.It was argued in Refs.[31,45] that the transverse electron-nucleus interaction is not negligible due to the interference term with the Coulomb contribution.This fact is illustrated in the last row of Table III, where we present the SE-NP corrections with the full electromagnetic electron-nucleus interaction taken into account by means of the formulas from Ref. [31].While the 0 + , 2 + , and 3 − contributions are not influenced much by the inclusion of the transverse part, the dominant 1 − correction is increased by as much as 35%.
In order to gain more insight into this difference, we plot in Fig. 3 the sum of the ladder and the cross contributions as a function of the intermediate electron energy for fixed 1 − (a) and 2 + (b) nuclear excitations with the energies of 12.8 MeV and 12.2 MeV, respectively.Both of these nuclear states represent the largest contributions from the giant resonances.The dashed lines are used for the Coulomb approximation, while the solid lines correspond to the full electron-nucleus interaction.Different colors represent the two allowed values of the relativistic angular quantum number κ of the intermediate electronic states in both cases.First, it is worth pointing out the large cancellations that take place between the contributions from the positive-and negative-energy parts of the Dirac spectrum, which makes the precise evaluation of the total SE-NP correction a challenging task.Then, by comparing Fig. 3(a) and Fig. 3(b) one notices a striking difference in how the inclusion of the transverse electronnucleus interaction affects these two cases.While it does not create much of a difference for the 2 + mode, sharp peaks appear in the low-energy region for the 1 − one.This unique behaviour is due to the fact that only in the dipole case there are transverse form factors containing the non-vanishing spherical Bessel function j 0 at zero momentum transfer [31].Here we note that the transverse interaction also leads to a small magnetic 1 + correction of approximately +1.5 meV for the ground state of 208 Pb 81+ , thus contributing with the opposite sign compared to the electric excitation modes.
We now present in Table IV the results of our full SE-NP calculations including the 0 + , 1 − , 2 + , 3 − , 4 + , 5 − , and 1 + excitation modes.The energy shifts are given for five hydrogenlike ions, namely 208 Pb 81+ , 120 Sn 49+ , 90 Zr 39+ , 60 Ni 27+ , and 40 Ca 19+ , in the 1s 1/2 , 2s 1/2 , and 2p 1/2 states.Again, in order to explore the nuclear model dependence, the values have been obtained using eleven different Skyrme parametrizations, except for the case of 60 Ni 27+ where the LNS calculation breaks down.In general, it can be seen that the SE-NP correction exhibits a rather steep decrease for lower charge numbers Z and loses about an order of magnitude when going to the next principal or angular quantum number.
A conspicuous feature of Table IV is that the inclusion of the transverse electron-nucleus interaction makes the SE-NP corrections to the 2p 1/2 state positive.As with the 1s 1/2 correction discussed above, the transverse interaction mainly affects the dominant 1 − contribution also in this case, which would otherwise be negative in the Coulomb approximation.Moreover, the sum of the ladder and the cross diagrams for the 2p 1/2 state is still negative in the full calculation, and it is the seagull term that plays the crucial role in making the total 1 − contribution positive.This behaviour has been confirmed by performing the calculations in both Coulomb and Feynman gauges, which lead to quite different contributions from the individual diagrams but eventually produce the same total value.To the best of our knowledge, this unusual result has not been reported in the literature so far, and its exact physical interpretation is rather elusive.
Finally, when it comes to high-precision tests of QED, one of the most important aspects of the calculations is to quantify the theoretical uncertainty.In the case of   IV the NRAPR parametrization represents significant outliers, while such a behaviour is not observed for muonic atoms [59].Even though there is no strong physical argument either in favor or against the reliability of a particular model, it seems nevertheless likely that the NRAPR parametrization tends to overestimate the magnitudes of the NP effect in electronic systems, thereby providing conservative upper bounds.

IV. CONCLUSIONS AND OUTLOOK
In this work we have presented the field-theoretical formalism for the NP effect as well as our full leading-order computational results involving a detailed microscopic nuclear description.In general, the formalism predicts two distinct NP contributions of order α 2 .One of them, which is represented by the effective vacuum-polarization diagram, has been somewhat of a mystery in the literature.On the one hand, it was suggested in Ref. [28] that a part of this correction could be distinguished from the ordinary FNS effect due to its different scaling with the nuclear charge radius.On the other hand, it still remained unclear how to extract such a contribution, and this VP-NP correction has been largely ignored in NP calculations, with the only exceptions of Refs.[29,30].It is especially important to resolve this ambiguity, as inclusion of the VP-NP term would modify the total NP corrections significantly, e.g., by canceling more than half of the usual SE-NP contribution for the 1s 1/2 state of 208 Pb 81+ .By considering the effect of the induced vacuum polarization on the nuclear size, we have demonstrated that the VP-NP contribution is already contained in the standard FNS correction.Importantly, we have checked that this result holds for both electronic and muonic systems, as the latter are used to experimentally measure nuclear charge radii.Consequently, we conclude that the VP-NP term should indeed be omitted in order to avoid double counting.
The other NP contribution of order α 2 is given by the effective self-energy diagram, which can also be interpreted as the two-photon exchange between the bound electron and the nucleus.This correction is often evaluated by neglecting the effects of the nuclear threecurrents, resulting in what we call the Coulomb approximation.In addition, the crucial contributions from the giant nuclear resonances are usually estimated by means of phenomenological EWSR.While we have found this approach to be in remarkably good agreement with our more detailed microscopic calculations, we have also demonstrated the importance of including the transverse part of the electron-nucleus interaction, which plays an essential role for an accurate evaluation of the dipole part of the SE-NP correction.For instance, the transverse interaction in 208 Pb 81+ increases the dipole contribution to the 1s 1/2 energy shift by as much as 35% and may even lead to a sign change from negative to positive in the case of the 2p 1/2 state as suggested by our results.Therefore, we emphasize that the commonly used Coulomb approximation as well as the corresponding effective potentials [60][61][62] significantly underestimate the magnitude of the SE-NP correction, and full calculations are needed in order to deliver satisfactory accuracy.
Finally, we have extensively analyzed the nuclear model dependence of our NP calculations, allowing us to reliably estimate the dominant theoretical uncertainty.The smallest spread of the results for eleven different models has been observed for the heaviest 208 Pb 81+ ion, with the uncertainties being below 10% for the 1s 1/2 and 2s 1/2 states.From one point of view, this level of precision for nuclear-structure effects makes 208 Pb 81+ an appealing platform for stringent tests of QED in extreme electromagnetic fields.Alternatively, a set of different predictions from various nuclear models in combination with high enough experimental precision may even potentially offer an opportunity to test and improve our understanding of the inner workings of the nucleus itself.

Figure 2 .
Figure 2. Effective self-energy (a) and vacuum-polarization (b) diagrams as a result of applying the NP insertion (shaded circle)to the photon propagator (wavy line).The former can be interpreted as a two-photon exchange between the bound atomic electron (double line) and the nucleus (single solid line), which includes the ladder, the cross and the seagull diagrams; while the latter can be understood as the interaction of the induced vacuum polarization with the bound electron via additional virtual nuclear excitations.

Figure 3 .
Figure 3. Contributions from different intermediate electronic states to the sum of the ladder and the cross parts of the SE-NP correction to the 1s 1/2 state of 208 Pb 81+ for single fixed 1 − (a) and 2 + (b) nuclear excitations with the energies of 12.8 MeV and 12.2 MeV, respectively.Dashed lines represent the Coulomb approximation, while the solid lines correspond to the full electron-nucleus interaction including the transverse part.Different colors are used for the two allowed values of the relativistic angular number κ of the electronic states in both cases.The nuclear part of the calculations is performed using the SLy5 Skyrme parametrization.

Table I .
Comparison between the VP-NP corrections and the increase of the FNS effects stemming from the enlargement of the nuclear RMS charge radius (first row) due to introduction of the Uehling potential into the mean field of the protons.The values are given for the 208 Pb nucleus with a single bound electron or muon in the 1s 1/2 , 2s 1/2 , and 2p 1/2 states.Three Skyrme parametrizations SKX, KDE0, and SkP are chosen as examples.

Table II .
[42]l dependence of contributions with different multipolarities to the Coulomb SE-NP shift (in meV) of the 1s 1/2 level of hydrogenlike 208 Pb 81+ .The corrections are obtained by means of the formulas from Ref.[27]using the same experimental and EWSR nuclear parameters (last row) as well as utilizing the B(EL) λ and ω λ values obtained from the skyrme_rpa program[42]with eleven different Skyrme parametrizations.
B(EL) λ and ω λ values are available only for low-lying nuclear states, while the crucial contributions from the giant resonances have to be estimated by resorting to phenomenological energy-weighted sum rules (EWSR).