Mon, 16 Dec 2024

Three-loop moments and spectral density of photonpolarization function in QED

A. I. Onishchenko*1,2,3

1 Bogoliubov Laboratory of Theoretical Physics, Joint Institute for Nuclear Research, Dubna, Russia
2 Budker Institute of Nuclear Physics, Novosibirsk, Russia
3 Skobeltsyn Institute of Nuclear Physics, Moscow State University, Moscow, Russia

DOI: 10.54546/NaturalSciRev.100102

Abstract

We calculate three-loop photon spectral density in QED with Ndifferent species of electrons. The obtained results were expressed in terms of iterated integrals, which can be either reduced to Goncharov’s polylogarithms or written in terms of one-fold integrals of harmonic polylogarithms and complete elliptic integrals. In addition, we provide threshold and high-energy asymptotics of the calculated spectral density. It is shown that the use of the obtained spectral density correctly reproduces separately calculated moments of corresponding photon polarization operator.

Keywords: photon polarization operator, radiative corrections, quantum electrodynamics


1. Introduction

Recent advances in the precision of low-energy \(e^+e^−\) experiments (VEPP-2M, DAFNE,BEPC, PEP-II and KEKB colliders) call for a comparable precision of theoretical predictions.In particular, one has to compute complete NNLO corrections to both leptonic (for example,\(μ^+μ^−\)) and hadronic (for example,\(π^+π^−\)) production. In the present paper we take a first step towards calculation of NNLO QED corrections to the \(e^+e^− \to μ^+μ^−\) total production cross section, which is a key process for the center-of-mass energy calibration at present and future \(e^+e^−\) colliders. Specifically, we will be interested in contribution related to photon vacuum polarization. At three loops the photon spectral density contains NNLO contribution to \(e^+e^− \to μ^+μ^−\), NLO contribution to \(e^+e^− \to μ^+μ^− + γ\), and LO contributions to \(e^+e^− \to μ^+μ^− + 2γ\) and \(e^+e^− \to μ^+μ^−μ^+μ^−\) total production cross sections. All these contributions can be separated from each other, but in the present paper we will restrict ourselves only to their sum. Such a restricted setup nevertheless allows us to test different approximate threshold and high-energy expansions of photon spectral density by comparing them with exact results and make conclusions on the applicability of similar expansions for the calculation of full cross sections1. The latter should greatly reduce the complexity of future full calculations. Moreover,the provided techniques can be further used for the calculation of NNLO corrections to the \(e^+e^− \to μ^+μ^−\) production cross section in the framework of scalar QED.

