The two-loop perturbative correction to the (g-2)$_{\mu}$ HLbL at short distances

The short-distance behaviour of the hadronic light-by-light (HLbL) contribution to $(g-2)_{\mu}$ has recently been studied by means of an operator product expansion in a background electromagnetic field. The leading term in this expansion has been shown to be given by the massless quark loop, and the non-perturbative corrections are numerically very suppressed. Here, we calculate the perturbative QCD correction to the massless quark loop. The correction is found to be fairly small compared to the quark loop as far as we study energy scales where the perturbative running for the QCD coupling is well-defined, i.e.~for scales $\mu\gtrsim 1\, \mathrm{GeV}$. This should allow to reduce the large systematic uncertainty associated to high-multiplicity hadronic states.


9)
PΠ 6 µ 1 µ 2 µ 3 µ 4 ν 4 = g µ 1 µ 3 g µ 2 ν 4 q 1,µ 4 , (2.10) (2.14) 23) which has been built in such a way that, combined with the crossing symmetries of the HLbL tensor, theΠ satisfy the following crossing symmetries (2.24) The operator C ij interchanges two momenta q i and q j . Notice how, from the knowledge of five of them, for exampleΠ 1,7,10,13,19 , one can easily infer the rest from these crossing symmetries. Thesẽ Π i are well-defined and are both ultraviolet and infrared finite to the order in α we are working. For our calculation we use two different sets ofΠ i , related by gauge invariance, which thus provides a cross-check of our results. An OPE is only valid for large Euclidean momenta [53]. As a consequence, it cannot be directly applied to the tensor in (2.1) for the (g − 2) µ kinematics, since by definition the external photon is soft, q 4 → 0, even though the other Euclidean momenta are large, −q 2 i ≡ Q 2 i Λ 2 QCD [54,57]. However, precisely the same fact allows one to connect the tensor in (2.1) to the OPE of the tensor operator with the background photon field The OPE in question holds for large photon virtualities Q 2 1 ∼ Q 2 2 ∼ Q 2 3 Λ 2 QCD . In this expansion, any local operator with the same quantum numbers as F µν , including F µν itself, can absorb the remaining soft static photon and, as a consequence, give a contribution [7,54,57].
Higher-dimensional operators are suppressed by extra powers of , providing a hierarchy of contributions with a systematic counting. A very detailed study of this OPE can be found in Ref. [57], where the different power corrections were computed and found to be small compared to the leading contribution 1 . The leading term comes from the F µν operator itself and is given by the massless quark loop at order α 0 s , and the leading mass effects are very small. In fact, this quark loop corresponds to the zero momentum limit of the derivative of the naive massless perturbative QCD tensor of (2.1), i.e. (2. 26) In this work, we compute the leading α s correction to the direct F µν contribution in the OPE of (2.25). This corresponds to a two-loop massless QCD calculation with three external legs offshell.
Before discussing the gluonic correction to the quark loop in the next section, we first write down the quark loop result for theΠ i basis. As will be remarked upon below, this basis was not used in Ref. [57]. The whole solution for the quark loop can be simply written as Here N c is the number of colours and C 123 (0) is a loop integral function that is defined in App. A. For the (g − 2) µ integration, it is convenient using the generic results of Refs. [24,28,57]. Following them, the HLbL tensor can be expanded in a basis of 54 scalar functions Π i weighted with Lorentz structures T µ 1 µ 2 µ 3 µ 4 (2.32) the 19Π i defined in (2.4) can be identified with the static q 4 → 0 limit of certain linear combinations of the Π i . Denoting these linear combinationsΠ i it can further be shown that for a HLbL µ only six Π i contribute, namelyΠ 1,4,7,17,39,54 . 2 In particular, the a HLbL µ may be written [24,28] a HLbL The integration variable τ is defined via Q 2 3 = Q 2 1 + Q 2 2 + 2τ Q 1 Q 2 , the T i (Q 1 , Q 2 , τ ) are functions and the Π i are functions of the sixΠ i . The latter set of functions is related to theΠ i through (2. 35) In summary, knowledge ofΠ 1,4,7,17,39,54 is enough to determine a HLbL µ from (2.34). TheΠ i can be obtained from the derivative of the HLbL tensor in the static limit with the projectors given in Ref. [57]. There we definedΠ with the projectors PΠ i µ 1 µ 2 µ 3 µ 4 ν 4 given in App. A of [57]. Using the definitions of theΠ in (2.4) the relation between theΠ and theΠ follows immediately. We have checked that this procedure reproduces the massless quark loop results as given in Ref. [57]. TheΠ i representation of the massless quark loop was given in (2.27)-(2.31), and it can be noted that it is much simpler than the expressions for theΠ i in Ref. [57].

