# Dimensional Regularization and Two–Loop Vacuum Polarization Operator:
Master Integrals, Analytic Results and Energy Shifts

Stefano Laporta  Dipartimento di Fisica e Astronomia, Università di Padova e
Istituto Nazionale di Fisica Nucleare, Sezione di Padova, Via Marzolo 8,
I-35131 Padova, Italy  Ulrich D. Jentschura  Department of Physics and LAMOR,
Missouri University of Science and Technology, Rolla, Missouri 65409, USA

###### Abstract

We present a complete reevaluation of the irreducible two-loop vacuum-
polarization correction to the photon propagator in quantum electrodynamics,
i.e. with an electron-positron pair in the fermion propagators. The
integration is carried out by reducing the integrations to a limited set of
master integrals, which are calculated using integration-by-parts identities.
Dimensional regularization is used in $D=4-2\epsilon$ dimensions, and on-mass
shell renormalization is employed. The one-loop effect is given to order
$\epsilon$, to be combined with the $1/\epsilon$ divergence of the two-loop
amplitude. Master integrals are given. Final evaluations of two-loop energy
shifts for $1S$, $2S$, and $2P$ states are done analytically, and results are
presented, with an emphasis on muonic hydrogen. For relativistic Dirac–Coulomb
reference states, higher-order coefficients are obtained for the
$Z\alpha$-expansion. We compare the results obtained to the existing
literature.

###### Contents

  1. I Introduction
  2. II One–Loop Vacuum Polarization
    1. II.1 Orientation
    2. II.2 Calculation
    3. II.3 Renormalization
  3. III Two–Loop Vacuum Polarization
    1. III.1 Orientation
    2. III.2 Reduction to Master Integrals
    3. III.3 Differential Equations and Master Integrals
    4. III.4 Results for the Master Integrals
    5. III.5 Renormalization
    6. III.6 Results
    7. III.7 Comparison with the Literature
  4. IV Contributions to Energy Shifts
    1. IV.1 Muonic Hydrogen
    2. IV.2 Relativistic Ordinary Hydrogen
  5. V Conclusions
  6. A Terms of Order $\maybebm{\epsilon^{2}}$
  7. B Coefficient Functions
  8. C Higher–Order Terms in the Master Integrals
  9. D Analytic Integrations
    1. D.1 One–Loop Diagram
    2. D.2 Two–Loop Reducible Diagram
    3. D.3 Two–Loop Irreducible Diagram
    4. D.4 Individual Energy Shifts

##  I Introduction

The two-loop vacuum-polarization correction to bound-state energy levels is an
important contribution to the Lamb shift in muonic bound systems.
Specifically, the irreducible two-loop effect lowers the 2S level energy in
muonic hydrogen by roughly 1.25 meV in comparison to the $2P$ level Pa1996mu ;
Je2011aop1 . This energy shift needs to be compared to the proton-size
discrepancy commonly referred to as the proton radius puzzle (0.88 fm versus
0.84 fm, for a recent brief summary see Ref. Je2022proton ). The proton radius
puzzle corresponds to a 0.3 meV shift of the $2S$ states, when expressed in
energy units Pa1996mu ; Je2011aop1 . The two-loop effect is thus
phenomenologically extremely relevant; any conceivable inaccuracy in its
evaluation could lead to at least a partial theoretical explanation of the
proton size discrepancy. Vacuum-polarization effects are drastically enhanced
in muonic as compared to electronic bound systems Pu1957 ; JeSoIvKa1997 ;
KaIvJeSo1998 ; KaJeIvSo1998 . Here, we present a reevaluation of the
irreducible two-loop diagrams, on the basis of modern integration-by-parts
techniques, employing dimensional regularization. We note that dimensional
regularization was not yet sufficiently developed in 1973 (despite Refs.
tHVe1972 ; BoGi1972  which appeared in the literature at the time) to make an
evaluation of the two-loop effect using dimensional regularization techniques
feasible BaRe1973 . For vacuum-polarization effects in particular, the modern
techniques lead to drastic simplifications of the calculations. We find both
the imaginary as well as the real part of the two-loop vacuum tensor in closed
analytic form.

We also derive analytic expressions for the expectation value of the two-loop
potential, evaluated on (nonrelativistic) Schrödinger–Coulomb eigenstates,
generalizing the treatment originally outlined in Ref. Pu1957  to two-loop
order. In general, for predominantly nonrelativistic bound systems, energy
shifts can be represented in terms of a semi-analytic expansion (in powers of
$Z\alpha$ and $\ln[(Z\alpha)^{-2}]$, where $\alpha$ is the fine-structure
constant and $Z$ is the nuclear charge number). As a byproduct, we derive a
few higher-order coefficients from the expectation value of two-loop
potential, evaluated on the relativistic Dirac–Coulomb eigenstates.

This paper is organized as follows. We start by briefly discussing the
evaluation of the one-loop vacuum-polarization effect in Sec. II. Throughout
this paper, we use dimensional regularization with the number $D$ of
dimensions expressed as

$D=4-2\epsilon\,.$ (1)

The one-loop effect needs to be evaluated to order $\epsilon$, because finite
contributions are generated at two-loop order when the one-loop terms of order
$\epsilon$ are multiplied by the $1/\epsilon$-terms from the two-loop
amplitudes. The irreducible two-loop diagrams are discussed in Sec. III,
including master integrals and their series expansions, as well as complete
results. A comparison to the available literature is also performed. In Sec.
IV, we demonstrate that energy shifts due to the two-loop effect can be
evaluated analytically, for muonic bound systems. Conclusions are reserved for
Sec. V.

##  II One–Loop Vacuum Polarization

###  II.1 Orientation

Our aim is to evaluate the one-loop vacuum-polarization effect to order
$\epsilon$, i.e., in a form suitable for the later two-loop calculation. For
simplicity we set the lepton mass (electron mass) in the vacuum-polarization
loops to be $m=1$. We use for a metric tensor
$g_{\mu\nu}=\text{diag}(1,-1,-1,-1)$ (see, for example, Refs. BjDr1965 ;
PeSc1995 ), where

$\displaystyle p^{\mu}=$ $\displaystyle\;(p_{0},\vec{p}),\qquad
p_{\mu}=(p_{0},-\vec{p}),\qquad$ (2)
$\displaystyle\\{\gamma^{\mu},\gamma^{\nu}\\}=$
$\displaystyle\;2g^{\mu\nu},\qquad
S(p)=\frac{1}{\cancel{p}-m+\mathrm{i}\epsilon}\,,\qquad m=1\,.$ (3)

Here, $m=1$ is the electron mass.

Due to current conservation, the vacuum polarization function $\Pi_{\mu\nu}$
has this tensorial structure

$\Pi_{\mu\nu}(q^{2})=(q_{\mu}q_{\nu}-q^{2}\,g_{\mu\nu})\,\Pi(q^{2})\,.$ (4)

We follow the convention of Ref. ItZu1980  so that

$\left.\frac{\mathrm{d}\Pi(q^{2})}{\mathrm{d}q^{2}}\right|_{q^{2}=0}>0$ (5)

at the one-loop and two-loop level. For completeness, we note that the
opposite conventions are used elsewhere in the literature, e.g., in Refs.
BeLiPi1982vol4 ; JeAd2022book . The conventions employed here are consistent
with those employed in Ref. ItZu1980 . Summation over the Lorentz indices in
$D$ dimensions leads to the result

${\Pi^{\mu}}_{\mu}(q^{2})=(-q^{2}\,{g^{\mu}}_{\mu}+q^{\mu}\,q_{\mu})\,\Pi(q^{2})=q^{2}\,(1-D)\,\Pi(q^{2})\,.$
(6)

So, we can extract the scalar vacuum polarization function $\Pi(q^{2})$ as

$\Pi(q^{2})=\frac{{\Pi^{\mu}}_{\mu}(q^{2})}{q^{2}\,(1-D)}\,.$ (7)

The threshold for pair production of $\Pi(q^{2}=s)$ is $q^{2}=s=4$. The
argument of the polarization tensor is

$q^{2}=(q^{0})^{2}-\vec{q}^{\,2}\,,$ (8)

which is the four-momentum square of the photon entering the vacuum-
polarization loop.

###  II.2 Calculation

For the renormalization of the two-loop diagram with a self-mass insertion we
need the one-loop vacuum polarization function $\Pi^{(1)}$, and the function
$\Pi^{(1a)}$ from the same diagram with iterated electron propagators on one
side of the loop. The unrenormalized tensors read as follows,

$\displaystyle\Pi^{(1)}_{\mu\nu}(q^{2})=$
$\displaystyle\;{-\mathrm{i}}e^{2}\int\dfrac{\mathrm{d}^{D}k}{(2\pi)^{D}}\left[-{\rm
Tr}\left(\gamma_{\mu}S(k)\gamma_{\nu}S(k-q)\right)\right]\,,$ (9)
$\displaystyle\Pi^{(1a)}_{\mu\nu}(q^{2})=$
$\displaystyle\;{\mathrm{i}}e^{2}\int\dfrac{\mathrm{d}^{D}k}{(2\pi)^{D}}\left[-{\rm
Tr}\left(\gamma_{\mu}[S(k)]^{2}\gamma_{\nu}S(k-q)\right)\right]\,,$ (10)

The needed quantities are expressed in terms of two one-loop master integrals,

$\displaystyle\Pi^{(1)}(q^{2})=$
$\displaystyle\;\left(\frac{e^{2}\,C(\epsilon)}{4\pi\alpha}\right)\,\left(\dfrac{\alpha}{\pi}\right)\,P^{(1)}(q^{2})\,,$
(11a) $\displaystyle\Pi^{(1a)}(q^{2})=$
$\displaystyle\;\left(\frac{e^{2}\,C(\epsilon)}{4\pi\alpha}\right)\,\left(\dfrac{\alpha}{\pi}\right)\,P^{(1a)}(q^{2})\,.$
(11b)

The prefactor in these results simplifies to unity provided we set

$e^{2}=\frac{4\pi\alpha}{C(\epsilon)}\,,\qquad
C(\epsilon)=\Gamma(1+\epsilon)\;(4\pi)^{\epsilon}\;{\mbox{\sl\textmu}}^{-2\epsilon}\,.$
(12)

One observes that the identification (12) is not completely canonical in
dimensional regularization. Namely [see, e.g., Eq. (10.173) of Ref.
JeAd2022book ], in the $\overline{\mathrm{MS}}$ scheme, one normally has the
relation
$e^{2}=(4\pi)^{1-\epsilon}\,\alpha\,\mbox{\sl\textmu}^{2\epsilon}\mathrm{e}^{\gamma_{E}\epsilon}$,
where $\gamma_{E}=0.57721\,56649\dots$ is the Euler–Mascheroni constant, and µ
is the renormalization scale induced by dimensional regularization. This is
equivalent to Eq. (12) up to terms linear in $\epsilon$. The prefactor
$1/C(\epsilon)$ helps to simplify results after integration.

The one-loop integrals are as follows,

$\displaystyle P^{(1)}(q^{2})=$
$\displaystyle\;\frac{(D-2)q^{2}+4}{2(D-1)q^{2}}M_{12}(q^{2})+\frac{(D-2)M_{2}(q^{2})}{(D-1)q^{2}}\,,$
(13) $\displaystyle P^{(1a)}(q^{2})=$
$\displaystyle\;\frac{(D-4)q^{2}+4}{q^{2}\left(q^{2}-4\right)}M_{12}(q^{2})+\frac{2(D-2)M_{2}(q^{2})}{q^{2}\left(q^{2}-4\right)}\,,$
(14)

where

$\displaystyle M_{12}(q^{2})=$
$\displaystyle\;N(\epsilon)\int\frac{[\mathrm{d}^{D}k_{1}]}{D_{1}D_{2}}\,,$
(15a) $\displaystyle M_{2}(q^{2})=$
$\displaystyle\;N(\epsilon)\int\frac{[\mathrm{d}^{D}k_{1}]}{D_{2}}\,,$ (15b)
$\displaystyle D_{j}=$
$\displaystyle\;-p_{j}^{2}+1-\mathrm{i}\epsilon\,,\qquad j=1,\ldots,2\,.$
(15c) $\displaystyle p_{1}=$ $\displaystyle\;k_{1}\,,\qquad
p_{2}=k_{1}-q\;\;\,,$ (15d) $\displaystyle[\mathrm{d}^{D}k]=$
$\displaystyle\mbox{\sl\textmu}^{-2\epsilon}\mathrm{d}^{D}k\,.$ (15e)

The above definition of $[\mathrm{d}^{D}k]$ contains the factor
$\mbox{\sl\textmu}^{-2\epsilon}$. It simplifies with the factor
$1/\mbox{\sl\textmu}^{-2\epsilon}$ of $N(\epsilon)$ so that the master
integrals are independent of µ. The normalization factor is

$N(\epsilon)=\frac{\mbox{\sl\textmu}^{2\epsilon}}{\pi^{2-\epsilon}\;\Gamma(1+\epsilon)}=\frac{4^{\epsilon}}{\pi^{2-2\epsilon}C(\epsilon)}=\frac{4(4\pi^{2})}{(2\pi){{}^{4-2\epsilon}}C(\epsilon)}\,.$
(16)

The latter form exhibits the dimension $D=4-2\epsilon$. With the normalization
factor $N(\epsilon)$, all master integrals $M$ occurring in our calculations
are dimensionless and not only uncluttered from logarithms of µ, but also,
from logarithms of $4\pi$.

One can show that the one-loop master integrals satisfy the following
differential equations,

$\displaystyle M_{2}^{\prime}\left(q^{2}\right)=0\,,$ (17a)
$\displaystyle\left[(D-4)q^{2}+4\right]\,M_{12}(q^{2})+2(D-2)M_{2}\left(q^{2}\right)$
$\displaystyle-2\left(q^{2}-4\right)q^{2}M_{12}^{\prime}\left(q^{2}\right)=0\
,$ (17b)

where the differentiation is with respect to $q^{2}$. One finds that $M_{2}$
is a momentum-independent vacuum integral,

$M_{2}(q^{2})=-\frac{1}{\epsilon(1-\epsilon)}\,.$ (18)

