Numerical calculation of 5-loop QED contributions to the electron anomalous magnetic moment

Preliminary numerical results of calculating the 5-loop QED contributions to the electron (and muon) anomalous magnetic moment are presented. The results include the total contribution of the 5-loop Feynman diagrams without lepton loops and the contributions of nine gauge-invariant classes that form that set. A discrepancy with known results is revealed. The contributions of the gauge-invariant classes are presented for the first time. The method of the computation is briefly described (with the corresponding references). The calculation is based on the following elements: 1) a subtraction procedure for removing infrared (IR) and ultraviolet (UV) divergences point by point in Feynman parametric representation before integration; 2) a nonadaptive Monte Carlo integration method that is founded on probability density functions (PDF) that are constructed for each Feynman diagram individually using its combinatorial structure; 3) a GPU-based numerical integration with the help of the supercomputer"Govorun"(JINR, Dubna, Russia).


Introduction
The electron anomalous magnetic moment a e is known with a very high precision. The most accurate measured value [1] is a e = 0.00115965218073(28).
The most precise theoretical prediction [2] was obtained using the following representation: (m e /m µ , m e /m τ ).
The contributions of different parts of this expression were obtained by different researchers. The corresponding value a e = 0.001159652181606(229) (11) (12) based on a relatively independent from a e measurement of α giving the value α −1 = 137.035999046(27); see [3] 1 . The first uncertainty came from the inaccuracy of α; the second one came from the numerical uncertainty of A (10) 1 due to the statistical error of the Monte Carlo integration; and the last one came from a e (hadronic) + a e (electroweak). The first uncertainty is far bigger than the last two ones. However, we should take into account that: • These calculations are used for improving α; see [2].
• Until recently, A (10) 1 had been calculated only by one researcher team [2]; an independent calculation is required for reliability.
Thus, the problem of calculating A (2n) 1 is still relevant. We note that the values for A (2n) 1 , n = 1, 2, 3, 4 are reliable: 1 = 0.5 was first calculated by J. Schwinger in 1948 [5,6]. A recalculation is contained in many textbooks of quantum field theory.
(ii) The value A (4) 1 = −0.32847... was first obtained by A. Petermann [7] and independently by C. Sommerfield [8] in 1957. There are some another independent calculations for this value. (iii) A (6) 1 was being calculated numerically during 1960's and 1970's by different scientific groups [9,10,11]; the most precise value A (6) 1 = 1.195 ± 0.026 for that time was obtained by T. Kinoshita and P. Cvitanović in 1974. Simultaneously, a work of analytical calculation was in progress. That work was contributed by many different researchers. The final value A (6) 1 = 1.181241... was obtained in 1996 by E. Remiddi and S. Laporta [12]. See also independent calculations in [13,14].
(iv) The most precise value A  [15]. That value was confirmed and improved by S. Laporta in 2017 [16] with the help of a semianalytical approach: A (8) 1 = −1.9122457.... See also independent (but not such precise) calculations in [13,17] and a calculation for diagrams without lepton loops in [18].

