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 (Qi2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {Q}_i^2 $$\end{document}) 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 Qi2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {Q}_i^2 $$\end{document} 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 figure 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 refs. [7][8][9][10][11][12][13][14].
The first full calculations of the HLbL were made using models in refs. [15][16][17]. More recently a dispersion theory based approach as in refs. [18,19] 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, examples are refs. [7,8,13,14,16,20]. One should of course be careful in comparing at the correct kinematics.
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 and in various models see e.g. refs. [21][22][23][24][25][26] and for the contributions from heavy quarks [27]. However, the naive OPE breaks down for the (g − 2) µ kinematics with q 4 → 0 [28]. The

JHEP10(2020)203
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. [28] (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 resulting OPE, originally formulated for baryon magnetic moment sum rules in refs. [29] and [30], and whose application to other hadronic (g −2) µ contributions was introduced in ref. [31], has the massless quark loop as the leading term and the di-quark magnetic susceptibility of the vacuum as the leading quark-mass suppressed contribution. In this work we extend the results of ref. [28] by computing the leading non-perturbative corrections not suppressed by masses. We also provide some more details about the calculations of ref. [28]. Our result should be useful to help reducing the error coming from the intermediate and short-distance regime [1].
In section 2 we briefly recapitulate the definitions of the four-point function of four electromagnetic currents, its decomposition in scalar functions and how it can be used for the muon g − 2 HLbL contribution. This follows the conventions of refs. [18,19]. It also gives the relation between the needed derivative of the four-point function and the three-point function in a constant field background that is used in the remainder of this paper.
In section 3 we give a complete description of the OPE in a constant background field and compare it to the usual vacuum OPE [32] and the one used in flavour-breaking transitions. We in addition comment on the physical meaning of the obtained matrix elements and build a renormalisation program. The renormalisation program is needed to systematically separate short-distance and long-distance effects while cancelling divergences. Both infrared and ultraviolet divergences are addressed. We also estimate the values of the matrix elements. The content of this section can in the future be used to obtain predictions for other Green functions in different kinematic regions.
Details on the calculation of the different non-perturbative pieces are provided in section 4, in particular we explain the different tools used, existing and newly developed, to obtain our analytic results and how the different infrared divergences systematically cancel. Finally, making use of the results of that section and the estimates of the matrix elements, in section 5 we calculate the numerical contribution of the different pieces for the (g − 2) µ . Final remarks and conclusions are made in section 6. Several intermediate derivations are relegated to the appendices as well as the full formulae.

The HLbL tensor
As can be seen in figure 1, the HLbL process involves a four-point correlation function of electromagnetic quark currents, i.e. Π µ 1 µ 2 µ 3 µ 4 (q 1 , q 2 , q 3 ) ≡ −i The currents are given by J µ (x) =q Q q γ µ q with the quark fields q = (u, d, s) and charge matrix Q q = diag(e q ) = diag(2/3, −1/3, −1/3). The tensor in (2.1) is the same as Π µνλσ in ref. [28], but with the internal notation slightly changed in the definition in order to make manifest some of its symmetry properties. This will be systematically exploited in the following sections. Moreover, the integral over q 4 is introduced to remove the δ-function of conservation of momentum, 1 i.e.

JHEP10(2020)203
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. [18,19] 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 [18,19] 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. [18,19]. In particular, the T µ 1 µ 2 µ 3 µ 4 i satisfy the same Ward identities as the HLbL tensor. As a consequence, in the static limit q 4 → 0 In the static limit q 4 → 0, the remaining tensor can thus 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 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. [18,19] and The exact definition of theΠ i functions is given in refs. [18,19]. 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.
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.

JHEP10(2020)203
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 appendix 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 appendix 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. [28], we focus on the kinematic region where the three loop momenta are large. As was shown in ref. [28], if one defines then in the static limit for the studied kinematic region, one can factor out the soft photon contributions according to 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 [28] lim 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 section 3.1 whereas the mixing of these operators is elaborated on in section 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 [34,35], which integrates out the heavy degrees of freedom leaving the low-energy (long-distance) dependence encoded in effective operators, in such a way JHEP10(2020)203 that the contributions from higher-dimension operators become suppressed by extra powers of the high-energy scale [36][37][38]. 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 [32], 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. [36]) 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 [29,30], 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 [28,31]. Higherdimensional operators are suppressed by extra powers of , providing a hierarchy of contributions with a systematic counting. Up to dimension 6 and order α s a basis of those operators is 3

