PHYSICAL REVIEW A 89, 042121 (2014)

Coordinate-space approach to vacuum polarization

Paul Indelicato,1,* Peter J. Mohr,2,† and J. Sapirstein3,‡

1Ecole Normale Sup  ́ erieure, CNRS, Universit  ́ e Pierre et Marie Curie - Paris 6, 4 Place Jussieu, 75252 Paris Cedex 05, France  ́ 2National Institute of Standards and Technology, Gaithersburg, Maryland 20899-8420, USA 3University of Notre Dame, Notre Dame, Indiana 46556, USA

(Received 29 January 2014; published 28 April 2014)

The vacuum-polarization correction for bound electrons or muons is examined. The objective is to formulate a

framework for calculating the correction from bound-state quantum electrodynamics entirely in coordinate space,

including the Uehling potential, which is usually isolated and treated separately. Pauli-Villars regularization is

applied to the coordinate-space calculation and the most singular terms are shown to be eliminated, leaving

the physical correction after charge renormalization. The conventional derivation of the Uehling potential in

momentum space is reviewed and compared to the coordinate-space derivation.

DOI: 10.1103/PhysRevA.89.042121 PACS number(s): 03.65.Pm, 31.30.J−

I. INTRODUCTION

In atomic hydrogen, an electron bound to a proton, the

largest radiative correction to the energy levels of S states

is the self-energy correction, which results from the electron

emitting and reabsorbing a photon. This is in contrast to

muonic hydrogen, a muon bound to a proton, in which the

vacuum-polarization correction dominates. The reason for the

difference is that for the muon, the lighter electron-positron

pair that can occur in a vacuum-polarization loop produces a

large effect, while for the electron, there is no lighter particle

to produce a correspondingly large correction.

There has been a renewed interest in the effects of vacuum

polarization due to the recent measurement of the Lamb shift

in muonic hydrogen [1–3]. As mentioned above, in this atom,

the Lamb shift is predominantly the effect of electron vacuum

polarization. However, the radius of the proton, which is

deduced by comparison of the measured transition frequencies

to the theoretical predictions, differs from the value obtained

from spectroscopy of hydrogen and deuterium and electron

scattering experiments by 7σ [4]. In view of this discrepancy,

a thorough review of the theory is warranted, and a number of

such investigations have been carried out, with recent reviews

given in Refs. [5–12].

In this paper, we reexamine the theory of vacuum polarization. The objective is to provide a formulation for

the calculation of the vacuum-polarization effect entirely in

coordinate space for any spherically symmetric binding field,

in parallel with an earlier analysis of the self-energy correction

in coordinate space [13]. This provides a framework for a direct

numerical evaluation of the correction for strong Coulomb

fields and for non-Coulombic binding fields for which a

perturbation expansion is not feasible. Vacuum polarization is

particularly problematic, because it contains the most severe

divergences of all the bound-state corrections.

The calculation of vacuum polarization done by Wichmann and Kroll [14] is also based on a coordinate-space

formulation, but it requires the explicit analytic expression

*paul.indelicato@lkb.ens.fr

†

mohr@nist.gov

‡

jsapirst@nd.edu

for the Green function for a point-charge nucleus. The

present calculation is not based on explicit solutions, so

it provides a framework for more general potentials. Here,

the vacuum polarization is examined in coordinate space

using Pauli-Villars regularization for an arbitrary spherically

symmetric charge distribution in the nucleus. An alternative

method of differential regularization has also been used to

treat divergences in coordinate space [15]. We reprise the

conventional derivation of the Uehling potential in momentum

space in order to compare it to the coordinate-space version.

II. VACUUM POLARIZATION

In bound-state quantum electrodynamics (QED), the

second-order electron vacuum-polarization correction for an

electron or muon in an external potential V (x) in the Feynman

gauge is given, in units where {c}={}={me} = 1 [16], by

(see, for example, [17,18])

E(2)

VP = 4πiα



d(t2 − t1)



dx2



dx1 DF(x2 − x1)

× Tr[γμSF(x2,x2)]φn(x1)γ μφn(x1). (1)

In Eq. (1), φn(x) is a four-component wave function given by

φn(x) = φn(x) e−iEnt, where φn(x) is an eigenfunction of the

Dirac equation

[−iα · ∇ + V (x) + βm − En] φn(x) = 0. (2)

Here m is the bound lepton mass, the Dirac matrices are

β =



I 0

0 −I



, α =



0 σ

σ 0



, (3)

with Pauli matrices σ, and the γ matrices are γ 0 = β and

γ i = βαi

, i = 1,2,3. The photon propagation function is given

by

DF(x2 − x1) = − i

(2π)4



d4

q e−iq(x2−x1)

q2 + iδ

= i

8π2

∞

−∞

dq0 e−iq0(t2−t1) e−b|x2−x1|

|x2 − x1| , (4)

1050-2947/2014/89(4)/042121(10) 042121-1 ©2014 American Physical Society
PAUL INDELICATO, PETER J. MOHR, AND J. SAPIRSTEIN PHYSICAL REVIEW A 89, 042121 (2014)

with b = −i(q02 + iδ)

1/2, Re(b) > 0. Integration over the

time difference yields

E(2)

VP = −α



dx2



dx1

1

|x2 − x1|

× Tr[γμSF(x2,x2)]φn(x1)γ μφn(x1). (5)

Evidently, the quantity

eTr[γ0SF(x2,x2)] (6)

takes the role of the vacuum-polarization charge density. The

electron propagation function can be written as

SF(x2,x1) = 1

2πi

∞

−∞

dz G(x2,x1,z(1 + iδ))γ 0e−iz(t2−t1),

(7)

where the Green function is the solution of the equation

[−iα · ∇2 + V (x2) + βm − z]G(x2,x1,z) = δ(x2 − x1),

(8)

