Short-distance HLbL contributions to the muon anomalous magnetic moment beyond perturbation theory

The hadronic light-by-light contribution to the muon anomalous magnetic moment depends on an integration over three off-shell momenta squared ($Q_i^2$) of the correlator of four electromagnetic currents and the fourth leg at zero momentum. We derive the short-distance expansion of this correlator in the limit where all three $Q_i^2$ are large and in the Euclidean domain in QCD. This is done via a systematic operator product expansion (OPE) in a background field which we construct. The leading order term in the expansion is the massless quark loop. We also compute the non-perturbative part of the next-to-leading contribution, which is suppressed by quark masses, and the chiral limit part of the next-to-next-to leading contributions to the OPE. We build a renormalisation program for the OPE. The numerical role of the higher-order contributions is estimated and found to be small.


Introduction
The Standard Model (SM) is the theoretical framework developed to describe particle physics at its most fundamental level, and is able to predict the anomalous magnetic moments of the leptons with a high number of significant digits. At the current level of precision, all the building blocks of the SM leave sizable numerical imprints for the muon anomalous magnetic moment, or, a µ = (g − 2) µ /2. Summing all the contributions, one finds [1] a SM µ = 116 591 810(43) × 10 −11 , showing a 3.7σ tension with the very precise experimental value [2,3], a exp µ = 116 592 091(63) × 10 −11 . (1. 2) The experimental value is expected to be significantly improved [4,5]. In case the discrepancy grows this could be a sign of new physics beyond the SM. The current uncertainties in the theoretical prediction are dominated by contributions from the hadronic sector. Since the relevant energy scale, i.e. the muon mass, is far below the region of applicability of perturbative QCD, the assessment of these contributions resorts to the use of non-perturbative tools. Further improvements are needed in order to find a SM value at the level of precision competing with that of the future experimental one. Decreasing the errors on the hadronic contributions would therefore shed some light on whether or not the current tension is a hint of new physics. An overview and assessment of the current theoretical situation is the white-paper [1].
In this paper, we focus on the hadronic light-by-light (HLbL) scattering contribution, represented by the diagram in Fig. 1. The calculation of the (g − 2) µ requires the integration of the HLbL tensor over q 1 , q 2 and q 3 , with the fourth photon in the static limit, i.e. q 4 → 0. Working with three Euclidean squared loop momenta q 2 i = −Q 2 i , this means one has to consider different kinematic regions of Q 2 i . We consider the short-distance regime with photon virtualities Q 2 1 ∼ Q 2 2 ∼ Q 2 3 Λ 2 QCD , and derive so-called short-distance constraints by means of an operator product expansion (OPE). The second important regime is with mixed virtualities, namely Q 2 3 Λ 2 QCD Q 2 1 ∼ Q 2 2 , and has been considered in Ref. [6]. There has been a lot of recent work in the latter regime, examples are Refs. [7][8][9][10][11][12][13].
The first full calculations of the HLbL were made using models in Refs. [14][15][16]. More recently a dispersion theory based approach as in Refs. [17,18] has allowed for better control of the low-energy region. In the latter approach one considers individual intermediate states, for which short-distance constraints such as those derived herein can be used [7,8,13,15,19].
The naive OPE in the vacuum for the HLbL tensor, which is valid for Q 2 1 ∼ Q 2 2 ∼ Q 3 3 ∼ Q 2 4 Λ 2 QCD , has the perturbative quark loop as its first term. The quark loop has always been used as an estimate for the whole contribution, using constituent quarks see e.g. Ref. [20], and for the contributions from heavy quarks. However, the naive OPE breaks down for the (g − 2) µ kinematics with q 4 → 0 [21]. The OPE of the HLbL tensor in this kinematics must be performed by taking into account that the static photon needs to be formulated as a soft degree of freedom. It was shown in Ref. [21] (see also Ref. [6]) how this could be done by factoring out the soft photon as a background field. The background field can originate either from the hard degrees of freedom (e.g. the massless quark loop) or the soft ones (e.g. the so-called di-quark magnetic susceptibility contribution). The Π µ 1 µ 2 µ 3 µ 4 (q 1 , q 2 , q 3 ) ≡ −i This defines q 4 as the negative of that in Ref. [21], which again is a choice to maximise the number of explicit symmetries. The conservation of the electromagnetic current implies that the HLbL tensor satisfies the following Ward identities, q i, µ i Π µ 1 µ 2 µ 3 µ 4 (q 1 , q 2 , q 3 ) = 0, ∀i ∈ [1,4] , (2.3) where q 4 must be rewritten in terms of the other three momenta through (2.2). Note that the last Ward identity implies that all the information on the HLbL tensor is contained in its derivative [26], In the (g−2) µ kinematics, the loop integral over the loop momenta q 1 , q 2 and q 3 can be rewritten as an integral over the Euclidean momenta Q 2 i ≡ −q 2 i > 0. The fourth photon, i.e. with momentum q 4 , is taken in the static limit. This corresponds to taking the q 4 = −q 1 − q 2 − q 3 −→ 0 limit after doing the derivative in (2.4). Notice how in this limit lim q 4 →0 ∂Π µ 1 µ 2 µ 3 ν 4 ∂q 4, µ 4 (q 1 , q 2 , q 3 ) = − lim ∂Π µ 1 µ 2 µ 3 µ 4 ∂q 4, ν 4 (q 1 , q 2 , q 3 ) , (2.5) i.e. it is anti-symmetric in the indices µ 4 ν 4 . This can be proven by multiplying both sides of (2.4) by q 4, µ 4 , then taking the derivative with respect to q 4,α and setting α = ν 4 . The resulting linear and anti-symmetric structure of the HLbL tensor is directly related to the fact that F µν ≡ ∂ µ A ν − ∂ ν A µ is the lowest dimension gauge invariant photon operator. Let us take advantage of the work of Refs. [17,18] to find general relations between the tensor ∂Π µ 1 µ 2 µ 3 ν 4 ∂q 4, µ 4 and its explicit contribution to the (g − 2) µ . Making use of the Ward identities above, one can rewrite in full generality the HLbL tensor as a linear combination of 54 scalar functions Π i (q 1 , q 2 , q 3 ) according to [17,18] Π µ 1 µ 2 µ 3 µ 4 (q 1 , q 2 , q 3 ) = 54 i=1 T µ 1 µ 2 µ 3 µ 4 i Π i (q 1 , q 2 , q 3 ) . (2.6) Expressions for the T µ 1 µ 2 µ 3 µ 4 i in terms of the Lorentz basis built with q 1 , q 2 , q 3 and the metric g µν can be found in Refs. [17,18]. In particular, the T µ 1 µ 2 µ 3 µ 4 i satisfy the same Ward identities as the HLbL tensor. As a consequence, ∂T µ 1 µ 2 µ 3 ν 4 i ∂q µ 4 4 Π i (q 1 , q 2 , q 3 ) . (2.7) Taking the static limit q 4 → 0, the remaining tensor can be written as a function of 19 linear combinationsΠ i of the original Π i (q 1 , q 2 , q 3 ). Only 6Π i functions contribute to the (g − 2) µ and one finds a HLbL where the integration variable τ is defined via Q 2 3 = Q 2 1 +Q 2 2 +2τ Q 1 Q 2 . The functions T i (Q 1 , Q 2 , τ ) can be found in Refs. [17,18] and (2.9) The exact definition of theΠ i functions is given in Refs. [17,18]. The C ij in (2.9) represent interchanging q i and q j . In order to find general Lorentz projectors from the derivative of the tensor in the static limit, i.e. lim q 4 →0 ∂Π µ 1 µ 2 µ 3 ν 4 ∂q µ 4 4 , (2.10) to theΠ i functions, we start by taking 19 independent projectors in the {q 1 , q 2 , q 3 , g} basis. Any other projector can be related to that set through the Ward identities given above. Applying them to (2.7) returns a system of 19 equations dependent on the 19Π i . A solution to that system of equations for the relevantΠ i is given in App. A. In practice, this means that for any contribution to the HLbL tensor in any basis, one can compute the associated (g − 2) µ contribution by taking the derivative with respect to q 4 , then taking the static limit q 4 → 0, applying the 6 Lorentz projectors given in App. A to find the associatedΠ i and finally using them in the integral of (2.8).
As explained above, the integral of (2.8) requires the knowledge of the HLbL tensor with three Euclidean momenta Q i at different kinematic regions and the fourth in the static limit q 4 → 0. In this work, which extends the results of Ref. [21], we focus on the kinematic region where the three loop momenta are large. As was shown in Ref. [21], if one defines Π µ 1 µ 2 µ 3 (q 1 , q 2 ) ≡ − 1 e d 4 q 3 (2π) 4 then in the static limit for the studied kinematic region, one can factor out the soft photon contributions according to Π µ 1 µ 2 µ 3 (q 1 , q 2 ) = Π µ 1 µ 2 µ 3 µ 4 ν 4 F (q 1 , q 2 ) 0|e q F ν 4 µ 4 |γ(q 4 ) . (2.12) As mentioned in the introduction, the soft background field F µν can originate from the hard degrees of freedom or the soft ones. One then finally has [21] lim . (2.13) For this OPE the massless quark loop is the leading term and the di-quark magnetic susceptibility of the vacuum is the leading, quark-mass suppressed, non-perturbative contribution. In this work we compute the leading non-perturbative (not mass-suppressed) corrections.