For $M_{12}(q^{2})$, one writes the result as a function of the variable
$v=v(q^{2})$,

$v=\sqrt{1-\frac{4}{q^{2}}}\,.$ (19)

The result is

$\displaystyle(1-$ $\displaystyle
2\epsilon)M_{12}(q^{2})=\frac{1}{\epsilon}+v\ln\left(\frac{v-1}{v+1}\right)-\frac{\epsilon}{6}v\biggl{[}\pi^{2}+$
$\displaystyle
12\text{Li}_{2}\left(\frac{1-v}{v+1}\right)+12\ln\left(\frac{2v}{v+1}\right)\ln\left(\frac{v-1}{v+1}\right)$
$\displaystyle-3\ln^{2}\left(\frac{v-1}{v+1}\right)\biggr{]}+\epsilon^{2}\frac{v}{6}N_{12}(v)+{\mathcal{O}}(\epsilon^{3})\,,$
(20)

where we include the terms of order $\epsilon$ and $\epsilon^{2}$ which will
be useful in the context of the two-loop calculation. The explicit form of
$N_{12}(q^{2})$ is given in Appendix A.

###  II.3 Renormalization

Substituting Eqs. (18) and (II.2) into (13), we can find the expansions of
$P^{(1)}(q^{2})$ and $P^{(1a)}(q^{2})$. But for the renormalization we need to
subtract the values at $q^{2}=0$,

$P^{(1)}(0)=\frac{1}{3\epsilon}\ ,\qquad P^{(1a)}(0)=\frac{1}{3}\,.$ (21)

In our formalism, the scaled functions $P^{(1)}$ and $P^{(2)}$ do not have
powers of µ, but the scalar vacuum polarization functions $\Pi^{(1)}$,
$\Pi^{(2)}$ depend on µ. Yet, after renormalization the dependence disappears
in the limit $\epsilon\to 0$. The renormalized expression for the one-loop
function $P^{(1)}_{\rm R}(q^{2})$ is as follows,

$P^{(1)}_{\rm R}(q^{2})=P^{(1)}(q^{2})-P^{(1)}(0)\\\
=\frac{1}{9}\left(8-3v^{2}\right)-\frac{v}{6}\left(v^{2}-3\right)\ln\left(\frac{v-1}{v+1}\right)\\\
+\epsilon\,Q^{(1)}_{1}(v)+\epsilon^{2}\,Q^{(1)}_{2}(v)+{\mathcal{O}}(\epsilon^{3})\,.$
(22)

The coefficient of the term linear in $\epsilon$ is

$Q^{(1)}_{1}(v)=\frac{v}{36}\left(v^{2}-3\right)\biggl{[}12\text{Li}_{2}\left(\frac{1-v}{v+1}\right)\\\
-3\ln^{2}\left(\frac{v-1}{v+1}\right)+12\ln\left(\frac{2v}{v+1}\right)\,\ln\left(\frac{v-1}{v+1}\right)+\pi^{2}\biggr{]}\\\
+\left(v-\frac{4v^{3}}{9}\right)\ln\left(\frac{v-1}{v+1}\right)-\frac{8v^{2}}{9}+\frac{52}{27}\,.$
(23)

For the explicit form of the term of order $\epsilon^{2}$, we refer to
Appendix A.

The renormalized value for the function $P^{(1a)}_{\rm R}(q^{2})$ with
iterated electron propagators is as follows,

$P^{(1a)}_{\rm R}(q^{2})=P^{(1a)}(q^{2})-P^{(1a)}(0)\\\
=\frac{1-v^{2}}{4v}\biggl{\\{}\frac{2v\left(3v^{2}-5\right)}{3\left(v^{2}-1\right)}+\left(v^{2}-1\right)\ln\left(\frac{v-1}{v+1}\right)\biggr{\\}}\\\
+\epsilon
Q^{(1a)}_{1}(v)+\epsilon^{2}Q^{(1a)}_{2}(v)+{\mathcal{O}}(\epsilon^{3})\,.$
(24)

The term linear in $\epsilon$ is given as follows,

$Q^{(1a)}_{1}(v)=\frac{1-v^{2}}{4v}\biggl{\\{}\biggl{[}\left(v^{2}-1\right)\biggl{(}-2\text{Li}_{2}\left(\frac{1-v}{v+1}\right)\\\
+\frac{1}{2}\ln^{2}\left(\frac{v-1}{v+1}\right)-\frac{\pi^{2}}{6}\biggr{)}+2\ln\left(\frac{v-1}{v+1}\right)\\\
\times\biggl{(}v^{2}-\left(v^{2}-1\right)\ln\left(\frac{2v}{v+1}\right)\biggr{)}+4v\biggr{]}\biggr{\\}}$
(25)

The quadratic term in $\epsilon$, denoted as $Q^{(2a)}_{1}(v)$, is given in
Appendix A. Actually, only the terms proportional to $\epsilon$ are necessary
for the renormalization of the two-loop vacuum polarization, but higher-order
terms are useful for higher-loop calculations, and are thus relegated to
Appendix A.

In the limit $\epsilon\to 0^{+}$, one recovers the known result for the one-
loop vacuum polarization,

$\Pi^{(1)}_{\rm R}(q^{2})=\left(\dfrac{\alpha}{\pi}\right)P^{(1)}_{\rm
R}(q^{2})\\\
=\left(\dfrac{\alpha}{\pi}\right)\left[\frac{1}{9}\left(8-3v^{2}\right)-\frac{1}{6}v\left(v^{2}-3\right)\ln\left(\frac{v-1}{v+1}\right)\right]\,.$
(26)

The dispersion relation is

$\Pi^{(1)}_{\rm
R}(q^{2})=\frac{q^{2}}{\pi}\int_{4}^{\infty}\mathrm{d}q^{\prime 2}\frac{{\rm
Im}\,\Pi^{(1)}_{\rm R}(q^{\prime 2}+\mathrm{i}\epsilon)}{q^{\prime
2}(q^{\prime 2}-q^{2})}\,.$ (27)

The imaginary part is positive infinitesimally above the cut, and negative
infinitesimally below the cut. One has the simple representation

${\rm Im}\,\Pi^{(1)}_{\rm
R}(q^{2}+\mathrm{i}\epsilon)=\frac{\alpha}{6}v(3-v^{2})\,\Theta(q^{2}-4)\,.$
(28)

On the cut, the $v$ variable runs from $v=0$ to $v=1$. Here, $\Theta$ denotes
the Heaviside step function. The above result given in Eq. (26) in the $v$
representation fulfills [see also Eq. (5)]

$\Pi^{(1)}_{\rm
R}(q^{2})=\frac{\alpha}{\pi}\frac{q^{2}}{15}+{\mathcal{O}}(q^{4})\,,$ (29)

for $q^{2}=(q^{0})^{2}-\vec{q}^{\,2}$.

##  III Two–Loop Vacuum Polarization

###  III.1 Orientation

We focus on the diagrams in Figs. 1. The first diagram 1 is the proper two-
loop vacuum polarization diagram, while diagram 2 is a one-loop vacuum
polarization diagram with a self-energy insertion. The expression
$\Pi^{(2;1)}_{\mu\nu}(q^{2})$ corresponds to diagram 1,

$\displaystyle\Pi^{(2;1)}_{\mu\nu}(q^{2})=e^{4}\int\dfrac{\mathrm{d}^{D}k_{1}}{(2\pi)^{D}}\dfrac{\mathrm{d}^{D}k_{2}}{(2\pi)^{D}}$
$\displaystyle\;\times\left[-{\rm
Tr}\left(\gamma_{\mu}S(p_{4})\,\gamma^{\rho}\,S(p_{2})\gamma_{\nu}S(-p_{1})\gamma^{\sigma}\,S(-p_{3})\,D_{\rho\sigma}(p_{5})\right)\right]\,,$
(30a) while $\Pi^{(2;2)}_{\mu\nu}(q^{2})$ corresponds to the self-energy
insertion into a fermionic line of the vacuum-polarization loop,
$\displaystyle\Pi^{(2;2)}_{\mu\nu}(q^{2})=e^{4}\int\dfrac{\mathrm{d}^{D}k_{1}}{(2\pi)^{D}}\dfrac{\mathrm{d}^{D}k_{2}}{(2\pi)^{D}}$
$\displaystyle\;\left[-{\rm
Tr}\left(\gamma_{\mu}\,S(p_{4})\,\gamma^{\rho}\,S(p_{2})\,\gamma^{\sigma}\,S(p_{4})\,\gamma_{\nu}(-p_{3})\,D_{\rho\sigma}(p_{5})\right)\right]\,.$
(30b)

Here, the photon propagator is given by [see Eq. (9.133) of Ref. JeAd2022book
]

$D_{\mu\nu}(k)=-\frac{1}{k^{2}+\mathrm{i}\epsilon}\left(g_{\rho\sigma}-\lambda\frac{k_{\rho}k_{\sigma}}{k^{2}}\right)\,,$
(31)

and the momenta and denominators are

$\displaystyle p_{1}=$ $\displaystyle\;k_{1}\,,\qquad p_{2}=q-k_{1}\,,\qquad
p_{3}=k_{2}\,,$ (32a) $\displaystyle p_{4}=$ $\displaystyle\;q-k_{2}\,,\qquad
p_{5}=k_{1}-k_{2}\,,$ (32b) $\displaystyle D_{j}=$
$\displaystyle\;-p_{j}^{2}+1-\mathrm{i}\epsilon\,,\qquad j=1,\ldots,4\,,$
(32c) $\displaystyle D_{5}=$ $\displaystyle\;-p_{5}^{2}-\mathrm{i}\epsilon\,.$
(32d)

These conventions are adapted to FORM FORM ; KuUeVeVo2013 ; in an East–Coast
metric with $g_{\mu\nu}=\mathrm{diag}(-1,-1,-1,1)$, one would replace
$p_{j}^{2}\to-p_{j}^{2}$.

Figure 1:  The diagrams concern the electronic two-loop vacuum-polarization
diagrams in muonic hydrogen (the negatively charged muon line is denoted by
$\mu^{-}$). These are are naturally divided into diagram (1), which is a
proper two-loop diagram, diagram (2), which is a one-loop diagram with a self-
energy insertion, and the loop-after-loop (reducible) diagram (3).

###  III.2 Reduction to Master Integrals

After the calculation of traces, performed with computer algebra (using FORM,
see Refs. FORM ; KuUeVeVo2013 ), and a Wick rotation, the scalar contribution
$\Pi^{(1)}(q^{2})$ and $\Pi^{(2)}(q^{2})$ of the two diagrams are reduced to a
combination of 24 and 41 Feynman integrals, respectively. In order to further
reduce the expressions, we use the integration by parts (IBP) identities
ChTk1981npb ; Tk1981plb . A system of IBP identities is built and solved with
the program $\tt{SYS}$, which is based on an algorithm described in Refs.
La2000ijmpa ; La2001dc . These identities allow to reduce $\Pi^{(1)}(q^{2})$
and $\Pi^{(2)}(q^{2})$ to a linear combination of five irreducible master
integrals.

The choice of the master integrals is somewhat arbitrary; we chose the
following ones:

$\displaystyle M_{1234}=$
$\displaystyle\;[N(\epsilon)]^{2}\int\frac{[\mathrm{d}^{D}k_{1}][\mathrm{d}^{D}k_{2}]}{D_{1}D_{2}D_{3}D_{4}}\,,$
(33a) $\displaystyle M_{234}=$
$\displaystyle\;[N(\epsilon)]^{2}\int\frac{[\mathrm{d}^{D}k_{1}][\mathrm{d}^{D}k_{2}]}{D_{2}D_{3}D_{4}}\,,$
(33b) $\displaystyle M_{235}=$
$\displaystyle\;[N(\epsilon)]^{2}\int\frac{[\mathrm{d}^{D}k_{1}][\mathrm{d}^{D}k_{2}]}{D_{2}D_{3}D_{5}}\,,$
(33c) $\displaystyle M_{235k}=$
$\displaystyle\;[N(\epsilon)]^{2}\int\frac{[\mathrm{d}^{D}k_{1}][\mathrm{d}^{D}k_{2}]\;q\cdot
k_{2}}{D_{2}D_{3}D_{5}}\,,$ (33d) $\displaystyle M_{24}=$
$\displaystyle\;[N(\epsilon)]^{2}\int\frac{[\mathrm{d}^{D}k_{1}][\mathrm{d}^{D}k_{2}]}{D_{2}D_{4}}\,.$
(33e)

The following master integrals are reducible and factorize into a product of
two one-loop master integrals,

$\displaystyle M_{24}(q^{2})=$ $\displaystyle\;[M_{2}(q^{2})]^{2}\,,$ (34a)
$\displaystyle M_{234}(q^{2})=$ $\displaystyle\;M_{2}(q^{2})M_{12}(q^{2})\,,$
(34b) $\displaystyle M_{1234}(q^{2})=$ $\displaystyle\;[M_{12}(q^{2})]^{2}\,.$
(34c)

The expressions for the contribution of diagrams 1 and 2 are given as follows,

$\displaystyle\Pi^{(2;1)}(q^{2})=$
$\displaystyle\;[C(\epsilon)]^{2}\left(\dfrac{\alpha}{\pi}\right)^{2}\left(P^{(2;1)}(q^{2})+\lambda
G^{(2;1)}(q^{2})\right)\,,$ (35a) $\displaystyle\Pi^{(2;2)}(q^{2})=$
$\displaystyle\;[C(\epsilon)]^{2}\left(\dfrac{\alpha}{\pi}\right)^{2}\left(P^{(2;2)}(q^{2})+\lambda
G^{(2;2)}(q^{2})\right)\,,$ (35b)

where $\lambda$ is the gauge parameter [see Eq. (31)], and

$P^{(2;1)}(q^{2})={\mathcal{P}}^{(2;1)}_{24}\,M_{24}(q^{2})+{\mathcal{P}}^{(2;1)}_{234}\,M_{234}(q^{2})\\\
+{\mathcal{P}}^{(2;1)}_{1234}\,M_{1234}(q^{2})+{\mathcal{P}}^{(2;1)}_{235}\,M_{235}(q^{2})\\\
+{\mathcal{P}}^{(2;1)}_{235k}\,M_{235k}(q^{2})\,.$ (36)