so that

E(2)

VP = iα

2π

∞

−∞

dz dx2



dx1

1

|x2 − x1|

×{Tr[G(x2,x2,z(1 + iδ))]φ†

n(x1) φn(x1)

− Tr[α G(x2,x2,z(1 + iδ))] · φ†

n(x1) α φn(x1)}. (9)

For a spherically symmetric potential V (x), the second term

vanishes [14], but only after a formal cancellation of infinite

terms (see Appendix A).

III. REGULARIZATION

The expression in Eq. (9) must be modified by regularization in order to produce a valid function for either analytical

or numerical evaluation. In particular, Tr[G(x2,x2,z)] is

meaningless, because for x2 ≈ x1 the trace of the Green

function has a limiting form given by

Tr [G (x2,x1,z)] = z

π |x2 − x1|

+··· (10)

and is undefined for equal coordinates. In addition, the

integration over z is divergent.

To obtain a physical prediction from Eq. (9), we carry out

Pauli-Villars regularization by replacing the Green function

with the regulated Green function given by [19]

GR (x2,x1,z) = 

i

CiGi (x2,x1,z) , (11)

where Gi is the Green function for a Dirac particle with mass

mi, the solution of

[−iα · ∇2 + V (x2) + βmi − z] Gi(x2,x1,z) = δ(x2 − x1),

(12)

and the Ci are functions of the masses. The leading term has

C0 = 1; m0 = m. (13)

The coefficients Ci for i > 0 are chosen to eliminate the

divergent terms from the level shift and will be seen to be

quotients of polynomials in the mi. Then for x2 = x1 we have

the limit

lim

mi→∞

i>0

GR (x2,x1,z) = G (x2,x1,z) , (14)

where G is the unregulated bound-electron Green function.

The level shift is calculated by replacing G by GR in Eq. (9)

and eventually taking the limit mi → ∞ for i > 0 after

renormalization.

IV. EXPANSION IN V

To separate the singularities from the finite physical contribution, it is useful to employ the power series in V for the Green

function, an approach that has been extensively discussed in

connection with the evaluation of QED corrections for bound

states (see, for example, [14,20,21]). The expansion is

G(x2,x1,z) = F(x2,x1,z) −



dw F(x2,w,z)

× V (w)F(w,x1,z) +··· , (15)

where F is the free Green function, given by

F(x2,x1,z) = [−iα · ∇2 + βm + z] e−c|x2−x1|

4π |x2 − x1|

, (16)

with c = √m2 − z2, Re(c) > 0. It satisfies the equation

[−iα · ∇2 + βm − z] F(x2,x1,z) = δ(x2 − x1). (17)

The expansion is illustrated in Fig. 1.

In this context, Furry’s theorem is seen by the following

consideration. In the expansion of the electron Green function

in powers of V , for the term with n powers of V , there will be

altogether n + 1 vertices and free Green functions F. Thus the

expression will consist of the trace of the sum of products of i

α matrices, jβ matrices, and k powers of z, where i + j + k =

n + 1. The trace will vanish unless both i and j are even, and

the integration over z will vanish unless k is even. Hence only

terms with n + 1 even may be nonzero, which is just Furry’s

theorem.

= + + + + ...

(a) (b) (c) (d) (e)

FIG. 1. Expansion of the vacuum-polarization correction in powers of the binding potential V . Diagram (a) represents the complete

vacuum polarization of order α. Diagrams (b) through (e) depict the

zero-, one-, two-, and three-potential contributions. Only diagrams (c)

and (e) and higher-order diagrams with an even number of verticies

in the loop are nonzero according to Furry’s theorem.

042121-2
COORDINATE-SPACE APPROACH TO VACUUM POLARIZATION PHYSICAL REVIEW A 89, 042121 (2014)

A. Zero potential

The leading term in the expansion is the free Green function,

with the zero-component trace given by

Tr [F(x2,x1,z)] = z e−c|x2−x1|

π |x2 − x1|

. (18)

The contribution of this term to the vacuum-polarization

correction formally vanishes because of Furry’s theorem.

However, the term is problematic when x2 = x1, which

is remedied by Pauli-Villars regularization. In view of the

analyticity of F(x2,x1,z) as a function of z and the fact

that its branch points are located at z = ±m(1 − iδ), we

can modify the contour of integration to be a straight line

along the imaginary axis plus contributions from two quarter

circles in the first and third quadrants of the complex z plane.

The contribution from the two quarter circles will vanish as

their radii increase, provided the integrand of the regulated

expression falls off for large |z| faster than 1/|z| if x2 = x1.

Letting z = iu, we have

ci =



m2

i − z2 = √

u2 + m2

i

2

√

u2 − m4

i

8u2

√

u2

+··· , (19)

and the leading terms in the expansion in 1/u are

e−ci |x2−x1|

|x2 − x1| = e−

√

u2 |x2−x1|

|x2 − x1| − m2

i e−

√

u2 |x2−x1|

2

√

u2

+··· . (20)

The two terms on the right-hand side of Eq. (20) lead

to divergent contributions, so we eliminate them with two

auxiliary mass propagators that fulfill the conditions

1 + C1 + C2 = 0, (21)

m2

0 + C1m21 + C2m22 = 0, (22)

satisfied by

C1 = m2

0 − m22

m2

2 − m21

, (23)

C2 = m2

1 − m20

m2

2 − m21

. (24)

The regularization is implemented by writing

FR(x2,x1,z) = 

2

i=0

CiFi(x2,x1,z), (25)

where Fi is given by Eq. (16) with m replaced by mi. We thus

have

Tr [FR(x2,x1,z)] = 

2

i=0

Ci

iue−ci |x2−x1|

π |x2 − x1|

, (26)

and

Tr[FR(x2,x2,z)] = lim

|x2−x1|→0