The operator product expansion: a theoretical description
In this section we describe the OPE and the associated renormalisation program. From this we systematically separate the long-distance effects from the short-distance ones. The general framework and operators involved are presented in Sec. 3.1 whereas the mixing of these operators is elaborated on in Sec. 3.2. Finally, the OPE developed as well as operators and corresponding matrix elements involved are discussed.

General framework
Perturbative calculations are known to provide a huge predictive power in the framework of Quantum Field Theory. However, when the calculation of a given observable involves the interplay of two (or more) very different scales, large logarithms between them slow down, if not spoil, the convergence of the series. These large logarithms can be avoided in many cases through the OPE [27,28], which integrates out the heavy degrees of freedom leaving the low-energy (longdistance) dependence encoded in effective operators, in such a way that the contributions from higher-dimension operators become suppressed by extra powers of the high-energy scale [29][30][31]. There are cases in which this logarithmic re-summation is not enough, since one (or several) relevant couplings of the theory diverge when its running is performed. This is the case for QCD, where a low-energy description in terms of approximately free quarks and gluons does not hold and the matrix elements between initial and final states cannot be computed within perturbative QCD. The contributions from the operators whose quantum numbers are compatible with the transition must be fitted to data, computed with effective field theories or other non-perturbative methods, such as lattice QCD, dispersion relations or model estimates.
In the OPE of two-point correlation functions [25], all the local operators with the same quantum numbers as the QCD vacuum, such as the identity orqq, can give a contribution to the Green functions 2 . In the OPE used in flavour-breaking transitions (e.g. see Ref. [29]) all the local operators with quantum numbers compatible with the studied transition among hadrons can give a contribution. In the OPE we are working with [22,23], applied to (2.11), any local operator with the same quantum numbers as F µν can absorb the remaining soft static photon and, as a consequence, give a contribution [21,24]. Higher-dimensional operators are suppressed by extra powers of Λ QCD Q i d , providing a hierarchy of contributions with a systematic counting. Up to dimension 6 and order α s a basis of those operators is 3 We use here the notation of Ref. [32]. In particular, G µν = ig S λ a G a µν ,Ḡ µν ≡ i 2 µνλρ G λρ , covariant derivatives act on all objects to their right and tr γ 5 γ µ γ ν γ α γ β = −4iε µναβ . For S 1...7,µν the quark field q refers to a given flavour and there are in principle different operators for different flavours q. Notice, however, that taking into account that the (massless) QCD vacuum preserves SU (3) V , the contributions of the operators to the studied transition depend, in the chiral limit (m u = m d = m s = 0), on a common constant multiplied by the corresponding quark electric charge. Note that with the conventions of Ref. [32], G µν is order g s . 4 The four-quark operators are only indicated generically in (3.8). A decomposition valid in the chiral limit into twelve operators is given in App. B. In fact, only two combinations of four-quark operators contribute as discussed in Sec. 4. The adopted notation is analogous to the one in Ref. [29].
In order to perform the OPE of the tensor in (2.11), one applies Wick's theorem with any number of needed (suppressed) extra (QCD or QED) vertices coming from the Dyson series. The uncontracted operators must then be Taylor expanded (e.g. see Ref. [32]), so that the resulting expression is of the form 5 where X S contains the magnetic susceptibilities of the operators, 6 S i,µν = ee q X i S F µν .
3 Notice how the order in α s depends on the short-distance structure of the studied Green function. For example in baryon sum rules, the four-quark operators S {8} do not enter α s -suppressed. In our calculation, at order g 3 s one may have a contribution from a three-gluon operator [22], but it enters suppressed by g s , loops and flavour (in the SU (3) V limit its contribution vanishes, since it transforms as a singlet and the photon field transforms as an octet). 4 In fact, the gluon tensor is named F µν in that reference. We rename it as G µν to avoid ambiguity with the electromagnetic field strength F µν . Similar renamings should be obvious. 5 Further technical details on the computation of the different pieces are given in Sec. 4. 6 We will define the susceptibilities more precisely later. The di-quark magnetic susceptibility. A complete separation between short-distance and long-distance effects should subtract possible divergent low-momenta contributions from the quark loop and possible divergent perturbative series arising from soft lines.
Even when this was a first step to achieve the separation between short-distance and longdistance effects, such a separation is not yet complete. Let us illustrate this with the simplest contributions: the quark loop and the (di-quark) magnetic susceptibility, represented in Fig. 2. In Fig. 2b one has short-distance contributions that arise from expanding the Dyson series and introducing vertices in the soft lines, represented by a blob. Since there is no momentum flowing, the resulting series, c n α n s (0), is manifestly divergent. This kind of effects must be subtracted, since they belong to the non-perturbative domain. A slightly different problem arises with Fig. 2a. The quark loop runs over all possible momenta. If the low-momenta contributions do not vanish, they must somehow be subtracted and reabsorbed into the long distance matrix elements. This is the case for the O(m 2 q ) mass corrections, whose (not regulated) Wilson coefficient, C m 2 q , leads to divergent series n c n α n s (Q 2 ) log n Q 2 /m 2 q , i.e. the coefficients are not well defined in the chiral limit. The m 2 q dependence in the coefficients originates from the low-energy domain and must be included in the operator expectation values as well.
The way of achieving these subtractions is by dressing and renormalising the S µν operators (normal-ordered operators in the notation of Ref. [33], tree-level operators in the notation of Ref. [29]). Following a notation close to the one in Ref. [29], the dressed operators Q µν 0 in terms of the tree-level ones S µν are obtained by reinserting them into the Dyson series, leading to a result of the form where the = − d−4 2 dependence appears in dimensional regularization as a consequence of ultraviolet divergences in the loop diagrams. The infrared divergences that can appear are regularised using the quark mass. The resulting redefinition of the matrix elements 7 involves two steps, calculating the relevant contributions leading to an expression at one-loop order of the form whereMˆ andM rem are perturbations either in e, in g s or in powers of mq Λ QCD and 1 = 1 − γ E + log(4π). a depends on the dimension of the operators when d = 4. The ultraviolet divergences are unphysical, and are removed via renormalisation. A convenient and simple renormalisation scheme is by performing the operator renormalisation in the M S scheme, which basically removes the 1 factors and takes out from the bare operators the non-canonical part of their dimension, proportional to 2a : (3.14) Putting this back into (3.9) one finds This equation defines the regularized and renormalised Wilson coefficients. The renormalised Wilson coefficients, become free from long-distance effects and infrared divergences, completing the desired separation.
For the matrix elements we define the magnetic susceptibilities X = (X 1 , . . .). 8 Q M S,µν (µ) = e X e q F µν (3.17) The contributions we calculate explicitly in Sec. 4 are: The leading contribution stems from Q µν 1 at d = 2 and corresponds to the massless quark loop. This operator receives mass corrections suppressed by powers of m 2 q Λ 2 . The first correction from a different operator comes from the d = 3 di-quark operator, Q µν 2 , which happens to enter suppressed by one power of the quark mass becoming effectively d = 4. This contribution using the mixing as defined in (3.11-3.14) removes the m 2 q log(m 2 q ) part of the quark loop and together with that forms the d = 4 contribution. The d = 5 contributions from Q µν 3−5 are also suppressed by one power of the quark mass. Q µν 6−8 give the first contributions that are not suppressed by the quark masses. The mixing here again removes contributions proportional to log(m 2 q ) and other infrared divergences. We therefore calculate with respect to the massless quark loop the corrections of orders g s through operator mixing as defined above and calculated in Sec. 3.2. This is analogous to the mixing between bi-linear quark condensate and gluon condensate in the usual vacuum OPE [34][35][36]. Figure 3: An example of a cancellation for contributions arising from attaching extra vertices to soft, at zero momentum, lines. The circle with "2" inside refers to the S 2 operator defined in (3.2), B indicates the extra factor compared to the left diagrams. All Lorentz indices have been suppressed for clarity.

