Numerical calculation of the full two-loop electroweak corrections to muon (g-2)

Numerical calculation of two-loop electroweak corrections to the muon anomalous magnetic moment ($g$-2) is done based on, on shell renormalization scheme (OS) and free quark model (FQM). The GRACE-FORM system is used to generate Feynman diagrams and corresponding amplitudes. Total 1780 two-loop diagrams and 70 one-loop diagrams composed of counter terms are calculated to get the renormalized quantity. As for the numerical calculation, we adopt trapezoidal rule with Double Exponential method (DE). Linear extrapolation method (LE) is introduced to regularize UV- and IR-divergences and to get finite values. The reliability of our result is guaranteed by several conditions. The sum of one and two loop electroweak corrections in this renormalization scheme becomes $a_\mu^{EW:OS}[1{\rm+}2{\rm -loop}]= 151.2 (\pm 1.0)\times 10^{-11}$, where the error is due to the numerical integration and the uncertainty of input mass parameters and of the hadronic corrections to electroweak loops. By taking the hadronic corrections into account, we get $a_\mu^{EW}[1{\rm+}2 {\rm -loop}]= 152.9 (\pm 1.0)\times 10^{-11}$. It is in agreement with the previous works given in PDG within errors.


I. INTRODUCTION
In order to get a sign of beyond the standard model physics from high precision experimental data, we need higher order radiative corrections within Standard Model (SM). For this purpose our group has been developing the automatic calculation system GRACE [2] since the late 1980's. The measurement of the muon anomalous magnetic moment a µ ≡ (g-2)/2 is the one of the most precise experiments to check the SM. QED correction was calculated by T.Kinoshita et al. [3] up to tenthorder. The two-loop electroweak (ELWK) correction to a µ was calculated approximately by Kukhto et al. [4] in 1992. Surprisingly, the two-loop correction is almost 20% of the one-loop correction [5][6][7][8]. We started to calculate the full two-loop corrections in 1995 and presented our formalism at Pisa conference [9]. We also showed that the two-loop QED value [10][11][12][13] was correctly reproduced within our general formalism. However, the number of diagrams is huge and the numerical integration requires the big CPU-power to achieve required accuracy, we must wait until various environments are improved.
During these days, the several groups did the approximate calculations [14][15][16][17][18][19] and the approximate value of the two-loop ELWK correction is widely accepted [1,20]. In 2001, BNL-Experiment 821 [21,22] announced that the precise experimental value deviates from that of SM around (2.2 ∼ 2.7)σ [23]. It brought much interest in the theoretical value. The main theoretical concern is now Although the two-loop ELWK correction is almost established, we try to get the value without any approximation to confirm the validity of the earlier studies 1 . This work is an important milestone to extend GRACEsystem from one-loop to two-loop calculation. In ELWK theory, there are so many fields, mass parameters and complex couplings that it is hard to get reliable higher order corrections to physical quantities, in general. It is desirable to construct the framework to calculate these higher order corrections as automatically as possible.
The key point is to perform Feynman integration numerically by using a sophisticated method with good convergence and high power CPU machine. We propose to call such framework as Perturbative Numerical Quantum Field Theory (PNQFT). The concepts of PNQFT are based on the following principles.
(a): It is essential to assume amplitudes as meromorphic functions of space time dimension n for regulariza-tion and getting gauge invariant renormalized values of physical quantities.
(b): The source program for numerical integration is automatically generated by GRACE together with a symbolic manipulation system such as FORM [30].
(c): A high precision numerical integration method should be adopted.
(d): Linear Expansion method (LE) (see subsection III B) [31,32] is crucial to extract both UV-and IR-divergences by taking advantage of the above analyticity. By LE method, we can expand the amplitude in any order of ε(= 2 − n/2) (Laurent expansion), so that it is a powerful tool for higher order calculation.
(e): To guarantee the validity of the calculation, several conditions must be cleared. An example is the cancellation of non-linear gauge (NLG) parameters.
(f ): It is crucial to reduce human intervention to avoid careless mistakes. We must minimize the handmade operations necessary for getting the physical quantities.
The following calculation is based on these principles. In section II and III, we briefly explain the flow and framework of our calculation. In section IV, we touch on our method of numerical calculation. We emphasize that the Linear Extrapolation (LE) method is simple and efficient method to regularize UV-and IR-divergences and also to get finite values. We also explain our consistency conditions to ensure the results. Some examples of calculations are explained. In section V, we give our results on a µ . In the last section, we give some comments to make extensive progress. In Appendices, we explain the technical parts of our calculation.

