$\mathcal{O}\left(m\alpha^{2}(Z\alpha)^{6}\right)$ contribution to Lamb shift from radiative corrections to the Wichmann-Kroll potential.

Recent advances in spectroscopy of ordinary hydrogen Bezginov2019 ; Grinin2020 and deuterium PhysRevLett.104.233001 and in their muonic analogs antognini2013proton ; pohl2016laser as well as in the electron-proton scattering xiong2019small provide new opportunities to perform subtle tests of the bound-state quantum electrodynamics (QED).

The hydrogen atom is a system in which the role of the relativistic, radiative, and recoil effects can be investigated with high precision, both experimentally and theoretically. The corresponding corrections to energy levels are small compared to the leading term; however, with a present accuracy achieved in measuring the frequency of specific hydrogenic transitions, those contributions are large compared with the experimental uncertainty.

One of the most important bound-state QED effects is the Lamb shift. It is responsible for the $2s_{1/2}$ – $2p_{1/2}$ energy splitting, which in the Dirac equation approximation would be zero. The calculations of various contributions to the Lamb shift have a long history starting from Refs. Bethe1947 ; Karplus1951 ; Kroll1951 ; Karplus1952 , see also review Eides2007 and references therein. For the $s$-states all corrections have been calculated up to the order $m\alpha^{2}(Z\alpha)^{6}\ln\left(1/(Z\alpha)^{2}\right)$. The $m\alpha^{2}(Z\alpha)^{6}$ contribution has not yet been calculated, although this correction may be important already in the next series of spectroscopic measurements.

In the present paper we calculate one of the previously unknown corrections to the Lamb shift of order $m\alpha^{2}(Z\alpha)^{6}$, which is connected with the radiative corrections to the Wichmann-Kroll (WK) potential. We use modern multiloop methods and obtain analytic result in terms of conventional polylogarithmic constants. Recently, we applied similar approach to the calculation of certain two-loop corrections to Lamb shift and hyperfine splitting in hydrogen, Ref. Krachkov:2023tly . The present calculation provides yet another example of the effectiveness of multiloop methods for obtaining analytic results in atomic physics.

Feynman diagrams for the radiative correction $\delta\tilde{\rho}(q)$ to the Wichmann-Kroll charge density are depicted in Fig. 1. On the same figure we also show the diagrams with one-loop conterterms. The black squares on solid lines correspond to the mass counter terms $i\delta m=im\frac{(4\pi\alpha)m^{-2\epsilon}(3-2\epsilon)\Gamma(\epsilon)}{(4\pi)^{2-\epsilon}(1-2\epsilon)}$, while the black squares on dashed lines correspond to $-i\delta Z_{A}(q^{2}g_{\mu\nu}-q_{\mu}q_{\nu})=i\frac{(16\pi\alpha)m^{-2\epsilon}\Gamma(\epsilon)}{3(4\pi)^{2-\epsilon}}(q^{2}g_{\mu\nu}-q_{\mu}q_{\nu})$. In principle, we should also account for the diagrams with one-loop fermion field renormalization and vertex counterterms, but their contributions cancel due to the Ward identity.

The corresponding correction to the potential $\delta\widetilde{V}(q)=\delta\tilde{\rho}(q)\cdot(e/q^{2})$ contributes to energy shifts. The characteristic atomic momenta $q$ are small compared to the electron mass. Therefore, we need the small-$q$ asymptotics of $\tilde{\rho}(q)$. This asymptotics can be analyzed using the expansion by regions approach. There are two regions which are relevant to this asymptotics. The hard region corresponds to all loop momenta $\sim m$. In this region we expand the integrand of $\delta\tilde{\rho}(q)$ in Taylor series in $q$. The zeroth term corresponds to the sum of diagrams in Fig. 1 at $q=0$. It is easy to see that, due to the identity

$$\sum_{k=1}^{N}\operatorname{Tr}\Big{[}\gamma^{\mu_{1}}(\hat{p}-\hat{l}_{1}-m)^{-1}\ldots\gamma^{\mu_{k}}(\hat{p}-\hat{l}_{k}-m)^{-1}\gamma^{\alpha}(\hat{p}-\hat{l}_{k}-m)^{-1}\ldots\gamma^{\mu_{N}}(\hat{p}-\hat{l}_{N}-m)^{-1}\Big{]}\\ =-\frac{\partial}{\partial p_{\alpha}}\operatorname{Tr}\left[\gamma^{\mu_{1}}(\hat{p}-\hat{l}_{1}-m)^{-1}\ldots\gamma^{\mu_{N}}(\hat{p}-\hat{l}_{N}-m)^{-1}\right]\,,$$

(1)

this sum can be written as the integral of total derivative which is zero in dimensional regularization. The sum in the left-hand side of the above identity corresponds to all possible insertions of the vertex $\gamma^{\alpha}$ in the fermion loop.

The contribution to $\delta\tilde{\rho}(q)$ linear in $\bm{q}$ is zero due to rotation symmetry. Therefore, the expansion in the hard region starts from $q^{2}$ term,