The operator mixing
In this section we calculate the mixing matrixÛ M S . However, there are some parts that can be ignored. These we discuss first. The original Wilson coefficients contain in principle infrared divergent parts that arise from attaching extra vertices to the soft zero-momentum lines. These long-distance contributions, explicitly independent of momenta, systematically cancel with analogous terms in the dressing procedure, which in addition are independent of the studied physical process, so we simply ignore them in both sides of the calculation. An example of this is sketched in Fig. 3. The top graphs show the contribution from S 2,µν to the correlator we calculate and the bottom line contributions to the mixing matrix. The contribution from the top right is via the diagram in the bottom right absorbed into the definition of the operator Q 2,µν .
As a consequence, only mixing terms with loops need to be considered for the calculation of U M S . This limits the type of terms that can show up. Then, at this stage, we have two possibilities. 1) If no more (QCD or QED) vertices are added, the original operator can only mix with lower dimensional ones, since loops imply connecting the fields of both operators with propagators. The only remaining scale to compensate dimensions is m q and as a consequenceM ( ) becomes a m n q perturbation, with n positive. The calculation forÛ 21 M S is of this type. 2) Alternatively, one can have mixing terms with operators with equal or higher dimensions that show up by adding fields through introducing extra vertices. However, these new vertices come with an extra cost (g s or e), and then they can also be regarded as perturbations, an example of this type isÛ 76 M S . Finally, when studying the mixing of lower dimensional operators with higher dimensional ones, one finds terms that go as m −n q . A well-known case of this in the usual vacuum OPE is the gluon condensate mixing with the quark condensate However this kind of mixing is simply absorbing in the renormalised operator unphysical (infrared) divergences contained in the perturbative Wilson coefficients. The quark mass is used as an infrared  regulator as well as for calculating genuine quark mass corrections. 9 The first step is thus the dressing of the tree level operators and the computation of the associated matrixM ( ). All the diagrams involved are shown in Fig. 4. The first diagram for each operator always corresponds to identity term in (3.11). The procedure should be done, as is also the case for the usual vacuum OPE, only including operators up to the dimension considered.
The only way Q µν 1 mixes with other operators is by introducing extra electromagnetic vertices. Since Q µν 1 is already O(e), the resulting corrections are O(e 2 ), and thus do not need to be considered. The first non-trivial case is that of the operator Q µν 2 . The non-zero element ofM (ε) and thusÛ M S are calculated in the next subsubsections. As in Sec. 4 the calculations in this section benefit very much from using the radial gauge.
The mixing with Q µν 1 = ee q F µν requires an extra electromagnetic vertex and closing the quark lines. This corresponds to the second Q 2 -diagram in Fig. 4. Let us sketch it in some detail to 9 The procedure is equivalent to what is used in the usual vacuum OPE [33][34][35][36][37].
illustrate the procedure: where we have used that we are working in the radial gauge, A µ (x) = 1 2 x ν F νµ [32,34]. 10 Taking into account (3.11) and (3.14), one findŝ (3.20) Through the M S renormalisation we have introduced a subtraction point, µ, which is the scale of separation between long-distance and short-distance effects. In particular this mixing term, when plugged into (3.16), removes the (infrared divergent) long-distance parts from the loops associated to the perturbative O(m 2 q ) corrections by effectively replacing log Q 2 i /m 2 q → log Q 2 i /µ 2 . 11 Notice how, in contrast with the usual vacuum OPE, where these divergences can only start at O(m 4 q ) (and be regulated by the quark condensate), the divergences in this OPE with respect to the quark loop start at O(m 2 q ). They become regulated by Q µν 2 . At the order we are working with, Q µν 2 also mixes with Q µν 6 through the third Q 2 -diagram in (3.21) The Wilson coefficient proportional to S µν 2 m 2 q combined with this matrix element cancels the power divergence of the (unregulated) Wilson coefficient associated to S µν 6 , while the mass correction to the Wilson coefficient associated to S µν 1 gives a finite contribution that needs to be included. Once again, this interplay is analogous to the corresponding one in the vacuum OPE [33][34][35][36].