II. OUTLINE OF OUR FRAME WORK
Our calculation is formulated under the following conditions.
1. The calculation is done within SM.

On mass shell renormalization scheme (OS) is
adopted [33][34][35]. We adopt α, M Z , M W , M H and fermion masses as physical parameters. Weinberg angle, Higgs fermion coupling and other quantities are expressed by these parameters.
3. Free quark model (FQM) is adopted and as for quarks, constituent masses are used.
5. Dimensional regularization is applied for both Ultra Violet (UV)-and Infrared (IR)-divergences. 6. Linear Extrapolation method (LE) is fully used for regularization and getting finite values.
Next, we briefly explain the flow of our calculation. 1. GRACE system generates all the diagrams we need in SM, automatically [36]. There are 1780 two-loop diagrams 2 and 70 one-loop diagrams composed of one-loop order counter term (CT). Two-loop order CT is not necessary in our case, because a µ is not related to the charge renormalization part.
2. These 1780 diagrams are classified into 14 types of topology. Types of the topology are displayed in Fig.1. Among these types, some of them give the same contribution because of symmetry. (an example: 5-a vs. 5-b) The diagrams including CT are classified into two types, namely, vertex and self-energy types.
FIG. 1. Types of topology 3. GRACE system generates the amplitude of each diagram in accordance with Feynman rules for ELWK theory with NLG [35]. 4. In order to express the amplitude as a function of Feynman parameters used to combine denominators, we define the following quantities for each topology in advance [37,38].
5. Using these tools, contribution of each diagram to a µ is expressed as function of Feynman parameters, according to the formulas given in the next section. We make use of a symbolic manipulation system FORM exhaustively.

III. LOGIC OF THE CALCULATION
A. Cvitanović-Kinoshita procedure In order to extract a µ factor from muon vertex function, we adopt Cvitanović-Kinoshita procedure [37,38]. We briefly explain the procedure in the case where there are six propagators. Starting formula is the two-loop muon vertex, where p j = 2 s=1 η s (j) s + q j , is the momentum on the internal line (j). The function η s (j)(= ±1, 0) defines the weight of loop momentum s on the internal line (j). The z j 's are the Feynman parameters to combine six propagators. F µ (D) is the numerator function and µ is the external photon polarization. Next we diagonalize the denominator function with respect to loop momenta 1 , 2 and perform integration. The result is, Where U is well known 2×2 matrix, composed of Feynman parameters (z j ). V (z j , m j , q j ) is the denominator function. The argument is easily extended to the case with five-propagators (diagrams with four-point coupling).
In order to generate the numerator function we use the following differential integral operator D µ j .
The operator D µ j generates momentum p µ j on the internal line (j). By operating D µ j to the denominator function V , we get the following expression.
Using the above formulas, we can write down the numerator functions in terms of B ij algebraically. The equivalence of the above method and the well known method of shifting loop momentum to diagonalize the denominator function is verified. The correspondence between two methods are symbolized as follows.
Next step is to extract the a µ factor by using projection operator, from the photon muon vertex Γ µ . The quantity a µ is given as follows. ( m 0 =muon mass ) where we set momentum of incoming µ − , outgoing µ − and incoming photon, as (p − q/2), (p + q/2) and q, respectively. The final expression for numerical integration is summarized in the following formula.
The numerators f 0 , f 2 represent the coefficient of 0 , 2 term, respectively, after projection operator is applied. They are also the function of dimension n.

B. Regularization Method
Next step is the regularization of UV-and IRdivergences. By adopting n-dimensional regularization method, any integrand F of Feynman parameter integration is regarded as a function of ε = 2 − n/2, F (ε). We adopt two methods for regularization.