Tr [FR(x2,x1,z)]

= −iu

π



2

i=0

Cici

= i

8πu√u2





2

i=0

Cim4

i + O

 1

u2



. (27)

Evidently, the contribution from the quarter circles vanishes,

and since the branches of the square root are specified to give

√

u2 = |u| for real values of u, the integrand is an odd function

of u, and we have

∞

−∞

du Tr [FR(x2,x2,i u)] = 0. (28)

B. One potential

The next term of the expansion in V is

G(1) (x2,x1,z) = − 

dw F(x2,w,z)V (w)F(w,x1,z). (29)

In Appendix B we recast Tr[G(1)(x2,x1,z)] into a form in

which the singularity at x2 = x1 is isolated in a relatively

simple term,

Tr[G(1)(x2,x1,z)]

= −Zα

2π



dw e−c|x2−w|

|x2 − w|

ρ(w)

e−c|w−x1|

|w − x1|

− z2

2π2



dw e−c|x2−w|

|x2 − w|

V (w)

e−c|w−x1|

|w − x1|

− 1

2π

[V (x2) + V (x1)]

e−c|x2−x1|

|x2 − x1|

, (30)

where Z is the charge of the nucleus and ρ is the nuclear

charge density normalized to 1. The limit as x2 → x1 of

Tr[G(1)(x2,x1,z)] is undefined due to the third term, so we

introduce a counterterm that cancels the singularity pointwise

and vanishes when integrated over z. The counterterm is

G(1)

A (x2,x1,z)

= − dw F(x2,w,z)

V (x2) + V (x1)

2

F(w,x1,z), (31)

and (see Appendix B)

Tr	G(1)

A (x2,x1,z)




= − z2

2π2



dw e−c|x2−w|

|x2 − w|

V (x2) + V (x1)

2

e−c|w−x1|

|w − x1|

− 1

2π

[V (x2) + V (x1)]

e−c|x2−x1|

|x2 − x1|

. (32)

The difference

G(1)

B (x2,x1,z) = G(1)(x2,x1,z) − G(1)A (x2,x1,z) (33)

has the trace

Tr	G(1)

B (x2,x1,z)




= −Zα

2π



dw e−c|x2−w|

|x2 − w|

ρ(w)

× e−c|w−x1|

|w − x1|

− z2

2π2



dw e−c|x2−w|

|x2 − w|

×



V (w) − V (x2) + V (x1)

2



e−c|w−x1|

|w − x1|

,

(34)

042121-3
PAUL INDELICATO, PETER J. MOHR, AND J. SAPIRSTEIN PHYSICAL REVIEW A 89, 042121 (2014)

with the limit

Tr	G(1)

B (x2,x2,z)




= −Zα

2π



dw ρ(w)

e−2c|x2−w|

|x2 − w|

2

− z2

2π2



dw[V (w) − V (x2)] e−2c|x2−w|

|x2 − w|

2 (35)

for x1 → x2.

To show that the integral over the counterterm vanishes, we

consider the expansion

F (x2,x1,z + δz)

= F (x2,x1,z) +



dw F (x2,w,z) δzF (w,x1,z)

+



dv



dw F (x2,v,z) δzF (v,w,z) δz

× F(w,x1,z) +··· . (36)

This yields

∂

∂ δz

F(x2,x1,z + δz)













δz=0

=



dw F(x2,w,z)F(w,x1,z),

(37)

so we can write

G(1)

A (x2,x1,z) = −V (x2) + V (x1)

2

∂

∂z

F(x2,x1,z).

(38)

With the contour modified as in the zero-potential case, the

leading surviving term in the regulated trace is

∂

∂z

Tr [FR(x2,x2,z)]

= ∂

∂u

1

8πu√u2





2

i=0

Cim4

i + O

 1

u2



, (39)

so that integration gives

∞

−∞

dz

∂

∂z Tr[FR(x2,x1,z)] = Tr[FR(x2,x1,iu)]|

u=∞

u=−∞ = 0

(40)

and the counterterm vanishes.

The subtracted one-potential contribution to the Green

function is regulated to be

G(1)

BR(x2,x1,z) = − 

dw



V (w) − V (x2) + V (x1)

2



× 

2

i=0

CiFi(x2,w,z) Fi(w,x1,z), (41)

and the trace of the corresponding equal coordinate propagation function is

Tr	γ 0

S(1)

FBR(x2,x2)




= 1

2πi

∞

−∞

dz Tr	G(1)

BR(x2,x2,z(1 + iδ))




=

∞

0

du



2

i=0

Ci



− Zα

2π2



dw ρ(w)

e−2ci |x2−w|

|x2 − w|

2

+ u2

2π3



dw [V (w) − V (x2)] e−2ci |x2−w|

|x2 − w|

2



. (42)

In this equation the second line follows from a contour rotation

to the imaginary z axis and a variable change z = iu. The

contribution from the quarter circles as |z|→∞ vanishes due

to the exponential falloff of the integrand. For the possibly

singular case of |x2 − w| ≈ 0, we note that



2

i=0

Ci



dw ρ(w)

e−2ci |x2−w|

|x2 − w|

2

→ ρ(x2)



2

i=0

Ci



dw e−2ci |x2−w|

|x2 − w|

2

= 3π

4 ρ(x2)



2

i=0

Ci

m4

i

u5 + O

 1

u7



(43)

and



2

i=0

Ci u2



dw [V (w) − V (x2)] e−2ci |x2−w|

|x2 − w|

2

→ 

2

i=0

Ci

u2

6

	

∇2

2V (x2)


 

dw e−2ci |x2−w|

= 5π2Zα

4 ρ(x2)



2

i=0

Ci

m4

i

u5 + O

 1

u7



, (44)

which provides sufficient convergence in the integration over

u for the contour rotation to be valid.

Integration by parts in Eq. (42) yields