The two-loop perturbative correction
In this section we present the calculation of the two-loop contribution. For the analytic calculation we use FORM [60]. The master integral reduction is done by means of Kira [61], which employs a Laporta algorithm to reach a minimal set of master integrals. Explicit analytic expressions of the master integrals can be found in the literature.
The gluonic corrections to the quark loop are obtained by including two quark-gluon vertices from the Dyson series expansion in (2.1), i.e.
Denoting colour indices with bars, the interaction Lagrangians above are of the form where B a i ν i is the gluon field, g S is the strong coupling and λ a i is an SU (3) c Gell-Mann matrix. The only nonzero topology at this order is obtained by connecting all the quarks to the same line. Two examples of the diagrams in question are shown in Fig. 2. As a consequence of the topology, both the quark electric charge e 4 q and the colour factor, Tr(λ a λ b )δ ab = 2(N 2 c − 1), can be factored out, allowing to re-express the total contribution as a sum of all possible hexagons where two of the external lines are to be contracted to form the second loop, i.e., where Here, S(p) = / p p 2 is the massless quark propagator and σ(1, 2, 4, 5, 6) the set of pairwise permutations of µ i and q i for i = 1, 2, 3, 5, 6. The correspondingΠ i for the two loops arẽ (3.5) After taking the derivative, using the limit q 4 → 0 and the projectors, we have for everyΠ i a large set of scalar two-loop integrals depending on two external momenta, q 1 , and q 2 , which can be parametrized as (3.7) Using KIRA [61] they can be reduced to the ones in Table 1, whose corresponding topologies are represented in Fig. 3. This reduction is done in d = 4 − 2 = 4 dimensions. This is a good point to discuss how we handle renomalization and regularization. Both ultraviolet and infrared divergences are regulated using dimensional regularization. We work to the lowest Table 1: List of master integrals M (i 1 , ..., i 7 ) needed for the massless fully off-shell triangle at two-loop order. The last one is not needed at this order.
order in α and to first order in α S in the massless quark limit. There are no counterterms needed to this order and infrared divergences must vanish since the three photon amplitude vanishes because of charge-conjugation. However, individual diagrams and master integral can be infrared and ultraviolet divergent. The quantitiesΠ i are finite and the cancellation of all divergences, up to 1/ 3 provides another good check on our calculations. Strong efforts have been made to successfully obtain compact analytical expressions for all those two-loop integrals. All of the appearing master integrals can be found in terms of classical polylogarithms in Refs. [62,63] up to the order that we need. They are collected together with their corresponding expansions in App. A. Using these we find, as expected, that all the intermediate divergences exactly cancel, leading to a result of the following form :  The remaining loop functions, C ijk (0), W ijk (0), F ijk (2) can also be found in App. A. The explicit numerical coefficients can be found in the file pitildes.txt of the supplementary material. In Table 2 we give numerical results for them in a benchmark point, (Q 2 1 , Q 2 2 , Q 2 3 ) = (1, 1.3, 1.7) GeV 2 , giving also the analogous quark loop ones for comparison. The loop corrections are found of the order of ∼ − αs π , i.e. they are found to be small as far as α s is not large. The scale at which α s should be set is similar to the scale Q 2 at which theΠ are evaluated. Otherwise, large logarithms ln µ Q i appearing at higher orders would break the perturbative series. As a consequence, the series are found to be reliable as far as we do not go below ∼ 1 GeV.
Taking the linear combinations of theΠ i which lead to theΠ, one finds analogous expressions for them, but with explicit negative powers of Källén functions λ = (Q 2 1 + Q 2 2 − Q 2 3 ) 2 − 4Q 2 1 Q 2 2 . They introduce singularities which are, however, spurious. When expanding around them, they cancel against the zeros of the polylogarithms, as explicitly checked in different kinematic limits. Details on these expansions can be found in App. A. Numerical values for theΠ i in the same benchmark point and comparison with the corresponding quark loop are given in Table 3. The analytical expressions are too long to be included here and are given in the file resultsgluon.txt of the supplementary material. The equivalent results for the massless quark loop are in the file resultsquark.txt. We have, however, included analytical expressions for both the quark loop and gluonic correction at the symmetric point Q 1 = Q 2 = Q 3 in App. B.
A possible check on our result is taking the limit where one of the virtualities is much smaller than the other two, i.e. the limit of Ref. [45], where it is argued that the leading term should have no corrections. The consequences of this limit for theΠ i has been analyzed in Ref. [  . Obviously, the identification of theΠ i with the ones obtained using the OPE only makes sense when such an expansion is valid, i.e., above some cut for the Euclidean momenta, Q 1,2,3 > Q min . We restrict ourselves to those integration regions, keeping in mind that the (dominant) contributions from the remaining regions, necessarily computed with non-perturbative methods, must be added to the ones computed here.
The numerical integration has been done with the VEGAS implementation in the CUBA library, as well as our own implementation of two deterministic algorithms. We have checked that the results agree. The general expressions for the quark loop and the gluonic corrections have large negative powers of λ and become numerically unstable whenever λ is small. We therefore use, as in our previous work for the quark loop [57], expansions whenever that happens. There are six different expansions that need to be done. This is explained in more detail in App. B. We have checked that the numerical results are not sensitive to changing the boundaries where we use the different expansions.
We perform the integrals of the 12Π i contributions at different Q min , both for the leading OPE contribution, the quark loop, and the gluonic corrections. They are displayed for Q min = 1 GeV in Table 4.
Consistently with the size of the gluonic corrections found for theΠ i in the previous section, we find that they are negative and of order αs π . Given the power fall-off of the contributions of theΠ with respect to the studied energies, the quantitative contribution above some energy cut Q min is saturated by the regions nearby such a cut. As a consequence, a natural scale to effectively avoid large logarithms in the corresponding perturbative series is µ ∼ Q min , however the exact choice of it is ambiguous. In order to estimate perturbative uncertainties we will vary the scale dependence, a consequence of cutting the series at two-loops or at the first α S correction, in the interval µ 2 ∈ ( 1 2 , 2)Q 2 min . At the studied order, the whole scale dependence comes from α s (µ). Taking α N f =5 s (M Z ) from Ref. [67], we run it at five loops to α N f =3 s (m τ ). 3 As a further conservative 3 We implement the running, in the conventional MS scheme, using Rundec [68].  Table 4: Leading contributions to the (g − 2) µ integration from Q min = 1 GeV in 10 −11 units.
estimates of perturbative uncertainty, we add quadratically the difference obtained by taking the one obtained running from α with the five loop running (which we take as our central result) with the one obtained keeping α N f =3 s with a fixed scale, µ = m τ , which at the order we are working with is also a legitimate choice. Finally we also add quadratically the subleading uncertainty coming from α The result, where we show the quark loop, the gluonic corrections and the obtained uncertainties is shown in Figure 4. While in general we consider our uncertainty estimates reliable, we notice that they may be slightly over-conservative in the region just below 1 GeV and over-optimistic just above it. This is a consequence of the sharp break down of the α s running at µ ∼ 0.7 GeV which makes our uncertainty strongly dependent on the exact scale interval chosen to estimate them. In essence, we find that the correction is small and negative and that the series are well-behaved, having a gluonic correction of around −10% above the perturbative breakdown.

Conclusions
One of the main sources of uncertainties entering in (g − 2) µ comes from the contributions of the short-distance regions of the HLbL tensor contributions. In this work, which can be regarded as a continuation of Refs. [54] and [57], we have culminated our task of giving a precise and systematic description of the contributions for three large loop momenta.
For years, it was assumed that some form of the quark loop, maybe with constituent quark masses, should be the leading order of some systematic expansion of the HLbL contribution tensor to the (g − 2) µ for large loop momenta. However, it was shown in Ref. [54] how applying an OPE directly to the HLbL tensor, where the massless quark loop is indeed the leading order, does not make sense for the (g − 2) µ kinematics. The correct expansion in this kinematic region was presented in that reference, where the massless quark loop was shown to be the leading order and the leading non-perturbative quark mass-suppressed correction was computed. A very comprehensive analysis to study the role of both the quark mass-suppressed and not suppressed non-perturbative corrections to the expansion was made in Ref. [57], where many formal aspects and subtleties of the expansion were developed and presented in full detail, showing that it is well founded. The obtained results showed how above 1 GeV the non-perturbative corrections, even when functionally more important than in other expansions, are still typically below 1%.
In view of that, the most important corrections to the leading massless quark loop, and the one that ultimately allows to understand from where the expansion is valid, is the pure gluonic correction, which has been the subject of this work.
While in principle a multi-scale four-loop integral could be regarded as a formidable task, it has become feasible through combining existing tools developed for generic contributions of the HLbL tensor to the (g − 2) µ , methods on finding compact expressions developed in Ref. [57], optimized software on reduction to master integrals, analytic reduction of those remaining master integrals and numerical integration routines.
Our final result brings good news. The size of the gluonic corrections are found small, typically of size −10% above the perturbative breakdown scale, and, as a consequence, the expansion is able to give a precise description of the (g − 2) µ contributions above it.
Taking all of this into account, we suggest as a legitimate method to compute the HLbL contribution to (g − 2) µ to use the results of this expansion from some point between Q min = 1 GeV and Q min = 2 GeV, which should give a more precise prediction than resonance models, and possible discontinuities in the matching should be incorporated as systematic model uncertainty.

A Master integrals
The aim of this appendix is to list the expressions of the master integrals needed in Sec. 3 (see also Fig. 3). They can be found in Refs. [62,63,69] 4 . All n-loop master integrals contain the overall factor S n D , where The functions C 123 and W 312 are finite in the limit d → 4. Their -expansions can be written as follows W 312 ≡ dp 1 dp 2 1 p 2 where q 1 + q 2 + q 3 = 0 and d = 4 − 2 . The integrations are defined as The remaining integrals appearing in the calculation contain some singularities which cancel out in the final amplitude. The -expansions of the integral B q and S q read and S q ≡ dp 1 dp 2 The integrals V 123 andV 123 can be expressed as functions of the finite integrals C 123 and W 312 . Their -expansions can be written aṡ V 312 ≡ dp 1 dp 2 1 and V 312 ≡ dp 1 dp 2 + log(−q 2 1 ) log(−q 2 3 ) − log(−q 2 2 ) + log(−q 2 2 ) log(−q 2 3 ) + log 2 (−q 2 3 ) − 10 log(−q 2 3 ) + 19 where C 123 (i) and W 312 (i) are the coefficients of the -expansion of their corresponding Master integral (c.f. (A.2) and (A.3)). Not all these coefficients survive in the gluonic corrections we are computing in section 3: C 123 (1), C 123 (2), W 312 (1) and W 312 (2) cancel in the final expression. The coefficients C 123 (0) and W 312 (0) that contribute, as well as the function F 312 (2) which appears in the expansion of V 312 , are given below and finally and z is given by The P i (z) are real (purely imaginary) functions over the complex plane when i is odd (real). They can be expressed using polylogarithms: The polylogarithms can be defined recursively 14) The P i satisfy a number of relations P 2 (z) = −P 2 (1/z) , P 3 (z) = P 3 (1/z) , P 4 (z) = −P 4 (1/z) , P 2 (z) = P 2 (1 − 1/z) = −P 2 (1 − z) = P 2 (1/(1 − z)) = −P 2 (z/(z − 1)) , which can be used to show that the master integrals have the required symmetries under interchange of momenta.