Linear Expansion method
First one is very simple and powerful if the accuracy of numerical integration is sufficiently guaranteed. We call it Linear Expansion method (LE). Just after the dimensional regularization method was introduced [39,40], the analyticity with respect to ε was discussed extensively. It is shown that the Feynman amplitude is a meromorphic function of ε [39,40,41]. This is a key point to utilize the LE-method to the Feynman amplitude. The followings are the steps to get the divergent and finite terms.
3. According to the analyticity, we can expand G(ε(i)) in Laurent series. In our case, it is evident that the expansion starts from (1/ε(i)) because of the lack of two-loop counter terms. We truncate the series The coefficients C −1 and C 0 correspond to the divergent and finite parts, respectively. In the case ε = 2 − n/2, C −1 expresses the UV-divergent part and if we set ε → ε R = (n/2 − 4), C −1 represents the IR-divergent part.
5. In order to improve the convergence, we set M ∼ 18 and α ∼ 1.1 by trial and error. Examples setting these parameters are shown in [31,32] 6. Various methods are known to extract C −1 , C 0 from G(ε(i)) [42], however, LE method is simple and appropriate in our case.
In order to get the reliable value of C −1 , C 0 up to 4 digits, we need the accuracy of the numerical integration at least 8 digits.

Subtraction Method
To complement the above calculation, we also adopt the well known subtraction method to separate divergent part and finite part. We extract 1/ε singularity from G(ε) when one of Feynman parameters approaches 0, (x → 0). The followings are the steps to extract the singularity.
3. In the case where there is self-energy type diagram, the factor x ε−2 appears in the head of integrant. If we expand it in ε by using analytic continuation the following formula is obtained.
We use this method partly to complement the LE method.

Counter terms
As for counter terms, GRACE has a library of renormalization constants at one-loop level based on OSrenormalization scheme. We make use of this library for 70-diagrams composed of counter terms. Generally speaking, it is necessary to expand one-loop renormalization constants up to order ε = 2 − n/2 at two-loop level. However, the divergent part of diagrams composed of CT does not contribute to a µ , the O(ε) term is unnecessary in our case. Here we comment on the wave function renormalization constant of goldstone fields χ, χ 3 . We keep the finite part of the constant in the form (−1/2){dΠ χχ (q 2 )/dq 2 } at q 2 = M 2 W . However, the final answer is independent of the finite part. For the renormalization of unphysical fields, the UV-divergent part is only useful to erase divergence.

IV. NUMERICAL CALCULATION
A. Double exponential method The final step to get the value a µ is the numerical integration over Feynman parameters. We employ trapezoidal rule with Double Exponential (DE) transformation method [43]. It is also called as tanh-sinh transformation method. It is very powerful if the integrand has singular behavior at the edge of the integration domain. Speed of convergence is accelerated by the DE transformation, The maximum dimension of multiple integration is five. We apply DE-method to any integration variable involved. As we need the accuracy greater than 8 digits to see the UV cancellation, the adaptive Monte Carlo method is not adopted in our two-loop calculation.

B. Criterion to ensure the validity of the result
In order to ensure the validity of our results, we impose several conditions given below. Originally non-linear gauge was introduced to reduce the number of diagrams, particularly containing bosonboson couplings [44][45][46][47][48][49]. Here, we adopt NLG to check the validity of our calculation. The gauge fixing Lagrangian is constructed as, Here,α,β,δ,ε andκ are non-linear gauge parameters specific to this formalism. The parameters s W and c W are the sine and cosine of Weinberg angle θ W . In our calculation we set ξ = ξ W = ξ Z = 1 to make the gauge boson propagators simple. NLG parameters are distributed among so many diagrams of different types of topologies.
So it is very powerful if we can verify the cancellation of these NLG parameters. We show the sample of cancellation in Appendix F.

Successive method
Diagrams with self energy type two-point function can be calculated by successive method using renormalized two point function. An example is diagram with (γ − γ) or (γ − Z) vacuum polarization type diagrams. We decompose the renormalization constants δZ 1/2 AA , δZ 1/2 ZA , δZ 1/2 AZ , δM 2 Z etc. into components according to the particles involved in the loop. By adding the counter term to corresponding one-loop unrenormalized two-point function, one-loop ( 1 ) integration is performed without divergence and we obtain the renormalized twopoint function Π R . By inserting Π R into the second loop( 2 ), we get finite value of a µ . We use this alternative method to reconfirm the results obtained by the methods given in section III. An example is shown in Appendix G.