$$\delta\tilde{\rho}(q)\approx(\delta\widetilde{V}(0)/e)\cdot q^{2}.$$

(2)

Note that the sum of the two diagrams on the last row of Fig. 1(b) is suppressed by an additional factor $q^{2}$, therefore, they can be neglected within our present accuracy.

There is also a soft region, corresponding to all momentum transfers to the nucleus $q_{1},q_{2},q_{3}=q-q_{1}-q_{2}$ being small. Due to the gauge invariance of the light-by-light block it is easy to see that the corresponding contribution starts from $q^{4-4\epsilon}$, which is also too small for our present accuracy.

Eq. (2) shows that the potential $V(\bm{r})=e\int e^{i\bm{q}\bm{r}}\delta\tilde{\rho}/q^{2}d\bm{q}/(2\pi)^{3}$ is proportional to the delta function and the Lamb shift contribution can be written as:

$$\delta E=|\psi_{n\ell}(0)|^{2}\delta\widetilde{V}(0)=m\frac{\alpha^{2}(Z\alpha)^{6}}{\pi^{2}n^{3}}B\delta_{\ell,0}\,,$$

(3)

where $B$ is a numerical coefficient to be calculated, $m$ is the electron mass, $n$ and $\ell$ are the principal and angular quantum numbers, respectively.

The small-$q$ expansion of the diagrams in Fig. 1 can be expressed in terms of the integrals of the family

$$j(n_{1},\cdots,n_{9})=\int\frac{dq_{1}dq_{2}dl_{1}dl_{2}}{\pi^{2d}}\prod_{k=1}^{12}D_{k}^{-n_{k}}\times\prod_{s=13}^{14}\frac{\delta^{(n_{s}-1)}\left(D_{s}\right)}{(n_{s}-1)!}\,,$$

(4)

where

	$\displaystyle D_{1}=1-l_{1}^{2}\,,\quad D_{2}=1-l_{2}^{2}\,,\quad D_{3}=1-(l_{2}-q_{2})^{2}\,,$
	$\displaystyle D_{4}=1-(l_{2}+q_{1})^{2}\,,\quad D_{5}=1-(l_{1}+q_{1})^{2}\,,\quad D_{6}=-(l_{1}-l_{2})^{2}\,,$
	$\displaystyle D_{7}=-q_{1}^{2}\,,\quad D_{8}=-q_{2}^{2}\,,\quad D_{9}=-(q_{1}+q_{2})^{2}\,,\quad D_{10}=(l_{1}n)\,,$
	$\displaystyle D_{11}=(l_{2}n)\,,\quad D_{12}=(l_{1}-q_{2})^{2}\,,\quad D_{13}=(q_{1}n)\,,\quad D_{14}=(q_{2}n)\,.$		(5)

Here $n=(1,\bm{0})$ is a time-like unit vector and we put $m=1$. The functions $D_{1-9}$ and $D_{13-14}$ correspond to the denominators of the topology depicted in Fig. 2.

The $\delta$-functions in Eq. (4) secure the energy conservation. Note that $n_{10-12}\leqslant 0$ and the prescription $-i0$ for $D_{1-9}$ is implied.

Making the IBP reduction Chetyrkin1981 ; Tkachov1981 with LiteRed Lee2014 , we reveal $14$ master integrals, see Fig. 3.

Note that the counter-term diagrams in the last line of Fig. 1 can also be expressed in terms of the four-loop master integrals in Fig. 1 although they have only three loops. To this end we multiply the corresponding integrals by $1=\frac{-1}{\Gamma[1-d/2]}\int\frac{d^{d}l_{2}}{i\pi^{d/2}D_{2}}$. Then the contribution of counter-terms is expressed via the master integrals with unit mass tadpole loop, namely, via $j_{1}$ and $j_{5}$.

Since the integral family (4) contains no dimensionless free parameter, the differential equations method can not help directly. Therefore, there is a temptation to calculate the master integrals with the DRA method Lee2010 . Unfortunately, there is a nontrivial $2\times 2$ diagonal block in the matrix of dimensional recurrence, corresponding to the integrals $j_{8}$ and $j_{9}$ which belong to one and the same sector. Although it is, in principle, possible to apply the DRA method also for the cases with nontrivial diagonal $2\times 2$ blocks as discussed in Ref. Lee:2017ftw ; Lee:2017yex , its application in this case is much more laborious than that for the triangular matrix. Fortunately, for the present task we may apply a combined approach. First, we use the DRA method for all integrals but $j_{8}$, $j_{9}$, and $j_{14}$ (the latter integral belongs to a super-sector of the sector of $j_{8,9}$). In order to obtain information about analytical properties necessary for fixing periodic functions in homogeneous parts of the solutions, we use the approach of Ref. Lee:2022art . Namely, we choose the integrals which are obviously finite on a sufficiently wide vertical stripe in the complex plane of $d$. For example, in order to reveal the analytic properties of the most complicated integral $j_{13}$ we use finiteness of the integral $j_{13f}=\raisebox{-5.69054pt}{\includegraphics[width=56.9055pt]{Fig/j13finite.pdf}}$ on the stripe $\operatorname{Re}d\in(3,5]$. Reducing $j_{13f}$ to the master integrals, we obtain a number of nontrivial constraints for the leading expansion terms of $j_{13}$ on the basic stripe $(3,5]$. The results of the DRA approach have the form of $n$-fold triangular sums with factorized summand and $n\leqslant 3$, We use SummerTime package Lee:2015eva to calculate the $\epsilon$-expansion of these sums with high precision and PSLQ algorithm ferguson1999analysis to recognize the result in terms of multiple zeta values.