The coefficient functions ${\mathcal{P}}^{(2;1)}_{24}$,
${\mathcal{P}}^{(2;1)}_{234}$, ${\mathcal{P}}^{(2;1)}_{1234}$,
${\mathcal{P}}^{(2;1)}_{235}$ and ${\mathcal{P}}^{(2;1)}_{235k}$ are given in
Appendix B. The gauge-dependent part of diagram 1 is given as

$G^{(2;1)}(q^{2})={\mathcal{G}}^{(2;1)}_{24}\,M_{24}(q^{2})+{\mathcal{G}}^{(2;1)}_{235}\,M_{235}(q^{2})\\\
+{\mathcal{G}}^{(2;1)}_{235k}\,M_{235k}(q^{2})\,,$ (37)

with the coefficient functions ${\mathcal{G}}^{(2;1)}_{24}$,
${\mathcal{G}}^{(2;1)}_{235}$ and ${\mathcal{G}}^{(2;1)}_{235k}$ also being
indicated in Appendix B. For diagram 2, one has the following reduction,

$P^{(2;2)}(q^{2})={\mathcal{P}}^{(2;2)}_{234}\,M_{234}(q^{2})+{\mathcal{P}}^{(2;2)}_{235k}\,M_{235k}(q^{2})\\\
+{\mathcal{P}}^{(2;2)}_{235}\,M_{235}(q^{2})+{\mathcal{P}}^{(2;2)}_{24}\,M_{24}(q^{2})\,.$
(38)

Again, the coefficient functions ${\mathcal{P}}^{(2;2)}_{24}$,
${\mathcal{P}}^{(2;2)}_{234}$, ${\mathcal{P}}^{(2;2)}_{1234}$,
${\mathcal{P}}^{(2;2)}_{235}$, and ${\mathcal{P}}^{(2;2)}_{235k}$ are given in
Appendix B. The gauge-dependent part of the second diagram cancels the gauge-
dependent part of the first diagram, in view of the relation
$G^{(2;2)}(q^{2})=-\tfrac{1}{2}G^{(2;1)}(q^{2})$.

###  III.3 Differential Equations and Master Integrals