V. RESULTS OF OUR CALCULATION
As the physical input parameters, we use the following fermion and boson masses (unit GeV).
As we mentioned before, we adopt OS renormalization, however, the expression in the preceding works is parametrized using Fermi constant G F = 1.1663787 × 10 −5 GeV −2 and α [18].
The difference of the 2-loop correction between our value and that in G F parametrization is due to the fact that one loop correction in G F parametrization partially includes the α 2 correction in our scheme.
So the comparison should be done to the sum of oneand two-loop. The one loop correction in our OS scheme is written down as follows.
We can carry out the numerical calculation without any approximation and get the value.
Here we add the error due to neglecting the uncertainty in electroweak loops involving hadrons.
When we compare our result with the value obtained by using G F parametrization, we need the naive free light quark model calculation with the same quark masses as ours. This is given 3 in ref. [16]. In this case the two loop correction becomes, −42.97(±1) × 10 −11 . In the G F parametrization, the one loop correction becomes 194.80(±0.01) × 10 −11 [18], so that we get, This is consistent with our value Eq. (22) We also add a comment on the relation between the well known PDG value [1] shown below and our value. If we include the hadronic correction to light quark contribution, by adding the difference of the following expression [18], a EW (2) µ (e, µ, u, c, d, s) = (−6.91 ± 0.20 ± 0.3) × 10 −11 (24) and the value quoted in the footnote below [16], our value becomes as follows.

VI. DISCUSSIONS AND COMMENTS
We developed the system to calculate the full ELWK two-loop corrections to a µ , by fully using GRACE and FORM on the basis of OS-scheme. The work we need beforehand is only to prepare several files which only depend on the type of the topology of diagrams as we explained in section III. We adopt the dimensional regularization to regularize UV-and IR-divergences and to get finite gauge invariant values of the physical quantity. To extract the (1/ε) terms, we use Linear Expansion Method explained in sectionIII-B. This method is very simple and attractive, compared with the conventional method to take out the (1/ε) term by extrapolating one of the Feynman parameters close to 0 . If we adopt the conventional method, it is crucial to introduce the most suitable transformations from (z j ) to [0,1] integration variables (x, y, u, v, w). Furthermore, we need rather complex operations including differentiation of the amplitude, etc. As a result, the necessary CPU-time increases extensively.
In the case of Linear Expansion method (LE), however, the choice of integration variables is not sensitive to get the reliable results and this method decreases the number of operation drastically. It is sufficient to define the quantity as function of ε(= 2 − n/2) . We only need to treat Dirac matrices and various vectors appeared in 3 In ref. [16], the contribution of light quarks in FQM is a EW (2) µ (e, µ, u, c, d, s) = −(4.0 + 4.65) × 10 −11 = −8.65 × 10 −11 the numerator, in n-dimension. This is easily done by using symbolic manipulation system such as FORM. The operation is simple and we can make use of the resultant short sources for both UV-(ε > 0) and IR-(ε R = −ε > 0) regularization and also to get finite results. We conclude that LE-method is the most simple and reliable method, at this moment. In order to get reliable physical value by this method, high precision numerical integration over Feynman parameters is inevitable. The DE-method introduced in section IV A is the suitable candidate.
Introduction of NLG-parameters makes the calculation very complex, however, it is very powerful to check the calculation of so called Boson contribution. The number of diagrams consisting of different types of topology are connected through NLG-parameters. The maximum number of diagrams mutually entangled reaches 864. So this is a very tough condition to be cleared.
By making use of these technical approaches mentioned above, we clear all the constraints given in section IV B . Namely, (i) reproduction of QED values, (ii)(iii) cancellation of UV-and IR-divergences, (iv) Independence of NLG-gauge parameters. We show some samples in Appendices how they are cleared.
The final value of the sum of one and two loop weak corrections is approximately the same as the one obtained by previous works using different parametrization.
Based on this work we can proceed to construct PN-QFT(Perturbative Numerical Quantum Field Theory), which we discussed in section I A. Wide range of application to ELWK higher loop expansion for several physical reactions will be opened. We expect that this work provides the fruitful foundation to formulate PNQFT.