3−6
Since Q µν 3 and Q µν 4 are already O(g s ) and two lines need to be closed to form a loop to give a non-zero contribution, they only mix with the gluon matrix element at the order we are working with (see Fig. 4), giving a finite contribution to its associated regulated Wilson coefficient. One finds, Since Q µν 5 is already order e, it can only mix with Q µν 1 and Q µν 6 at the order we are working with. Dimensional analysis shows that the mixing with Q µν 1 only modifies the very tiny and safely neglected O(m 4 q ) contribution. The mixing with Q µν 6 is the same as the mixing of the quark condensate with the gluon condensate, as in (3.18), since at the order in e we are working with the photon does not see strong interactions [25,[33][34][35][36], The last operator to treat is Q µν 7 which through the topologies shown in Fig. 4 only mixes with Q µν 6 . One has

Full matrixÛ M S (µ)
Putting all the elements together one finds for N c = 3

Values of the matrix elements
Apart from providing model-independent information on the kinematic dependence of the nonperturbative corrections for the short-distance HLbL, in principle this OPE can be used to study different Green functions with all its (Euclidean) momenta large except for one soft photon. This might also allow to obtain more information on the expectation values (or the susceptibilities). However in absence of the latter we need to determine the values in a different way. We can find values for all of them with a number of assumptions that should at least give the correct order of magnitude.
The matrix elements associated to Q µν 5 and Q µν 6 are directly related to the quark and the gluon condensates, respectively. The former one is well-known, since it is the order parameter of the spontaneous chiral symmetry breaking of QCD, both from chiral perturbation theory and the lattice. Updated lattice values can be found in Ref. [38]. The gluon condensate is not so wellknown, since separating its effect from those of the perturbative series is non-trivial. However its order of magnitude is known, namely X 6 ∼ 0.02 GeV 4 [25].
The matrix elements 0|Q µν i |γ(q 4 ) = X i 0|e e q F µν |γ(q 4 ) are directly related to the values of QCD two-point functions at zero momentum, Π V Q i (q 2 ). In order to see that, let us take a generic dressed operator Q µν i,0 . Its contribution to a matrix element with a final static photon can only arise through an extra electromagnetic vertex. Then, in the static limit, These two-point functions can be computed at large Euclidean momenta through the OPE in the vacuum [25]. One could compute these matrix elements in chiral perturbation theory. In fact, promoting the global SU (3) V symmetries to local ones lead to trivial transformations for the operators Q µν 2 , Q µν 3 , Q µν 4 and Q µν 8 . The resulting effective low-energy Lagrangians are given in Ref. [39]. For X 2 , X 3 and X 4 we obtain an educated guess by making use of (3.27) and extrapolating the argument for X 2 from Ref. [40]. First of all, in the large-N c limit the QCD spectrum is made of an infinite number of free, stable meson states [41][42][43]. The two-point functions in that limit are then saturated by the exchange of resonances. Owing to the quantum numbers of the studied twopoint functions, the corresponding resonances must be vector mesons [44]. The low-energy part of the N c = 3 QCD spectrum is actually close to the sum of narrow width resonances predicted by the large-N c limit, while at higher energies a transition towards the flat perturbative QCD spectrum is observed. Taking all this into account, and that in the chiral limit the V T two-point function vanishes in the perturbative regime, it is reasonable to assume that the physical spectrum is saturated by the contribution of the lowest vector meson, i.e. the ρ meson. Using the formalism developed in Ref. [45] and adding a tensor source [39,40], one can write the two-point functions of (3.28) on the form (3.29) Here, the C T i are constants. It then follows that (3.30) 12 Since the symmetry transformation of Q µν 3 , Q µν In order to estimate the C T i , we can match the ansatz of (3.29) with the short-distance OPE of (3.28). This is sometimes referred to as a vector-meson-dominance (VMD) estimate. We find This leads to The obtained value for X 2 is, in fact, in very good agreement with a precise lattice determination [46,47]. 13 No precise modern evaluation of m 2 0 can be found in the literature. However, once again, some numerical estimates are available [48]. The estimates of X 3 and X 4 are new to the best of our knowledge.
Notice how, once again, the mass correction to (3.31) leads to an infrared divergent term which is regularised by the mixing of the tensor current with F µν of (3.16), Our largest uncertainty comes from Q 7 , where we simply perform a dimensional guess inspired in its mixing with the gluon matrix element. One has Notice how the derivative term of this operator makes non-trivial its low-energy effective realization when trying to promote invariance under the global symmetry to a local one.

8,2
For the four-quark operators there are only two combinations that contribute. These are These do not mix with other operators at the order we work so the Q µν 8,i are the same. Note that the two operators only differ in the way the quark charges appear.
For both operators we define a generalised magnetic susceptibility In the massless limit we can also use which is a definition more similar to (3.17). This is not possible for Q µν 8,2 . The operators Q µν 8,1 and Q µν 8,2 can be decomposed in a basis of 12 four-quark operators containing different flavour and Dirac matrices (see App. B for details on the reduction). However, in the large-N c limit not all of these survive. In particular, one finds that only two are non-vanishing due to the factorisation of two colour singlet currents in this limit. These are In the first equality we have performed the sum over flavour indices for the quark condensate. In the second step we have projected out the flavour matrix λ 8 which is traced with the quark charge matrix. Here one sees why only these two matrix elements can survive. First of all, the quark charge matrix is a linear combination of λ 3 and λ 8 , so the relation Tr(λ a λ b ) = 2 δ ab implies that only matrix elements with λ 3 or λ 8 are non-zero. In addition, the only non-vanishing two-quark condensates are qq and the di-quark matrix elementqσ µν q. To leading order in N c , one therefore finds Alternatively one can directly evaluate the large-N c limit of the two matrix elements needed by using λ aαβ λ aγδ = 2δ αδ δ γβ . Then use Fierzing and the charge matrix equivalent of (3.44). The result agrees with (3.45). This way one sees also directly that q A = q B in the non-zero part of the matrix elements at large N c . The entire contribution to HLbL is then proportional to A e 4 q A .

Calculation of the HLbL contributions
In this section we consider the analytic calculations of the various contributions in the OPE discussed above. They are the fully connected quark loop, diagram topologies with one quark line non-contracted (related to two-quark operator matrix elements), diagram topologies with two quark lines non-contracted (giving rise to four-quark operator matrix elements) as well as the gluon matrix contribution. One check we have performed on all contributions is that This work relied heavily on FORM [49,50]. The Feynman integral reduction for the quark loop and the gluon matrix element contributions was done with Reduze 2 [51] and Kira [52]. In the supplementary material we provide the analytic results as FORM output in the file results.txt.

The quark loop
The quark loop contribution arises from allowing for a soft emission from one hard vertex, which is equivalent to modifying a quark propagator by an external background field [22,34], as shown explicitly in Ref. [21]. Starting from (2.11), the needed background field, F ν 4 µ 4 , is simply obtained by Taylor expanding the photon field appearing in the Dyson series. In the static limit in the radial gauge one has Define the quark propagator for a quark of mass m q as One may then write Π µ 1 µ 2 µ 3 µ 4 ν 4 F in a very compact way. After having contracted all the quark fields in the definition of the tensor in question one finds where σ(i, j, k) denotes a member of the permutation group acting on the set {i, j, k} = {(q , µ )} ∈{i,j,k} . In other words, σ(i, j, k) simply states that we sum over all permutations of momentum and Lorentz index pairs. Using iteratively the relation allows for a systematic computation of the quark loop. Applying the projectors given in App. A and reducing the (ultraviolet finite) integrals, the result is left in terms of scalar tadpole, self-energy and triangle integrals. Expanding these in the masses, 14 one findŝ (4.8) (4.10) The different coefficients (c, d, f, g, h, p, q, r, s) can be found in App. C. Explicit analytic formulas for F (Q 2 1 , Q 2 2 , Q 2 3 ) in terms of Clausen, Glaisher and L functions can be found in Ref. [53]. Spurious singularities in the λ → 0 limit cancel against the zeros of the triangle function F for the different Π i .
As explained above, one finds logarithmic infrared divergences forΠ . Rearranging the logarithms, they can be expressed on the form log As a consequence, while the massless quark loop corresponds to the leading term in the shortdistance regime, the naive mass correction does not. Infrared divergent logarithms must be substracted first.

Contributions from diagrams with one cut quark line
Several kinds of contributions need to be taken into account up to the computed order from topologies in which, starting from (2.11), one quark line is left uncontracted: see Fig. 5. There 14 Taking into account that the mass dependence goes as ∼ A + B m 2 q log(m 2 q ) + C m 2 q + O(m 4 q ), one needs to be careful in order to obtain the correct coefficients as the naive Taylor expansion does not hold. are several expansions involved. First, the uncontracted quark fields must be Taylor expanded. Working in the radial gauge both for the gluon and for the photon, the Taylor expansion [32,34] of the quark bilinears can be written as: Since our computation goes up to dimension D = 6, we need to expand up to three derivatives. A lower number of derivatives in that expansion can give contributions (apart from the di-quark magnetic susceptibility, which, as shown in Ref. [21], gives a contribution already at D = 4 when combined with masses) when combined with mass terms from the hard propagators or from soft gluons or photons coming from hard propagators. A first simplification consists in realizing that one can put the gluons and photons together with covariant derivative terms. Extending the results of Refs. [32,34], one finds (4.13) In fact, it can be shown that for a given flavour the sum of all the contributions that enter into our computation can be reduced to a compact form. Define Γ A to be an element in the set of Clifford matrices according to (4.14) The compact expression for Π µ 1 µ 2 µ 3 is then where σ(1, 2, 3) again denotes a pairwise permutation over q i and µ i . The coefficients c A are defined as such that one in a standard fashion can decompose in a spinor basis according tō Here, all dependence on the other quantum numbers such as colour or flavour has been suppressed. The proof of (4.15) up to the order that we need, i.e. up to p = 3, can be found in App. D. Note that already for p ≤ 3 the proof involves a very large cancellation of contributions, and the compact form allows for a much simplified calculation of the diagram topologies with one cut quark line.
The reduction of the matrix elements 0|qD ν 1 . . . D νn Γ A q|γ(q 4 ) into the matrix elements of Sec. 3 is rather involved. One needs to recursively exploit spinor algebra relations, symmetry transformations under Lorentz, parity and charge conjugation as well as the equations of motion of the quarks and the gluons. The resulting non-zero matrix elements are of eight types. With zero derivatives we have 1 With one derivative one has and with two they are 1 (4.22) Here, we have defined the linear combinations

23)
(4.24) For three derivatives there are two contributions. These are 1 Here,F µν ≡ i 2 µναβ F αβ and the A i are given by