The photon spectral density can be conveniently defined as a discontinuity of photon polarization operator \(Π(s)\). The latter is given by an extra factor \((1 + Π(s))\) in the denominator of renormalized photon propagator and is one of the several fundamental quantities arising ina study of quantum electrodynamics. At present we have exact results for one- and two-loop contributions [1,2]. However, starting from three loops there are only approximate results [36]. For example, the derivation of Baikov and Broadhurst [3] employs a simple Pad`e ap-proximation for three-loop contribution using a few terms of the asymptotic expansions near 3 special points: \(s= 0,4m^2,∞\) (\(m\) is the electron mass).

In the present paper we will provide for the first time exact as well as approximate threshold and high-energy results for three-loop photon spectral density. To check the obtained exact and asymptotic expressions for spectral density, the latter are used for the calculation of the moments of photon polarization operator. The comparison is performed with similar moments calculated from generalized Frobenius power series expansion of photon polarization operator at \(s= 0\). Also, the knowledge of exact photon spectral density already allows us to check the importance of the missed threshold \((s= 16m^2)\) in the reconstruction analysis of three-loop photon polarization operator performed by Baikov and Broadhurst.

2. Spectral density calculation

To calculate three-loop QED photon spectral density \(ρ(s)\), we followed standard procedure:

Figure 1. Diagrams contributing to three-loop photon self-energy in QED.

This way, considering QED with \(N\) electron flavors\(^4\) and using on-shell renormalization scheme\(^5\),we get

\begin{equation} \rho(s) = \rho^{(1)}(s) \frac{\alpha}{4\pi} + \rho^{(2)}(s) \left(\frac{\alpha}{4\pi}\right)^2 + \rho^{(3)}(s) \left(\frac{\alpha}{4\pi}\right)^3 + \dots, \tag{1} \end{equation}

where\(^6\) \( \left( \bar{s} = \frac{s}{m^2}, \quad \beta = \sqrt{1 - \frac{4}{\bar{s}}} \right) \),

\begin{align} \rho^{(1)}(s) &= \frac{4 N\pi (2+\bar{s})\beta }{3\bar{s}}\, ,\tag{2} \\[1pt] \rho^{(2)}(s) &= \frac{16 N\pi}{3}\Big\{ \beta^2 \left(2 \mathrm{II}_{l_1,l_2} + 2\mathrm{II}_{l_2,l_1}+\mathrm{II}_{l_1,l_0}+\mathrm{II}_{l_0,l_1}\right) -\frac{\beta (2+\bar{s})}{\bar{s}}\left( \mathrm{II}_{l_0}+2\mathrm{II}_{l_2}\right) \nonumber \\[1pt] & \quad\, - \frac{7+8\log 2 + \bar{s} (2-(3+\log 4)\bar{s})}{\bar{s}^2}\mathrm{II}_{l_1} + \frac{4\pi^2 (\bar{s}^2-4)-\beta\bar{s} (8 (2+\bar{s})\log 2 - 3 (6+\bar{s}))}{4\bar{s}^2} \Big\}. \tag{3} \end{align}

The three-loop contribution \(\rho^{(3)}(s)\) can be naturally separated into two pieces corresponding to \(2m\) and \(4m\) cuts, so that we have

\begin{equation} \rho^{(3)}(s) = \rho^{(3)}_{2m}(s) + \theta (s - 16 m^2)\rho^{(3)}_{4m}(s)\, , \tag{4} \end{equation}

where

\begin{align} \rho^{(3)}_{2m}(s) &= \frac{16\pi N^2}{3}\Big\{-\frac{4 c_{25}}{27}+\frac{2}{27} c_{30} \text{II}_{l_1}+\frac{4}{9} c_{27} \text{II}_{l_1,l_1}-\frac{2}{3} c_9 \text{II}_{l_1,l_1,l_1}\Big\} +\frac{16\pi N}{3}\Big\{\frac{c_{45}}{12}+\frac{1}{2} c_{36} \text{II}_{l_0}\nonumber \\[1pt] & \quad\, -\frac{1}{24} c_{46} \text{II}_{l_1}-c_{35} \text{II}_{l_2}+\frac{1}{6} c_{37} \text{II}_{l_0,l_1}-\frac{1}{4} c_{39} \text{II}_{l_1,l_0}-c_{34} \text{II}_{l_1,l_1}+\frac{1}{2} c_{38} \text{II}_{l_1,l_2}+4 c_{43} \text{II}_{l_2,l_1} \nonumber \\[1pt] &\quad\, +8 c_2 \left(\text{II}_{l_0,l_0}+2 \left(\text{II}_{l_0,l_2}+\text{II}_{l_2,l_0}+2 \text{II}_{l_2,l_2}\right)\right)+2 c_{15} \text{II}_{l_0,l_0,l_1}-2 c_5 \text{II}_{l_0,l_1,l_1} \nonumber \\[1pt] &\quad\, +4 c_3 \left(\text{II}_{l_0,l_1,l_0}+2 \text{II}_{l_0,l_1,l_2}\right)-2 c_{42} \text{II}_{l_1,l_0,l_1}+c_{44} \text{II}_{l_1,l_1,l_1}+4 c_{40} \left(\text{II}_{l_1,l_1,l_0}+2 \text{II}_{l_1,l_1,l_2}\right)\nonumber \\[1pt] &\quad\, -4 c_{41} \text{II}_{l_1,l_2,l_1}-4 c_6 \text{II}_{l_2,l_1,l_1}-8 c_{12} \left(\text{II}_{l_1,l_0,l_0}+2 \text{II}_{l_1,l_0,l_2}+2 \text{II}_{l_1,l_2,l_0}+4 \text{II}_{l_1,l_2,l_2} \right. \nonumber \\[1pt] &\quad\, +\left. \text{II}_{l_2,l_0,l_1} -\text{II}_{l_2,l_1,l_0}-2 \text{II}_{l_2,l_1,l_2}-2 \text{II}_{l_2,l_2,l_1}\right) -8 c_{17} \left(\text{II}_{l_0,l_2,l_1}+\text{II}_{l_0,l_1,l_1,l_1}\right) \nonumber \\[1pt] &\quad\, +c_{21} \text{II}_{l_1,l_1,l_0,l_1} +c_{11} \left(3 \left(2 \log (2) \text{II}_{l_4,l_1}+\text{II}_{l_4,l_1,l_0}+2 \text{II}_{l_4,l_1,l_2}\right)-\text{II}_{l_4,l_0,l_1}\right) \nonumber \\[1pt] &\quad\, +c_{10} (3 \log (2) \text{II}_{l_0,l_4,l_1}-\frac{1}{2} \text{II}_{l_0,l_4,l_0,l_1}+\frac{3}{2} \text{II}_{l_0,l_4,l_1,l_0}+3 \text{II}_{l_0,l_4,l_1,l_2})-2 c_{22} \text{II}_{l_1,l_2,l_1,l_1}\nonumber \\[1pt] &\quad\, +2 c_4 \left(\text{II}_{l_1,l_1,l_1,l_0}+2 \text{II}_{l_1,l_1,l_1,l_2}\right) +c_{26} \left(-18 \log (2) \text{II}_{l_1,l_4,l_1}-\text{II}_{l_1,l_0,l_0,l_1}-2 \text{II}_{l_1,l_0,l_1,l_0}\right.\nonumber \\[1pt] &\quad\, -4\text{II}_{l_1,l_0,l_1,l_2}+4 \text{II}_{l_1,l_0,l_2,l_1}+4 \text{II}_{l_1,l_1,l_0,l_0}+8 \text{II}_{l_1,l_1,l_0,l_2}+4 \text{II}_{l_1,l_1,l_1,l_1}+8 \text{II}_{l_1,l_1,l_2,l_0}\nonumber \\[1pt] &\quad\, +16\text{II}_{l_1,l_1,l_2,l_2}+4 \text{II}_{l_1,l_2,l_0,l_1}-4 \text{II}_{l_1,l_2,l_1,l_0}-8 \text{II}_{l_1,l_2,l_1,l_2}-8 \text{II}_{l_1,l_2,l_2,l_1}+3 \text{II}_{l_1,l_4,l_0,l_1}\nonumber \\[1pt] &\quad\, -\left. 9\text{II}_{l_1,l_4,l_1,l_0}-18 \text{II}_{l_1,l_4,l_1,l_2}\right)-c_{23} \text{II}_{l_1,l_0,l_1,l_1}+2 c_{20} \text{II}_{l_1,l_1,l_2,l_1}\Big\} \tag{5} \end{align}

and

\begin{align} \rho^{(3)}_{4m}(s) &= \frac{2\pi N^2}{9} \Big\{-\frac{32}{27} c_{24} f''\left(\bar{s}\right)+\frac{4}{27} c_{32} f'\left(\bar{s}\right)-\frac{1}{9} c_{33} f\left(\bar{s}\right)+\frac{1}{3} c_9 \left(2 \text{II}_{r_2}-\text{II}_{l_1,\tilde{r}_3}\right)-\frac{2}{9} c_{28} \text{II}_{\tilde{r}_3}\Big\}\nonumber \\[1pt] & \quad\, +\frac{2\pi N}{9} \Big\{-48 c_1 f''\left(\bar{s}\right)+2 c_{29} f'\left(\bar{s}\right)-\frac{1}{2} c_{31} f\left(\bar{s}\right)+2 c_7 \left(\text{II}_{r_1}-\text{II}_{l_0,\tilde{r}_3}\right) +3 c_{14} \text{II}_{\tilde{r}_3} \nonumber \\[1pt] &\quad\, +c_{19} \left(\frac{3}{2} \text{II}_{l_1,\tilde{r}_3}-3 \text{II}_{r_2}\right)+c_8 \left(2 \text{II}_{l_2,\tilde{r}_3}-4 \text{II}_{r_3}\right)+c_{13} \left(2 \text{II}_{l_0,l_1,\tilde{r}_3}-4 \text{II}_{l_0,r_2}\right) +2 c_{16} \text{II}_{r_0} \nonumber \\[1pt] &\quad\, +c_{26} \left(-2 \text{II}_{l_1,l_1,\tilde{r}_3}-\text{II}_{l_1,r_0}+4 \text{II}_{l_1,r_2}\right)+2 c_{18} \left(-\text{II}_{l_1,l_0,\tilde{r}_3}+\text{II}_{l_1,l_2,\tilde{r}_3}+\text{II}_{l_1,r_1}- 2 \text{II}_{l_1,r_3}\right)\Big\}. \tag{6} \end{align}

Here, \(f(\bar{s})\) is the function given in Appendix B and defined in terms of products of complete elliptic integrals of first kind. The expressions for \(c_i\) coefficients can be found in Appendix C. The iterated integrals with \(l_i\) weights only can be straightforwardly rewritten in terms of Goncharov’s multiple polylogarithms, while those with elliptic kernels (\(r_i\) or \( r̃_i\) weights) in terms of one-fold integrals of harmonic polylogarithms and complete elliptic integrals as shown in [10]. The latter representations are already well suited for numerical evaluations. However, as we will see in the next section, much more faster numerics in the whole range of \(s\) values can be obtained with threshold and high-energy expansions of spectral densities.

3. Asymptotics and checks

The asymptotic expansions of obtained spectral densities can be done through the asymptotic expansions of iterated integrals as described in [10]. While it is straightforward to do asymptotic expansions in the threshold cases, the high-energy expansions are more involved and may require finding PLSQ relations [12] for polylogarithmic constants at unity argument.To avoid this, one can obtain required asymptotic expansions for master integrals themselves through the Frobenius solution of corresponding differential equations\(^7\). The calculation of spectral densities with the asymptotic expressions for cut master integrals then gives

\begin{equation} \rho^{(1)}_{\text{thr}}(s) = \frac{2}{3} N \pi \beta (3 - \beta^2), \tag{7} \end{equation}
\begin{equation} \rho^{(1)}_{\rm high} (\bar{s}) = \frac{4N\pi}{3} \left\{ 1 - \frac{6}{\bar{s}^2} - \frac{8}{\bar{s}^3} - \frac{18}{\bar{s}^4} \right\} + \mathcal{O}\left(\frac{1}{\bar{s}^5}\right), \tag{8} \end{equation}
\begin{equation}\rho^{(2)}_{{\rm thr}}(\bar{s}) = 4N\pi\Big\{ \pi^2 - 8\beta +\frac{2\pi^2\beta^2}{3}-\frac{\pi^2\beta^4}{3}+\frac{4}{9}\beta^3 (-37+36\log 2+24\log\beta) \Big\} + \mathcal{O}\left(\beta^5\right) , \tag{9}\end{equation}
\begin{equation} \begin{aligned} \rho^{(2)}_{\rm high}(\bar{s}) &= 4N\pi \Big\{ 1 + \frac{12}{\bar{s}} + \frac{2(5 + 12\log\bar{s})}{\bar{s}^2} \\ & \quad\, + \frac{16(-47 + 87\log\bar{s})}{27\bar{s}^3} + \frac{-983 + 1218\log\bar{s}}{9\bar{s}^4} \Big\} + \mathcal{O}\left(\frac{1}{\bar{s}^5}\right), \end{aligned} \tag{10} \end{equation}
\begin{equation} \begin{aligned} \rho^{(3)}_{2m, {\rm thr}}(\bar{s}) &= \frac{32N^2\pi}{9} \Big\{ 4(11 - \pi^2)\beta - \frac{1}{9}(245 - 24\pi^2)\beta^3 \Big\} \\ &\quad + \frac{8N\pi}{9} \Big\{ \frac{3\pi^4}{\beta} - 72\pi^2 + (351 - 70\pi^2 + 5\pi^4 + 48\pi^2\log 2 - 24\pi^2\log\beta - 36\zeta_3)\beta \\ &\quad\quad + 2(-43\pi^2 + 24\pi^2\log 2 + 48\pi^2\log\beta)\beta^2 \\ &\quad\quad + (1411 + 51\pi^2 + \pi^4 - 1152\log 2 - 42\pi^2\log 2 \\ &\quad\quad\quad - (768 - 20\pi^2)\log\beta + 117\zeta_3)\beta^3 \\ &\quad\quad + \frac{8}{75}(-947\pi^2 + 480\pi^2\log 2 + 960\pi^2\log\beta)\beta^4 \Big\} + \mathcal{O}(\beta^5), \end{aligned} \tag{11} \end{equation}
\begin {align}\rho^{(3)}_{4m, {\rm thr}}(s) &= \frac{\pi^2 N^2 (\bar{s}-16)^{9/2}}{516096}\Big\{ 1 - \frac{629 (\bar{s}-16)}{2640} + \frac{10243 (\bar{s}-16)^2}{274560} -\frac{7973 (\bar{s}-16)^3}{1647360}\Big\}\nonumber\\[1pt] &\quad\, + \frac{\pi^2 N (\bar{s}-16)^{9/2}}{5160960}\Big\{ 1 - \frac{7(\bar{s}-16)}{48} + \frac{307 (\bar{s}-16)^2}{27456} - \frac{193 (\bar{s}-16)^3}{658944} \Big\}\nonumber\\[1pt] &\quad\, + \mathcal{O} \left(\bar{s}-16\right)^{17/2} , \tag{12} \end{align}
\begin{align} \rho^{(3)}_{2m, {\rm high}}(s) &= \frac{16N^2\pi}{81} \Big\{ 766 - 66\pi^2 - 265\log\bar{s} + 12\pi^2\log\bar{s} + 57\log^2\bar{s} - 6\log^3\bar{s} \nonumber \\[1pt] &\quad\, + \frac{2}{\bar{s}}(-65+24\pi^2-216\log\bar{s} + 108\log^2\bar{s}) \Big\} + \frac{2N\pi}{135}\Big\{ 16065 - 900\pi^2 + 76\pi^4 \nonumber \\[1pt] &\quad\, - 1440\pi^2\log 2 - 1440\zeta_3 + \log\bar{s} (-2340+360\pi^2-1440\zeta_3) + \frac{1}{\bar{s}} (-35100 + 3600\pi^2 \nonumber \\[1pt] & \quad\, - 114\pi^4 - 11520\zeta_3 + \log\bar{s} (-7200+240\pi^2+1440\zeta^3) + (-7200+60\pi^2)\log^2\bar{s} \nonumber \\[1pt] &\quad\, + 240\log^3\bar{s} - 30\log^4\bar{s}) \Big\} + \mathcal{O} \left(\frac{1}{\bar{s}^2}\right) , \tag{13}\end{align}
\begin{align} \\[1pt] \rho^{(3)}_{4m,{\rm high}}(s) &= \frac{8N^2\pi}{81}\Big\{ -1829 + 132\pi^2 + 216\zeta_3 + (584-24\pi^2)\log\bar{s} - 114\log^2\bar{s} + 12\log^3\bar{s} \nonumber \\[1pt] & \quad\, + \frac{1}{\bar{s}} (-1144-96\pi^2+1512\log\bar{s} - 432\log^2\bar{s}) \Big\} + \frac{4N\pi}{135}\Big\{ -8100+450\pi^2\nonumber \\[1pt] & \quad\, -38\pi^4+720\pi^2\log 2 + 720\zeta_3 + (1170-180\pi^2+720\zeta_3)\log\bar{s} + \frac{1}{\bar{s}}( 18360 \nonumber \\[1pt] & \quad\, -1800\pi^2 + 57\pi^4 + 5760\zeta_3 - (6120+120\pi^2+720\zeta_3)\log\bar{s} \nonumber \\[1pt] &\quad\, + (3600-30\pi^2)\log^2\bar{s} - 120\log^3\bar{s} + 15\log^4\bar{s}) \Big\} + \mathcal{O} \left(\frac{1}{\bar{s}^2}\right), \tag{14}\end{align}

and results with more terms in the expansions can be found in accompanying Mathematica notebook. Using the latter, it is easy to get a very accurate representation of the above spectral densities for the whole range of \(\bar s\) values. The corresponding plot of different contributions can be found in Fig. 2. It is also interesting to compare the ratio of \(4m\) and \(2m\) cuts contributions to three-loop spectral density in the region of \(\bar s\) values where mass effects become important. The corresponding plot can be found in Fig. 3. From the latter we may conclude that the account of second threshold \(\bar s = 16\) missed in the reconstruction of photon polarization operator performed in [3] is not actually important.

Figure 2. Values of one-, two- and three-loop spectral densities. Here \( \rho^{(i), N} \) and \( \rho^{(i), N^2} \) are coefficients in front of \(N\) and \(N^2\) contributions to full spectral density \( \rho^{(i)} = \rho^{(i), N} N + \rho^{(i), N^2} N^2\).

To check the obtained results for the spectral densities, we first performed the whole calculation in an arbitrary gauge and made sure that the gauge dependence of spectral densities

Figure 3. Ratio of \(4m\) and \(2m\) cuts contributions to three-loop spectral density. Here \( \rho_*^{(i), N} \) and \( \rho_*^{(i), N^2} \) denote \(N\) and \( N^2 \) contributions to spectral density \( \rho_*^{(i)}:\rho_*^{(i)}= \rho_*^{(i), N} N + \rho_*^{(i), N^2} N^2 \).

cancels. Second, we have numerically checked that the first three moments\(^8\) \(M_n\) of polarization operator computed in [3] agree with a very high precision (at least 15 digits) with moments obtained with the use of our spectral densities and dispersion relation

\begin{equation}\Pi(s) = \frac{1}{\pi} \int_{4m^2}^{\infty} \frac{\rho(s')}{s' - s} . \tag {15}\end{equation}

Moreover, we performed a comparison with additionally calculated 97 moments, see Appendix D and accompanying Mathematica notebook. This latter calculation is similar to the one we did for spectral densities except that in this case we used uncut master integrals. The Frobenius solutions for the latter can be easily obtained from corresponding differential equations using standard techniques. Next, with the computed spectral densities, we can deduce approximate analytical expressions for moments with large n values. Indeed, making variable change in the dispersion relation from \(s'\) to \(\beta' = 1 - \sqrt{{4m^2}/{s'}}\), it is easy to see that

\begin{equation} M_n^{(i)} \simeq \int\limits_0^1 d\beta' \frac{2\beta'}{\pi} (1-\beta'^2)^{n-1} \rho^{(i)}_{\rm thr}(\beta')\,\tag {16}\end{equation}

where for large n values the largest contribution to the integral comes from small \(\beta'\) values and for three-loop spectral density it is sufficient to consider only \(2m\) cut contribution. At one loop due to finite \(\beta'\) expansion of spectral density, this expression is actually exact and we have

\begin{equation} M_n^{(1)} = \frac{(1 + n)\sqrt{\pi} \, \Gamma(n)}{\Gamma\left(\frac{5}{2} + n\right)} \, N. \tag{17}\end{equation}

At two and three loops the formulae are approximate and are given by

\begin{equation} M_n^{(2)} \simeq \sum_{i=0}^{\text{order}_\beta} \frac{\Gamma\left(\frac{1+i}{2}\right)\Gamma(n)}{4\Gamma\left(\frac{1+i+2n}{2}\right)} \left\{ 2 d_{1,i}^{(2)} + d_{2,i}^{(2)} \left[ \Psi^{(0)}\left(\frac{1+i}{2}\right) - \Psi^{(0)}\left(\frac{1+i+2n}{2}\right) \right] \right\} \tag{18} \end{equation}
\begin{align} M_n^{(3)}\simeq &\sum_{i=0}^{\text{order}_\beta}\frac{\Gamma (\frac{1+i}{2})\Gamma (n)}{8 \Gamma (\frac{1+i+2n}{2})} \Big\{ 4 d_{1,i}^{(3)} + 2 d_{2,i}^{(3)} \left[ \Psi^{(0)} \left(\frac{1+i}{2}\right) - \Psi^{(0)}\left( \frac{1+i+2n}{2} \right) \right] \nonumber \\[2pt] & + d_{3,i}^{(3)} \Big[ \Big( \Psi^{(0)}\Big( \frac{1+i}{2} \Big) - \Psi^{(0)}\Big( \frac{1+i+2n}{2} \Big) \Big)^2 + \Psi^{(1)}\Big(\frac{1+i}{2}\Big) - \Psi^{(1)}\Big( \frac{1+i+2n}{2} \Big) \Big] \Big\} , \label{eqn:moment3loop-largen} \tag{19} \end{align}

where \(d\) coefficients are defined as

\begin{equation} \frac{2\beta}{\pi} \, \rho_\mathrm{thr}^{(2)}(\beta) \simeq \sum_{i=0}^{\text{order}_\beta} d^{(2)}_{1,i} \beta^i + d^{(2)}_{2,i} \beta^i \log \beta, \tag{20}\end{equation}
\begin{equation} \frac{2\beta}{\pi} \, \rho^{(3)}_{\mathrm{thr}}(\beta) \simeq \sum_{i=0}^{\text{order}_\beta} d^{(3)}_{1,i} \beta^i + d^{(3)}_{2,i} \beta^i \log \beta + d^{(3)}_{3,i} \beta^i \log (\beta^2). \tag{21}\end{equation}

Here, order\(_\beta\) is the expansion order in \(\beta\) of corresponding spectral densities, and

\begin{equation}\Psi^{(\nu)}(z)=\frac{d^{v+1}\log \Gamma(z)}{dz^{n+1}}\tag{22}\end{equation}

is the usual polygamma function. With order\(_\beta\) = 80 the approximate expression for \(M_1^{(2)}\) is accurate with 0.04 percent level, \(M_{10}^{(2)}\) with \(10^{−10}\) percent level and \(M_{50}^{(2)}\) with \(10^{−26}\) percent level. For order\(_\beta\) = 10 the corresponding errors are \(2,10^{(−3)}\) and \(10^{−6}\) percent. With order\(_\beta\) = 120 the accuracy for approximate expression of \(M_1^{(3)}\) is 0.4 percent, for \(M_{10}^{(3)}\) it is \(10^{−10}\) percent, and for \(M_{50}^{(3)}\) it is \(10^{−32}\) percent. For lower value of order\(_\beta\) = 10 the corresponding values are \(13, 10^{−30}\) and \(10^{−7}\) percent. So, we see that for sufficiently large values of moments our approximate formulae are very accurate except for a few first moments.

4. Conclusion

In the present paper we have performed calculation of three-loop photon spectral density in the framework of QED with \(N\) different electron flavors. The obtained results contain both exact and asymptotic expressions. The exact results were written in terms of iterated integrals, which reduce either to Goncharov’s polylogarithms or to one-fold integrals of harmonic polylogarithms and complete elliptic integrals. The asymptotic expressions, on the other hand, allow us to have very accurate expressions for the spectral density in the whole s range. It is shown that the obtained spectral density correctly reproduces first hundred moments of the photon polarization operator. In addition, we supply approximate analytical formulae for the moments of the photon polarization operator. The latter are very accurate for almost all moments except maybe the first few.

Finally, we would like to note that the performed calculation showed that for practical purposes it is sufficient to know required master integrals in terms of their asymptotic threshold and high-energy expansions. The latter can be easily obtained from generalized Frobenius solutions to corresponding differential equations. This gives us a hope that similar approximate solutions will also be applicable to other master integrals required to obtain full NNLO contribution for \(e^+e^− \to μ^+μ^−\) total production cross section.

Acknowledgments

The author is grateful to R. N. Lee for interest in this work and valuable discussions. This work was supported by the Russian Science Foundation, grant 20-12-00205.

Conflict of Interest

The author declares no conflict of interest.

Appendix A. On-shell renormalization constants

The required on-shell renormalization constants for QED with N similar electron species can be extracted from renormalization of two-loop photon and electron self-energies. For photon wave function renormalization \(A_0=Z_AA\) we have \( \left( \alpha = \frac{e^2}{4\pi} \right) \)

\begin{equation} Z_A = 1 + Z_{A}^{(1)} \cdot \left(\frac{\alpha \mu^{2\varepsilon} e^{-\varepsilon \gamma_E}}{(4\pi)^{1 - \varepsilon} m^{2\varepsilon}}\right) + Z_{A}^{(2)} \cdot \left(\frac{\alpha \mu^{2\varepsilon} e^{-\varepsilon \gamma_E}}{(4\pi)^{1 - \varepsilon} m^{2\varepsilon}}\right)^2 + \dots, \tag{A.1} \end{equation}

where

\begin{align} Z_A^{(1)} &= N\Big\{ -\frac{4}{3\varepsilon} -\frac{\pi^2\varepsilon}{9} + \frac{4}{9}\varepsilon^2\zeta_3- \frac{\pi^4\varepsilon^3}{120} + \frac{\varepsilon^4}{135} (5\pi^2\zeta_3 + 36\zeta_5) \Big\}+ \ldots \, , \tag{A.2} \end{align}
\begin{equation} Z_A^{(2)} = N \left( -\frac{2}{\varepsilon} - 15 - \frac{1}{6}(93 + 2\pi^2)\varepsilon + \varepsilon^2 \left( -\frac{223}{4} - \frac{5\pi^2}{2} + \frac{4\zeta_3}{3} \right) \right) + \ldots \tag{A.3} \end{equation}

Here \(m\) is the electron pole mass, \(\mu\) is the dimensional parameter entering dimensional regularization, and \(e\) is the on-shell electron charge at \(q^2 = 0\). The charge renormalization constant \( \left( \alpha_0 = Z_\alpha \alpha \mu^{2\varepsilon} \right) \) is then obtained using Ward identity as

\begin{equation} Z_{\alpha} = Z_{A}^{-1} \tag{A.4} \end{equation}

Similarly, from electron self-energy for on-shell mass renormalization constant \( \left( m_0 = Z_m m \right) \) we get

\begin{equation} Z_m^{(1)} = 1 + Z_m^{(1)} \cdot \left( \frac{\alpha \mu^{2\varepsilon} e^{-\varepsilon \gamma_E}}{(4\pi)^{1 - \varepsilon} m^{2\varepsilon}} \right) + Z_m^{(2)} \cdot \left( \frac{\alpha \mu^{2\varepsilon} e^{-\varepsilon \gamma_E}}{(4\pi)^{1 - \varepsilon} m^{2\varepsilon}} \right) ^2 + \dots , \tag{A.5} \end{equation}

where

\begin{align} Z_m^{(1)} ={} & -\frac{3}{\varepsilon} - 4 -\left(8+\frac{\pi^2}{4}\right)\varepsilon + \varepsilon^2 \left(-16-\frac{\pi^2}{3}+\zeta_3\right) + \varepsilon^3 \left(-32-\frac{2\pi^2}{3}-\frac{3\pi^4}{160}+\frac{4\zeta_3}{3}\right) \\ & + \varepsilon^4 \left(-64 - \frac{4\pi^2}{3} - \frac{\pi^4}{40} + \frac{8\zeta_3}{3} + \frac{1}{12}\pi^2\zeta_3 + \frac{3\zeta_5}{5}\right) + \ldots, \tag{A.6} \end{align}
\begin{align} Z_m^{(2)} ={} & \frac{9}{2\varepsilon^2} + \frac{45}{4\varepsilon} + \frac{199}{8} - \frac{17\pi^2}{4} + 8\pi^2\log 2 - 12\zeta_3 \notag \\ & + \varepsilon \left( \frac{677}{16} - \frac{205\pi^2}{8} + \frac{14\pi^4}{5} + 48\pi^2\log 2 - 16\pi^2\log^2 2 - 8\log^4 2 - 192 a_4 - 135\zeta_3 \right) \notag \\ & + \varepsilon^2 \left( -\frac{110857}{128} + \frac{7039\pi^2}{192} + \frac{209\pi^4}{64} - 54\pi^2\log 2 + 8\pi^2\log^2 2 + 4\log^4 2 + 96a_4 + \frac{4915}{24}\zeta_3 \right) \notag \\ & + N \left( -\frac{2}{\varepsilon^2} + \frac{5}{3\varepsilon} + \frac{143}{6} - 3\pi^2 \right. \notag \\ & \quad + \varepsilon \left. \left( \frac{1133}{12} - \frac{235\pi^2}{18} + 16\pi^2\log 2 - \frac{164}{3}\zeta_3 \right) \right) + \ldots \tag{A.7} \end{align}

and \(a_4 = \text{Li}_4\left(\frac{1}{2}\right) = \sum_{n=1}^{\infty} \frac{1}{2^n n^4}\). Note that at \(N = 1\) these renormalization constants reduce to already known QED values [1315].

Appendix B. Notation for iterated integrals

The iterated integrals present in the current paper involve the following set of integration kernels or weights:

\begin{equation}l_0(\bar{s}) = \frac{1}{\bar{s}}, \quad l_1(\bar{s}) = \frac{1}{\bar{s}\beta}, \quad l_2(\bar{s}) = \frac{1}{\bar{s} - 4}, \quad l_4(\bar{s}) = \frac{1}{\bar{s} - 1}, \tag{B.1}\end{equation}
\begin{equation} r_k(\bar{s}) = \frac{f(\bar{s})}{\bar{s}\beta^k} \theta(\bar{s} - 16), \quad (k = 0, 1, 2, 3), \tag{B.2}\end{equation}
\begin{equation} r_3(\bar{s}) = \frac{8\beta(\bar{s} + 2) f(\bar{s})}{(\bar{s} - 16)(\bar{s} - 4)^2} \theta(\bar{s} - 16), \tag{B.3}\end{equation}

where

\begin{gather}\label{eq:fdef} f(\bar{s})=\frac{16 (\bar{s}-16)}{\bar{s}}\left[ \mathrm{K}(1-k_{-})\mathrm{K}(k_{+}) - \mathrm{K}(k_{-})\mathrm{K}(1-k_{+}) \right], \tag{B.4} \end{gather}
\begin{equation} \\[1pt] k_{\pm}=\frac12\left[1\pm\left(1-\frac{8}{\bar{s}}\right) \sqrt{1-\frac{16}{\bar{s}}}+\frac{16}{\bar{s}} \sqrt{1-\frac{4}{\bar{s}}}\right], \tag{B.5}\end{equation}

and \(K\) is complete elliptic integral of first kind. The required iterated integrals are then defined as

\begin{equation} \Pi_{w_n,\ldots w_1} = I(w_n,\ldots w_1|\bar{s})= \int\limits_{\bar{s}>s_n>\ldots >s_1>4} \prod_{k=1}^{n} ds_k w_k(s_k)\, , \tag{B.6}\end{equation}

where weights \(w_k(s)\) take values from the set above.

Appendix C. Expressions for \(c_i\) coefficients

For \(c_i\) coefficients present in the expression for three-loop photon spectral density, we have the following expressions:

\begin{align} c_1 &= -\frac{29}{8 (\beta +1)}+\frac{15}{4 (\beta +1)^2}+\frac{29}{8 (\beta -1)}+\frac{15}{4 (\beta -1)^2}-\frac{1}{4}\, , \; c_2 = \frac{3 \beta }{4}-\frac{\beta ^3}{4}\, ,\tag{C.1}\end{align}
\begin{equation} c_3 = -\frac{\beta^4}{4} + \frac{\beta^2}{2} - \frac{1}{4} \, ,c_4 = -\frac{3 \beta ^4}{4}+\frac{7 \beta ^2}{2}-\frac{11}{4}\, , \; c_5 = -\frac{69 \beta ^3}{4}+\frac{113 \beta }{4}+\frac{3}{\beta }\, , \\[1pt]\tag{C.2}\end{equation}
\begin{equation} c_6 = -\frac{15 \beta ^3}{2}+\frac{35 \beta }{4}+\frac{3}{4 \beta }\, , \; c_7 = -\frac{15 \beta ^3}{4}+8 \beta +\frac{3}{4 \beta }\, , \; c_8 = -\frac{9 \beta ^3}{4}+\frac{7 \beta }{2}+\frac{3}{4 \beta }\, , \\[1pt]\tag{C.3}\end{equation}
\begin{equation} c_9 = -\frac{5 \beta ^4}{4}+\frac{5 \beta ^2}{2}+\frac{3}{4}\, , \; c_{10} = -\frac{3 \beta ^4}{4}+\frac{\beta ^2}{2}+\frac{17}{4}\, , \; c_{11} = -\frac{\beta ^4}{2}+\frac{11 \beta ^2}{4}+\frac{23}{4}\, , \\[1pt]\tag{C.4}\end{equation}
\begin{equation} c_{12} = -\frac{\beta ^4}{4}+\frac{\beta ^2}{2}+\frac{3}{4}\, , \; c_{13} = -\frac{\beta ^4}{4}+\frac{\beta ^2}{2}+\frac{7}{4}\, , \; c_{14} = -\frac{\beta ^3}{4}+\frac{3 \beta }{4}+\frac{2}{\beta }\, ,\tag{C.5}\end{equation}
\begin{equation} c_{15} = -\frac{\beta ^4}{4}+\beta ^2+\frac{5}{4}\, , \; c_{16} = \frac{9}{4}-\frac{\beta ^4}{4}\, ,\; c_{17} = \frac{5}{4}-\frac{\beta ^4}{4}\, , \; c_{18} = -\frac{\beta ^4}{4}+\frac{3 \beta ^2}{4}-\frac{3}{2}\, , \\ \tag{C.6}\end{equation}
\begin{equation} c_{19} = -\frac{5 \beta ^4}{4}+\frac{9 \beta ^2}{2}-\frac{5}{4}\, , \; c_{20} = -\frac{7 \beta ^4}{4}+\frac{13 \beta ^2}{2}-\frac{35}{4}\, , \;c_{21} = -\frac{9 \beta ^4}{4}+8 \beta ^2-\frac{47}{4}\, ,\tag{C.7}\end{equation}
\begin{equation} \\[1pt] c_{22} = -\frac{11 \beta ^4}{4}+12 \beta ^2-\frac{45}{4}\, , \; c_{23} = -\frac{17 \beta ^4}{4}+\frac{33 \beta ^2}{2}-\frac{81}{4}\, , \\[1pt]\tag{C.8}\end{equation}
\begin{equation} c_{24} = -\frac{407}{16 (\beta +1)}+\frac{1829}{16 (\beta +1)^2}+\frac{407}{16 (\beta -1)}+\frac{1829}{16 (\beta -1)^2}-\frac{771}{4}\, ,\\[1pt]\tag{C.9}\end{equation}
\begin{equation} c_{25} = -\frac{1}{4} \left(15 \pi ^2-236\right) \beta ^3+\frac{3}{4} \left(27 \pi ^2-340\right) \beta +\frac{9}{2 \beta }\, , \; c_{26} = -\frac{\beta ^5}{4}+\frac{\beta ^3}{4}+\frac{5 \beta }{4}+\frac{3}{4 \beta }\, , \\[1pt]\tag{C.10}\end{equation}
\begin{equation} c_{27} = -\frac{15 \beta ^3}{4}-\frac{3}{4 \beta ^3}+4 \beta +\frac{10}{\beta }\, , \; c_{28} = -\frac{39 \beta ^3}{4}-\frac{3}{4 \beta ^3}-\frac{19 \beta }{2}-\frac{7}{2 \beta }\, , \\[1pt]\tag{C.11}\end{equation}
\begin{equation} c_{29} = -\frac{39 \beta ^2}{2}+\frac{6507}{16 \left(4 \beta ^2-3\right)}+\frac{77}{2 (\beta -1)}-\frac{77}{2 (\beta +1)}+\frac{1997}{16}\, , \\[1pt]\tag{C.12}\end{equation}
\begin{equation} c_{30} = -\frac{1}{4} \left(15 \pi ^2-236\right) \beta ^4+\frac{1}{2} \left(15 \pi ^2-167\right) \beta ^2+\frac{9}{\beta ^2}+\frac{9}{4} \left(\pi ^2-52\right) , \\[1pt]\tag{C.13}\end{equation}
\begin{equation} c_{31} = -2 \beta ^4-\frac{5 \beta ^2}{8}-\frac{567}{4 \left(4 \beta ^2-3\right)}+\frac{6507}{64 \left(4 \beta ^2-3\right)^2}+\frac{24}{\beta ^2}+\frac{1197}{64}\, , \\[1pt]\tag{C.14}\end{equation}
\begin{equation} c_{32} = -\frac{967 \beta ^2}{2}+\frac{72063}{16 \left(4 \beta ^2-3\right)}+\frac{18}{\beta ^2}+\frac{1245}{4 (\beta -1)}-\frac{1245}{4 (\beta +1)}+\frac{21117}{16}\, , \\[1pt]\tag{C.15}\end{equation}
\begin{equation} c_{33} = -4 \beta ^4+\frac{3}{\beta ^4}+\frac{481 \beta ^2}{8}-\frac{1693}{4 \left(4 \beta ^2-3\right)}+\frac{24021}{64 \left(4 \beta ^2-3\right)^2}+\frac{11}{\beta ^2}+\frac{1011}{64}\, , \\[1pt]\tag{C.16}\end{equation}
\begin{align} c_{34} &= -\frac{1}{2} \beta ^5 \left(-6+3 \pi ^2-4 \log ^2(2)\right)-\frac{7 \pi ^2 \beta ^4}{4}+\frac{1}{4} \beta ^3 \left(191+6 \pi ^2-8 \log ^2(2)+24 \log (2)\right)\nonumber \\[1pt] &\quad\, +\frac{13 \pi ^2 \beta ^2}{2}+\frac{\beta}{2} \left(-306+15 \pi ^2-20 \log ^2(2)+20 \log (2)\right)\nonumber \\[1pt] &\quad\, +\frac{3 \left(1+6 \pi ^2-8 \log ^2(2)\right)}{4 \beta }-\frac{35 \pi ^2}{4}\, , \tag{C.17}\end{align}
\begin{align} c_{35} &= -2 \pi ^2 \beta ^4+\frac{1}{4} \beta ^3 \left(9+14 \pi ^2+32 \log (2)\right)+4 \pi ^2 \beta ^2\nonumber \\[1pt] &\quad\, +\frac{24 \beta }{\beta ^2+3}-\frac{32 \beta }{\left(\beta ^2+3\right)^2}-\frac{1}{4} \beta \left(-9+14 \pi ^2+96 \log (2)\right)+6 \pi ^2\, , \tag{C.18}\end{align}
\begin{align} c_{36} &= -2 \pi ^2 \beta ^4-\frac{1}{4} \beta ^3 \left(9+66 \pi ^2+32 \log (2)\right)-\frac{24 \beta }{\beta ^2+3}+\frac{32 \beta }{\left(\beta ^2+3\right)^2}\nonumber \\[1pt] & \quad\, +\frac{3 \pi ^2}{\beta }+\frac{1}{4} \beta \left(-9+118 \pi ^2+96 \log (2)\right)+10 \pi ^2\, , \tag{C.19}\end{align}
\begin{align} c_{37} &= -\frac{1}{4} \beta ^4 \left(-213+4 \pi ^2+48 \log (2)\right)-\frac{132}{\beta ^2+3}+\frac{432}{\left(\beta ^2+3\right)^2}\nonumber \\[1pt] &\quad\, -\frac{384}{\left(\beta ^2+3\right)^3}-2 \beta ^2 \left(39+\pi ^2-12 \log (2)\right)+\frac{3}{4} \left(-295+4 \pi ^2-16 \log (2)\right),\tag{C.20}\end{align}
\begin{align} c_{38} &= -2 \pi ^2 \beta ^5-\frac{1}{4} \beta ^4 \left(-41+6 \pi ^2-64 \log (2)\right)+2 \pi ^2 \beta ^3+\frac{264}{\beta ^2+3} -\frac{864}{\left(\beta ^2+3\right)^2}+\frac{768}{\left(\beta ^2+3\right)^3} \nonumber \\[1pt] & \quad\, +\frac{1}{2} \beta ^2 \left(3+14 \pi ^2-64 \log (2)\right)+10 \pi ^2 \beta +\frac{6 \pi ^2}{\beta }+\frac{1}{4} \left(-247-22 \pi ^2-192 \log (2)\right) , \tag{C.21}\end{align}
\begin{align} c_{39} &= -2 \pi ^2 \beta ^5+\frac{1}{4} \beta ^4 \left(-41+14 \pi ^2-64 \log (2)\right)+2 \pi ^2 \beta ^3-\frac{264}{\beta ^2+3}+\frac{864}{\left(\beta ^2+3\right)^2}-\frac{768}{\left(\beta ^2+3\right)^3}\nonumber \\[1pt] & \quad\, -\frac{1}{2} \beta ^2 \left(3+26 \pi ^2-64 \log (2)\right)+10 \pi ^2 \beta +\frac{6 \pi ^2}{\beta }+\frac{1}{4} \left(247+70 \pi ^2+192 \log (2)\right) , \tag{C.22}\end{align}
\begin{equation} c_{40} = -\frac{1}{2} \beta ^5 \log (2)+\frac{1}{4} \beta ^3 (2 \log (2)-3)+\frac{5}{4} \beta (2 \log (2)-1)+\frac{3 \log (2)}{2 \beta }\, ,\tag{C.23}\end{equation}
\begin{equation} c_{41} = -\frac{1}{2} \beta ^5 \log (2)+\frac{1}{4} \beta ^3 (33+2 \log (2))+\frac{1}{4} \beta (10 \log (2)-59)+\frac{6 \log (2)-6}{4 \beta }\, , \\[1pt] \tag{C.24}\end{equation}
\begin{equation} c_{42} = -\frac{1}{2} \beta ^5 \log (2)+\frac{1}{4} \beta ^3 (48+2 \log (2))+\frac{1}{4} \beta (10 \log (2)-91)+\frac{6 \log (2)-9}{4 \beta }\, , \\[1pt] \tag{C.25}\end{equation}
\begin{equation} c_{43} = -\frac{1}{2} \beta ^4 (2 \log (2)-19)+\frac{1}{4} \beta ^2 (8 \log (2)-101)-3+3 \log (2)\, , \\[1pt]\tag{C.26}\end{equation}
\begin{equation} c_{44} = -\frac{1}{4} \beta ^4 (31+12 \log (2))+\frac{7}{2} \beta ^2 (1+4 \log (2))-\frac{11}{4} (5+4 \log (2))\, , \\[1pt]\tag{C.27}\end{equation}
\begin{align} c_{45} &= -2 \pi ^2 \left(\pi ^2-3\right) \beta ^5-12 \pi ^2 \beta ^4 (19+2 \log (2))+6 \pi ^2 \beta ^2 (101+8 \log (2))\nonumber \\[1pt] & \quad\, -\frac{288 \beta \log (2)}{\beta ^2+3} +\frac{384 \beta \log (2)}{\left(\beta ^2+3\right)^2}+\frac{6 \pi ^4}{\beta }+\frac{1}{4} \beta ^3 \left(-1764 \zeta (3)-615+628 \pi ^2+8 \pi ^4\right.\nonumber \\[1pt] & \quad\,-\left. 192 \log ^2(2)-108 \log (2)-840 \pi ^2 \log (2)\right)+\frac{1}{4} \beta \left(1764 \zeta (3)+2757-1872 \pi ^2\right. \nonumber\\[1pt] & \quad\, +\left.40 \pi ^4+576 \log ^2(2) -108 \log (2)+1416 \pi ^2 \log (2)\right)+72 \pi ^2 (1+\log (2))\, , \tag{C.28}\end{align}
\begin{align} c_{46} &= -24 \pi ^2 \beta ^5 \log (2)+12 \pi ^2 \beta ^3 (2 \log (2)-33)-\frac{96 (33 \log (2)-4)}{\beta ^2+3}+\frac{384 (27 \log (2)-1)}{\left(\beta ^2+3\right)^2}\nonumber \\[1pt] & \quad\, -\frac{9216 \log (2)}{\left(\beta ^2+3\right)^3}+12 \pi ^2 \beta (59+10 \log (2))+\frac{72 \pi ^2 (1+\log (2))}{\beta }+\frac{1}{4} \beta ^4 \left(756 \zeta (3)+597\right. \nonumber\\[1pt] & \quad\, +\left.104 \pi ^2-384 \log ^2(2)-492 \log (2)+168 \pi ^2 \log (2)\right)-\frac{1}{2} \beta ^2 \left(1692 \zeta (3)+621+152 \pi ^2\right. \nonumber\\[1pt] & \quad\, -\left.384 \log ^2(2)+36 \log (2)+312 \pi ^2 \log (2)\right)+\frac{1}{4} \left(2052 \zeta (3)-243-1144 \pi ^2\right.\nonumber\\[1pt] & \quad\, + \left. 1152 \log ^2(2)+2964 \log (2)+840 \pi ^2 \log (2)\right) .\tag{C.29}\end{align}

Appendix D. Moments of photon correlation function

The moments of photon correlation function are given by

\begin{equation} \Pi (\bar{s}) = \sum_{n>0} M_n \left(\frac{\bar{s}}{4}\right)^n\, ,\tag{D.1}\end{equation}

where

\begin{equation} M_n = M^{(1)}_n \left(\frac{\alpha}{4\pi}\right) + M^{(2)}_n \left(\frac{\alpha}{4\pi}\right)^2 + M^{(3)}_n \left(\frac{\alpha}{4\pi}\right)^3 + \ldots\tag{D.2}\end{equation}

and the first five moments at one, two and three loops are given by

\begin{align} M^{(1)}_1 = \frac{16 N}{15}, \quad M^{(2)}_1 = \frac{1312 N}{81}, \tag{D.3}\end{align}
\begin{align} M^{(1)}_2 = \frac{16 N}{35}, \quad M^{(2)}_2 = \frac{7184 N}{675}, \tag{D.4}\end{align}
\begin{align} M^{(1)}_3 = \frac{256 N}{945}, \quad M^{(2)}_3 = \frac{3998656 N}{496125}, \tag{D.5}\end{align}
\begin{align} M^{(1)}_4 = \frac{128 N}{693}, \quad M^{(2)}_4 = \frac{831776 N}{127575}, \tag{D.6}\end{align}
\begin{align} M^{(1)}_5 = \frac{2048 N}{15015}, \quad M^{(2)}_5 = \frac{6918163456 N}{1260653625},\tag{D.7}\end{align}
\begin{equation} M^{(3)}_1 = N^2 \left(\frac{406 \zeta_3}{27}-\frac{45628}{729}+\frac{256 \pi ^2}{45}\right)+N \left(\frac{22781 \zeta_3}{108}-\frac{8687}{54}+\frac{32 \pi ^2}{3}-\frac{256}{15} \pi ^2 \log (2)\right) , \\[1pt]\tag{D.8}\end{equation}
\begin{align} M^{(3)}_2 &= N^2 \left(\frac{14203 \zeta_3}{1152}-\frac{1520789}{25920}+\frac{512 \pi ^2}{105}\right)\nonumber \\[1pt] &\quad\, +N \left(\frac{4857587 \zeta_3}{2880}-\frac{223404289}{116640}+\frac{64 \pi ^2}{7}-\frac{512}{35} \pi ^2 \log (2)\right) , \tag{D.9}\end{align}
\begin{align} M^{(3)}_3 &= N^2 \left(\frac{12355 \zeta_3}{864}-\frac{83936527}{1458000}+\frac{4096 \pi ^2}{945}\right)\nonumber \\[1pt] & \quad\, +N \left(\frac{33067024499 \zeta_3}{3225600}-\frac{885937890461}{72576000}+\frac{512 \pi ^2}{63}-\frac{4096}{315} \pi ^2 \log (2)\right) , \tag{D.10}\end{align}
\begin{align} M^{(3)}_4 &= N^2 \left(\frac{2522821 \zeta_3}{147456}-\frac{129586264289}{2239488000}+\frac{8192 \pi^2}{2079}\right)\nonumber \\[1pt] & \quad\, +N \left(\frac{1507351507033 \zeta_3}{25804800}-\frac{269240669884818833}{3840721920000}+\frac{5120 \pi^2}{693}-\frac{8192}{693} \pi ^2 \log (2)\right) , \tag{D.11}\end{align}
\begin{align} M^{(3)}_5 &= N^2 \left(\frac{1239683 \zeta_3}{61440}-\frac{512847330943}{8692992000}+\frac{32768 \pi^2}{9009}\right)\nonumber \\[1pt] &\quad\, +N \left(\frac{939939943788973 \zeta_3}{2980454400}-\frac{360248170450504167133}{950578675200000}+\frac{20480 \pi^2}{3003}-\frac{32768 \pi ^2 \log (2)}{3003}\right) .\tag{D.12}\end{align}

The additional moments up to \(n\) = 100 can be found in accompanying Mathematica notebook.

V. B. Berestetskii, E. M. Lifshitz, L. P. Pitaevskii, Quantum Electrodynamics: Volume 4, Vol. 4, Butterworth-Heinemann, 1982.
G. Källén, A. Sabry, Fourth order vacuum polarization, Matematisk-Fysiske Meddelelser Kongelige Danske Videnskabernes Selskab 29 (17) (1955) 556.
P. A. Baikov, D. J. Broadhurst, Three loop QED vacuum polarization and the four loop muon anomalous magnetic moment, in: Proceedings of the 4th International Workshop on Software Engineering and Artificial Intelligence for High-Energy and Nuclear Physics, 1995, pp. 0167–172. arXiv:hep-ph/9504398
T. Kinoshita, M. Nio, Sixth-order vacuum-polarization contribution to the lamb shift of muonic hydrogen, Physical Review Letters 82 (16) (1999) 3240.
T. Kinoshita, M. Nio, Accuracy of calculations involving \(α^3\) vacuum-polarization diagrams: Muonic hydrogen lamb shift and muon g – 2, Physical Review D 60 (5) (1999) 053008.
T. Kinoshita, W. Lindquist, Parametric formula for the sixth-order vacuum polarization contri- bution in quantum electrodynamics, Physical Review D 27 (4) (1983) 853.
T. Hahn, Generating Feynman diagrams and amplitudes with FeynArts 3, Computer Physics Communications 140 (2001) 418–431. doi:10.1016/S0010-4655(01) 00290-9 arXiv:hep-ph/0012260
K. G. Chetyrkin, F. V. Tkachov, Integration by parts: The algorithm to calculate β-functions in 4 loops, Nuclear Physics B 192 (1981) 159.
F. V. Tkachov, A theorem on analytical calculability of 4-loop renormalization group functions, Physics Letters B 100 (1) (1981) 65–68, available at http://www.sciencedirect.com/science/article/B6TVN-46YSNNV-109/1/886c25adc81acf1d171b80d7b1e7cb7f
R. N. Lee, A. I. Onishchenko, Master integrals for bipartite cuts of three-loop photon self energy, Journal of High Energy Physics 04 (2021) 177. doi:10.1007/JHEP04(2021)177 arXiv:2012.04230
R. N. Lee, A. I. Onishchenko, \(ε-regular\) basis for non-polylogarithmic multiloop integrals and total cross section of the process \(e^+e^− \to 2(Q\bar{Q})\), Journal of High Energy Physics 12 (2019) 084. doi:10.1007/JHEP12(2019)084 arXiv:1909.07710
H. Ferguson, D. Bailey, A polynomial time, numerically stable integer relation algorithm, RNR Technical Report RNR-91-032 (1992).
N. Gray, D. J. Broadhurst, W. Grafe, K. Schilcher, Three-loop relation of quark \(\bar{\text{MS}}\) and pole masses, Zeitschrift f ̈ur Physik C 48 (1990) 673–679.doi:10.1007/BF01614703
D. J. Broadhurst, N. Gray, K. Schilcher, Gauge-invariant on-shell \(Z_2\) in QED, QCD and the effective field theory of a static quark, Zeitschrift f ̈ur Physik C 52 (1991) 111–122. doi:10.1007/BF01412333
K. Melnikov, T. van Ritbergen, The three loop on-shell renormalization of QCD and QED, Nuclear Physics B 591 (2000) 515–546. doi:10.1016/S0550-3213(00)00526-5 arXiv:hep-ph/0005131
A. I. Onishchenko
Corresponding author e-mail address: onish@theor.jinr.ru