ACKNOWLEDGMENTS
We would like to thank Prof.T.Kaneko for his important contribution to construct the framework of calculation at the early stage of this work. We also wish to thank Prof.K.Kato, Prof.F.Yuasa and Prof.M.Kuroda for discussions. Last but not least, we express our deep appreciation to late Prof.Y.Shimizu for his continual encouragement. This research is partially supported by Grant-in-Aid for Scientific Research (15H03668,15H03602) of JSPS and Grand-in-Aid for High Performance Computing with General Purpose Computers (Research and development in the next-generation area) of MEXT. In Fig.1, we show types of topology to formulate the two-loop contributions. However, in order to check the consistency of numerical values, it is useful to classify diagrams by distinguishing fermion and boson lines in each diagrams in Fig.1. We briefly figure out the classification method in Fig.2.

FIG. 2. Classification of Diagrams
In the figure, the straight lines and wavy lines represent fermion and boson, respectively. The circle indicates one-loop diagram. We classify diagrams depending on the place where one-loop diagram is inserted. It is summarized in the following Table I. You can easily see which one of Fig.1 is classified into which category.   Table I we add two types of topology having no divergence, namely, Fig.1-(1) and (6).

Appendix C: UV-cancellation
In this Appendix, we show the cancellation of UVdivergence in linear gauge ('t Hooft-Feynman gauge). Examples of a group of diagrams are shown in Fig.3. The diagrams in Fig.3 belong to several groups in Table I.
In this case, total 13-diagrams and 1-counter term (e) make a group to cancel UV-divergence. In the Table II, the coefficient of C (2) U V = (1/ε−2γ+2 ln(4π)), corresponding to each diagram is shown. As you can read from the Table II, the cancellation is marvelous, up to almost 15 digits. U V of a set of diagrams. Here we show the set given in Fig.3 as a sample. The 15-digit cancellation is realized in the sum.

Diagram Particles on
Value (unit 10 −11 ) in Fig.3 line Appendix D: IR-cancellation LE method is applied to check IR-cancellation. Among two-loop diagrams, the 8 diagrams in Fig.(4) have IRdivergence.
The diagrams with CT also have IRdivergence through δZ 1/2 W (19 diagrams) and δZ 1/2 µ (28diagrams). The IR-divergence at two-loop level is proportional to C (2) IR = (−1/ε R − 2γ + 2 ln(4π)), ε R = (n/2 − 2) > 0. It is easily shown that the IR-divergence coming from δZ 1/2 W cancels among the 19 CT-diagrams. As for the diagrams with δZ 1/2 µ , IR-divergence is cancelled by the corresponding two-loop diagrams. We show coefficients of C (2) IR in Table III. The correspondence between small photon mass (λ) method and LE-method for IR-regularization is checked in the case of QED ladder diagram. Analytic value of a coefficient of ln(λ 2 /m 2 µ ) in unit of (α/π) 2 is (1/4) [10]. It is 0.249999998 by our calculation using small photon mass. In LE-method, the coefficient of C (2) IR becomes 0.249999999. We understand the correspondence between ln(λ 2 /m 2 µ ) and C (2) IR is established. As we show in Table III, no IR-divergence remains in the final expression. W are taken into account. The gray circles and white circles correspond to Self CT and Vertex CT in Table III , respectively.

Appendix E: Sample calculation of the two-loop diagram
As a sample, we show the calculation of the two-loop diagram which contains both UV-and IR-divergences. The diagram is shown in Fig.4-(b). The line numbers are given in Fig.6. Following Eq.(9) given in section III-A, the essential part of expression of two-loop diagram contribution is written as the following form. Quantities in Eq.(E1),(E2) are expressed by Feynman parameters z j , (j = 1, · · · , 6), ( j z j = 1). Masses are made dimensionless using muon mass.
The f 0 and f 2 part give IR-and UV-divergences, respectively. By adopting DE and LE methods for Feynman parameter integration, we get the following expansion.
We multiply the Γ functions arising from n-dimensional integration and factor 1/(4π) n , to the above quantities, We must pay attention that the term (2C −1 ) appears, because Γ(6 − n) contains the term (1 − 2ε R ). In order to get the correct finite value, we must calculate the counter terms using the same regularization method for both IRand UV-divergences as in two-loop case. In this case, we need muon wave function renormalization constant δZ µ . The δZ µ is obtained by calculating muon self energy diagrams. The δZ µ (γ) represents photon exchange part and δZ µ (Z) represents Z-exchange part. In Fig.6-(b), CT Z represents δZ µ (Z), and in Fig.6-(c), CT 1 ∼ CT 5 represent only IR-divergent part of δZ µ (γ). These IR-divergent part cancels that of Fig.6-(a) as is shown in the Table III. We will show the numerical results to see the situation clearly. (unit=10 −11 ) The term CT Z has both UV and finite parts.
UV-cancellation is excellent and cancels 17 digits.
To cancel IR part of Fig.6-(a) , we adopt IR-part of δZ µ (γ), CT 1 ∼ CT 5 . It is explained in the Table III, row (Z boson). The finite value is obtained as follows.
We can also regularize the IR divergent part by employing the small photon mass λ 2 , to see log(λ 2 ) term.
As we notice that we must use the same regularization method for both two-loop and CT diagrams. At first sight, the finite part changes compared with dimensional regularization method, however, the sum of two-loop and CT contribution is the same in both regularization methods. We confirm it by numerical calculation. As we can see from Eq.(E9), we need very careful treatment to discuss order of (a few) ×10 −11 quantity of a µ .
Appendix F: NLG parameter cancellation As an example of cancellation of NLG parameters, we show the UV-divergent part and finite part of diagrams which containα n (n = 1, 2, 3) terms. Total number of diagrams depending onα amount to 864. In order to see how the NLG-parameters are cancelled, we classify diagrams according to Fig.2 and Table I. We summarize the result of numerical calculation in Table IV. In the last row in Table IV, we show the maximum absolute value among the individual terms in the column and its type of diagram, to indicate the degree of the cancellation. First we see the cancellation of the parameterα n in the UV-divergent part, namely, the coefficients of C  Table IV. We can see the cancellation works very well and it strongly guarantees the validity of our numerical calculation. In the right half of Table IV, we also show the finite contribution to a µ . The columnα 1 ∼α 3 show that how the NLGcancellation works well also in finite part.

Appendix G: Successive method
In some cases, in order to check our calculation, we perform two-loop integrations successively. By integrating the first loop ( 1 ), we make the renormalized effective function and insert it to the second loop( 2 ) integration. As an example, we calculate the diagrams consisting of photon-Z meson mixing vacuum polarization. There are two types of topology 4 and 10 in Fig.1. In general, unrenormalized two point γ − Z function is written as follows.
Each one-loop diagram contributes to Π un T (q 2 ), Π un L (q 2 ) in the following form.
The variables A 1 ∼ A 4 , B, D Q , m depend on the specific one-loop diagram. The A 4 -term comes from the 4-point boson coupling diagram. However, this term drops out by the renormalization process.
In the followings, the renormalized quantity Π are expressed asΠ.
We notice that the renormalization conditions are fixed by Π un T (q 2 ), there is no freedom to renormalize Π un L (q 2 ). By using these renormalization conditions, we can write down renormalized a R (q 2 ). It has no divergence, we can perform loop( 1 ) integration. The result is given as Eq.(G5),(G6). where M represents the mass of a particle circulating the loop. As forb(q 2 ), we must check whether it is really finite or not. The C U V -part ofb(q 2 ) disappears after summing all the one-loop diagrams and integration of Feynman parameter x. We also notice that theb(q 2 )-part does not contribute to a µ , because it is proportional to q µ q ν . Next step is to insert the renormalized two point functionâ(q 2 )g µν into the triangle diagram and integrate over second loop momentum 2 . The integration has logarithmic divergent part, however, it drops out by the projection operator to a µ , Eq. In unit of 10 −11 we get, −6.1048×10 −3 by this method. On the other hand, we get −6.104 × 10 −3 by our twoloop formalism with CT-terms. The coincidence is quite good. Notice that, in this case the two-loop formalism takes huge cpu-time, especially for (W −W ) diagram. Its contribution is around 8.8×10 6 in unit of 10 −11 so that we need almost 15 digits number to cancel UV-divergence. So, the effective method is not only important to check the reliability of our general formalism but also is useful to get the numerical result. In the case of self-energy type diagrams, we can construct effective method in several cases, however, for vertex type diagrams, to construct effective method is complex.