31)
The operator enters from using the gluon equation of motion. Its susceptibility is defined as Using the above decompositions in (4.15) and rewriting the S µν i into the Q µν i and replacing the X i S with the corresponding X i , one findŝ The numerical coefficients c m,n,p i,j,k can be found in App. C.

Contributions from four-quark operators
The four-quark operator contributions arise from cutting two quark lines. 15 The resulting diagrams become split in two parts and, as a consequence, introduce flavour mixing. An extra gluon Figure 6: Contributions from four-quark operators obtained cutting two quark lines. All possible ways to connect the gluon to the quark lines must be considered.
propagator needs to be included from the Dyson series expansion to connect the quark lines. The resulting contribution, up to permutations of σ(1, 2, 3), is shown in Fig. 6 where the gluon can be connected in three different positions in the quark line above and in two different positions in the one below. These diagram contributions can be compactly written as where Γ ωP ∈ {γ ω , γ ω γ 5 } and Γ ω P ∈ {γ ω , −γ ω γ 5 }. In fact charge conjugation requires that A and B must be different to get a non-zero matrix element and one of the remaining two possible contributions vanishes then taking the traces. Recalling the definition of X 8,2 in (3.42) we find A reduction of all possible four-quark matrix elements into a basis of 12 independent ones is given in App. B for the chiral limit, i.e. m u = m d = m s = 0.