Tr	γ 0

S(1)

FBR(x2,x2)




=

∞

0

du



2

i=0

Ci



−Zαu2

π2ci



dw ρ(w)

e−2ci |x2−w|

|x2 − w|

+ u4

3π3ci



dw [V (w) − V (x2)] e−2ci |x2−w|

|x2 − w|



, (45)

where estimates analogous to those in Eqs. (43) and (44) show

that the integration over u converges and that the surface terms

at u = ∞ from the partial integration vanish even for |x2 −

w| ≈ 0. For the second term we have



dw [V (w) − V (x2)] e−2ci |x2−w|

|x2 − w|

= −Zαdw



dr ρ(r)

 1

|w − r|

− 1

|x2 − r|

e−2ci |x2−w|

|x2 − w|

= πZα

c2

i



dr ρ(r)

e−2ci |x2−r|

|x2 − r| (46)

042121-4
COORDINATE-SPACE APPROACH TO VACUUM POLARIZATION PHYSICAL REVIEW A 89, 042121 (2014)

so that

Tr	γ 0

S(1)

FBR(x2,x2)




= Zα

π2

∞

0

du



dr ρ(r)



2

i=0

Ci



−u2

ci

+ u4

3c3

i

e−2ci |x2−r|

|x2 − r| .

(47)

Thus the one-potential level shift

E(2,1)

VP = −α



dx2



dx1

1

|x2 − x1|

Tr	γ 0

S(1)

FBR(x2,x2)




× φ†

n(x1) φn(x1) (48)

is the expectation value of a vacuum-polarization potential

given by

V (2,1)

VP (x1) = Zα2

π2



dx2

1

|x2 − x1|

∞

0

du

×



dr ρ(r)



2

i=0

Ci

u2

ci

− u4

3c3

i

e−2ci |x2−r|

|x2 − r| .

(49)

Integration over x2,



dx2

1

|x2 − x1|

e−2ci |x2−r|

|x2 − r| = π

c2

i |x1 − r|

(1 − e−2ci |x1−r|

),

(50)

yields

V (2,1)

VP (x1) = Zα2

π

∞

0

du



dr ρ(r)



2

i=0

Ci

u2

c3

i

− u4

3c5

i



× 1 − e−2ci |x1−r|

|x1 − r| . (51)

We have

∞

0

du



2

i=0

Ci

u2

c3

i

− u4

3c5

i



= −1

3



2

i=0

Ci ln m2

i , (52)

which produces a potential corresponding to a mass-dependent

charge proportional to the charge distribution ρ. This is

ultimately eliminated by charge renormalization. For the

remaining part of Eq. (51), the contribution from each i is

separately finite because of the exponential factor, and there is

no contribution from the terms with i > 0 in the limit of large

auxiliary masses. This can be seen from the nonrelativistic

estimate for the level shift from each of these terms

−Zα2

π

|φn(0)|

2



dx1

∞

0

du

×



dr ρ(r)Ci

u2

c3

i

− u4

3c5

i

e−2ci |x1−r|

|x1 − r|

= −Zα2 |φn(0)|

2

∞

0

duCi

u2

c5

i

− u4

3c7

i



= −4Zα2

15 |φn(0)|

2 Ci

m2

i

. (53)

The surviving term of Eq. (51) with i = 0 is just the Uehling

potential VU(x1) [22,23]. We thus have

V (2,1)

VP (x1) = −Zα2

3π



dr

ρ(r)

|x1 − r|



2

i=0

Ci ln m2

i + VU(x1),

(54)

where

VU(x1) = −Zα2

π



dr ρ(r)

∞

0

du

u2

c3

0

− u4

3c5

0

e−2c0|x1−r|

|x1 − r|

= −Zα2

3π



dr ρ(r)

∞

1

dt



t 2 − 1

 2

t 2 +

1

t 4



× e−2tm0|x1−r|

|x1 − r| , (55)

and in the case of a point charge

VU(x1) = −Zα2

3π

∞

1

dt



t 2 − 1

 2

t 2 +

1

t 4

e−2tm0x1

x1

, (56)

where x1 = |x1|.

To summarize, the Green function in the expression for

the vacuum polarization is not defined for equal coordinates as it appears formally, so a counterterm that removes

this singularity is subtracted, the Pauli-Villars regularization

sum is made for unequal coordinates, and the regulated

expression is taken to be the limit as the coordinates become

equal. Then, integration over the energy parameter in the Green

function is carried out, the charge is renormalized, after which

the auxiliary masses are taken to the infinite limit. The result

is just the Uehling potential.

C. Two potential

The next term in the expansion of the Green function in

powers of V is

G(2)(x2,x1,z) =



dv



dw F(x2,v,z)V (v)F(v,w,z)

× V (w)F(w,x1,z). (57)

This gives a divergent contribution to the level shift, which is

expected to vanish when regulated because of Furry’s theorem.

To confirm that the regulated contribution vanishes in our

framework, we consider the approximation given by

G(2)

A (x2,x1,z) =



dv



dw F(x2,v,z)V (x2)F(v,w,z)

× V (x1)F(w,x1,z). (58)

From Eq. (36), we have

1

2

∂2

∂ δz2 F(x2,x1,z + δz)













δz=0

=



dv



dw F(x2,v,z)F(v,w,z)F(w,x1,z), (59)

so we can write

G(2)

A (x2,x1,z) = 1

2

V (x2)V (x1)

∂2

∂z2 F(x2,x1,z). (60)

042121-5
PAUL INDELICATO, PETER J. MOHR, AND J. SAPIRSTEIN PHYSICAL REVIEW A 89, 042121 (2014)

The zero-component trace is

∂2

∂z2 Tr[Fi(x2,x1,z)] = ∂2∂z2

z e−ci |x2−x1|

π|x2 − x1|

= − ∂2