A system of differential equations in $q^{2}$ is being created; we need to
compute the 5 master integrals. Three of them are reducible into products of
one-loop integrals, as outlined in Eq.(34. We obtain a system of differential
equations in $q^{2}$ satisfied by the irreducible master integrals Ko1991plb ;
Re1997ny ; GeRe2000 . The first of these equations is of second order,

$\displaystyle
q^{2}\left[-3\left((D-4)q^{2}+4\right)M_{235}^{\prime}(q^{2})\right.$
$\displaystyle\left.+2q^{2}\left(q^{2}-4\right)M_{235}^{\prime\prime}(q^{2})\right]+(D-2)^{2}M_{24}(q^{2})$
$\displaystyle+(D-3)\left[D\left(q^{2}+2\right)-4\left(q^{2}+1\right)\right]\,M_{235}(q^{2})=0\,,$
(39a) while the remaining one is a first-order differential equation,
$\displaystyle\left[D\left(q^{2}+2\right)-q^{2}-6\right]\,M_{235}(q^{2})+6(D-2)\,M_{235k}(q^{2})$
$\displaystyle+(D-2)\,M_{24}(q^{2})+q^{2}\,\left(q^{2}-4\right)\,M_{235}^{\prime}(q^{2})=0\,,$
(39b)

where the differentiation is with respect to $q^{2}$. The initial conditions
at $q^{2}=0$ for these diagrams are simple vacuum integrals.

###  III.4 Results for the Master Integrals

Due to the appearance of the factor $D-4$ in the denominators of Eqs.
(36)—(38) (for the coefficients, see Appendix B), the master integrals have to
be evaluated up to order ${\mathcal{O}}(\epsilon)$. It is not sufficient to
keep only the finite part. The results are as follows. For $M_{235}$, one
finds the following relation, where we note the prefactor $(1-v^{2})$,

$\displaystyle(1$
$\displaystyle-v^{2})M_{235}(q^{2})=\frac{v^{2}-1}{\epsilon^{2}}+\frac{3v^{2}-2}{\epsilon}+\biggl{[}6v^{2}+\frac{1}{2}$
$\displaystyle+\frac{1}{4}\left(-v^{4}-2v^{2}+3\right)\ln^{2}\left(\frac{v-1}{v+1}\right)$
$\displaystyle-\left(v^{2}-3\right)v\ln\left(\frac{v-1}{v+1}\right)\biggr{]}+\epsilon
N_{235}(q^{2})+{\mathcal{O}}\left(\epsilon^{2}\right)\,.$ (40)

The order-$\epsilon$ coefficient $N_{235}(q^{2})$ is found in Appendix C. For
$M_{235k}$, one finds the relation

$\displaystyle(1$
$\displaystyle-v^{2})^{2}M_{235k}(q^{2})=\frac{1-v^{2}}{\epsilon^{2}}+\frac{1}{\epsilon}\biggl{[}\frac{13}{6}-\frac{7v^{2}}{2}\biggr{]}$
$\displaystyle+\frac{1}{8}\left(-v^{6}+3v^{4}+5v^{2}-7\right)\ln^{2}\left(\frac{v-1}{v+1}\right)$
$\displaystyle-\frac{1}{6}v\left(3v^{4}-8v^{2}+21\right)\ln\left(\frac{v-1}{v+1}\right)$
$\displaystyle+\frac{1}{12}\left(-6v^{4}-91v^{2}-7\right)+\epsilon
N_{235k}(q^{2})+{\mathcal{O}}\left(\epsilon^{2}\right)\,.$ (41)

Again, $N_{235k}(q^{2})$ can be found in Appendix C. Similar master integrals
have recently been considered in Refs. ScTa1993 ; AnEtAl2007 .

###  III.5 Renormalization

We now substitute the master integrals (34a), (34c), (34b), (III.4) and
(III.4) into the reduction formulas given in Eqs. (36) and (38). We can thus
find the series expansions of the vacuum-polarization functions
$\Pi^{(2;1)}(q^{2})$ and $\Pi^{(2;2)}(q^{2})$. For the renormalization, we
need the values at $q^{2}=0$. The behavior, for small $q^{2}$, is the
following:

$\displaystyle P^{(2;1)}(q^{2})$
$\displaystyle=-\frac{12\epsilon^{5}+4\epsilon^{4}-23\epsilon^{3}-14\epsilon^{2}-21\epsilon-6}{12(\epsilon-2)(\epsilon-1)\epsilon^{2}(2\epsilon+1)(2\epsilon+3)}$
(42)
$\displaystyle-\frac{1}{q^{2}}\frac{1}{(1-\epsilon)(1+2\epsilon)(1-2\epsilon)}+{\mathcal{O}}(q^{2})\
,$ (43) $\displaystyle P^{(2;2)}(0)$
$\displaystyle=-\frac{(\epsilon+1)\left(4\epsilon^{4}-8\epsilon^{3}-8\epsilon^{2}+21\epsilon+3\right)}{12(\epsilon-2)(\epsilon-1)\epsilon^{2}(2\epsilon+1)(2\epsilon+3)}$
(44)
$\displaystyle+\frac{1}{2q^{2}}\frac{1}{(1-\epsilon)(1+2\epsilon)(1-2\epsilon)}+{\mathcal{O}}(q^{2})\
.$ (45)

Terms containing $1/q^{2}$ appear, due to the fact that the single diagrams
are not gauge invariant. Taking the gauge invariant sum
$P^{(2;1)}(q^{2})+2P^{(2;2)}(q^{2})$, the divergent terms cancel out, so that,
taking the limit $q^{2}\to 0$ one gets

$P^{(2;1)}(0)+2P^{(2;2)}(0)=\frac{10\epsilon^{2}-7\epsilon-9}{12(\epsilon-2)\epsilon(2\epsilon+1)}\,.$
(46)

In order to carry out the renormalization, we need the renormalization
constants

$Z_{J}=1+\sum_{i}Z_{J}^{(L)}\left(\frac{\alpha}{\pi}\right)^{L}\left(\frac{e^{2}\,C(\epsilon)}{4\pi\alpha}\right)^{L}\,,$
(47)

where $J$ stands for the vertex renormalization ($J=1$), the wave function
renormalization ($J=2$), and the mass renormalization ($J=m$), and $L$
enumerates the loop order. We define the one-loop quantities

$F=-Z_{1}^{(1)}\,,\qquad S_{0}=\delta m=-Z_{m}^{(1)}\,,\qquad
S_{1}=-Z_{2}^{(1)}\,,$ (48)

and obtain the results

$F=\frac{3-2\epsilon}{4\epsilon-8\epsilon^{2}}\,,\quad
S_{0}=\frac{3-2\epsilon}{4\epsilon-8\epsilon^{2}}\,,\quad S_{1}=S_{0}\,.$ (49)

We note that $Z_{1}^{(1)}=Z_{2}^{(1)}$ because of the Ward–Takahashi identity.
The identity $Z_{2}^{(1)}=Z_{m}^{(1)}$ is a simple coincidence, which occurs
only at one-loop order in dimensional regularization, but to all orders in
$\epsilon$. The constants $F$, $S_{0}$ and $S_{1}$ in Eq. (48) are defined
with the prefactor $(-1)$ so that the expansion of the the self-energy
insertion into the fermion line about the mass shell acquires positive terms,
for $\cancel{q}\approx m=1$:

$\displaystyle\frac{1}{\cancel{q}-m}\Sigma(q)\frac{1}{\cancel{q}-m}=$
$\displaystyle\;\frac{1}{\cancel{q}-m}[S_{0}+(\cancel{q}-m)S_{1}]\frac{1}{\cancel{q}-m}$
$\displaystyle=$
$\displaystyle\;\frac{S_{1}}{\cancel{q}-m}+\frac{S_{0}}{(\cancel{q}-m)^{2}}\,.$
(50)

Finally, the renormalized contributions of the diagrams are found as follows,

$\displaystyle P^{(2;1)}_{\rm R}(q^{2})=$
$\displaystyle\;P^{(2;1)}(q^{2})-P^{(2;1)}(0)-2FP^{(1)}_{\rm R}(q^{2})\,,$
(51) $\displaystyle P^{(2;2)}_{\rm R}(q^{2})=$
$\displaystyle\;P^{(2;2)}(q^{2})-P^{(2;2)}(0)$
$\displaystyle-(S_{1}P^{(1)}_{\rm R}(q^{2})+S_{0}P^{(1a)}_{\rm
R}(q^{2}))(-1)\,,$ (52)

and the renormalized two-loop function $P^{(2)}_{\rm R}(q^{2})$ is obtained as

$P^{(2)}_{\rm R}(q^{2})=P^{(2;1)}_{\rm R}(q^{2})+2P^{(2;2)}_{\rm R}(q^{2})\,.$
(53)

###  III.6 Results

We recall the definition of $v$ in Eq. (19). In the limit $\epsilon\to 0^{+}$,
one obtains the two-loop vacuum-polarization function

$\Pi^{(2)}_{\rm R}(q^{2})=\left(\dfrac{\alpha}{\pi}\right)^{2}\,P^{(2)}_{\rm
R}(q^{2})\,.$ (54)

We find a compact representation (see also Ref. LaJe2024suppl )

$\Pi^{(2)}_{\rm
R}(q^{2})=\left(\dfrac{\alpha}{\pi}\right)^{2}\biggl{\\{}\frac{v^{4}-2v^{2}-3}{12}\biggl{[}\Phi_{1}\left(\frac{v-1}{v+1}\right)\\\
\times\ln^{2}\left(\frac{v-1}{v+1}\right)-4\Phi_{2}\left(\frac{v-1}{v+1}\right)\ln\left(\frac{v-1}{v+1}\right)+6\Phi_{3}\left(\frac{v-1}{v+1}\right)\\\
+3\zeta(3)\biggr{]}+\frac{3v-v^{3}}{12}\biggl{[}4\Phi_{1}\left(\frac{v-1}{v+1}\right)\ln\left(\frac{v-1}{v+1}\right)\\\
-4\Phi_{2}\left(\frac{v-1}{v+1}\right)+3\ln^{2}\left(\frac{v-1}{v+1}\right)\biggr{]}+\frac{5v-3v^{3}}{8}\ln\left(\frac{v-1}{v+1}\right)\\\
+\frac{1}{96}\left(7v^{4}-22v^{2}-33\right)\ln^{2}\left(\frac{v-1}{v+1}\right)-\frac{13(v^{2}-1)}{24}+\frac{5}{24}\biggr{\\}}\,,$
(55)

where

$\Phi_{n}(x)={\rm Li}_{n}(x)+2{\rm Li}_{n}(-x)\ ,\quad{\rm
Li}_{1}(x)=-\ln(1-x)\ .$ (56)

For small positive $q^{2}$, one obtains the expression

$\Pi^{(2)}_{\rm
R}(q^{2})=\left(\dfrac{\alpha}{\pi}\right)^{2}\,\frac{41q^{2}}{162}+{\ldots}\,.$
(57)

The representation for the imaginary part just above the cut is even more
compact,

${\rm Im}\,\Pi^{(2)}_{\rm
R}(q^{2}+\mathrm{i}\epsilon)=\frac{\alpha^{2}}{\pi}\Theta(q^{2}-4)\\\
\times\biggl{\\{}\frac{v^{4}-2v^{2}-3}{6}\biggl{[}\Phi_{1}\left(\frac{v-1}{v+1}\right)\ln\left(\frac{1-v}{1+v}\right)\\\
-2\Phi_{2}\left(\frac{v-1}{v+1}\right)\biggr{]}+\frac{3v-v^{3}}{6}\biggl{[}2\Phi_{1}\left(\frac{v-1}{v+1}\right)+3\ln\left(\frac{1-v}{1+v}\right)\biggr{]}\\\
+\frac{1}{48}\left(7v^{4}-22v^{2}-33\right)\ln\left(\frac{1-v}{1+v}\right)+\frac{5v-3v^{3}}{8}\biggr{\\}}\,.$
(58)

###  III.7 Comparison with the Literature

An essential ingredient of our calculations is the master integrals, which
should be compared to results communicated in Refs. Da2000 ; DaKa2001 ;
DaKa2004 . In Refs. Da2000 ; DaKa2001 , the authors define a master integral
$J_{011}$ which is equivalent to our $M_{12}$ up a factor

$J_{011}=\mathrm{i}\,\pi^{2-\epsilon}\,\Gamma(1+\epsilon)M_{12}\,.$ (59)

In Eq.(1) of Ref. Da2000 , the general term of the $\epsilon$-expansion of
$J_{011}$ is given in terms of log-sine integrals which can be written in
terms of Nielsen polylogarithms; the analytic continuation is shown in Eqs.
(2.9)—(2.23) of Ref. DaKa2001 . Moreover, an closed analytical expression for
$J_{011}$ is given in Eqs. (2.10)—(2.14) of Ref. DaKa2001 , containing
hypergeometric ${}_{2}F_{1}$ functions.

Changing the variable from $t=q^{2}$ to $v$ and using the transformations for
Nielsen polylogarithm $S_{n,p}(-1/z)\to S_{n,p}(-z)$, we found that the
results in our Eqs. (II.2) and (A) agree perfectly with the corresponding
terms of Eq. (1) of Ref. Da2000  and the expansion given in Eqs. (2.10)—(2.14)
of Ref. DaKa2001 .

In Eqs. (4.9), (4.10), and (E.4) of Ref. DaKa2004 , analytical results for
master integrals equivalent to our $M_{235}$ and $M_{235k}$ were presented.
These master integrals are related to ours by the following relations: one
finds for the ultraviolet-convergent integral (4.9) of Ref. DaKa2004 ,

$4q^{2}(q^{2}-4)^{2}J_{011}(1,2,2)=-24(D-3)(D-2)\left(q^{2}+2\right)\\\ \times
M_{235k}-(D-2)\left((D-2)q^{4}-2Dq^{2}+16D-40\right)M_{24}\\\
-4(D-3)\left(2(D-2)q^{4}-(D-6)q^{2}+8D-20\right)M_{235}\ ,$ (60)

while, for the ultraviolet-convergent linear combination of integrals given in
Eq. (4.10) of Ref. DaKa2004 , one has

$4(D-4)q^{2}(q^{2}-4)^{3}\left(J_{011}(1,2,2)+2J_{011}(2,1,2)\right)=\\\
-(D-2)\biggl{[}(D-4)(D-2)q^{6}+2(D-6)(D-4)q^{4}\\\
+16(D(5D-31)+47)q^{2}-32(2D-7)(2D-5)\biggr{]}M_{24}\\\
-4(D-3)\biggl{[}4(D-4)(D-2)q^{6}+3D(3D-10)q^{4}\\\
+24(D-4)(D-3)q^{2}-16(2D-7)(2D-5)\biggr{]}M_{235}\\\
-48(D-3)(D-2)\left(q^{2}-1\right)\\\
\times\left(D\left(q^{2}+8\right)-4\left(q^{2}+7\right)\right)M_{235k}\,,$
(61)

and for the ultraviolet divergent integral (E.4) of Ref. DaKa2004 ,

$J_{011}(1,1,1)=M_{235}\,.$ (62)

The $\epsilon-$expansion obtained from those of $M_{235}$ and $M_{235k}$ agree
perfectly with those presented in Ref. DaKa2004 .

Our results for the complete two-loop function given in Eq. (55) and for its
imaginary part, as given in Eq. (58), agree with the calculation of Ref.
KaSa1955 . One notes that Eq. (49) of Ref. KaSa1955  contains the imaginary
part, while Eq. (57) of Ref. KaSa1955  also contains the real part of the two-
loop vacuum polarization. One should note, though, that some integrals are
left unevaluated in Ref. KaSa1955 . Furthermore, it needs to be pointed out
that the expressions in Ref. KaSa1955  also contain the reducible diagram (see
diagram 3 in Fig. 1) with two iterated one-loop vacuum polarizations on the
same line. Our result for the imaginary part also agrees with the result given
in Eq. (5-4.200) of Ref. Sc1970vol3 . Our results are also in agreement with
the paper BrFlTa1993 , specifically, with the results given in Eqs. (71),
(78), (79), (80), and (81) of Ref. BrFlTa1993 . Furthermore, we can refer to
related investigations on the vector part of the quantum chromodynamic vacuum-
polarization tensor Kn1990 ; DjGa1994 ; ChKuSt1996 .

In Ref. BaRe1973 , the two-loop vacuum polarization was also calculated. The
result given in Eq. (49) of Ref. BaRe1973  is expressed in terms of the
variable $\theta$, which in our notation reads as $\theta=-\frac{1-v}{1+v}$.
The last-but-one line of Eq. (49) of Ref. BaRe1973  contains two terms, which,
in the notation of Ref. BaRe1973 , read as $2\ln(\theta){\rm Li}_{2}(\theta)$
and $2\ln(\theta){\rm Li}_{2}(-\theta)$. The term $2\ln(\theta){\rm
Li}_{2}(-\theta)$ contains a misprint; it should have read as
$4\ln(\theta){\rm Li}_{2}(-\theta)$. Formula (49) of Ref. BaRe1973  was
referred to in Ref. Bo1975zpa , and it was stated in Ref. Bo1975zpa  that Eq.
(49) of Ref. BaRe1973  has a misprint. Indeed, Eq. (49) of Ref. BaRe1973  is
rewritten as Eq. (6) of Ref. Bo1975zpa , with a change of notation: The
variable $\delta$ used in Ref. Bo1975zpa  is related to the variable $\theta$
used in Eq. (49) of Ref. BaRe1973  by the relation
$\delta=(1+\theta)/(1-\theta)=v$. In Eq. (6) of Ref. Bo1975zpa , the wrong
term $2\ln(\theta){\rm Li}_{2}(-\theta)$ is corrected to the right
$4\ln(\theta){\rm Li}_{2}(-\theta)$, but unfortunately new misprints are
inserted. The terms (in the notation adapted in Ref. BaRe1973 )

$\displaystyle T_{1}=$ $\displaystyle\;-(\delta(5-3\delta^{2})/8)\ln\theta\,,$
(63a) $\displaystyle T_{2}=$
$\displaystyle\;+\tfrac{1}{4}\ln^{2}\theta\left[\frac{33+22\delta^{2}-7\delta^{4}}{24}-\delta(3-\delta^{2})\right]\,,$
(63b)

have the wrong signs and should have read

$\displaystyle T_{1}\to$
$\displaystyle\;{\widetilde{T}}_{1}=+(\delta(5-3\delta^{2})/8)\ln\theta\,,$
(64a) $\displaystyle T_{2}\to$
$\displaystyle\;{\widetilde{T}}_{2}=-\tfrac{1}{4}\ln^{2}\theta\left[\frac{33+22\delta^{2}-7\delta^{4}}{24}-\delta(3-\delta^{2})\right]\,.$
(64b)

In Ref. BoRi1982 , Eq. (6) of Ref. Bo1975zpa  is rewritten as Eq. (140) of
Ref. BoRi1982 , but the two terms involving $\ln(\theta)$’s keep the wrong
sign.

In Ref. Pa1993pra , the work from Ref. KaSa1955  is being cited for the
imaginary part of the proper two-loop vacuum polarization, and it is expressed
in terms of the variable $d$, which equals our $v$. However, Eq. (16) of Ref.
Pa1993pra  contains a misprint. Namely, the term $+3/2\ln((1+d)/(1-d))$ (in
Ref. Pa1993pra ) has a logarithm missing and should be read
$+3/2\ln((1+d)/(1-d))\ln((1+d)/2)$. We should clarify that, despite the
typographical error in Eq. (16) of Ref. Pa1993pra , final results for the
energy shifts due to the two-loop vacuum polarization obtained in Ref.
Pa1993pra  are in agreement with those reported here.

##  IV Contributions to Energy Shifts

###  IV.1 Muonic Hydrogen

Armed with the compact expressions given in Eqs. (55) and (58) for the
imaginary part of the two-loop vacuum-polarization function, we can evaluate
energy shifts in hydrogenlike ions. In Chap. 10 of Ref. JeAd2022book , the
calculation of the vacuum-polarization energy shift is outlined in detail;
specifically, we refer to Eqs. (10.244) and Eq. (10.245) of Ref. JeAd2022book
. (We note the different sign conventions for the sign of the imaginary part
of the renormalized scalar function $\Pi_{\rm R}$, as compared to Ref.
JeAd2022book .) From the one-loop and two-loop scalar functions, one infers
the one-loop and two-loop irreducible vacuum-polarization potentials
$V^{(1)}_{\rm R}(r)$ and $V^{(2)}_{\rm R}(r)$, respectively,

$\displaystyle V^{(1)}_{\rm R}(r)=$
$\displaystyle\;-\frac{Z\alpha}{\pi}\int\limits_{4}^{\infty}\frac{\mathrm{d}(q^{2})}{q^{2}}\;\frac{\mathrm{e}^{-qr}}{r}\;\mathrm{Im}\\!\left[\Pi^{(1)}_{\rm
R}(q^{2}+\mathrm{i}\epsilon)\right]\,.$ (65) $\displaystyle V^{(2)}_{\rm
R}(r)=$
$\displaystyle\;-\frac{Z\alpha}{\pi}\int\limits_{4}^{\infty}\frac{\mathrm{d}(q^{2})}{q^{2}}\;\frac{\mathrm{e}^{-qr}}{r}\;\mathrm{Im}\\!\left[\Pi^{(2)}_{\rm
R}(q^{2}+\mathrm{i}\epsilon)\right]\,.$ (66)

Of particular phenomenological importance for the proton radius puzzle
Je2022proton  is the contribution of the irreducible two-loop vacuum-
polarization effect on the $2P$–$2S$ energy interval in muonic hydrogen. We
use the following results,

$\displaystyle\langle 2S|\frac{\exp(-q\,r)}{r}|2S\rangle=$
$\displaystyle\;\frac{1+2(\beta q)^{2}}{4\beta\,(1+\beta q)^{4}}\,,$ (67a)
$\displaystyle\langle 2P|\frac{\exp(-q\,r)}{r}|2P\rangle=$
$\displaystyle\;\frac{1}{4\beta\,(1+\beta q)^{4}}\,.$ (67b)

Here, the $\beta$ parameter is

$\beta=\frac{1}{Z\alpha\mu}\,,$ (68)

where $Z$ is the nuclear charge number, $\alpha$ is the fine-structure
constant, and $\mu$ is the reduced mass of the system. For muonic hydrogen,
one has $Z=1$. In units where the electron mass is $m=1$ and the muon mass is
$m_{\mu}$, one has $\mu=m_{p}/(1+m_{p})$ for ordinary hydrogen and
$\mu=m_{\mu}m_{p}/(m_{\mu}+m_{p})$ for muonic hydrogen.

The contribution of the irreducible two-loop vacuum polarization diagrams to
the $2P$–$2S$ splitting in muonic hydrogen is obtained as

$E^{(2)}(2P-2S)=\left(\dfrac{\alpha}{\pi}\right)^{2}(Z\alpha)^{2}\mu
f^{(2)}(\beta)\,,\\\
f^{(2)}(\beta)=\int_{4}^{{\infty}}\frac{\mathrm{d}(q^{2})}{q^{2}}\frac{{\rm
Im}\,\Pi^{(2)}_{\rm R}(q^{2}+\mathrm{i}\epsilon)}{\alpha^{2}/\pi}\frac{(\beta
q)^{2}}{2(1+\beta q)^{4}}\,.$ (69)

Quite surprisingly, the two-loop energy shift can be expressed analytically.
We define

$\beta=\frac{\sqrt{1-w^{2}}}{2}\ ,\qquad w=\sqrt{1-4\beta^{2}}\ ,$ (70)

and find the result

$\displaystyle f^{(2)}$
$\displaystyle(w)=\frac{1}{288w^{6}}\biggl{\\{}-480\left(w^{2}-1\right)^{2}w^{6}I_{1}\left(\sqrt{\frac{1-w}{1+w}}\right)$
$\displaystyle+128\left(15w^{6}-25w^{4}+8w^{2}+3\right)w^{3}I_{2}\left(\sqrt{\frac{1-w}{1+w}}\right)$
$\displaystyle+404w^{8}-312w^{6}-176w^{4}+96w^{2}$
$\displaystyle-32\pi\sqrt{1-w^{2}}\left(15w^{4}-19w^{2}+3\right)w^{4}+4\left(w^{2}-1\right)\times$
$\displaystyle\left(41w^{6}+58w^{4}-10w^{2}-60\right)w\ln\left(\frac{1-w}{1+w}\right)$
$\displaystyle-16\pi^{2}\left(1-w^{2}\right)^{3/2}\left(15w^{6}-5w^{4}-2w^{2}+6\right)$
$\displaystyle-\left(199w^{10}-302w^{8}+51w^{6}-4w^{4}+176w^{2}-96\right)\times$
$\displaystyle\left[\ln^{2}\left(\frac{1-w}{1+w}\right)+\pi^{2}\right]\biggr{\\}}\,,$
(71)

where the expressions $I_{j}(a)$ with $j=1,2$ read as follows,

$\displaystyle I_{1}$
$\displaystyle(a)=\frac{1}{4}\biggl{[}6\Phi_{3}(-a^{2})-4\Phi_{2}(-a^{2})\ln(a^{2})$
$\displaystyle+\Phi_{1}(-a^{2})\ln^{2}(a^{2})-4\mathrm{i}\pi\left({\rm
Li}_{2}(\mathrm{i}a)-{\rm Li}_{2}(-\mathrm{i}a)\right)$
$\displaystyle-\pi^{2}\ln(1+a^{2})-4\pi^{2}\ln(1+a)+3\zeta(3)\biggr{]}\,,$
(72) $\displaystyle I_{2}$
$\displaystyle(a)=\frac{1}{16}\biggl{[}-4\Phi_{2}(-a^{2})+4\Phi_{1}(-a^{2})\ln(a^{2})$
$\displaystyle+3\left(\ln^{2}(a^{2})+\pi^{2}\right)-8\pi\arctan(a)\biggr{]}\,.$
(73)

Some of the terms could also be expressed in terms of the Legendre $\chi$
function, $\chi_{\nu}(z)=\tfrac{1}{2}\,[{\rm Li}_{\nu}(z)-{\rm
Li}_{\nu}(-z)]$.

For muonic hydrogen ($Z=1$), one obtains for the nonrelativistic energy shift
due to the irreducible two-loop diagrams, using CODATA values for the relevant
physical constants MoTaNe2008 ; TiMoNeTa2021  (see also Ref. LaJe2024suppl ),

$\left.E^{(2)}(2P-2S)\right|_{\mbox{$\mu$H}}=1.25298\,{\rm meV}\,.$ (74)

This is consistent with the literature (e.g., Ref. Je2011aop1 ). The result
(74) confirms the numerical value of a contribution to the proton radius
puzzle in muonic hydrogen whose numerical magnitude could have potentially
contributed to an explanation the discrepancy. The confirmation is obtained
based on modern quantum-field theoretical methods.

The reducible diagram, according to Eqs. (114) and (115), leads to an energy
shift of

$\left.E^{({\rm(r)})}(2P-2S)\right|_{\mbox{$\mu$H}}=0.25495\,{\rm meV}\,.$
(75)

For selected individual low-lying atomic reference states, an analytic
integration of the one-loop and two-loop energy shifts, for arbitrary $\beta$,
still in the nonrelativistic approximation, is presented in Appendix D.

Table 1: Coefficients from the reducible diagram are indicated for $S$ and $P$ states. | 1S | 2S | $2P_{1/2}$ | $2P_{3/2}$
---|---|---|---|---
$B^{{\rm(r)}}_{50}$ | $-\dfrac{23\pi}{567}$ | $-\dfrac{23\pi}{567}$ | $0$ | $0$
$B^{{\rm(r)}}_{60}$ | $\dfrac{289}{4050}$ | $\dfrac{181}{4050}$ | $\dfrac{2}{225}$ | $\dfrac{2}{225}$
Table 2: Coefficients from the irreducible diagrams are indicated for $S$ and $P$ states. | 1S | 2S |  $2P_{1/2}$ | $2P_{3/2}$
---|---|---|---|---
$B_{40}^{(2)}$ |  $-\dfrac{82}{81}$ |  $-\dfrac{82}{81}$ | 0 | 0
$B_{50}^{(2)}$ |  $\dfrac{15647\pi}{13230}-\dfrac{25\pi^{2}}{63}+\dfrac{52}{63}\pi\ln(2)$ |  $\dfrac{15647\pi}{13230}-\dfrac{25\pi^{2}}{63}+\dfrac{52}{63}\pi\ln(2)$ |  $0$ | $0$
$B_{60}^{(2)}$ |  $-\dfrac{82}{81}\ln\left(\dfrac{1}{2}\right)+\dfrac{77\zeta(3)}{72}-\dfrac{358}{75}$ |  $\dfrac{77\zeta(3)}{72}-\dfrac{12727}{2700}$ |  $-\dfrac{737}{2700}$ | $-\dfrac{449}{5400}$
$B_{61}^{(2)}$ |  $-\dfrac{82}{162}$ |  $-\dfrac{82}{162}$ | 0 | 0

###  IV.2 Relativistic Ordinary Hydrogen

Given our analytic approach, we may also also study the energy shifts to low-
lying $S$ and $P$ states for ordinary (electronic) hydrogen, with a
relativistic reference state. We focus on higher-order terms in the semi-
analytic expansion in powers and logarithms of $Z\alpha$, for relativistic
Dirac-Coulomb reference states (see Chap. 8 of Ref. JeAd2022book ).

The potential $V^{{\rm(r)}}_{\rm R}$ due to the reducible diagram is

$V^{{\rm(r)}}_{\rm
R}(r)=\frac{Z\alpha}{\pi}\int\limits_{4}^{\infty}\frac{\mathrm{d}(q^{2})}{q^{2}}\;\frac{\mathrm{e}^{-qr}}{r}\;\mathrm{Im}\left\\{\\!\left[\Pi^{(1)}_{\rm
R}(q^{2}+\mathrm{i}\epsilon)\right]^{2}\right\\}\,.$ (76)

We study the energy shifts

$\Delta E^{{\rm(r)}}(nS)=\langle\psi|V^{{\rm(r)}}_{\rm R}|\psi\rangle\,,$ (77)

where $|\psi\rangle$ is the relativistic reference-state wave function, and
$\langle\psi$ stands for the Hermitian adjoint (not for the Dirac adjoint
$\overline{\psi}=\psi^{\dagger}\,\gamma^{0}$). The relevant semi-analytic
expansion is as follows,

$\Delta E^{{\rm(r)}}_{\rm
NR}=\left(\dfrac{\alpha}{\pi}\right)^{2}\frac{(Z\alpha)^{4}}{n^{3}}\biggl{\\{}(Z\alpha)B_{50}^{{\rm(r)}}+(Z\alpha)^{2}B_{60}^{(r)}\\\
+(Z\alpha)^{3}\left[B_{71}^{(r)}\,\ln[(Z\alpha)^{-2}]+B_{70}^{(r)}\right]\biggr{\\}}\,.$
(78)

For the $nS$ levels, the dependence on $n$ is given, for the reducible
diagram, by the following formulas (see also Table 1),

$\displaystyle B_{50}^{{\rm(r)}}(nS)=$ $\displaystyle\;-\frac{23\pi}{567}\,,$
(79a) $\displaystyle B_{60}^{{\rm(r)}}(nS)=$
$\displaystyle\;\frac{29}{810}+\frac{8}{225n^{2}}\,,$ (79b) $\displaystyle
B_{71}^{{\rm(r)}}(nS)=$ $\displaystyle\;\tfrac{1}{2}B_{50}^{{\rm(r)}}(nS)\,,$
(79c) $\displaystyle B_{70}^{{\rm(r)}}(nS)=$
$\displaystyle\;-\frac{23\pi}{567}W_{n}-\frac{570029\pi}{7858620}-\frac{46}{567}\pi\ln(2)$
$\displaystyle+\frac{1294\pi}{18711n^{2}}\,,$ (79d) $\displaystyle W_{n}=$
$\displaystyle\;\ln\left(\frac{n}{2}\right)-\Psi(n)-\gamma_{E}+\frac{1}{n}\,.$
(79e)

Here, $\Psi(x)=\Gamma^{\prime}(x)/\Gamma(x)$ is the logarithmic derivative of
the Gamma function. For the $nP_{1/2}$ levels, one has
$B^{\rm(r)}_{50}(nP)=B^{\rm(r)}_{71}(nP)=0$ and the following results for the
other coefficients for the reducible diagram,

$\displaystyle B_{60}^{(r)}(nP_{1/2})=$
$\displaystyle\frac{8}{675}\frac{n^{2}-1}{n^{2}}\,,$ (80) $\displaystyle
B_{70}^{(r)}(nP_{1/2})=$
$\displaystyle-\frac{983\pi}{93555}\frac{n^{2}-1}{n^{2}}\,.$ (81)

For the $nP_{3/2}$ levels, one has the results

$\displaystyle B_{60}^{(r)}(nP_{3/2})=$
$\displaystyle\;B_{60}^{(r)}(nP_{1/2})\,,$ (82a) $\displaystyle
B_{70}^{(r)}(nP_{3/2})=$
$\displaystyle\;-\frac{127\pi}{17820}\frac{n^{2}-1}{n^{2}}\,.$ (82b)

The two-loop energy shift from the irreducible diagrams, evaluated on the
relativistic wave function, is

$\Delta E^{(2)}=\langle\psi|V^{(2)}_{\rm R}|\psi\rangle\,,$ (83)

where $V^{(2)}_{\rm R}$ is given in Eq. (66). It gives rise to the following
semi-analytic expansion,

$\Delta
E^{(2)}=\left(\dfrac{\alpha}{\pi}\right)^{2}\frac{(Z\alpha)^{4}}{n^{3}}\biggl{\\{}B_{40}^{(2)}+(Z\alpha)B_{50}^{(2)}\\\
+(Z\alpha)^{2}\left[B_{61}^{(2)}\,\ln[(Z\alpha)^{-2}]+B_{60}^{(2)}\right]\biggr{\\}}\\\
+(Z\alpha)^{3}\left[B_{71}^{(2)}\,\ln[(Z\alpha)^{-2}]+B_{70}^{(2)}\right]\biggr{\\}}\,.$
(84)