The gluon matrix element
The gluon matrix element contribution arises from all the possible combinations in which, starting from (2.11), one extra QED and two extra QCD vertices are added (see Fig. 7). The gauge boson fields F µν , G a µν and G b µν are then Taylor expanded according to (4.2). Since all the quark fields are connected, the colour chain always leads to the same colour trace, namely Once the colour and the space-time terms from the Taylor expansions have been factored out, the remaining six-point function is fully symmetric under the exchange of indices. Taking advantage of this symmetry, one can rewrite the whole contribution as a sum of permutations according to (1,2,4,5,6) Tr γ µ 3 S(p + q 1 + q 2 + q 4 + q 5 + q 6 )γ µ 1 S(p + q 2 + q 4 + q 5 + q 6 ) × γ µ 2 S(p + q 4 + q 5 + q 6 )γ µ 4 S(p + q 5 + q 6 )γ µ 5 S(p + q 6 )γ µ 6 S(p) . Here, σ(1, 2, 4, 5, 6) is the set of pairwise permutations of µ i and q i for i = 1, 2, 3, 5, 6. Also, the equation has been written in terms of d = 4 − 2 for renormalisation purposes. Using (4.4) iteratively in the above equation, then taking the momentum limits, calculating the Dirac trace and projecting the results into theΠ i , one finds the results in terms of ultraviolet finite integrals. We do this by reducing the loop integrals to combinations of triangle, self-energy and tadpole integrals with the help of the package KIRA, all the time carefully performing both the expansions in and in the quark masses. Without including operator mixing the result takes the form 42) where c, f, g and h are numerical coefficients given in App. C.
When operator mixing is taken into account, the found divergences, which scale as 1 , exactly cancel respectively with the X 1 and X 7 contributions. The dependence on the triangle integral also cancels, leading to a fully analytic gluon matrix element contribution. One should note that the final expressions for this contribution are very simple, even compared to the quark loop, and there are substantial cancellations along the way that lead to this simple form.