∂z2

z ci

π

+ O(|x2 − x1|), (61)

which is finite for equal coordinates. We thus have

∂2

∂z2 Tr[Fi(x2,x2,z)] = ∂2∂u2

iu

π



m2

i + u2

= ∂

∂u

i

√

u2

π



2 + m4

i

4u4 +···

(62)

and the regulated Green function has the factor

∂2

∂z2 Tr[FR(x2,x2,z)]

= ∂

∂u

i

4πu2

√

u2





2

i=0

Cim4

i + O

 1

u2



, (63)

which vanishes when integrated over u. The trace of the full

regulated two-potential Green function Tr[G(2)

R (x2,x1,z)] will

have a similar analytic behavior, and it is an odd function of z,

so it also vanishes.

D. Three potential

The three-potential Green function is

G(3)(x2,x1,z) = −

ds



dv



dw F(x2,s,z)V (s)F(s,v,z)

× V (v)F(v,w,z)V (w)F(w,x1,z). (64)

This expression gives the leading term in powers of Zα in the

all-order calculation by Wichmann and Kroll [14] for a point

nucleus. For the effect of the finite size of the nucleus on this

and higher-order terms, see Soff and Mohr [24] and papers

cited therein. A first approximation for this function is given

by

G(3)

A (x2,x1,z) = −

ds



dv



dw F(x2,s,z)V (x2)F(s,v,z)

× V (x2)F(v,w,z)V (x2)F(w,x1,z)

= −[V (x2)]3

6

∂3

∂z3 F(x2,x1,z). (65)

The zero-component trace of the ith mass term is proportional

to

∂3

∂z3 Tr[Fi(x2,x1,z)] = ∂3∂z3

z e−ci |x2−x1|

π|x2 − x1|

= − ∂3

∂z3

z ci

π

+ O(|x2 − x1|), (66)

which is finite for equal coordinates, and

∂3

∂z3 Tr[Fi(x2,x2,z)] = ∂3∂u3

u

π



m2

i + u2

= ∂

∂u

√

u2

π

2

u − 3m4

i

4u5 +···

. (67)

Evidently, the integral over u of this expression is finite

and nonzero, because the relevant branch of the square root is

positive at u = ±∞. This well-known property of the “lightby-light” Feynman diagram yields a spurious finite gaugenoninvariant part [25]. On the other hand, the integral of the

regulated approximate expression vanishes with no ambiguity

from the quarter circles. It is of interest to note that the spurious

contribution is not present if the correction is calculated from

an expansion of the Green function in angular momentum

eigenfunctions [24].

E. All-order generalization

The higher-order terms in the potential expansion of

the vacuum polarization are finite and unambiguous, and

according to Furry’s theorem, only closed loops with an even

number of vertices are nonzero. In the case considered here,

there is one vertex from the interaction with the bound electron

or muon and an odd number from the expansion of the bound

Green function in powers of the external potential.

V. MOMENTUM-SPACE APPROACH

Here, the conventional momentum-space derivation of the

Uehling potential for a point-charge nucleus is reprised in order

to provide a comparison to the coordinate-space approach.

We consider the Fourier transform of the trace of the

propagation function in Eq. (5) starting from the unequal

coordinate case Tr[γ 0SF(x2,x1)]. From the transforms

F(x2,x1,z) = − 1

(2π)3



d p ei p·x2 1

γ · p − m γ 0 e−i p·x1 ,

(68)

where p0 = z, and

V (x) = − Zα

2π2



d p ei p·x 1

p2 (69)

for the point-nucleus potential, we have from Eq. (29), after

integration over w,

G(1)(x2,x1,z) = Zα

16π5



dk



d p eik·x2 ei p·(x2−x1)

× 1

γ · (k + p) − m γ 0 1k2

1

γ · p − m γ 0

,

(70)

where k0 = 0. The singularity for x2 → x1, as seen in Eq. (30),

is here manifested in the divergence of the integral over p for

large momenta when x2 = x1. This singularity is removed by

regulating with the Pauli-Villars summation in the integrand,

which gives

G(1)

R (x2,x2,z)

= Zα

16π5



dk eik·x2

k2



d p



2

i=0

Ci

1

γ · (k + p) − mi

γ 0

× 1

γ · p − mi

γ 0. (71)

042121-6
COORDINATE-SPACE APPROACH TO VACUUM POLARIZATION PHYSICAL REVIEW A 89, 042121 (2014)

We thus have

Tr	γ 0

S(1)

FR(x2,x2)




= Zα

32π6i



dk eik·x2

k2



d4

p 

2

i=0

Ci

× Tr 1

γ · (k + p) − mi

γ 0 1

γ · p − mi

γ 0



, (72)

where it is understood that the variable p0 in the integrand

includes a factor of (1 + iδ), which is equivalent to specification of the Feynman contour. This expression can be written

as

Tr	γ 0

S(1)

FR(x2,x2)




= Zα

32π6i



dk eik·x2

k2 I 00(k), (73)

where

Iμν (k) =



d4p 

2

i=0

Ci Trγ μ 1

γ · (k + p) − mi

γ ν 1

γ · p − mi



. (74)

Rotation of the contour of the variable p0 to the imaginary axis and application of standard Feynman integral evaluation methods

yields (see Appendix C for details)

Iμν (k) = 16π2 i

∞

0

dr r3

2

i=0

Ci

1

0

dy

(gμν k2 − kμkν )y(1 − y)

	

r2 + k2y(1 − y) + m2

i


2

+ 2π2 i

∞

0

dr d

dr r4

2

i=0

Ci

1

0

dy

gμν

r2 + k2y(1 − y) + m2

i

. (75)

The second term is not formally gauge invariant, but the regulated integrand falls off sufficiently rapidly that the integral of the

derivative vanishes. A single mass counterterm regularization would yield a nonzero result for this term. For the first term, the