For the irreducible diagrams, the dependence on $n$ is the following (see also
Table 2):

$\displaystyle B^{(2)}_{40}(nS)=$ $\displaystyle\;-\frac{82}{81}\,,$ (85a)
$\displaystyle B^{(2)}_{50}(nS)=$
$\displaystyle\;\frac{15647\pi}{13230}-\frac{25\pi^{2}}{63}+\frac{52}{63}\pi\ln(2)\,,$
(85b) $\displaystyle B^{(2)}_{60}(nS)=$ $\displaystyle\;B_{40}(nS)\,W_{n}\
+\frac{77\zeta(3)}{72}-\frac{2311}{405}+\frac{1313}{675n^{2}}\,,$ (85c)
$\displaystyle B^{(2)}_{61}(nS)=$
$\displaystyle\;\tfrac{1}{2}B_{40}(nS)\,,\qquad
B^{(2)}_{71}(nS)=\tfrac{1}{2}B_{50}(nS)\,.$ (85d) The Catalan constant
$G=0.915\,965\dots$ enters the $B^{(2)}_{70}$ coefficient, $\displaystyle
B^{(2)}_{70}(nS)=$
$\displaystyle\;C_{50}\,W_{n}+C^{(a)}_{70}+\frac{C^{(b)}_{70}}{n^{2}}\,,$
$\displaystyle C^{(a)}_{70}=$
$\displaystyle\frac{104}{63}\pi\ln^{2}(2)-\frac{25}{63}\pi^{2}\ln(2)-\frac{200\,\pi}{63}\,G+\frac{13\pi^{3}}{189}$
$\displaystyle+\frac{79318\pi\ln(2)}{19845}-\frac{82361\pi^{2}}{79380}+\frac{89515961\pi}{20003760}\,,$
(85e) $\displaystyle C^{(b)}_{70}=$
$\displaystyle-\frac{3304831\pi}{1428840}+\frac{1891\pi^{2}}{2268}-\frac{1076}{567}\pi\ln(2)\,.$
(85f)

For the $nP_{1/2}$ levels, the coefficients are as follows,

$\displaystyle B^{(2)}_{60}(nP_{1/2})=$
$\displaystyle-\frac{737}{2025}\frac{n^{2}-1}{n^{2}}\,,$ (86a) $\displaystyle
B^{(2)}_{70}(nP_{1/2})=$
$\displaystyle-\pi\biggl{[}\frac{1764797}{7144200}-\frac{809\pi}{11340}$
$\displaystyle+\frac{328\ln(2)}{2835}\biggr{]}\frac{n^{2}-1}{n^{2}}\,.$ (86b)

For the $nP_{3/2}$ levels, the coefficients are as follows,

$\displaystyle B^{(2)}_{60}(nP_{3/2})=$
$\displaystyle\;-\frac{449}{4050}\,\frac{n^{2}-1}{n^{2}}\ ,$ (87a)
$\displaystyle B^{(2)}_{70}(nP_{3/2})=$
$\displaystyle-\pi\left[\frac{75763}{510300}-\frac{31\pi}{810}+\frac{19\ln(2)}{405}\right]\frac{n^{2}-1}{n^{2}}\,.$
(87b)

The fine-structure difference

$B^{(2)}_{60}(nP_{3/2})-B^{(2)}_{60}(nP_{1/2})=\frac{41}{162}\frac{n^{2}-1}{n^{2}}$
(88)

is consistent with the result communicated in Eq. (7.7) of Ref. JeCzPa2005 .

##  V Conclusions

We have completed a complete reevaluation of the two-loop vacuum-polarization
tensor in dimensional regularization, on the basis on integration-by-parts
identities. The two-loop vacuum-polarization tensor constitutes a numerically
significant contribution to the Lamb shift of muonic hydrogen which influences
the determination of the proton radius from muonic hydrogen spectroscopy
Je2011aop1 .

In Sec. II, we have discussed the evaluation of the one-loop vacuum
polarization insertion into the photon propagator, evaluated to order
$\epsilon$, and $\epsilon^{2}$, and thus, in a form suitable for inclusion
into higher-order loop calculation where knowledge of the terms of higher
orders in $\epsilon$ is indispensable. We note that we use somewhat
nonstandard conventions for the $\overline{\mathrm{MS}}$-charge, as given in
Eq. (12).

In Sec. III, the irreducible two-loop vacuum polarization insertion has been
evaluated, by expressing it in terms of master integrals. The renormalization
has been carried out, and final results have been presented for the real and
imaginary parts, in Eqs. (55) and (58). A comparison to the existing
literature is being performed as well (Sec. III.7).

In Sec. IV and Appendix D, we demonstrate that, for arbitrary reduced mass of
a two-body bound system, the two-loop vacuum-polarization corrections to the
energy can be evaluated analytically (for nonrelativistic reference states)
and expressed in terms of dilogarithmic, and trilogarithmic, functions. This
applies both to the $2P$–$2S$ difference (see Sec. IV.1) as well as to
individual hydrogenic levels (see Appendix D). Higher-order coefficients for
the semi-analytic expansion of the two-loop vacuum-polarization energy shift
could be evaluated with the help of the fully relativistic hydrogen wave
function. Results have been presented in Sec. IV.2.

## Acknowledgments

The authors acknowledge insightful conversations with Professor Ettore
Remiddi. Support from the National Science Foundation (grant PHY–2110294) is
gratefully acknowledged. S.L. also acknowledges support from Italian Ministry
of University and Research (MUR) via the PRIN 2022 project n. 20225X52RA —
MUS4GM2 funded by the European Union — Next Generation EU.