4)
S 5, µν ≡qq e e q F µν , (3.5) We use here the notation of ref. [39]. In particular, 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 2 The words n-point function, correlation function and Green function are all used but have the same meaning. 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 [29], but it enters suppressed by gs, 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).

JHEP10(2020)203
multiplied by the corresponding quark electric charge. Note that with the conventions of ref. [39], 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 appendix B. In fact, only two combinations of four-quark operators contribute as discussed in section 4. The adopted notation is analogous to the one in ref. [36].
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. [39]), so that the resulting expression is of the form 5 (3.9) where X S contains the magnetic susceptibilities of the operators, 6 S i,µν = ee q X i S F µν . Even when this was a first step to achieve the separation between short-distance and long-distance 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 figure 2. In figure 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 figure 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. [40], tree-level operators in the notation of ref. [36]). Following a notation close to the one in ref. [36], 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 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 section 4. 6 We will define the susceptibilities more precisely later. 7 We need here the condensates induced by the external field. We refer to those as matrix elements to distinguish them from the usual (vacuum) condensates. involves two steps, calculating the relevant contributions leading to an expression at oneloop order of the formM 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 noncanonical part of their dimension, proportional to 2a : (3.14) Putting this back into (3.9) one finds  The contributions we calculate explicitly in section 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 The computation of the last three is needed, since they give contributions to the g 2 s Λ 4 QCD Q 4 through operator mixing as defined above and calculated in section 3.2. This is analogous to the mixing between bi-linear quark condensate and gluon condensate in the usual vacuum OPE [41][42][43].

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 figure 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Û 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 U 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

JHEP10(2020)203
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 figure 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 section 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 figure 4. Let us sketch it in 9 The procedure is equivalent to what is used in the usual vacuum OPE [40][41][42][43][44]. some detail to illustrate the procedure: where we have used that we are working in the radial gauge, A µ (x) = 1 2 x ν F νµ [39,41]. 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 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),

JHEP10(2020)203
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 2diagram in figure 4, leading toÛ (3.21) 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 [40][41][42][43].

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 figure 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 [32,[40][41][42][43], The last operator to treat is Q µν 7 which through the topologies shown in figure 4 only mixes with Q µν 6 . One has JHEP10(2020)203

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 non-perturbative 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. [45]. The gluon condensate is not so well-known, 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 [32].
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 [32].

JHEP10(2020)203
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. [46]. 12 At the lowest order the X M S i (µ) are directly proportional to the low-energy However, this gives no extra information by itself in the SU(3) V limit, since the Λ i,M S 1 (µ) are not known. 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 refs. [47][48][49][50]. First of all, in the large-N c limit the QCD spectrum is made of an infinite number of free, stable meson states [51][52][53]. The two-point functions in that limit are then saturated by the exchange of resonances. Owing to the quantum numbers of the studied two-point functions, the corresponding resonances must be vector mesons [54]. 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. [55] and adding a tensor source [46,50], 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 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 Since the symmetry transformation of Q µν 3 , Q µν 4 and Q µν 8 are identical, their Lagrangians are functionally equivalent and only the low-energy couplings are different.

JHEP10(2020)203
The obtained value for X 2 is, in fact, in very good agreement with a precise lattice determination [56,57]. 13 No precise modern evaluation of m 2 0 can be found in the literature. However, once again, some numerical estimates are available [58]. 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.20), 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 .

JHEP10(2020)203
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 appendix 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 (3. 44) 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 element qσ µν 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 nonzero 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 [59,60]. The Feynman integral reduction for the quark loop and the gluon matrix element contributions was done with Reduze 2 [61] and Kira [62]. 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 [29, JHEP10(2020)203 41], as shown explicitly in ref. [28]. 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 appendix 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) 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. Here, λ = λ(q 2 1 , q 2 2 , q 2 3 ) is the Källén function defined through (4.10) The different coefficients (c, d, f, g, h, p, q, r, s) can be found in appendix 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. [63]. 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 short-distance 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 figure 5.
There 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 [39,41] of the quark bilinears can be written as: (4.12)