change of variable y = (t − √t 2 − 1)/2t and integration by parts gives

1

0

dy

y(1 − y)

	

r2 + k2y(1 − y) + m2

i


2 =

1/2

0

dy

2y(1 − y)

	

r2 + k2y(1 − y) + m2

i


2 =

∞

1

dt

1

4t 4

√t 2 − 1

1



r2 + k2/4t 2 + m2

i

2

= 1

6

1



r2 + m2

i

2 −

∞

1

dt

√t 2 − 1

12  2t 4 +

1

t 6

 k2



r2 + k2/4t 2 + m2

i

3 , (76)

and hence

∞

0

dr r3

2

i=0

Ci

1

0

dy

y(1 − y)

	

r2 + k2y(1 − y) + m2

i


2 = − 1

12



2

i=0

Ci



ln m2

i +

∞

1

dt



t 2 − 1

 2

t 2 +

1

t 4

 k2

k2 + 4t 2m2

i



= − 1

12



2

i=0

Ci ln m2

i +

1

3

∞

1

dt



t 2 − 1



2 +

1

t 2



2

i=0

Cim2

i

k2 + 4t 2m2

i

. (77)

Integration over k in Eq. (73) yields

Tr	γ 0

S(1)

FR(x2,x2)




= Zα

3π



2

i=0

Ci ln m2

i δ(x2) − Zα

3π2

∞

1

dt



t 2 − 1



2 +

1

t 2



2

i=0

Cim2

i

e−2tmix2

x2

. (78)

According to Eq. (48), the level shift corresponds to a potential energy given by

V (2,1)

VP (x1) = −α



dx2

1

|x2 − x1|

Tr	γ 0

S(1)

FR(x2,x2)




= −Zα2

3π





2

i=0

Ci

ln m2

i

x1

−

∞

1

dt



t 2 − 1

 2

t 2 +

1

t 4



2

i=0

Ci

1 − e−2tmix1

x1



= −Zα2

3π





2

i=0

Ci

ln m2

i

x1

+

∞

1

dt



t 2 − 1

 2

t 2 +

1

t 4

e−2tm0x1

x1



. (79)

The term on the second line that is independent of mass vanishes from the condition on Ci, and the terms with the exponential

factor and with i > 0 make no contribution in the large mass limit, as shown by an estimate analogous to that in Eq. (53). The

first term on the third line is eliminated by charge renormalization, and the second term is just the Uehling potential, in agreement

with Eq. (56).

042121-7
PAUL INDELICATO, PETER J. MOHR, AND J. SAPIRSTEIN PHYSICAL REVIEW A 89, 042121 (2014)

The results of the coordinate-space and momentum-space

calculations agree, as they must, since Pauli-Villars regularization makes the result finite, but the charge renormalization

terms arise in different ways. This is seen by comparing the

point-charge special case of Eq. (47),

Tr	γ 0

S(1)

FBR(x2,x2)




= Zα

π2

∞

0

du



2

i=0

Ci



−u2

ci

+ u4

3c3

i



× e−2cix2

x2

, (80)

to Eq. (78). For x2 = 0, only the second term of Eq. (78)

is nonzero. In Eq. (80), the exponential factor provides

convergence for large u, so that the integrals in the sum over i

are separately finite. Thus, the mass-dependent variable change

t = ci/mi may be made for each i and the result is the same

as the second term in Eq. (78). However, if x2 = 0, then the

integrals over u are not separately finite and a mass-dependent

variable change is not valid. The effect of including the point

x2 = 0 can be checked by integrating either expression over

x2. The result of integration of Eq. (78) is



dx2 Tr	

γ 0S(1)

FR(x2,x2)




= Zα

3π



2

i=0

Ci ln m2

i , (81)

since the integral over the second term vanishes. Integration

of Eq. (80) gives



dx2 Tr	

γ 0S(1)

FBR(x2,x2)




= Zα

π

∞

0

du



2

i=0

Ci



−u2

c3

i

+ u4

3c5

i



= Zα

3π



2

i=0

Ci ln m2

i . (82)

These results are in agreement, although the evolution of the

logarithmic terms is quite different.

It is worth pointing out that if the (incorrect) variable

change u → miu were made in each term in the sum over

i in Eq. (82), the mass dependence would drop out and the

sum would vanish. This illustrates the fact that the order of

summation and integration is crucial and particular care is

needed in dealing with such potentially divergent expressions.

VI. CONCLUSION

Vacuum polarization for a spherically symmetric potential

is examined in coordinate space in this work. For the

contribution of first order in the expansion in powers of the

potential, a counterterm is introduced to remove the equal

coordinate singularity, which could be problematic for a purely

numerical calculation. Although this singularity is removed

in principle by Pauli-Villars regularization, in practice, the

subtraction can be expected to improve the numerical convergence by providing a pointwise removal of the singularity

before numerical integrations are carried out. Moreover, this

approach can be expected to be useful for more general

calculations.

The coordinate-space calculation is compared to the

momentum-space calculation. With Pauli-Villars regularization, the logarithmic charge renormalization term is seen to

be the same in either case, despite the fact that it arises

in a completely different way in the two approaches. As

with the coordinate-space calculation, the momentum-space

calculation is explicitly based on the use of two Pauli-Villars

auxiliary mass subtraction terms, and it can be seen that the

result is not well defined unless both terms are included.

APPENDIX A: VECTOR VACUUM POLARIZATION

For a spherically symmetric binding potential, the vector

contribution to the vacuum polarization, in the last line of

Eq. (9), is finite and vanishes when Pauli-Villars regularization

is applied. The fact that it formally vanishes follows from the

spin-angular momentum expansion of the Green function for

a spherically symmetric potential. In particular, in this case

the wave function can be written in terms of the Dirac spinangle functions χ as (see, for example, [26] and references

therein)