Numerical results
In this section we present numerical results for the full a HLbL µ integral in (2.8). Using the relations in (2.9) between the set of functionsΠ i and theΠ i given in the appendices as well as Sec. 4.3, we evaluate a HLbL µ for the matrix element as well as loop contributions. We use the matrix elements as estimated in Sec. 3.3. For this purpose, we use the Cuba library [54], in particular the Vegas integrator building on Monte Carlo sampling the three-dimensional integral. The results have been checked as well with an adaptive deterministic integrator implemented by us. Care has to be taken in the numerics since λ can vanish or get very small and appears with rather high negative powers in some expressions. Those areas in the integration have to be treated by expanding the loop functions around the λ = 0 points analytically, some of these limits are given explicitly in App. C.2.
We investigate the various contributions first at two benchmark values of the lower momentum cut-off, i.e. Q 1,2,3 > Q min ∈ {1, 2} GeV. Then we proceed to investigate how the different pieces scale with Q min and compare the respective sizes. For notational convenience, we refer to the contributions with respect to the corresponding X i .
At the sought precision level it is sufficient to assess the order of magnitude of these corrections. Therefore, in want of precise input we resort to simplified input as discussed in the previous section. For this purpose, we use Given the smallness of the matrix element and quark mass correction contributions we did not take into account the running with µ of the various inputs but kept them fixed. The benchmark points for Q min ∈ {1, 2} GeV are presented in Table 1. Since the contributions from the matrix element X 2 come in both suppressed at order m q and at order m 3 q , we here present the respective contributions, labelled X 2,m and X 2,m 3 , of these. The table shows that power correction are suppressed by at least two orders of magnitude with respect to the quark loop. This is also visible in Figs. 8-10 where we consider the scaling with Q min . The T i (Q 1 , Q 2 , τ ) in (2.8), when expanded for large Q i , are of order m 2 µ , except for T 1 which is m 4 µ . The variation with Q min from dimensions is thus 1/Q 2 min for the massless quark loop, 1/Q 4 min for the d = 4 contributions and 1/Q 6 min for the d = 6 contributions. The scaling is found to agree with naive dimensional counting.
The power corrections not suppressed by quark masses, i.e. X 6 , X 7 , X 8,1 and X 8,2 , are found to be numerically suppressed. This is partially explained by their extra suppression in powers of Λ QCD . Their numerical impact is similar to the one of the di-quark magnetic susceptibility, X 2 , and clearly more important than the perturbative mass contributions. We find that even at 1 GeV, all these power corrections are suppressed by at least two orders of magnitude with respect to the massless quark loop. Even though this result motivates studying whether the smallness of the corrections pinpoint a trend also for the purely perturbative ones, no strong conclusions should be derived from it.
A comprehensive description of the resulting OPE has been provided, including a detailed explanation on how to achieve a complete and systematic separation of short and long distance effects while cancelling internal divergences. The physical meaning of some of the resulting matrix elements is also given, and we have used those to present estimates of all of them. The obtained results could in the future be used to analyse other Green functions and their phenomenological applications.
The resulting OPE is applied to the HLbL in the (g−2) µ kinematics for Q 1 , Q 2 , Q 3 Λ QCD . As a consequence of setting one of the momenta to zero, the long distance effects become functionally more important. The quark loop is still found to be the dominant contribution, but the first nonperturbative correction becomes suppressed by just one power of Λ QCD (plus one power of m q ), in contrast with the Λ 3 QCD (qq) suppression in the OPE applicable when all the Euclidean momenta are large.
However, no operators allowed by the symmetries are found to enter without quark mass suppression below Λ 4 QCD with respect to the quark loop contribution. These Λ 4 QCD contributions are computed and their role for the (g − 2) µ is estimated. They are found to be very small. Whether or not this may be indicating that the quark loop gives a precise description of the HLbL at relatively low momenta (i.e. Q i ∼ 1 GeV) will only be known once the two-loop perturbative corrections have been computed. This calculation is already under way and will be presented in a future publication.  X 1,0 X 3 X 4 X 5 Figure 9: Numerical contributions from the X 3−5 pieces in absolute value using the inputs from Table 1. The massless quark loop is shown for comparison.