##  Appendix A Terms of Order $\maybebm{\epsilon^{2}}$

We first give the terms quadratic in $\epsilon$ for Eq. (II.2), for the master
integral $M_{12}$,

$\displaystyle N_{12}(q^{2})$
$\displaystyle=-12\text{Li}_{3}\left(\frac{1-v}{v+1}\right)-24\text{Li}_{3}\left(\frac{2v}{v+1}\right)+\ln^{3}\left(\frac{v-1}{v+1}\right)-6\ln\left(\frac{2v}{v+1}\right)\ln^{2}\left(\frac{v-1}{v+1}\right)$
$\displaystyle+12\ln^{2}\left(\frac{2v}{v+1}\right)\left(\ln\left(\frac{v-1}{v+1}\right)-\ln\left(-\frac{v-1}{v+1}\right)\right)-\pi^{2}\ln\left(\frac{v-1}{v+1}\right)+6\pi^{2}\ln\left(\frac{2v}{v+1}\right)+12\,\zeta(3)\,.$
(89)

We now give the expressions for the terms of order $\epsilon^{2}$ in the one-
loop effect, as discussed in Eq. (22). For the quadratic term $Q^{(1)}_{2}(v)$
in Eq. (22), one obtains the result

$\displaystyle Q^{(1)}_{2}(v)=$
$\displaystyle\;\left(v-\frac{4v^{3}}{9}\right)\biggl{[}-2\text{Li}_{2}\left(\frac{1-v}{v+1}\right)+\frac{1}{2}\ln^{2}\left(\frac{v-1}{v+1}\right)-2\ln\left(\frac{2v}{v+1}\right)\ln\left(\frac{v-1}{v+1}\right)-\frac{\pi^{2}}{6}\biggr{]}$
$\displaystyle+\frac{1}{36}v\left(v^{2}-3\right)\biggl{[}12\text{Li}_{3}\left(\frac{1-v}{v+1}\right)+24\text{Li}_{3}\left(\frac{2v}{v+1}\right)-\ln^{3}\left(\frac{v-1}{v+1}\right)+6\ln\left(\frac{2v}{v+1}\right)\ln^{2}\left(\frac{v-1}{v+1}\right)$
$\displaystyle+12\left(\ln\left(\frac{1-v}{v+1}\right)-\ln\left(\frac{v-1}{v+1}\right)\right)\ln^{2}\left(\frac{2v}{v+1}\right)+\pi^{2}\ln\left(\frac{v-1}{v+1}\right)-6\pi^{2}\ln\left(\frac{2v}{v+1}\right)-12\zeta(3)\biggr{]}$
$\displaystyle-\frac{2}{81}\left(\left(39v^{3}-81v\right)\ln\left(\frac{v-1}{v+1}\right)+78v^{2}-160\right)\,.$
(90)

For the term $Q_{1a}(r)$ from the integral with two fermion propagators, one
finds for the term quadratic in $\epsilon$ the expression

$\displaystyle Q^{(1a)}_{2}(v)=$
$\displaystyle\;\frac{1-v^{2}}{4v}\,\biggl{\\{}-\frac{1}{6}\left(v^{2}-1\right)\biggl{[}12\text{Li}_{3}\left(\frac{1-v}{v+1}\right)+24\text{Li}_{3}\left(\frac{2v}{v+1}\right)-\ln^{3}\left(\frac{v-1}{v+1}\right)+6\ln\left(\frac{2v}{v+1}\right)\biggl{(}\ln^{2}\left(\frac{v-1}{v+1}\right)-\pi^{2}\biggr{)}$
$\displaystyle-24\ln^{2}\left(\frac{2v}{v+1}\right)\ln\left(\frac{v-1}{v+1}\right)+12\left(\ln\left(\frac{1-v}{v+1}\right)+\ln\left(\frac{v-1}{v+1}\right)\right)\ln^{2}\left(\frac{2v}{v+1}\right)+\pi^{2}\ln\left(\frac{v-1}{v+1}\right)$
$\displaystyle-12\zeta(3)\biggr{]}+v^{2}\biggl{[}-4\text{Li}_{2}\left(\frac{1-v}{v+1}\right)+\ln^{2}\left(\frac{v-1}{v+1}\right)-4\left(\ln\left(\frac{2v}{v+1}\right)-1\right)\ln\left(\frac{v-1}{v+1}\right)-\frac{\pi^{2}}{3}\biggr{]}+8v\biggr{\\}}\,.$
(91)

##  Appendix B Coefficient Functions

The coefficients in Eq. (36) are as follows,

$\displaystyle{\mathcal{P}}^{(2;1)}_{234}=$
$\displaystyle\;-\frac{(2-D)\left(D^{2}\left(-q^{2}-4\right)-D\left(q^{4}-5q^{2}-28\right)+2\left(q^{4}-2q^{2}-20\right)\right)}{2(D-4)(D-3)(D-1)q^{2}\left(4-q^{2}\right)}\,,$
(92) $\displaystyle{\mathcal{P}}^{(2;1)}_{24}=$
$\displaystyle\;-\frac{(D-2)\left(3D^{2}q^{4}-4D\left(3q^{4}-2q^{2}+8\right)+4\left(3q^{4}-8q^{2}+20\right)\right)}{2(D-4)(D-3)(D-1)q^{4}\left(4-q^{2}\right)^{2}}\,,$
(93) $\displaystyle{\mathcal{P}}^{(2;1)}_{1234}=$
$\displaystyle\;\frac{-\left(D^{3}-9D^{2}+30D-32\right)q^{4}+4\left(D^{3}-11D^{2}+38D-40\right)q^{2}+32}{8(D-4)(D-1)q^{2}\left(4-q^{2}\right)}\,,$
(94) $\displaystyle{\mathcal{P}}^{(2;1)}_{235k}=$
$\displaystyle\;-\frac{1}{(D-4)(D-3)(D-1)q^{4}\left(4-q^{2}\right)^{2}}\biggl{[}3\left(-4D^{3}q^{2}\left(1-q^{2}\right)-D^{2}\left(-q^{6}+30q^{4}-36q^{2}+16\right)\right.$
$\displaystyle\left.+4D\left(-q^{6}+18q^{4}-28q^{2}+20\right)-4\left(-q^{6}+14q^{4}-28q^{2}+24\right)\right)\biggr{]}\,,$
(95) $\displaystyle{\mathcal{P}}^{(2;1)}_{235}=$
$\displaystyle\;-\frac{1}{4(D-4)(D-3)(D-1)q^{4}\left(4-q^{2}\right)^{2}}\biggl{[}-D^{4}q^{2}\left(4-q^{2}\right)^{2}+D^{3}q^{2}\left(23q^{4}-100q^{2}+176\right)$
$\displaystyle+4D^{2}\left(q^{8}-35q^{6}+119q^{4}-176q^{2}-32\right)-8D\left(2q^{8}-40q^{6}+127q^{4}-148q^{2}-88\right)$
$\displaystyle+8\left(2q^{8}-31q^{6}+98q^{4}-84q^{2}-120\right)\biggr{]}\,.$
(96)

The coefficients in Eq. (37) are as follows,

$\displaystyle{\mathcal{G}}^{(2;1)}_{235k}=$
$\displaystyle\;-\frac{3(D-2)\left(-D^{2}\left(-q^{2}-2\right)q^{2}-D\left(4q^{4}+6q^{2}+8\right)+4\left(q^{4}-q^{2}+6\right)\right)}{(D-4)(D-1)q^{4}\left(4-q^{2}\right)^{2}}\,,$
(97) $\displaystyle{\mathcal{G}}^{(2;1)}_{24}=$
$\displaystyle\;\frac{(2-D)\left(3D^{2}q^{4}-4D\left(3q^{4}-2q^{2}+8\right)+4\left(3q^{4}-8q^{2}+20\right)\right)}{4(D-4)(D-1)q^{4}\left(4-q^{2}\right)^{2}}\,,$
(98) $\displaystyle{\mathcal{G}}^{(2;1)}_{235}=$
$\displaystyle\;\frac{1}{4(D-4)(D-1)q^{4}\left(4-q^{2}\right)^{2}}\biggl{[}-3D^{3}\left(q^{2}+2\right)q^{4}+D^{2}\left(18q^{6}+26q^{4}+64\right)$
$\displaystyle+4D\left(-9q^{6}-4q^{4}+8q^{2}-88\right)+8\left(3q^{6}-2q^{4}-10q^{2}+60\right)\biggr{]}\,.$
(99)

The coefficients in Eq. (38) are as follows,

$\displaystyle{\mathcal{P}}^{(2;2)}_{234}=$
$\displaystyle\;\frac{(D-2)\left(-2D^{2}q^{2}-D\left(q^{4}-14q^{2}+8\right)+2\left(q^{4}-10q^{2}+12\right)\right)}{16(D-3)q^{2}\left(4-q^{2}\right)}\,,$
(100) $\displaystyle{\mathcal{P}}^{(2;2)}_{235k}=$
$\displaystyle\;\frac{3(D-2)}{8(D-4)(D-3)(D-1)q^{4}\left(4-q^{2}\right)^{2}}\biggl{[}-D^{3}q^{2}\left(4-q^{2}\right)^{2}+D^{2}q^{2}\left(7q^{4}-48q^{2}+128\right)$
$\displaystyle+2D\left(-7q^{6}+40q^{4}-136q^{2}-32\right)-8\left(-q^{6}+4q^{4}-12q^{2}-24\right)\biggr{]}\,,$
(101) $\displaystyle{\mathcal{P}}^{(2;2)}_{235}=$
$\displaystyle\;-\frac{1}{8(D-4)(D-3)(D-1)q^{4}\left(q^{2}-4\right)^{2}}\biggl{[}D^{4}q^{2}\left(q^{2}-4\right)^{2}\left(q^{2}+1\right)-3D^{3}q^{2}\left(3q^{6}-19q^{4}\right.$
$\displaystyle\left.+28q^{2}+48\right)+4D^{2}\left(7q^{8}-40q^{6}+69q^{4}+112q^{2}+32\right)-4D\left(9q^{8}-45q^{6}+80q^{4}\right.$
$\displaystyle\left.+128q^{2}+176\right)+16\left(q^{8}-4q^{6}+6q^{4}+6q^{2}+60\right)\biggr{]}\,.$
(102) $\displaystyle{\mathcal{P}}^{(2;2)}_{24}=$
$\displaystyle\;\frac{D-2}{4(D-4)\left(D^{2}-4D+3\right)q^{4}\left(q^{2}-4\right)^{2}}\biggl{[}D^{4}q^{2}\left(q^{2}-4\right)-8D^{3}q^{2}\left(q^{2}-4\right)$
$\displaystyle+12D^{2}q^{2}\left(2q^{2}-7\right)+D\left(-34q^{4}+96q^{2}-32\right)+20q^{4}-64q^{2}+80\biggr{]}\,.$
(103)

##  Appendix C Higher–Order Terms in the Master Integrals

We have the following results for the higher-order terms from the master
integrals, with reference to Eqs. (III.4) and (III.4). For $N_{235}(q^{2})$,
we obtain

$\displaystyle N_{235}(q^{2})$
$\displaystyle=\frac{1}{6}v\left(v^{2}-3\right)\biggl{[}18\Phi_{2}\left(\frac{v-1}{v+1}\right)-18\ln\left(\frac{v-1}{v+1}\right)\Phi_{1}\left(\frac{v-1}{v+1}\right)-6\text{Li}_{2}\left(\frac{v-1}{v+1}\right)-6\ln\left(\frac{2}{v+1}\right)\ln\left(\frac{v-1}{v+1}\right)$
(104)
$\displaystyle-39\ln\left(\frac{v-1}{v+1}\right)+\pi^{2}\biggr{]}+\frac{1}{12}\left(v^{4}+2v^{2}-3\right)\biggl{[}36\Phi_{3}\left(\frac{v-1}{v+1}\right)-18\ln\left(\frac{v-1}{v+1}\right)\Phi_{2}\left(\frac{v-1}{v+1}\right)$
$\displaystyle-6\text{Li}_{2}\left(\frac{v-1}{v+1}\right)\ln\left(\frac{v-1}{v+1}\right)-3\ln^{3}\left(\frac{v-1}{v+1}\right)-6\ln\left(\frac{2}{v+1}\right)\ln^{2}\left(\frac{v-1}{v+1}\right)+\pi^{2}\ln\left(\frac{v-1}{v+1}\right)+18\zeta(3)\biggr{]}$
$\displaystyle+\frac{1}{8}\left(60v^{2}+\left(-5v^{4}-16v^{3}-14v^{2}+48v+3\right)\ln^{2}\left(\frac{v-1}{v+1}\right)+170\right)\,.$

For $N_{235k}(q^{2})$, we obtain

$\displaystyle N_{235k}(v)$
$\displaystyle=\frac{1}{36}v\left(3v^{4}-8v^{2}+21\right)\biggl{[}18\Phi_{2}\left(\frac{v-1}{v+1}\right)-18\ln\left(\frac{v-1}{v+1}\right)\Phi_{1}\left(\frac{v-1}{v+1}\right)-6\text{Li}_{2}\left(\frac{v-1}{v+1}\right)$
(105)
$\displaystyle-6\ln\left(\frac{2}{v+1}\right)\ln\left(\frac{v-1}{v+1}\right)+\pi^{2}\biggr{]}+\frac{1}{24}\left(v^{6}-3v^{4}-5v^{2}+7\right)\biggl{[}36\Phi_{3}\left(\frac{v-1}{v+1}\right)-18\ln\left(\frac{v-1}{v+1}\right)\Phi_{2}\left(\frac{v-1}{v+1}\right)$
$\displaystyle-6\text{Li}_{2}\left(\frac{v-1}{v+1}\right)\ln\left(\frac{v-1}{v+1}\right)-3\ln^{3}\left(\frac{v-1}{v+1}\right)-6\ln\left(\frac{2}{v+1}\right)\ln^{2}\left(\frac{v-1}{v+1}\right)+\pi^{2}\ln\left(\frac{v-1}{v+1}\right)+18\zeta(3)\biggr{]}$
$\displaystyle+\frac{1}{48}\biggl{[}-2\left(86v^{4}+233v^{2}+601\right)-4v\left(37v^{4}-108v^{2}+279\right)\ln\left(\frac{v-1}{v+1}\right)$
$\displaystyle+\biggl{(}-13v^{6}-48v^{5}+63v^{4}+128v^{3}+105v^{2}-336v-27\biggr{)}\ln^{2}\left(\frac{v-1}{v+1}\right)\biggr{]}\,.$
(106)