φn(x) =

 f1(x)χμ

κ (xˆ)

f2(x)χμ

−κ (xˆ)



(A1)

and the Green function is given by

G(x2,x1,z) = 

κμ  G11

κ (x2,x1,z)χμκ (xˆ 2) χμ†κ (xˆ 1) −i G12κ (x2,x1,z)χμκ (xˆ 2) χμ†

−κ (xˆ 1)

i G21

κ (x2,x1,z)χμ

−κ (xˆ 2) χμ†

κ (xˆ 1) G22κ (x2,x1,z)χμ

−κ (xˆ 2) χμ†

−κ (xˆ 1)



, (A2)

where

Gij

κ (x2,x1,z) = 

n

fi(x2) fj (x1)

En − z

(A3)

and



μ

χμ

κ (xˆ 2) χμ†κ (xˆ 1) = |κ|

4π



IPκ+ (ξ ) +

i

κ

σ · (xˆ 2 × xˆ 1)P 

κ+ (ξ )



, (A4)



μ

χμ

−κ (xˆ 2) χμ†

κ (xˆ 1) = 1

4π

κ

|κ|

{σ · xˆ 2 P 

κ− (ξ ) − σ · xˆ 1 P κ+ (ξ )}, (A5)

042121-8
COORDINATE-SPACE APPROACH TO VACUUM POLARIZATION PHYSICAL REVIEW A 89, 042121 (2014)

with ξ = xˆ 2 · xˆ 1 and κ± = |κ ± 1/2| − 1/2. The vector trace

in Eq. (9) is

Tr[α G(x2,x1,z)]

= i

2π



κ

κ

|κ|



G21

κ (x2,x1,z)[xˆ 2 P κ− (ξ ) − xˆ 1 P κ+ (ξ )]

+ G12

κ (x2,x1,z)[xˆ 2 P κ+ (ξ ) − xˆ 1 P κ− (ξ )]

. (A6)

The equal coordinate limit is ambiguous due to the discontinuity in the radial factor at x2 = x1, but the angular factors are

just

[xˆ 2 P 

κ± (ξ ) − xˆ 1 P κ∓ (ξ )]xˆ 1=xˆ 2 = ±κ xˆ 2, (A7)

and integration of the level-shift expression in Eq. (9) includes



d2

xˆ 2

|x2 − x1|

= 4π

3

xˆ 1

x<

x2

>

, (A8)

so the vector term is proportional to

φ†

n(x1) α · xˆ 1 φn(x1) = 0. (A9)

Although the vector term vanishes as shown above, the

integrals giving the level shift from the first few terms in the

expansion of the Green function in powers of the potential are

not convergent. However, these terms are well defined when

Pauli-Villars regularization is applied. The analysis is similar

to that employed in evaluating the nonvanishing terms of the

vacuum polarization and is not repeated in this case.

APPENDIX B: ONE-POTENTIAL GREEN FUNCTION

From Eqs. (16) and (29), we have

Tr[G(1)(x2,x1,z)] = − 1

4π2 [∇2 · ∇1 + m2 + z2

]

×



dw e−c|x2−w|

|x2 − w|

V (w)

e−c|w−x1|

|w − x1|

.

(B1)

Since 2 ∇2 · ∇1 = (∇2 + ∇1)

2 − ∇2

2 − ∇21, Eq. (30) follows

from

(∇2 + ∇1)

2



dw e−c|x2−w|

|x2 − w|

V (w)

e−c|w−x1|

|w − x1|

=



dw e−c|x2−w|

|x2 − w|

	

∇2

w V (w)


 e−c|w−x1|

|w − x1| (B2)

together with

∇2

w V (w) = 4πZα ρ(w) (B3)

and



∇2

i − c2 e−c|xi−w|

|xi − w| = −4πδ(xi − w). (B4)

For the counterterm in Eq. (31), we have

Tr	G(1)

A (x2,x1,z)




= − 1

4π2



[∇2 · ∇1 + m2 + z2

]

×



dw e−c|x2−w|

|x2 − w|

e−c|w−x1|

|w − x1|

V (x2) + V (x1)

2 (B5)

and Eq. (32) follows from

∇2 · ∇1



dw e−c|x2−w|

|x2 − w|

e−c|w−x1|

|w − x1|

= 2π

c ∇2 · ∇1 e−c|x2−x1|

= 2π

c

 2c

|x2 − x1|

− c2



e−c|x2−x1|

= 4π

|x2 − x1|

e−c|x2−x1| + (z2 − m2

)

×



dw e−c|x2−w|

|x2 − w|

e−c|w−x1|

|w − x1|

. (B6)

APPENDIX C: MOMENTUM-SPACE INTEGRATION

Here we give some details of the evaluation of the function

Iμν (k) =



d4p 

2

i=0

Ci

× Trγ μ 1

γ · (k + p) − mi

γ ν 1

γ · p − mi



(C1)

that appears in Sec. V. Rationalization of the propagation functions and application of the Feynman denominator formula

1

AB =

1

0

dy

1

[Ay + B(1 − y)]2 (C2)

yields

Iμν (k) =



d4p 

2

i=0

Ci

1

0

dy

× Tr[γ μ(γ · p + γ · k + mi)γ ν (γ · p + mi)]

	

(p + ky)2 + k2y(1 − y) − m2

i + i


2 ,

(C3)

where  = 2δp02 and δ2p02 is dropped. After the translation

p → p − ky, the numerator is

Tr{γ μ[γ · p + γ · k(1 − y) + mi]γ ν (γ · p − γ · ky + mi)}

→ 8(gμν k2 − kμkν )y(1 − y) − 4gμν

× 	p2 + k2y(1 − y) − m2

i




+ 8 pμpν , (C4)

where terms odd in p are not included, and the denominator

is