JHEP10(2020)203
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. [28], 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. [39,41], one finds 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 16) 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 appendix 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 section 3 is rather involved. One needs to recursively exploit spinor algebra JHEP10(2020)203 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 With one derivative one has and with two they are Here, we have defined the linear combinations , (4.23) (4.24) For three derivatives there are two contributions. These are (4.26)

JHEP10(2020)203
Here,F µν ≡ i 2 µναβ F αβ and the A i are given by (4.28) (4.32) The operator enters from using the gluon equation of motion. Its susceptibility is defined as (4.34) 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ŝ (4.35) The numerical coefficients c m,n,p i,j,k can be found in appendix 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 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 figure 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 when taking the traces. Recalling the definition of X 8,2 in (3.42) we findΠ (4.37) A reduction of all possible four-quark matrix elements into a basis of 12 independent ones is given in appendix 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 figure 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 Tr λ a 2 (4.40) 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.
(4.42) where c, f, g and h are numerical coefficients given in appendix C.
When operator mixing is taken into account, the found divergences, which scale as 1 m 2 q and log Q 2 m 2 q , 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 section 4.3, we evaluate a HLbL µ for the matrix element as well as loop contributions. We use the matrix elements as estimated in section 3.3. For this purpose, we use the Cuba library [64], in particular the Vegas integrator building on Monte Carlo sampling the threedimensional 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 appendix C.2.

JHEP10(2020)203
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 figures 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.

Conclusions and prospects
The leading asymptotic behaviour of the HLbL for the (g − 2) µ kinematics has been confirmed to be given by the massless quark loop contribution. Although this had been commonly assumed previously, this is by no means obvious due to the static limit associated to the (g − 2) µ definition with the external photon leg at zero momentum. The main result of this work and of ref. [28] is to show how a proper short-distance expansion can be done in the limit of Q 2 1 ∼ Q 2 2 ∼ Q 2 3 Λ 2 QCD . In order to show that the quark loop is the first order of a well-defined expansion, the soft photon has been formulated as a long-distance, or, background, degree of freedom, following previous works of refs. [29,30]. We stress that using the vacuum OPE valid for HLbL when all the four Euclidean momenta are large, one would for the (g −2) µ kinematics have obtained a divergent expansion [28].

JHEP10(2020)203
Contribution Inputs (GeV units) Q min = 1 GeV Q min = 2 GeV X 1,0 1.73 · 10 −10 4.35 · 10 −11 X 1,m 2 −5.7 · 10 −14 −3.6 · 10 −15  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.  X 1,0 X 6 X 7 X 8,1 X 8,2 Figure 10. Numerical contributions from the X 6−8 pieces in absolute value using the inputs from table 1. Even when they are not suppressed by the quark mass size, they are found to be small compared with the massless quark loop contribution.

JHEP10(2020)203
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 non-perturbative 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.

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. Obtaining the same results with both projectors is one of the checks we did. Below we only present one of the sets of projectors.

JHEP10(2020)203 B Four-quark reduction
In this section we reduce the number of four-quark matrix elements from the basis q iAᾱ q jBβq kCγ 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

JHEP10(2020)203
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 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 JHEP10(2020)203 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Ô ijkl Γ ijB Γ klA is trivially related toÔ ijkl Γ ijA Γ klB . Taking advantage of this and requiring that the reduced matrix elements should be odd under charge conjugation one finds 1 2 (q i λ 1 q jqk λ 2 q l −q i λ 2 q jqk λ 1 q l ) = 1 64 γ µ ji γ ν lk −γ ν ji γ µ lk qλ 1 γ µ qqλ 2 γ ν q−qλ 2 γ µ qqλ 1 γ ν q 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 section 4.3.

C.1 The quark loop
Recall The two terms above are given bŷ

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,

JHEP10(2020)203
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)
For the topologies with one gauge boson, the only change with respect to (D.2) is an extra vertex in the quark chain, which can be allocated in two different positions plus the boson field itself (see figure 5b): (−1) p−q+1 pu ω 1 · · · u ωp (p+1)q!(p−q)! ×qD ν 1 · · · D νn D ω 1 · · · D ωq D D ω q+1 · · · D ωp D ν 1 · · · D ν m q, This simplification occurs for every Dirac structure Γ A and therefore also for their sum. This completes the needed derivation. We conjecture that the duality holds at all dimensions 16 and that its trivial generalization holds for any number of external legs, greatly simplifying calculations for this kind of topology.
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.