##  Appendix D Analytic Integrations

###  D.1 One–Loop Diagram

We recall the imaginary part of the one-loop effect from Eq. (28),

$\mathrm{Im}\\!\left[\Pi^{(1)}_{\rm
R}(q^{2}+\mathrm{i}\epsilon)\right]=\frac{\alpha}{3}\sqrt{1-\frac{4}{q^{2}}}\,\left(1+\frac{2}{q^{2}}\right)\,,$
(107)

for $q^{2}\geq 4$. We use Eq. (10.245) of Ref. JeAd2022book  and write

$V^{(1)}_{R}(r)=-\frac{Z\alpha}{\pi}\int\limits_{4}^{\infty}\frac{\mathrm{d}(q^{2})}{q^{2}}\;\frac{\mathrm{e}^{-qr}}{r}\;\mathrm{Im}\\!\left[\Pi^{(1)}_{\rm
R}(q^{2}+\mathrm{i}\epsilon)\right]\,.$ (108)

One uses the well-known formula for the generalized Bohr radius in a
hydrogenlike ion with nuclear charge number $Z$, namely, $\beta\equiv
a_{0}=1/(Z\alpha\mu)$, where $\mu$ is the reduced mass of the two-body system,
measured in units of the electron mass. For the $2S$ and $2P$ states of atomic
hydrogen, one finds, respectively,

$\displaystyle\langle 2S|\frac{\exp(-q\,r)}{r}|2S\rangle=$
$\displaystyle\;\frac{1+2(\beta q)^{2}}{4\beta\,(1+\beta q)^{4}}\,,$ (109a)
$\displaystyle\langle 2P|\frac{\exp(-q\,r)}{r}|2P\rangle=$
$\displaystyle\;\frac{1}{4\beta\,(1+\beta q)^{4}}\,.$ (109b)

Hence,

$\langle 2P|\frac{\exp(-q\,r)}{r}|2P\rangle-\langle
2S|\frac{\exp(-q\,r)}{r}|2S\rangle\\\ =-\frac{(\beta
q)^{2}}{2\,\beta\,(1+\beta q)^{4}}\,,$ (110)

The contribution to the $2P$–$2S$ energy shift is

$E_{1}(2P-2S)=\langle 2P|V^{(1)}_{R}|2P\rangle-\langle
2S|V^{(1)}_{R}|2S\rangle\\\
=\frac{(Z\alpha)^{2}\mu}{2\pi}\int\limits_{4m^{4}}^{\infty}\frac{\mathrm{d}(q^{2})}{q^{2}}\;\mathrm{Im}\\!\left[\Pi_{\rm
R}^{(1)}(q^{2}+\mathrm{i}\epsilon)\right]\,\frac{(\beta q)^{2}}{(1+\beta
q)^{4}}\\\ =\left(\dfrac{\alpha}{\pi}\right)\,(Z\alpha)^{2}\mu
f_{1}(\beta)\,,\\\
f_{1}(\beta)=\int_{4}^{{\infty}}\frac{\mathrm{d}(q^{2})}{q^{2}}\frac{{\rm
Im}\Pi_{\rm R}^{(1)}(q^{2}+\mathrm{i}\epsilon)}{\alpha}\frac{(\beta
q)^{2}}{2(1+\beta q)^{4}}\,.$ (111)

We choose the $2P$–$2S$ interval in order to demonstrate that energy shifts
due to the one-loop and two-loop vacuum-polarization effects can be evaluated
analytically. For muonic hydrogen, one has the result
$\mu=m_{\mu}m_{p}/(m_{\mu}+m_{p})$, where $m_{\mu}$ and $m_{p}$ are the muon
and proton masses, respectively.

It turns out to be advantageous to express the result in terms of the variable
$w$ defined in Eq. (70), with the result

$\displaystyle f_{1}(w)=$
$\displaystyle\frac{\left(w^{2}-1\right)^{2}}{72w^{5}}\biggl{[}3\left(8w^{4}+4w^{2}+3\right)\ln\left(\frac{1-w}{1+w}\right)$
$\displaystyle+\frac{24\pi
w^{5}}{\sqrt{1-w^{2}}}+\frac{2\left(24w^{6}-28w^{4}-3w^{2}+9\right)w}{\left(w^{2}-1\right)^{2}}\biggr{]}\,.$
(112)

We can get the expansions for small and large $\beta$,

$\displaystyle f_{1}(\beta)=$
$\displaystyle\;\frac{1}{18}-\beta^{2}+\frac{8\pi\beta^{3}}{3}+\beta^{4}\left(10\ln\left(\beta^{2}\right)+\frac{32}{3}\right)$
$\displaystyle\;+\beta^{6}\left(\frac{140\ln\left(\beta^{2}\right)}{3}+\frac{236}{3}\right)+{\mathcal{O}}\left(\beta^{8}\right)\,,$
(113a) $\displaystyle f_{1}(\beta)=$
$\displaystyle\;\frac{1}{30\beta^{2}}-\frac{5\pi}{384\beta^{3}}+\frac{1}{28\beta^{4}}-\frac{35\pi}{4096\beta^{5}}+\frac{1}{54\beta^{6}}$
$\displaystyle-\frac{63\pi}{16384\beta^{7}}+{\mathcal{O}}\left(\frac{1}{\beta^{8}}\right)\,.$
(113b)

For muonic hydrogen, one has $\beta=0.737384\dots$, while, for ordinary
hydrogen, one has $\beta=137.110\dots$. Hence, the second of the above
equations is relevant for ordinary hydrogen, while none of the expansions can
be used with good accuracy for muonic hydrogen. It should be possible, though,
to generalize the results reported in this Appendix to reference states other
than $2S$ and $2P$ if desired. Relativistic corrections to the one-loop
vacuum-polarization shift have been analyzed in Ref. Je2011pra , with an
emphasis on muonic hydrogen.

###  D.2 Two–Loop Reducible Diagram

Next, we consider the contribution of the reducible diagram 3 in Fig. 1,

$E^{{\rm(r)}}(2P-2S)=\left(\dfrac{\alpha}{\pi}\right)^{2}(Z\alpha)^{2}\mu\,f^{{\rm(r)}}(\beta)\,,\\\
f^{{\rm(r)}}(\beta)=-\int_{4}^{\infty}\frac{\mathrm{d}(q^{2})}{q^{2}}\frac{{\rm
Im}\left[\left(\Pi^{(1)}_{\rm
R}(q^{2}+\mathrm{i}\epsilon)\right)^{2}\right]}{\alpha^{2}/\pi}\,\frac{(\beta
q)^{2}}{2(1+\beta q)^{4}}\,.$ (114)

We use the conventions outlined Eq. (70). An analytic evaluation leads to the
formula

$f^{{\rm(r)}}(w)=-\frac{1}{648w^{5}}\biggl{\\{}4w\biggl{(}-315w^{8}+1005w^{6}-836w^{4}\\\
+84w^{2}+72\biggr{)}+12\biggl{(}-105w^{10}+370w^{8}-381w^{6}\\\
+96w^{4}+10w^{2}+12\biggr{)}\ln\left(\frac{1-w}{1+w}\right)\\\
-45\left(w^{2}-1\right)^{2}\left(7w^{2}-13\right)w^{5}\biggl{[}\ln^{2}\left(\frac{1-w}{1+w}\right)+\pi^{2}\biggr{]}\biggr{\\}}\,.$
(115)

The expansion for small $\beta$ reads as follows,

$f^{{\rm(r)}}(\beta)=-\frac{\ln\left(\beta^{2}\right)}{27}-\frac{5}{81}+\beta^{2}\left(\frac{2}{3}\ln\left(\beta^{2}\right)+\frac{32}{9}\right)\\\
-\frac{\beta^{4}}{27}\left(180\ln^{2}\left(\beta^{2}\right)+774\ln\left(\beta^{2}\right)+180\pi^{2}-73\right)+{\mathcal{O}}(\beta^{6})\,.$
(116)

For large $\beta$, one obtains the result

$\displaystyle f^{{\rm(r)}}(\beta)=$
$\displaystyle\;\frac{23\pi}{4536\beta^{3}}-\frac{1}{45\beta^{4}}+\frac{127\pi}{19008\beta^{5}}-\frac{1}{60\beta^{6}}$
$\displaystyle\;+\frac{1675\pi}{439296\beta^{7}}+{\mathcal{O}}\left(\frac{1}{\beta^{8}}\right)\,.$
(117)

###  D.3 Two–Loop Irreducible Diagram

We continue with the contribution of the irreducible two-loop vacuum
polarization diagrams,

$E_{2}(2P-2S)=\left(\dfrac{\alpha}{\pi}\right)^{2}(Z\alpha)^{2}\mu
f_{2}(\beta)\,,\\\
f_{2}(\beta)=\int_{4}^{{\infty}}\frac{\mathrm{d}(q^{2})}{q^{2}}\frac{{\rm
Im}\,\Pi^{(2)}_{\rm R}(q^{2}+\mathrm{i}\epsilon)}{\alpha^{2}/\pi}\frac{(\beta
q)^{2}}{2(1+\beta q)^{4}}\,.$ (118)

The result is a bit more complex and given in Eq. (IV.1). The expansion for
small $\beta$ is

$\displaystyle f_{2}$
$\displaystyle(\beta)=\frac{1}{24}-\beta^{2}\left(3\ln(\beta)+\frac{11}{2}\right)+\beta^{3}\left(\frac{320\pi}{27}-\frac{56\pi^{2}}{9}\right)$
$\displaystyle+\beta^{4}\left(60\ln^{2}(\beta)+69\ln(\beta)-20\zeta(3)+15\pi^{2}+\frac{43}{12}\right)$
$\displaystyle+{\mathcal{O}}\left(\beta^{5}\right)\,.$ (119)

For large $\beta$, one finds the result

$\displaystyle f_{2}$
$\displaystyle(\beta)=\frac{41}{324\beta^{2}}+\frac{1}{\beta^{3}}\left(-\frac{15647\pi}{105840}+\frac{25\pi^{2}}{504}-\frac{13}{126}\pi\ln(2)\right)$
$\displaystyle+\frac{449}{2160\beta^{4}}+\frac{1}{\beta^{5}}\left(-\frac{75763\pi}{544320}+\frac{31\pi^{2}}{864}-\frac{19}{432}\pi\ln(2)\right)$
$\displaystyle+\frac{62479}{453600\beta^{6}}+\frac{1}{\beta^{7}}\biggl{[}-\frac{71317517\pi}{936714240}+\frac{37\pi^{2}}{2112}$
$\displaystyle-\frac{125}{8448}\pi\ln(2)\biggr{]}+{\mathcal{O}}\left(\frac{1}{\beta^{8}}\right)\,.$
(120)

###  D.4 Individual Energy Shifts

In Sec. IV.1, we had concentrated on the $2P$–$2S$ energy difference, for
muonic hydrogen. It is instructive to separately give the nonrelativistic
contributions to selected energy levels with principal quantum number $n$ and
orbital angular momentum quantum number $\ell$ (in the nonrelativistic
approximation) from the one-loop diagram, the two-loop irreducible and the
two-loop reducible diagrams,

$\displaystyle E^{(1)}(n\ell)=$
$\displaystyle\;\left(\dfrac{\alpha}{\pi}\right)\,(Z\alpha)^{2}\mu
f^{(1)}_{n\ell}(\beta)\,,$ (121a) $\displaystyle E^{{\rm(r)}}(n\ell)=$
$\displaystyle\;\left(\dfrac{\alpha}{\pi}\right)^{2}\,(Z\alpha)^{2}\mu
f^{{\rm(r)}}_{n\ell}(\beta)\,,$ (121b) $\displaystyle E^{(2)}(n\ell)=$
$\displaystyle\;\left(\dfrac{\alpha}{\pi}\right)^{2}\,(Z\alpha)^{2}\mu
f^{(2)}_{n\ell}(\beta)\,.$ (121c)

We consider the cases $|n\ell\rangle=|1S\rangle,|2S\rangle,|2P\rangle$. The
variable $w$ has been defined in Eq. (70). We define the variable $u$ as
follows,

$u=\sqrt{1-\beta^{2}}\,,\quad\beta=\sqrt{1-u^{2}}\,.$ (122)

For the ground state, we obtain

$f^{(1)}_{1S}(\beta)=\frac{1}{18u}\biggl{[}3\pi\sqrt{1-u^{2}}\left(4u^{2}-7\right)u-24u^{3}+46u\\\
-3\left(4u^{4}-9u^{2}+3\right)\ln\left(\frac{1-u}{1+u}\right)\biggr{]}\,.$
(123)

For the $2S$ state, one has the result

$f^{(1)}_{2S}(\beta)=\frac{1}{144w^{5}}\biggl{[}-336w^{7}+464w^{5}-18w^{3}-54w\\\
+12\pi w^{5}\left(14w^{2}-17\right)\sqrt{1-w^{2}}\\\
-3\left(56w^{8}-96w^{6}+27w^{4}+9\right)\ln\left(\frac{1-w}{1+w}\right)\biggr{]}\,.$
(124)

For the $2P$ level, we obtain

$f^{(1)}_{2P}(\beta)=\frac{1}{144w^{5}}\biggl{[}-240w^{7}+352w^{5}-30w^{3}-18w\\\
+12\pi w^{5}\left(10w^{2}-13\right)\sqrt{1-w^{2}}\\\
-3\left(40w^{8}-72w^{6}+21w^{4}+4w^{2}+3\right)\ln\left(\frac{1-w}{1+w}\right)\biggr{]}\,.$
(125)