p2 + k2

y(1 − y) − m2

i + i

= p02 − p2 − k2

y(1 − y) − m2

i + i. (C5)

Poles of the integrand are located at

p0 = ±	

p2 + k2y(1 − y) + m2

i − i


1/2 (C6)

in the second and fourth quadrants of the complex p0 plane,

so the contour of the p0 integration may be rotated to the

imaginary axis, and p may be replaced by a Cartesian vector

q, where p0 = iq0 and pi = qi for i = 1,2,3. The integrals

over the four-vector p are thus expressed as integrals over

042121-9
PAUL INDELICATO, PETER J. MOHR, AND J. SAPIRSTEIN PHYSICAL REVIEW A 89, 042121 (2014)

r = (q2

0 + q21 + q22 + q23 )

1/2, the magnitude of the Cartesian

four-vector q, where



d4

p f (p2) = 2π2 i

∞

0

dr r3f (−r2) , (C7)



d4

p pμpνf (p2) = −π2 i

2

∞

0

dr r5gμνf (−r2). (C8)

We thus have

Iμν (k) = 16π2 i

∞

0

dr r3

2

i=0

Ci

×

1

0

dy

(gμν k2 − kμkν )y(1 − y)

	

r2 + k2y(1 − y) + m2

i


2

+ 2π2 i

∞

0

dr d

dr r4

2

i=0

Ci

×

1

0

dy

gμν

r2 + k2y(1 − y) + m2

i

. (C9)

It is evident that the second term vanishes with Pauli-Villars

regularization from the identity



2

i=0

Ci

R2 + m2

i

= C0m2

1m22 + C1m20m22 + C2m20m21



R2 + m2

0

R2 + m2

1

R2 + m2

2

 , (C10)

where R2 = r2 + k2y(1 − y). Note that in the numerator on

the right-hand side of Eq. (C10), a possible term proportional

to R4 has the vanishing coefficient 2

i=0 Ci and a possible

term proportional to R2 has the coefficient

C0



m2

1 + m22



+ C1



m2

0 + m22



+ C2



m2

0 + m21



= C0



m2

1 + m22



+ C1



m2

0 + m22



+ C2



m2

0 + m21



+ 

2

i=0

Cim2

i =





2

i=0

Ci

 ⎛

⎝

2

j=0

m2

j

⎞

⎠ , (C11)

which also vanishes.

[1] R. Pohl, A. Antognini, F. Nez, F. D. Amaro, F. Biraben, J. a. M.

R. Cardoso, D. S. Covita, A. Dax, S. Dhawan, L. M. P. Fernandes

et al., Nature (London) 466, 213 (2010).

[2] R. Pohl, F. D. Amaro, A. Antognini, F. Biraben, J. a. M. R.

Cardoso, D. S. Covita, A. Dax, S. Dhawan, L. M. P. Fernandes,

A. Giesen et al., J. Phys. Conf. Ser. 264, 012008 (2011).

[3] A. Antognini, F. Nez, K. Schuhmann, F. D. Amaro, F. Biraben,

J. a. M. R. Cardoso, D. S. Covita, A. Dax, S. Dhawan, M.

Diepold et al., Science 339, 417 (2013).

[4] P. J. Mohr, B. N. Taylor, and D. B. Newell, Rev. Mod. Phys. 84,

1527 (2012).

[5] R. Pohl, R. Gilman, G. A. Miller, and K. Pachucki, Annu. Rev.

Nucl. Part. Sci. 63, 175 (2013).

[6] A. Antognini, F. Kottmann, F. Biraben, P. Indelicato, F. Nez, and

R. Pohl, Ann. Phys. (NY) 331, 127 (2013).

[7] P. Indelicato, Phys. Rev. A 87, 022501 (2013).

[8] E. Borie, Ann. Phys. (NY) 327, 733 (2012); arXiv:1103.1772v6.

[9] S. G. Karshenboim, V. G. Ivanov, and E. Y. Korzinin, Phys. Rev.

A 85, 032509 (2012).

[10] U. D. Jentschura, Ann. Phys. (NY) 326, 500 (2011).

[11] U. D. Jentschura, Ann. Phys. (NY) 326, 516 (2011).

[12] J. D. Carroll, A. W. Thomas, J. Rafelski, and G. A. Miller, Phys.

Rev. A 84, 012506 (2011).

[13] P. Indelicato and P. J. Mohr, Phys. Rev. A 46, 172 (1992).

[14] E. H. Wichmann and N. M. Kroll, Phys. Rev. 101, 843 (1956).

[15] D. Z. Freedman, K. Johnson, and J. I. Latorre, Nucl. Phys. B

371, 353 (1992).

[16] The curly brackets {} denote the numerical coefficient of the

value of the quantity. For example, in SI units{c} = 299 792 458.

Symbols for mass are understood to be dimensionless numbers

for the mass relative to the electron mass me and symbols for

energy to be relative to mec2, etc.

[17] P. J. Mohr, Phys. Rev. A 32, 1949 (1985).

[18] P. J. Mohr, G. Plunien, and G. Soff, Phys. Rep. 293, 227 (1998).

[19] W. Pauli and F. Villars, Rev. Mod. Phys. 21, 434 (1949).

[20] H. M. Fried and D. R. Yennie, Phys. Rev. 112, 1391 (1958).

[21] A. J. Layzer, J. Math. Phys. 2, 292 (1961).

[22] E. A. Uehling, Phys. Rev. 48, 55 (1935).

[23] R. Serber, Phys. Rev. 48, 49 (1935).

[24] G. Soff and P. J. Mohr, Phys. Rev. A 38, 5066 (1988).

[25] R. Karplus and M. Neuman, Phys. Rev. 80, 380 (1950).

[26] P. J. Mohr, Ann. Phys. (NY) 88, 26 (1974).

042121-10