Note that all three remaining integrals, $j_{8}$, $j_{9}$, and $j_{14}$, contain two disjoint massive loops. In order to calculate those integrals, we consider a family of integrals with different masses in these two loops. This family is defined by Eqs. (4) and (5) where one should replace $D_{1}\to\widetilde{D}_{1}=m^{2}-l_{1}^{2}$, $D_{5}\to\widetilde{D}_{5}=m^{2}-(l_{1}+q_{1})^{2}$ and assume that $n_{6,7}\leqslant 0$ (in addition to $n_{10-12}\leqslant 0$). This family corresponds to the denominators of $j_{14}$, where the unit mass in the left-most fermion loop is replaced by $m$. Performing the IBP reduction, we reveal 14 master integrals depicted in Fig. 4. We obtain the differential system

$$\partial_{m}\tilde{j}=M(\epsilon,m)\tilde{j}$$

(6)

for the column $\tilde{j}=(\tilde{j}_{1},\ldots,\tilde{j}_{14})^{\intercal}$ and reduce it to $\epsilon$-form using Libra Lee:2020zfb . The general solution of Eq. (6) is expressed in terms of harmonic polylogarithms Remiddi:1999ew . The boundary conditions are put at $m\to 0$. We use expansion by regions Beneke:1997zp ; Pak:2010pt to fix the required coefficients in the asymptotic expansion of integrals. The original integrals $j_{8,9,14}$ are recovered as

$$j_{8}=\left.\tilde{j}_{8}\right|_{m=1},\quad j_{9}=\left.\tilde{j}_{9}\right|_{m=1},\quad j_{14}=\left.\tilde{j}_{14}\right|_{m=1}.$$

(7)

Note that other integrals in Fig. 3 which contain two disjoint fermion loops are also expressed in terms of the master integrals in Fig. 4:

$$j_{1}=\left.\tilde{j}_{1}\right|_{m=1},\quad j_{3}=\left.\tilde{j}_{2}\right|_{m=1},\quad j_{4}=\left.\tilde{j}_{3}\right|_{m=1}=\left.\tilde{j}_{4}\right|_{m=1},\quad j_{7}=\left.\tilde{j}_{5}\right|_{m=1},\quad j_{11}=\left.\tilde{j}_{12}\right|_{m=1}.$$

(8)

This provides a number of nontrivial cross checks of the obtained results for the master integrals.

Finally, we expand the diagrams in Fig. 1 in $q$ and perform the Dirac algebra using FeynCalc Shtabovenko:2020gxv . After the IBP reduction and substitution of the results for the master integrals, we obtain our final result for the coefficient $B$ in Eq. (3). We present the contributions of diagrams in Figs. 1(a) and 1(b) separately:

	$\displaystyle B_{\ref{fig:ph}}$	$\displaystyle=\frac{1456}{45}\mathrm{Li}_{4}\big{(}\tfrac{1}{2}\big{)}-\frac{4511\pi^{4}}{16200}+\frac{182\ln^{4}{2}}{135}+\frac{274}{135}\pi^{2}\ln^{2}{2}-\frac{2387\pi^{2}\ln{2}}{1080}+\frac{199\zeta_{3}}{45}$
		$\displaystyle+\frac{13057}{3240}+\frac{3703\pi^{2}}{5760}=0.125181281880322\ldots\,,$		(9)
	$\displaystyle B_{\ref{fig:po}}$	$\displaystyle=\frac{71\zeta_{3}}{56}-\frac{479}{756}+\frac{38401\pi^{2}}{217728}-\frac{283}{756}\pi^{2}\ln{2}=0.070271202837585\ldots\,,$		(10)
	$\displaystyle B$	$\displaystyle=B_{\ref{fig:ph}}+B_{\ref{fig:po}}=0.195452484717907\ldots\,.$		(11)

In the present paper we obtain the contributions of order $\alpha^{2}(Z\alpha)^{6}m$ to the Lamb shift from radiative corrections to the Wichmann-Kroll potential depicted in Fig. 1. Numerically our result (11) appears to be rather small and compatible with the heuristic estimate $B=0.13\pm 0.13$ of Ref. Karshenboim2019 . For the calculation of master integrals we use a combination of the DRA method and the approach based on the differential equations. This calculation provides yet another example of the effectiveness of multiloop methods for obtaining analytic results in atomic physics.

The work has been supported by Russian Science Foundation under grant 20-12-00205.