For the reducible diagram, the results are as follows. Let us first indicate
the result for the ground state,

$f^{{\rm(r)}}_{1S}(\beta)=-\frac{1}{1620u}\biggl{[}-288\pi u\sqrt{1-u^{2}}\\\
+20u\left(63u^{4}-282u^{2}+304\right)\\\
+45u\left(7u^{6}-36u^{4}+51u^{2}-18\right)\left[\ln^{2}\left(\frac{1-u}{1+u}\right)+\pi^{2}\right]\\\
+60\left(21u^{6}-101u^{4}+126u^{2}-24\right)\ln\left(\frac{1-u}{1+u}\right)\biggr{]}\,.$
(126)

For the $2S$ state, one has the result

$f^{{\rm(r)}}_{2S}(\beta)=-\frac{1}{6480w^{5}}\biggl{[}-576\pi
w^{5}\sqrt{1-w^{2}}\\\
+20\left(1386w^{8}-4719w^{6}+4132w^{4}-396w^{2}-216\right)w\\\
+45w^{5}\left(154w^{6}-627w^{4}+792w^{2}-315\right)\left[\ln^{2}\left(\frac{1-w}{1+w}\right)+\pi^{2}\right]\\\
+60\left(462w^{10}-1727w^{8}+1851w^{6}-468w^{4}-54w^{2}-36\right)\times\\\
\ln\left(\frac{1-w}{1+w}\right)\biggr{]}\,,$ (127)

while the $2P$ result is

$f^{{\rm(r)}}_{2P}(\beta)=-\frac{1}{6480w^{5}}\biggl{[}-576\pi
w^{5}\sqrt{1-w^{2}}\\\
+60\left(252w^{8}-903w^{6}+820w^{4}-76w^{2}-24\right)w\\\
+45w^{5}\left(84w^{6}-357w^{4}+462w^{2}-185\right)\times\\\
\left[\ln^{2}\left(\frac{1-w}{1+w}\right)+\pi^{2}\right]+60\biggl{(}252w^{10}-987w^{8}\\\
+1089w^{6}-276w^{4}-34w^{2}-12\biggr{)}\ln\left(\frac{1-w}{1+w}\right)\biggr{]}\,.$
(128)

The difference of the results given in Eqs. (128) and (127) confirms Eq.
(115).

The results for the irreducible two-loop diagram are more complicated and read
as follows [we refer to the definitions of $I_{1}$ and $I_{2}$ in Eqs. (IV.1)
and (IV.1)]. For the ground state, one obtains

$f^{(2)}_{1S}(\beta)=\frac{1}{96u^{2}}\biggl{[}128\left(11-5u^{2}\right)u^{3}I_{2}\left(\sqrt{\frac{1-u}{1+u}}\right)\\\
+\left(43u^{6}-94u^{4}-45u^{2}-48\right)\left(\ln^{2}\left(\frac{1-u}{1+u}\right)+\pi^{2}\right)\\\
+4\left(43-37u^{2}\right)u^{3}\ln\left(\frac{1-u}{1+u}\right)\\\
+4u^{2}\left(8\left(5u^{4}-10u^{2}+1\right)I_{1}\left(\sqrt{\frac{1-u}{1+u}}\right)-57u^{2}+74\right)\\\
+\sqrt{1-u^{2}}\biggl{[}\pi\left(\frac{32}{9}\left(45u^{2}-113\right)u^{2}+\frac{1024}{3}u^{2}\ln(2)\right)\\\
+\pi^{2}\left(-80u^{4}+\frac{320u^{2}}{3}+48\right)\biggr{]}\biggr{]}\,.$
(129)

For the metastable $2S$ level, the analytic integration leads to the result

$\biggl{[}f^{(2)}_{2S}(\beta)=\frac{1}{1152w^{6}}\biggl{[}3\biggl{(}\biggl{(}653w^{10}-1058w^{8}+93w^{6}\\\
+56w^{4}+256w^{2}-192\biggr{)}\left(\ln^{2}\left(\frac{1-w}{1+w}\right)+\pi^{2}\right)\\\
-4w^{2}\left(447w^{6}-430w^{4}-88w^{2}+48\right)\biggr{)}\\\
-384\left(55w^{6}-97w^{4}+28w^{2}+6\right)w^{3}I_{2}\left(\sqrt{\frac{1-w}{1+w}}\right)\\\
+96\left(55w^{4}-110w^{2}+51\right)w^{6}I_{1}\left(\sqrt{\frac{1-w}{1+w}}\right)\\\
+12\left(-227w^{8}+77w^{6}+224w^{4}+52w^{2}-120\right)w\times\\\
\ln\left(\frac{1-w}{1+w}\right)+\sqrt{1-w^{2}}\biggl{(}\pi\biggl{(}2048w^{6}\ln(2)\\\
+\frac{32}{3}\left(495w^{4}-757w^{2}+54\right)w^{4}\biggr{)}-16\pi^{2}\\\
\times\left(165w^{8}-220w^{6}+21w^{4}+30w^{2}-36\right)\biggr{)}\biggr{]}\,.$
(130)

Finally, the irreducible two-loop diagram, for the $2P$ level, leads to the
energy shift

$f^{(2)}_{2P}(\beta)=\frac{1}{1152w^{6}}\biggl{[}\biggl{(}1163w^{10}-1966w^{8}+75w^{6}\\\
+184w^{4}+64w^{2}-192\biggr{)}\left(\ln^{2}\left(\frac{1-w}{1+w}\right)+\pi^{2}\right)\\\
-128\left(105w^{6}-191w^{4}+52w^{2}+6\right)w^{3}I_{2}\left(\sqrt{\frac{1-w}{1+w}}\right)\\\
+4w^{2}\biggl{(}24\left(35w^{4}-70w^{2}+31\right)w^{4}I_{1}\left(\sqrt{\frac{1-w}{1+w}}\right)\\\
-937w^{6}+978w^{4}+88w^{2}-48\biggr{)}\\\
-4\left(517w^{8}-299w^{6}-400w^{4}+44w^{2}+120\right)w\times\\\
\ln\left(\frac{1-w}{1+w}\right)+\sqrt{1-w^{2}}\biggl{(}\pi\biggl{(}2048w^{6}\ln(2)\\\
+\frac{32}{3}\left(315w^{4}-529w^{2}+18\right)w^{4}\biggr{)}\\\
-16\pi^{2}\left(105w^{8}-140w^{6}+9w^{4}-2w^{2}-12\right)\biggr{)}\biggr{]}\,.$
(131)

The difference of the $2P$ and $2S$ energy shifts confirms Eq. (74) for muonic
hydrogen.

## References

  * (1) K. Pachucki, Theory of the Lamb shift in muonic hydrogen, Phys. Rev. A 53, 2092–2100 (1996).
  * (2) U. D. Jentschura, Lamb Shift in Muonic Hydrogen. —I. Verification and Update of Theoretical Predictions, Ann. Phys. (N.Y.) 326, 500–515 (2011).
  * (3) U. D. Jentschura, Proton Radius: A Puzzle or a Solution!?, J. Phys. Conf. Ser. 2391, 012017 (2022).
  * (4) G. Pustovalov, Zh. Éksp. Teor. Fiz. 32, 1519 (1957), [JETP 5, 1234 (1957)].
  * (5) U. D. Jentschura, G. Soff, V. G. Ivanov, and S. G. Karshenboim, Bound $\mu^{+}\mu^{-}$ system, Phys. Rev. A 56, 4483–4495 (1997).
  * (6) S. G. Karshenboim, V. G. Ivanov, U. D. Jentschura, and G. Soff, Bound states of the muon-antimuon system: lifetimes and hyperfine splitting, JETP 86, 226–236 (1998), [ZhETF JETP 113, 409 (1998)].
  * (7) S. G. Karshenboim, U. D. Jentschura, V. G. Ivanov, and G. Soff, Next-to-leading and higher-order corrections to the decay rate of dimuonium, Phys. Lett. B 424, 397–404 (1998).
  * (8) G. ’t Hooft and M. Veltman, Regularization and Renormalization of Gauge Fields, Nucl. Phys. B 44, 189–213 (1972).
  * (9) C. G. Bollini and J. J. Giambiagi, Dimensional Renormalization: The Number of Dimensions is a Regularizing Parameter, Nuovo Cimento 12B, 20–26 (1972).
  * (10) R. Barbieri and E. Remiddi, Infra-red divergences and adiabatic switching. Fourth-order vacuum polarization, Nuovo Cim. A 13, 99–119 (1973).
  * (11) J. D. Bjorken and S. D. Drell, Relativistic Quantum Fields (McGraw-Hill, New York, 1965).
  * (12) M. E. Peskin and D. V. Schroeder, An Introduction to Quantum Field Theory (Perseus, Cambridge, Massachusetts, 1995).
  * (13) C. Itzykson and J.-B. Zuber, Quantum Field Theory (McGraw-Hill, New York, 1980).
  * (14) V. B. Berestetskii, E. M. Lifshitz, and L. P. Pitaevskii, Quantum Electrodynamics, Volume 4 of the Course on Theoretical Physics, 2 ed. (Pergamon Press, Oxford, UK, 1982).
  * (15) U. D. Jentschura and G. S. Adkins, Quantum Electrodynamics: Atoms, Lasers and Gravity (World Scientific, Singapore, 2022).
  * (16) J. A. M. Vermaseren, arXiv:math-ph/0010025.
  * (17) J. Kuipers, T. Ueda, J. A. M. Vermaseren, and J. Vollinga, FORM version 4.0, Comput. Phys. Commun. 184, 1453 (2013).
  * (18) K. G. Chetyrkin and F. V. Tkachov, Integration by Parts: The Algorithm to Calculate beta Functions in 4 Loops, Nucl. Phys. B 192, 159 (1981).
  * (19) F. V. Tkachov, A Theorem on Analytical Calculability of Four Loop Renormalization Group Functions, Phys. Lett. B 100, 65 (1981).
  * (20) S. Laporta, High precision calculation of multiloop Feynman integrals by difference equations, Int. J. Mod. Phys. A 15, 5087 (2000).
  * (21) S. Laporta, Calculation of master integrals by difference equations, Phys. Lett. B 504, 188 (2001).
  * (22) A. V. Kotikov, Differential equations method: New technique for massive Feynman diagrams calculation, Phys. Lett. B 254, 158 (1991).
  * (23) E. Remiddi, Differential equations for Feynman graph amplitudes, Nuovo Cim. A 254, 158 (1991).
  * (24) T. Gehrmann and E. Remiddi, Differential equations for two loop four point functions, Nucl. Phys. B 580, 458 (2000).
  * (25) R. Scharf and J. B. Tausk, Scalar two loop integrals for gauge boson selfenergy diagrams with a massless fermion loop, Nucl. Phys. B 412, 523–552 (1994).
  * (26) C. Anastasiou, S. Beerli, S. Bucherer, A. Daleo, and Z. Kunszt, Two-loop amplitudes and master integrals for the production of a Higgs boson via a massive quark and a scalar-quark loop, J. High Energy Phys. 0701, 082 (2007).
  * (27) S. Laporta and U. D. Jentschura, Dimensional Regularization and Two–Loop Vacuum Polarization Operator: Master Integrals, Analytic Results and Energy Shifts (Supplementary Material), submitted (2024).
  * (28) A. I. Davydychev, Explicit results for all orders of the epsilon expansion of certain massive and massless diagrams, Phys. Rev. D 61, 087701 (2000).
  * (29) A. I. Davydychev and M. Y. Kalmykov, New results for the epsilon expansion of certain one, two and three loop Feynman diagrams, Nucl. Phys. B 605, 266–318 (2001).
  * (30) A. I. Davydychev and M. Y. Kalmykov, Massive Feynman diagrams and inverse binomial sums, Nucl. Phys. B 699, 3–64 (2004).
  * (31) G. Källén and A. Sabry, Two–Loop Vacuum Polarization, Dan. Mat. Fys. Medd. 29, no. 17 (1955).
  * (32) J. Schwinger, Particles, Sources and Fields (Volume III) (Addison-Wesley, Reading, MA, 1970).
  * (33) D. J. Broadhurst, J. Fleischer, and O. V. Tarasov, Two-loop two-point functions with masses: asymptotic expansions and Taylor series, in any dimension, Z. Phys. C 60, 287–301 (1993).
  * (34) B. A. Kniehl, Two-loop corrections to the vacuum polarizations in perturbative QCD, Nucl. Phys. B 347, 86–104 (1990).
  * (35) A. Djouadi and P. Gambino, Electroweak gauge boson self-energies: Complete QCD corrections, Phys. Rev. D 49, 3499–3511 (1994).
  * (36) K. G. Chetyrkin, J. H. Kühn, and M. Steinhauser, Three-loop polarization function and $O(\alpha_{s}^{2})$ corrections to the production of heavy quarks, Nucl. Phys. B 482, 213–240 (1996).
  * (37) E. Borie, Lamb Shift in Light Muonic Atoms, Z. Phys. A 275, 347–349 (1975).
  * (38) E. Borie and G. A. Rinker, The energy levels of muonic atoms, Rev. Mod. Phys. 54, 67–118 (1982).
  * (39) K. Pachucki, Contributions to the binding, two-loop correction to the Lamb shift, Phys. Rev. A 48, 2609–2614 (1993).
  * (40) P. J. Mohr, B. N. Taylor, and D. B. Newell, CODATA recommended values of the fundamental physical constants: 2006, Rev. Mod. Phys. 80, 633–730 (2008).
  * (41) E. Tiesinga, P. J. Mohr, D. B. Newell, and B. N. Taylor, CODATA recommended values of the fundamental physical constants: 2018, Rev. Mod. Phys. 93, 025010 (2021).
  * (42) U. D. Jentschura, A. Czarnecki, and K. Pachucki, Nonrelativistic QED approach to the Lamb shift, Phys. Rev. A 72, 062102 (2005).
  * (43) U. D. Jentschura, Relativistic Reduced–Mass and Recoil Corrections to Vacuum Polarization in Muonic Hydrogen, Muonic Deuterium and Muonic Helium Ions, Phys. Rev. A 84, 012505 (2011).