B Analytical formulae
In this section we present analytical formulae for the scalar functions entering into the calculation of a HLbL µ . We in particular discuss the momentum expansions of the master integrals needed to make spurious singularities cancel numerically. As an explicit example, we also give the expressions for theΠ at the symmetric point Q 1 = Q 2 = Q 3 for the quark loop and gluonic correction in App. B.2.

B.1 Expansions
In the numerical evaluation of a HLbL µ there are certain limits of the kinematics requiring particular care. The integration domain can be divided into several regions as in Fig. 5. We there see the so-called side, corner and inside regions together with their boundaries. Also the cut-off µ has been indicated. Unless the side and corner regions are properly taken care of, the numerical integration will diverge as one obtains zeros in denominators that numerically do not cancel the zeros in numerators. Below we discuss the two types of problematic regions.
The precise definition of the regions is: Q i ≥ µ = Q min . The corners are defined by Q i /(Q 1 + Q 2 + Q 3 ) ≤ 1 for i = 1, 2, 3. The sides are the part of the remaining region that satisfy (2Q i /(Q 1 + Q 2 + Q 3 ) − 1 ≤ 2 for i = 1, 2, 3. The inside is the remaining allowed region.

B.1.1 Side regions
The side regions are defined as the kinematical limit where one Q i is close to Q j + Q k , or, in other words when Side Region S i : where δ is a small parameter. The inverse powers of the Källén function in theΠ i diverge in the side regions. These apparent singularities do, however, cancel when all the kinematical factors, the master integrals C 123 (0), W 312 (0), W 213 (0), W 123 (0), F 312 (2), F 213 (2) and F 123 (2) as well as the Källén function itself are expanded in δ. For a finite result we have to expand the master integrals up to order δ 9 . The analytical forms of these expansions are very long and we here therefore only give the first two orders for one case, S 3 . In the supplementary file sideexpansions.txt, however, we provide the full expansions needed for all S i . In region S +3 we have To obtain these one has to expand the relevant P i functions around a general z. Note that one obtains e.g. P 3 (1 + Q 2 /Q1), which, from the definition of the function in (A.13) gives rise to log (−Q 2 /Q 1 ). This lies on the branch-cut of the (poly-)logarithm. However, the P i are wellbehaved, single-valued functions without branch-cuts and one can safely neglect these issues. When theΠ are expanded, the negative powers of δ cancel. The expressions are not displayed here due to their length, but they can be found in the supplementary file resultsgluon.txt. Equivalent expressions are provided for the quark loop in resultsquark.txt.

B.1.2 Corner regions
In the corner regions the situation is different. There one has two small parameters instead of one Corner Region C i : Q i Q j , Q k and δ ≡ Q j − Q k Q i ≡ Q j + Q k . (B.9) Below, we list the expansions in the region C 3 5 . The expansions of the master integrals are given by . While all these expressions are finite when the small parameter δ tends to zero, some of them diverge when Q 2 3 → 0. However, this divergence has no physical meaning, since that limit lies outside the region of validity of the OPE. The expansion of theΠ in the other corner regions can be found in the supplementary file resultsgluon.txt. Equivalent expressions are provided for the quark loop in resultsquark.txt.

B.2 Symmetric Point
In this section, we write the expressions for theΠ at the symmetric point Q 1 = Q 2 = Q 3 = Q. 6 For the quark loop these are These expressions agree with those given in Ref. [36].