The method
For reducing the computation time in high-order QED calculations and for making that calculations practically feasible it is very important to remove all IR and UV divergences in Feynman diagrams before integration. It is known that Bogoliubov's R-operation removes all UV divergences in all individual diagrams in Schwinger parametric space [19,20] and in momentum space with noncovariant regularization [21] as well as the Zimmermann forest formula 2 does. The electron anomalous magnetic moment is free from IR divergences, because the on-shell renormalization removes them as well as it removes the UV ones. However, individual diagrams remain IR divergent after applying the direct subtraction on the mass shell. The structure of IR divergences in Feynman parametric space is quite complicated: they can not be recognized by a direct power counting based on a set of Feynman parameters [24,14]. Moreover, sometimes IR and UV divergences can not be separated from each other; see notes in [25]. For removing both IR and UV divergences different authors developed different subtraction procedures [2,9,10,24]. We use another one that is based on linear operators that are applied to Feynman amplitudes of UV divergent subdiagrams. A detailed explanation of the procedure is contained in [14]. Let us recapitulate the advantages of the developed subtraction procedure: (i) It is fully automated for any order of the perturbation series. is the sum of these contributions.
(v) Feynman parameters can be used directly, without any additional tricks.
After applying the subtraction we have a finite integrand for each diagram. Each n-loop integrand has 3n − 2 variables 3 . However, the shape of these integrands remains unconvenient for integration. The 5-loop integrands have 13 variables; we evaluate these integrals by Monte Carlo. The convergence speed of Monte Carlo integration depends on the choice of the PDF. Adaptive algorithms like VEGAS can fit the PDF to the shape of the integrand. However, these algorithms can adjust only relatively small number of parameters; this can be inefficient for large dimensions. In contrast, for diagrams without lepton loops we propose to use a nonadaptive algorithm that is based on some theoretical considerations about Feynman parametric integrands behavior. The problem of constructing the PDF is very similar to the problem of constructing a good upper bound for the absolute value of the integrand and to the problem of proving the finiteness of renormalized Feynman amplitudes [25]. Thus, we use some ideas from that area.
Let f (z 1 , . . . , z M ) be the Feynman parametric integrand. We split all the integration area into the Hepp sectors [19]: The PDF is defined as where Deg(s) are positive numbers 4 that are defined on all subsets of the set of all internal lines of the diagram except the empty set and the full set. In the hypothetical case when the diagram does not contain UV divergent subdiagrams and we do not take the infrared limit, it can be proved using Speer's lemma [26] that (1) is an upper bound for |f (z 1 , . . . , z M )|, if Deg(s) = −ω(s), where ω(s) is the ultraviolet degree of divergence of the set s. When we take the infrared limit, Deg(s) should be decreased. For considering the infrared limit we apply ω to so-called I-closures of sets [25]. We define the I-closure of a set s as s∪s ′ , where s ′ is the set of all photon lines, for which the electron path connecting the ends of this line is contained in s. The complete definition of Deg is quite complicated and is not fully justified mathematically; two versions of this definition are presented in [25,18]. A comparison of the developed algorithm and known adaptive algorithms with respect to the Monte Carlo convergence speed is presented in [25].

Preliminary numerical results
Despite the good convergence, the Monte Carlo integration for the 5-loop case still requires a lot of computer time. For 5-loop QED diagrams without lepton loops this computation was performed using the GPUs NVidia Tesla V100 as a part of the supercomputer "Govorun". That computation led to the preliminary result that corresponds to 1σ limits. For eliminating round-off errors that occur due to numerical subtraction of divergences we used the interval arithmetic with different precisions and speeds; see the details in [25]. The code for all 3213 integrands with all arithmetics was generated in C++ with CUDA and required 500 GB of disk space in the compiled form. For reliability, the Monte Carlo integration was performed with two different pseudorandom generators from the NVidia CURAND library: (i) MRG32k3a (Calc 1; 19515 GPU-hours); (ii) Philox 4x32 10 (Calc 2; 6282 GPU-hours).
That calculations gave the results 6.739(132) and 6.905(220) respectively. The results were statistically combined in (2). In total, 1.9 · 10 14 Monte Carlo samples were generated and processed. The total calculation time amounted 25797 GPU-hours and the calculation is still in progress. Let us remark that there is a significant discrepancy between (2) and the result from [2]:   Different QED calculations of a e performed by different researchers showed that relatively small  Table 1 clearly demonstrate this fact. If a j is the contribution of a diagram j, then the contribution of a class is obtained as a j . For the classes (k, m, n) it is easy to see that |a j | and max |a j | are much greater in absolute value than the contribution. Here N diag means the quantity of diagrams in the class 5 .
The developed method also provides a possibility to check the results by parts: the whole set of diagrams can be split into 809 sets for comparison with the direct subtraction on the mass shell in Feynman gauge. The results for these sets will be published in the further papers; see analogous 4-loop results in [18] and a comparison of the 2-loop and 3-loop results with the known analytical results in [14,18].