A A set of Lorentz projectors for theΠ i
In this appendix we present a set of projectors useful for projecting to the set ofΠ i . Note that due to gauge invariance this set is not unique. We have derived and used a second set of projectors.

B Four-quark reduction
In this section we reduce the number of four-quark matrix elements from the basisq iAᾱ q jBβqkCγ q lDδ , where the barred Greek indices denote colour, the capital ones flavour and the latin ones spinor, into one with only twelve non-zero elements.
First of all, due to confinement only colour singlet operators can give contributions. From this one has Next, taking into account that the QCD vacuum preserves SU (3) V in the flavour sector up to small quark mass corrections, the contributing four-quark operators must break SU (3) V in the same direction as the octet charge operator. There are therefore four independent flavour structures which contribute and can be taken to be The charge operator can be expressed with Gell-Mann matrices through and the third and fourth operators O 3 and O 4 can be rewritten using a Fierzlike identity: Orthogonal operatorsÕ i satisfyingÕ n = P n ABCDq iA q jBqkC q lD , with P n ABCD P m BADC = δ mn , can be obtained from linear combinations of the operators O i (i = 1, . . . , 4) as This leads toq iA q jBqkC q lD = n P n BADCÕ n . (B.13) Finally, the fact that the vacuum expectation values of the operators must be proportional to Q implies (B.16) q i λ 1 q jqk λ 1 q l = − q i λ 8 q jqk λ 8 q l = q i λ n q jqk λ n q l , for n = 2, . . . , 5 , (B.17) q i λ 1 q jqk λ 1 q l = − 1 2 q i λ n q jqk λ n q l , for n = 6, 7 , (B.18) q i λ 1 q jqk λ 1 q l = 1 2 √ 3 q i λ 3 q jqk λ 8 q l + q i λ 8 q jqk λ 3 q l , (B. 19) and the final flavour decomposition reads Only the spinor reduction remains. As usual, one can decompose each operator above according toÔ where Γ A is an element in the spinor basis of (4.14) and c A the corresponding normalisation defined in (4.16). However, many restrictions apply. Proportionality with F µν leaves a small number of independent Lorentz structures possible. Moreover, sinceÔ ijkl are by construction either symmetric or anti-symmetric under the exchange (ij) ↔ (kl), the reduced matrix element Taking advantage of this and requiring that the reduced matrix elements should be odd under charge conjugation one finds q i λ 8 q jqk q l ±q i q jqk λ 8 q l = 1 32 σ µν ji δ lk ± δ ji σ µν lk qσ µν λ 8 qqq ±qσ µν qqλ 8 q × qλ 8 γ µ 2 γ 5 qqγ ν 2 q ±qγ µ 2 γ 5 qqλ 8 γ ν 2 q + 1 32 σ µν ji γ 5 lk ± γ 5 ji σ µν lk qσ µν λ 8 qqγ 5 q ±qσ µν qqλ 8 γ 5 q .

(B.24)
This reduces the original set of 1679616 matrix elements to a basis of 12 non-zero ones.

C Explicit expressions for theΠ i
In this appendix, we seperately list the form factorsΠ i for the contributions from the quark loop, one cut quark line topologies and the gluon matrix element. Note that those for the two-cut quark line topologies were given in Sec. 4.3.

C.2 Some massless quark loop limits
Here we give the explicit limits for the massless quark loop in the regimes where two momenta become much larger than the other, this is a subset of the regions where λ becomes small. In order to find them, we change the variables to one of the large momenta, the small one and the angle between them. We then expand in the small over the large momentum. All negative powers of λ cancel and the leading contribution in the ratio of small over large momentum is always independent of the angle and is given below. The choice of third variable is not unique, however the results to the order given are always the same. (C.72)

C.3 Contributions from diagrams with one-cut quark lines
The explicit expressions of theΠ i for the one-cut quark line diagrams can be written aŝ where the non-zero coefficients are .
(C.132) D Derivation of (4.15) up to n = 3 Up to the order that we need, we have only contributions either without explicit gauge bosons or with one of them, which can be put together owing to (4.13). Let us start by writing down the expansion coming from the contributions without gauge bosons (see Fig. 5a). Starting from (2.11) and using the decomposition (4.17), one trivially finds (up to permutations of the set P = {1, 2, 3}): × Tr γ µ 3 Γ A γ µ 1 iS(x 1 − x 2 )γ µ 2 iS(x 2 − x 3 ) .