Effects of Longitudinal Short-Distance Constraints on the Hadronic Light-by-Light Contribution to the Muon $g-2$

We present a model-independent method to estimate the effects of short-distance constraints (SDCs) on the hadronic light-by-light contribution to the muon anomalous magnetic moment $a_\mu^\text{HLbL}$. The relevant loop integral is evaluated using multi-parameter families of interpolation functions, which satisfy by construction all constraints derived from general principles and smoothly connect the low-energy region with those where either two or all three independent photon virtualities become large. In agreement with other recent model-based analyses, we find that the SDCs and thus the infinite towers of heavy intermediate states that are responsible for saturating them have a rather small effect on $a_\mu^\text{HLbL}$. Taking as input the known ground-state pseudoscalar pole contributions, we obtain that the longitudinal SDCs increase $a_\mu^\text{HLbL}$ by $(9.1\pm 5.0) \times 10^{-11}$, where the isovector channel is responsible for $(2.6\pm 1.5) \times 10^{-11}$. More precise estimates can be obtained with our method as soon as further accurate, model-independent information about important low-energy contributions from hadronic states with masses up to 1-2 GeV become available.


Introduction
The persisting discrepancy between the Standard Model evaluation and the experimental determination [1] of the muon anomalous magnetic moment a µ is one of the outstanding open problems in particle physics and is traditionally considered a harbinger of New Physics. Moreover, the forthcoming results from the Fermilab E989 experiment, which aim to improve the present accuracy by a factor of 4 to reach an uncertainty of about 16×10 −11 (i.e. 0.14 ppm) [2], make it even more crucial and timely to further scrutinize and improve control over theory predictions.
In this framework a HLbL µ is evaluated via a two-loop integral of dispersively reconstructed scalar functions against analytically known kernels. At sufficiently small space-like photon virtualities, contributions from low-mass states accessible to a dispersive treatment are enhanced. At higher virtualities such an enhancement does not occur leading to important effects from higher intermediate states, which are constrained by operator product expansions (OPEs) and perturbative QCD (pQCD).
More specifically, there are two relevant kinematic regimes concerning short-distance constraints (SDCs) on a HLbL µ for asymptotic values of (subsets of) the photon virtualities. Since one of the photons corresponds to the static electromagnetic source in the definition of g − 2, one asymptotic regime is realized when the remaining three space-like photon virtualities are comparable and much larger than Λ 2 QCD , and the other when two space-like photon virtualities are much larger than both the third and Λ 2 QCD . The latter SDC was first derived by Melnikov and Vainshtein (MV) [25] using an OPE that leads to relations involving longitudinal and transversal amplitudes of the correlator of two vector and one axial current (VVA) in the chiral limit. The former SDC was also discussed in Ref. [25] based on the quark-loop calculation at leading order in pQCD and its derivation was recently put on a firmer theoretical ground by means of an OPE in an external electromagnetic background field [26].
Tree-level resonance exchanges cannot make a HLbL µ comply with all SDCs unless an infinite number of states is included. This is due to the fact that the transition form factors (TFFs) describing the resonance couplings to off-shell photons are subject themselves to asymptotic QCD constraints [27][28][29], which make the full HLbL four-point function decay too fast at high virtualities. 1 MV proposed a model to satisfy the longitudinal and transversal OPE SDCs through a modification of the pion pole contribution [25,31], which affects also the low-energy region. Recently, alternative model-dependent solutions have been investigated to fulfill both OPE and pQCD SDCs by instead adding degrees of freedom to the ground-state pseudoscalars. In this context, Refs. [32,33] proposed the inclusion of infinite towers of excited pseudoscalar poles in large-N c -inspired Regge models to satisfy longitudinal SDCs away from the chiral limit, 2 while the effect of summing over axial-vector contributions in holographic QCD was the subject of Refs. [40,41]. 3 Through the explicit summation of intermediate states, these models provide specific interpolations between the low-energy region and the asymptotic regimes for the scalar functions that determine a HLbL µ . The goal of this paper is complementary to these studies. We introduce an approach based on more general interpolating scalar functions, independent of the physical mechanism that is ultimately responsible for their actual form outside the low-energy region. The multi-parameter families of functions studied here satisfy all constraints rigorously derived from general principles: unitarity, analyticity and crossing in the low-energy domain, OPE and pQCD constraints in the mixed and high-energy regions. Here we focus on longitudinal SDCs since these are tightly related to the pseudoscalar poles for which accurate low-energy input is available and their implementation does not involve any mixing of OPE constraints among different scalar functions [33]. Error estimates as well as the role played by the various parameters and assumptions, can be easily and transparently addressed in our approach and are investigated in detail in our numerical study.
Crucial input for our analysis is provided by an accurate low-energy representation of the scalar functions. In the following we will mostly assume that this is given by the ground-state pseudoscalar poles. In this context, an important role is played by the lightest state (π 0 ), whose contribution is under firm theoretical control thanks to a dispersive evaluation [12][13][14]. Improved determinations of the effects of SDCs can be obtained in a straightforward way within our approach once similarly precise, model-independent information about further relevant intermediate states in the energy region up to 1-2 GeV become available. In order to illustrate this aspect and compare against a different way to estimate the contribution from SDCs, we have applied our method also to the case where the lightest axial-vector meson is included in the low-energy region using input from holographic QCD [40,41] and neglecting issues related to intrinsic model dependence.
The paper is structured as follows. In Sec. 2 we review the relevant constraints on HLbL and the assumptions made in their derivations. Sec. 3 describes our interpolation between the OPE and pQCD asymptotic constraints while in Sec. 4 we discuss its smooth connection with the low-energy region. In sec. 5 we present our numerical analysis with particular emphasis on the error estimation. Conclusions are drawn in Sec. 6. App. A is devoted to the analysis of the convergence properties of our interpolants.
2 Longitudinal short-distance constraints on HLbL 2.1 Master formula for a HLbL µ and pseudoscalar pole contributions In order to set up the notation, we start by summarizing the relevant definitions and results from Refs. [9,11]. The HLbL contribution to a µ is governed by the fourth-rank vacuum polarization tensor for fully off-shell photon-photon scattering in pure QCD, 2 For other calculations based on large-Nc arguments to satisfy long-and short-distance constraints on QCD correlators using finite or infinite sets of narrow resonances, see Refs. [30,[34][35][36][37][38][39].
3 See also Ref. [42] for a discussion of the role of axial-vector mesons in the saturation of the SDCs.
with momenta assigned as q 1 + q 2 + q 3 = q 4 . In this expression, the electromagnetic current for the light quark triplet is given by By generalizing the procedure introduced by Bardeen and Tung [43] and Tarrach [44] in their studies of doubly-virtual Compton scattering, it is possible to derive a generating redundant "BTT" set of 54 Lorentz structures, which is manifestly gauge invariant, closed with respect to crossing relations and such that the scalar functions Π i are free of kinematic singularities. The HLbL contribution to a µ can be derived from the tensor Π µνλσ by using projection operator methods [5,45,46] in the static limit q 4 → 0. After performing a Wick rotation to Euclidean momenta, angular averages [47,48] lead to the master formula [11] a HLbL where Q 1,2 denote the magnitudes of the Euclidean loop four-momenta, Q 1,2 = −q 2 1,2 , and τ is the cosine of the angle between these vectors. The scalar functionsΠ i are linear combinations of the previous Π i for q 4 → 0. The analytic expressions of the integration kernels T i are given in Ref. [11].
Parameterizing the three-dimensional integration domain by the coordinates [49] Σ which are related to the non-vanishing photon virtualities by will prove very useful in the following discussion about asymptotic constraints on HLbL. The master formula in Eq. (4) then takes the form In terms of the Q 2 i coordinates, the integration domain amounts to a cone with tip at the origin. In terms of (Σ, r, φ), a given point in this cone is specified by the distance Σ to the tip of the point's projection on the symmetry axis (Σ = Q 2 1 + Q 2 2 + Q 2 3 ), and by the polar coordinates r and φ on the plane containing the point and orthogonal to the symmetry axis, normalized such that r = 1 corresponds to the surface of the cone.
In the master formula, a special role is played by the scalar functionsΠ 1,2 , which fulfill where the crossing operator C i,j exchanges momenta and Lorentz indices of the photons i and j. These functions are the only ones describing the effects of pseudoscalar tree-level exchanges. For small values of Σ, the pion pole dominates yielding the largest contribution to a HLbL µ and also η/η poles yield sizable effects. Furthermore, distinctively, the OPE SDCs onΠ 1,2 do not involve other scalar functions [33].
The functional form ofΠ 1,2 in specific kinematic regimes is constrained according to analytic QCD results, which we will fully exploit to estimate the impact of (longitudinal) SDCs on where the shift in the variable φ in T 2 corresponds to the crossing operation onΠ 1 . Thus for our analysis, we only need to study one BTT scalar function in the g − 2 kinematics. For the purpose of later discussion, we stress here that a pole term inΠ 1 due to a single-particle intermediate state of mass M yielding the denominator Q 2 3 + M 2 leads to a hierarchy among contributions in the space-like momentum region relevant for a HLbL µ . For small values of Q 2 3 larger masses get suppressed, while for Q 2 3 comparable to the squared mass of the heavier state or larger, no suppression is expected. 4 This effect is of course modified by the numerator inΠ 1 , which encodes information on the strength of the coupling to two (off-shell) photons, but it still helps us identify which states can be relevant at specific energy scales and which not, independent of the values of Q 2 1,2 . The lightest state contributing toΠ 1 is π 0 . The unitarity relation for a single pseudoscalar intermediate state yields where the numerator is given by the product of a doubly-virtual and a singly-virtual transition form factor (TFFs) for an on-shell pseudoscalar meson (PS), which is defined by the matrix element with 0123 = +1. If we setΠ 1 =Π PS-pole 1 , then a long µ amounts to the pseudoscalar pole contribution a PS-pole µ . In the π 0 case, this has been evaluated within a few percent accuracy via a data-driven dispersive approach [12][13][14], This result agrees with other recent determinations based on lattice QCD [50], Canterbury approximants [51], Dyson-Schwinger equations [52,53] and AdS/QCD models [54]. While and Q 2 1 = Q 2 3 , respectively. The colored regions denote where SDCs onΠ 1 hold at large Σ. The blue domains yield contributions toΠ 1 from the OPE expansion of the VVA correlator that are sub-leading compared to the green one, while the orange region corresponds to the pQCD constraint. a dispersive analysis of the doubly-virtual η/η TFFs has not been completed yet, 5 the method of Canterbury approximants in Ref. [51] provides data-driven determinations and associated uncertainties also for the η/η TFFs. In our numerical analysis of SDCs, we have used as input the dispersive π 0 TFF from Refs. [13,14] and compared our final results against those with form factors from Canterbury and Dyson-Schwinger approaches, while for η/η we have used the TFFs in Ref. [51] and compared against Ref. [52].
The asymptotic constraints on a HLbL µ [25,26] (see also Refs. [42,58]) that we are going to discuss in the next sections have been translated into the BTT framework in Refs. [26,32,33]. In this context, there are two distinct relevant kinematic regimes. The first (asymmetric) one is realized when one of the photon virtualities is much smaller than the other two, which are large and comparable, e.g. Q 2 1 ∼ Q 2 2 Q 2 3 . The second (symmetric) limit occurs when all the Euclidean non-vanishing photon virtualities are large and comparable in size . Both asymptotic limits correspond to Σ → ∞ but for different values of r and φ: the asymmetric limit Q 2 1 ∼ Q 2 2 Q 2 3 corresponds to r = 1 and φ = π while the symmetric configuration occurs in a neighborhood of r = 0 (see Fig. 1).
In the following we will review the relevant constraints onΠ 1 at large Σ and describe in detail our method to provide general families of interpolants forΠ 1 (Σ, r, φ) between lowand high-energy regions in the g − 2 integral.

The asymmetric asymptotic region: OPE constraints
For large Euclidean values ofq ≡ (q 1 − q 2 )/2, one can expand the time-ordered product of two electromagnetic currents, which defines the tensor into a series of local operators. At leading order in α s , by matching connected single-quark matrix elements, one obtains for q 1 + q 2 = 0 (see e.g. Refs. [25,33]) where the axial current j µ 5 is defined by j µ 5 (x) =ψ(x)Q 2 γ µ γ 5 ψ(x) with charge matrix given in Eq. (2). The ellipsis denotes sub-leading terms suppressed by powers of {|q 1 + q 2 |/|q|, Λ QCD /|q|}. This result implies that, at leading order in the OPE and at leading order in α s , the HLbL tensor can be expressed in terms of the correlator of two vector currents with an axial current, This three-point function, which also appears in the calculation of fermion loop electroweak contributions to a µ [59,60], can be decomposed into Lorentz structures that are longitudinal and transversal with respect to the Lorentz index of the axial current (see e.g. Ref. [61]). The corresponding longitudinal scalar function determines the asymptotic behavior ofΠ 1 in the asymmetric region and is fixed by the axial Adler-Bell-Jackiw anomaly. For massless quarks, this translates into the following constraint [33] for the singlet and octet flavor components ofΠ 1 , defined by the decomposition of the axial current, where C a = Tr(Q 2 λ a )/2 in terms of the charge matrix Q and Gell-Mann matrices λ a . In particular, Eq. (16) holds in the kinematic limit Q 2 For the non-singlet components (a = 3, 8), since perturbative [62] as well as non-perturbative [63] corrections are absent in the chiral limit, the hierarchy between Λ 2 QCD and Q 2 3 can be dropped. This can be done only at leading order in the large-N c expansion in the singlet case, which is affected by the gluonic U (1) anomaly. Furthermore, in the crossed kinematics ( , the leading-order OPE contributions toΠ 1 vanish. Let us now compareΠ (a), OPE 1 against the pseudoscalar pole contributions, focusing on the pion pole first. In the chiral limit and using the fact that lim 3 F π at leading order in α s [64,65], one finds This expression has a pole at Q 2 3 = 0 since F πγ * γ * (0, 0) = 3C 3 /(2π 2 F π ). The location of this pole as well as its residue agree withΠ (3), OPE 1 in Eq. (16), which is consistent with the pion being the only massless isovector state in the chiral limit.
For finite quark masses, the pole inΠ π 0 -pole 1 is shifted from Q 2 3 = 0 to Q 2 3 = −m 2 π , which lies outside the integration domain for a HLbL µ . The closest point in the integration region for fixed asymptotic Σ (see Fig. 1 This is still close to the actual pole, which leads to the enhancement by m −2 π . Since no other contribution receives the same enhancement, the last expression is expected to provide an excellent approximation to the trueΠ (3) 1 in the specified limit. 6 We observe that the OPE result, which is derived in the chiral limit, reproduces Eq. (19) if the pole position is shifted by the pion mass as dictated by the pion pole contribution This is also consistent with the OPE result in Eq.
where chiral corrections are sub-leading. We extend the last relation to the η/η channels and write Here C π = C 3 but C η/η cannot be directly identified with C 0/8 due to η-η -mixing. In analogy to the pion channel, we assume that ground-state single-pseudoscalar exchanges dominateΠ , despite the fact that the η/η poles are further away from Q 2 3 = 0. This assumption implies that C η/η can be read off from the pole contributions One can show that in the chiral limit and at leading order in α s [33] lim At this point we note that, besides the α s corrections to the TFFs and OPE coefficient discussed in Sec. 2.3, which affect all ground-state pseudoscalars in the same way, the U (1) anomaly induces a running of the flavor singlet decay constant [66][67][68]. This running leads to an incomplete cancellation between the decay constants in the symmetric asymptotic and the real photon limits, which has a sizable impact due to the large scale separation [69,70]. SinceΠ can be expressed in terms of TFFs according to Eq. (10), assuming that corrections due to non-vanishing meson masses are negligible both in the real photon limit and in the symmetric asymptotic limit of the η/η TFFs, Eqs. (21-23) together imply 6 At variance with the MV model of Ref. [25], we do not neglect the momentum dependence of the singly-virtual TFF and we allow for the contribution from other states besides the pion at finite Q 2 3 .
up to the above-mentioned anomaly-induced scale-dependence, which leads to a violation of this equality (cf. Sec. 5.2). For Q 2 3 Λ 2 QCD , the additional Q 2 3 -suppression of the singly-virtual TFF leads to a mismatch between the pseudoscalar pole contributions and the OPE constraint. In Ref. [25] MV proposed to solve this issue by setting the singly-virtual TFF equal to a constant. This prescription is not compatible with the dispersive definition of the pole contributions in the framework summarized in Sec. 2.1, according to which, instead, an infinite tower of heavier intermediate states is needed to saturate the constraint (see e.g. Ref. [33]). For this purpose, summations of series of contributions from excited pseudoscalars [32,33] and axials [40,41] have been recently performed in the context of hadronic models. In Secs. 5.3 and 5.4, we will compare the outcome of our analysis against these estimates of the effects of longitudinal SDCs.

α s corrections to the OPE
The derivation of Eq. (14) has been performed at leading order in α s . Since no other operator of dimension 3 can appear in the OPE of two vector currents, α s corrections only affect the OPE coefficient of the axial-vector current. Thus, the next-to-leading order (NLO) version of Eq. (14) is where C is a purely numerical factor. C can be easily computed by observing that the two-current operator not only enters the HLbL tensor, but also the pion TFF (see Eq. (11)). Thus, any perturbative correction to the OPE Wilson coefficient automatically implies the same perturbative correction to the symmetric limit of the pion TFF and vice versa. Therefore, C can be determined from the pion TFF asymptotics, which has been calculated to NLO in Refs. [14,71] From this, we read off C = −1. 8 It follows that the NLO version of Eq. (21) reads

The symmetric asymptotic limit: perturbative QCD constraints
In Ref. [26] it has been shown that the pQCD quark loop is the leading term of an OPE in the kinematic limit where the fourth (external) photon has vanishing momentum in (g − 2)-kinematics. At leading order in this OPE and at leading order in α s [33], (28) In the symmetric limit, neglecting terms that are suppressed by powers of m 2 q /Q 2 , where we have chosen to adopt the same flavor decomposition as for the asymmetric OPE case, Eq. (16). If higher-order perturbative corrections are small, the leading-order result above is expected to be a good approximation also away from the fully symmetric configuration as long as large logarithms of ratios of momenta are absent. SinceΠ PS-pole 1 decays like Q −6 , (towers of) hadronic contributions beyond ground-state pseudoscalar poles have to be responsible for the behavior shown by Eq. (29). Following the MV prescription in Ref. [25], the parametric dependence on Q can be reproduced but with an incorrect coefficient.
In order to saturate the pQCD result in the isosinglet channels, we need coefficients C pQCD Since Eq. (24) is violated, we define where the parameter δ 0 is chosen such that Eq. (30) holds.

Interpolating between asymptotic constraints
We approximate the trueΠ 1 (Σ, r, φ) following a two-step procedure. We first select functional forms that are valid for asymptotic Σ and are compatible with the constraints discussed in the previous section. We then interpolate between this set of functions and various representations ofΠ 1 at small Σ determined by single-particle intermediate states.
Here we work at leading order in α s . Perturbative corrections will be discussed in our numerical analysis in Sec. 5. The relevant constraints onΠ 1 at large Σ are given by Eq.
if C PS = C pQCD PS . Thus, Eq. (32) interpolates between symmetric and asymmetric asymptotic limits. According to Sec. 2.4, δ 0 parameterizes the anomaly corrections to the singlet VVA correlator and the resulting shift in C pQCD PS with respect to C PS . Since a term proportional to Σ −2 and independent of (r, φ) does not change the leading behavior at Q 2 3 = 0 and thus does not spoil compatibility with the OPE constraint, we subtract 36δ 0 C 2 η/η /(π 2 Σ 2 ) from Eq. (32) in the case of η/η .
Obviously, the choice made in Eq. (32) and the exact form of the singlet correction are not unique and we are free to add a generic function such that the interpolant still satisfies the constraints. In order to have a non-negligible effect at asymptotic values of Σ, this additional function should also scale as Σ −2 and we demand it to be finite and analytic for all r ≤ 1. 9 Therefore it can be approximated by a Taylor series in r cos φ and r sin φ truncated after order M , where for integer j, due to the pQCD constraint and crossing symmetry.
Up to now, we have applied the quark-loop result only at r = 0. However, the fact that Eq. (28) holds also in a neighborhood of this point can be used to fix the coefficients a i,j . To this end, we fitted Eq. (33) at fixed asymptotic Σ with M = 2 to Eq. (28) for r < 0.9. 10 We chose a grid of equally separated points in this fitting region and minimized the sum of the relative squared differences between our interpolant and the leading-order quark-loop expression. The resulting 5 dimensionless fit parameters read a 1,0 = −0.170 , a 2,0 = 0.094 , a 0,2 = −0.554 , a 1,2 = −0.169 , a 2,2 = −0.756 (35) and are all at most O(1), as expected since in Eq. (33) they parameterize relative corrections. In Sec. 5.1.2 we will discuss uncertainties due to the chosen fitting range, the number of parameters in the fit and α s corrections to the asymptotic constraints.

Interpolating between low and high energies
The next step is to smoothly connect our representation ofΠ 1 for Σ {Λ 2 QCD , M 2 PS , . . . } given by Eq. (33) to an accurate low-energy description. We achieve this by adding suitable 9Π 1 cannot decay more slowly than Σ −2 for any (r, φ) region in order for the aµ integral in Eq. (9) to be finite. 10 Since the maximal ratio of two squared momenta for r = 0.9 is 14.7 and ln(14.7) ≈ 2.7, large logarithms do not occur in this region.
terms toΠ (a), asymp 1 that are sub-leading at large Σ. For each choice of r and φ, the coefficients of these terms are then matched onto an input low-energy representation ofΠ 1 at a suitably defined surface Σ match (r, φ). In Sec. 4.2 we will discuss how Σ match is related to the mass scale at which intermediate states beyond the ones explicitly considered start to affectΠ 1 .

Interpolation functions and matching procedure
For Σ > Σ match we consider the following two interpolation functions Π (a), int 1 1 whose leading terms at asymptotic Σ are given in Eq. (33), whereas below the matching surface we setΠ In App. A we will show that these functions converge to the trueΠ (a) 1 in the limit N → ∞ when matched to exact low-energy input using the convergence property of a Taylor series. The two different forms given in Eq. (36) will be used to estimate the sensitivity of our numerical results on the specific choice of interpolation between low and high energies. 11 The coefficients b i (r, φ) are fixed from the requirement that theΠ (a), int i 1 have the same value and the same N − 1 Σ-derivatives as the low-energy representation if evaluated at Σ = Σ match (r, φ) for each (r, φ). No expansion is performed in r and φ, which is crucial to obtain a smooth transition to the low-energy regime.
Determining the optimal value of N is a non-trivial issue. On the one hand, larger values of N seem to be preferable since the trueΠ 1 is analytic for space-like momenta and thus all derivatives are continuous. On the other hand, matching many derivatives leads to a function that is almost saturated by the low-energy input contribution up to considerably higher energies than Σ match (r, φ). Since it is desirable to match at least one derivative in order to haveΠ 1 differentiable at the matching point, we will use N ∈ {2, 3} in order to estimate the dependence on N .
Interpolation functions with a logarithmic dependence on Σ are not forbidden. This can stem, for example, from non-perturbative corrections leading to terms like ln (Q 2 i /M 2 ), where M is some non-perturbative mass scale. In fact, the Regge model considered in Refs. [32,33] leads to interpolants containing terms like Q −4 ln (Q 2 /σ 2 ) for Q 2 i = Q 2 → ∞, where σ 2 could e.g. be the Regge slope of the excited pseudoscalar masses. In order to allow for such a logarithmic approach of the asymptotic expression, we additionally consider the alternative interpolant Π (a), int 3 1 (37) and use again N ∈ {2, 3}.

The matching surface Σ match
The remaining crucial ingredient in our procedure is the function Σ match (r, φ), which determines the value of Σ at which the matching is performed for given r and φ. Choosing it too low leads to important modifications ofΠ For small values of Q 2 3 , the π 0 , η, η poles are assumed to dominate independently of Q 2 1,2 , due to the pole at Q 2 3 = −m 2 PS (see Secs. 2.1 and 2.2). This implies that no matching is needed in this regime, i.e. Σ match (1, π) = ∞. The most general function that is analytic for all (r, φ) except for a (first-order) pole at (r, φ) = (1, π) can be written as where m 2 determines the matching scale at r = 0 and P is a polynomial with two arguments and no constant term. The transformation property ofΠ 1 under crossing specified in Eq. (8) restricts P to contain only even powers of r sin φ.
The parameter m 2 sets the absolute mass scale of Σ match and should thus be related to the masses of the states affectingΠ 1 beyond the ones explicitly included, namely π 0 , η, η here. In the following, we will assume that contributions toΠ 1 in the g − 2 kinematics stemming from multi-particle intermediate states are dominated by narrow resonances while non-resonant effects lead to negligible corrections to the matching procedure and can be simply added to our final results. 13 This is realized for example in the large-N c limit of pure QCD: since the short-distance expressions forΠ 1 in both symmetric and asymmetric limits scale like N c , these can indeed be saturated by single-meson exchanges (see Ref. [73]). Non-resonant contributions from multi-hadron intermediate states (like 2π, 2K, πη, 3π, . . . ) are sub-leading for large N c and thus cannot contribute to the SDCs. Since scalar mesons have no impact onΠ 1 , the lightest states beyond the ground-state pseudoscalars that are the most relevant at small Q 2 3 (see Sec. 2.1) are the axial mesons like a 1 (1260) and the tensor mesons like f 2 (1270), with masses in the 1-2 GeV region, whose effects on a HLbL µ can presently be estimated only using hadronic models. For P (x, y) = 0, Σ = Σ match corresponds to Q 2 3 = m 2 . Since a state of mass M ceases to be suppressed by the denominator (Q 2 3 + M 2 ) compared to lighter states when Q 2 3 approaches M 2 , m 2 should be chosen well below M 2 . At the same time, it should not be taken too small, because we do not expect any large contribution toΠ 1 at Q 2 3 M 2 . We thus regard m 2 = 0.5 GeV 2 as a good starting point for our analysis. In Sec. 5.1.4 we will discuss a range of choices for this parameter as well as the effects of the polynomial which we have estimated by means of a Monte Carlo sampling over the coefficients p i,j . 12 We denote by a 1 . 13 We have checked the effects of the inclusion of the pion-loop contribution toΠ1 [11] in the low-energy representation. Since the two-pion state contains a five-dimensional representation of the isospin group, a full decomposition intoΠ is not possible. However, even if its complete contribution is added to the isovector channel, we find that at the current level of accuracy it is irrelevant whether the pion loop is included in the matching procedure or not.

The isovector channel
The isovector channel is best suited to our method since it is characterized by a large contribution from the low-energy region dominated by the well-known pion pole, which does not mix (strongly) into the other flavor channels. The lightest one-particle intermediate state beyond the π 0 in this channel is the a 1 (1260), whose effect at low energies is suppressed by the large mass gap. The numerical dominance of this channel at low energies, however, does not imply that the same holds true at intermediate and high energies. In fact, the values of C a in Eq. (17) make the flavor singlet channel the numerically most important one in the asymptotic region where meson masses can be neglected, i.e. for Q 2 i Λ 2 QCD . In Sec. 5.2 we will discuss the inclusion of η/η and in Sec. 5.4 also the case of the isovector ground-state axial, which is however affected by a larger degree of model dependence.
We start by selecting a "reference" set of assumptions and input parameters. The impact of their modifications will be assessed in the next sections and will define the range of our predictions in the form of an uncertainty band. This procedure allows us also to examine how the estimate of the effects of SDCs would be improved by more precise information on the pion pole, the contributions from states with masses around 1 GeV and the asymptotic regime.
As low-energy reference input, we took the leading-order dispersive π 0 singly-and doubly-virtual TFFs [13,14], while the corresponding O(α s ) correction is included in the uncertainty. As reference interpolating function, we usedΠ , is shown in Fig. 2 for r = 0 together with the uncertainty band for the interpolants that we are going to discuss in the next sections. Our reference outcome for the contribution to a long µ due to the longitudinal SDCs in the isovector channel is where a in the master formula Eq. (9), and a π 0 -pole µ,disp is given in Eq. (12) according to Refs. [13,14]. In Sec. 5.1.5 we will argue that our final result does not strongly depend on the choice of the reference set of parameters.

Pion TFF uncertainties
We shall now describe the effects of modifying the different ingredients of the reference configuration, one by one, starting from the pion TFF. By propagating the errors quoted in Refs. [13,14] for the dispersive determination of the pion TFF and by summing the different contributions in quadrature, taking as well into account that a modification of the TFF affects both terms in Eq. (40), we obtained an asymmetric error band around the reference result with boundary values which correspond to the asymmetric error for the dispersive π 0 TFF. Given the smallness of these uncertainties, the (negative) correlation between them and the uncertainties of a π 0 -pole µ,disp can be safely neglected.
In order to study the impact of different pion TFF parameterizations, we compared the previous results against the ones obtained using, both for the construction of the interpolant and the evaluation of a π 0 -pole µ , the C 1 2 Canterbury approximant with a π;1,1 = 2b 2 π of Ref. [51] and the Dyson-Schwinger TFF from Ref. [52]. We obtained ∆a which are both compatible with the reference result within the range given above. We conclude that the outcome of our analysis is very robust against changes in the TFF input and that the present knowledge of the pion TFF is sufficient for our purposes.

Asymptotic uncertainties
Here we focus on the uncertainties inΠ • α s corrections to the OPE constraint as given by Eq. (27); • α s corrections to the quark loop.
We start by discussing the fit to the quark-loop result. In Sec. 3, we chose M = 2, which leads to a strongly improved fit quality compared to M = 1. Considering a larger value of M gives an estimate of the errors made by approximating the pQCD result by a polynomial at r < r max at fixed asymptotic Σ and by extrapolating to the regime r > r max , which is unknown except for the OPE constraint. Choosing M = 3 shifts the result for ∆a (3) µ by only 0.02 × 10 −11 indicating that the truncation at M = 2 is sufficient. We also studied the effects of a substantial reduction of the radius, namely from r max = 0.9 down to r max = 0.5, where no logarithm of ratios of squared momenta is larger than 1. We found that this leads to a small shift (0.07 × 10 −11 ). We did not consider r max > 0.9 since fixed-order pQCD is not expected to converge for r close to 1 due to large logarithms. Combining linearly the uncertainties from the choice of M and the fitting radius gives with respect to the reference contribution of longitudinal SDCs to the pion pole input in Eq. (40).
Let us now focus on the estimate of the separate perturbative corrections to either the OPE or the pQCD result. Since those concerning the OPE should not be extrapolated into the domain of validity of pQCD, for asymptotic Σ we write (cf. Eq. (27) Here the Heaviside step function θ ensures that the perturbative correction only affects a region around Q 2 3 = 0, whose size can be varied via the free parameter A. By setting A = 1/29, this region does not extend into the r < 0.9 domain. The choice of the renormalization scale µ 2 is the same as in Ref. [14]. In our numerical analysis, for the running of α s we used the three-flavor one-loop beta function and matched to α s (µ 2 = M 2 τ ) = 0.35. According to the discussion in Sec. 2, the OPE constraint in the chiral limit is saturated by the pion pole at Q 2 3 = 0 to all orders in perturbation theory. For this reason in a consistent analysis the OPE coefficient and the pion TFF in the symmetric limit should be taken at the same perturbative accuracy. Hence we replaced in Eq. (36)Π  to the pion-pole contribution with TFFs including O(α s ) effects [14,71]. 14 Using this interpolant and this pion-pole result, our outcome for ∆a (3) µ is larger than the reference result Eq. (40) by For A = 1/3, the domain where the correction applies extends down to r = 0.25, but nevertheless the shift of ∆a (3) µ turns out to be −0.05 × 10 −11 and thus still negligible. The smallness of these shifts can be understood from the large values of Σ match in the region where these perturbative corrections apply. For this reason, the effect is almost completely included in the pion pole contribution, where it also has a small impact [14].
Since the NLO calculation of the quark loop has not been performed yet, we can only provide a rough estimate of the uncertainty related to unknown O(α s ) corrections. We assumed in analogy with Eq. (44), and as in the leading-order quark loop fit, we set r max = 0.9. Using this expression in Eq. (36) for the matching to the pion-pole with leading-order dispersive TFF, we obtained a shift of −0.18 × 10 −11 compared to the reference result. Even when inflating this uncertainty by a factor of 2, δ NLO pQCD ∆a (3) µ = 0.36 × 10 −11 (47) this effect is still sufficiently small compared to the current precision goal. We stress that once NLO calculations become available,Π should be constructed to analytically interpolate between the NLO expressions for the OPE and the quark loop. The discontinuous functions employed here only serve to provide a ballpark estimate of NLO effects.

Choice of interpolation functions
In Eqs. (36) and (37) we have introduced three different interpolation functions, characterized by two or three free parameters to be matched to the low-energy representation. The corresponding results for the contribution from longitudinal SCDs are where ∆a given by Eq. (40) has been included for completeness. We observe that the slower logarithmic approach to the asymptotic limits in the interpolant 3 leads to smaller results, especially when compared to the similar interpolant 1.
all values listed above are within the range ∆a

Choice of Σ match (r, φ)
The function Σ match (r, φ) in Eq. (38) contains the mass parameter m and the polynomial P (x, y), which has been set equal to zero so far. We have argued in Sec. 4.2 that m should be chosen considerably smaller than the mass M of the lightest resonances contributing tō Π 1 in addition to the ground-state pseudoscalar mesons. For this reason, for the reference interpolant we set m 2 = 0.5 GeV 2 . Here we discuss the effects of alternative choices for this parameter within a range between m min and m max . Since according to Sec. 4.2 a conservative choice for the upper end of the range is m max M , we set m 2 max = 1 GeV 2 . In order to determine an appropriate value for m min , one has to estimate isovector contributions beyond the π 0 -pole. Following our argument in Sec. 4.2, it is sufficient to restrict ourselves to single-particle intermediate states and focus on the one giving the largest effect at energies around m. We assumed this to be given by the pseudoscalar π(1300) for the following reasons. Models for tensor mesons around 1 GeV give similar or smaller contributions to a HLbL µ [16,74,75]. For ground-state axials, recent studies based on different approximations and hadronic models yield quite different numerical results, see e.g. Refs. [40,41,74,76], leading to large uncertainties. If future model-independent analyses show that axial-meson exchanges are responsible for significant effects inΠ (3) 1 also at relatively small momenta, then these contributions should be added to the pion pole before the matching is performed since our procedure relies on a sufficiently precise knowledge ofΠ min . The green band shows the sum of the π 0 -and π(1300)-pole contributions, where the latter has been calculated using input from RχT and phenomenology, including errors. dependence, in Sec. 5.4 we will discuss the inclusion in our procedure of information from holographic QCD on the lightest axial meson.
As for the light pseudoscalars, the interaction of π(1300) with two photons can be described by a TFF, which determines the contribution toΠ 1 as in Eq. (10). In our analysis we used as input the π(1300) TFF derived in Ref. [77] in the framework of Resonance Chiral Theory (RχT) [78]. We fixed the free parameters in Eq. (68) of Ref. [77] by requiring that (i) the π 0 -TFF satisfies the Brodsky-Lepage condition, i.e. F π 0 γ * γ * (−Q 2 , 0) = 2F π /Q 2 + O(Q −4 ) and (ii) the two-real-photon limit of the excited pion TFF is in the range F π(1300)γ * γ * (0, 0) ∈ [0, 0.0544]GeV −1 , argued for in Ref. [33] based on experimental results [79,80]. The upper boundary of this interval leads to the most conservative error estimate in our analysis, and is used in the following.
Our procedure to determine m min can be illustrated by means of Fig. 3. Given a value of m min , at large enough Σ, the range of interpolants (orange band) spanned by m ∈ [m min , m max ] safely includes the green band representing the sum of the contributions from π 0 and π(1300), including errors on the latter due to the range for F π(1300)γ * γ * (0, 0). Since this is not the case at small Σ, the contribution to a long µ from this region is underestimated in our approach. To gauge this effect, we calculated the integral in Eq. (9) with using the maximal π(1300) contribution and restricting the Σ-domain to the region below the point where the bands start to fully overlap (as a function of r and φ). This integral gives the missed contribution a  In the following we will use δ + m,1 ∆a µ as an alternative, even more conservative uncertainty.
We also considered a different parameterization for the π(1300) TFF, namely the one given by the Regge model in Refs. [32,33]. Using the empirical m π(1300) = 1.30 GeV instead of the Regge-model value of 1.36 GeV used in these references and following the same procedure discussed above, we obtained m 2 min, 1 = 0.53 GeV 2 and m 2 min, 2 = 0.20 GeV 2 . This leads to where we again added a missed µ to the uncertainties in the upward direction. For our final result we use the more conservative uncertainty estimates given in Eq. (51).
In order to study the effects of the polynomial in Σ match (r, φ), Eq. (39), we set M = 2 and sampled the free parameters according to a standard normal distribution. Since the pion pole gives an excellent approximation ofΠ 1 for very small Σ at any (r, φ), we only allowed for parameters giving Σ match (r, φ) > Σ t for all (r, φ), where Σ t is defined as the smallest value of Σ such thatΠ π(1300)-pole 1 (Σ, r, φ) holds for some (r, φ). WithΠ π(1300)-pole 1 calculated using RχT, we obtained Σ t = 0.57 GeV 2 . This condition ensures that there are no large contributions from our interpolation at points where RχT predicts a very small excited pion contribution. From this we calculated a distribution of results for ∆a (3) µ , which features a Gaussian-like peak close to the reference result and asymmetric tails, and read off the 16 % quantiles from both sides corresponding to the 1σ errors for a Gaussian. This gives We have checked that this result is stable against the inclusion of terms of order 3 in the polynomial and moderate changes in the value of the ratio in Eq. (53).

Estimate of the effects of longitudinal SDCs in the isovector channel
The π 0 -column of Tab. 1 collects all uncertainties in our estimate of the effects of longitudinal SDCs in the isovector channel, as described in the previous subsections. By combining them in quadrature we get of longitudinal SDCs assuming that the low-energy region is dominated by ground-state pseudoscalar poles, whose contributions are taken as input. In each flavor channel the results are presented as the shifts ∆a µ,ref with respect to the pole contributions for a specific reference set of parameters and a list of uncertainties corresponding to different choices for each of these parameters. In the last two rows, these uncertainties are added in quadrature and the final range is symmetrized. See main text for details. Since we do not regard the reference parameterization as the central value, we symmetrized the uncertainty to finally obtain the range Using instead δ + m,2 ∆a Eq. (51), the final result would be (3.8 ± 2.7) × 10 −11 . Notice that, despite the fact that it likely overestimates the range of longitudinal short-distance effects, this interval is still definitely compatible with the current precision goal. Fig. 4 shows the contributions to the quadratic error from the different sources discussed above. The vastly dominant effect stems from the interpolation between low and high energies, with an especially crucial role played by the choice of m, the scale at which the matching between the low-energy representation ofΠ 1 and the interpolant is performed. The uncertainties δ int , δ m and δ P (x,y) could be reduced by additional low-energy input concerning further intermediate states and higher-order terms in the symmetric and asymmetric OPEs, which would help constrain the coefficients b i (r, φ) in the interpolants in Eqs. (36) and (37). The uncertainties related to the perturbative corrections are considerably smaller. While we do not expect that calculations of α s corrections will crucially improve the final estimate, these perturbative results will definitely be important to better assess the regime of validity of the asymptotic constraints and thereby verify and sharpen some of our assumptions.
We have also checked that our results are robust against the choice of different reference sets of parameters.
where all other sources of uncertainty are included. We obtained similar results by selecting as reference different interpolants or different values of the number of free parameters N contained therein.

The isoscalar contributions
In this section the procedure presented above for the isovector case is applied to the isoscalar channels with η/η -poles as low-energy input. In our analysis, we employed the Canterbury TFFs from Ref. [51] in the reference solution. 15 We determined the parameters encoding η − η -mixing as explained in Secs. 2.2 and 2.4 and obtained which shows that δ 0 is indeed sizable. Following the same procedure for the construction of the reference interpolant as in Sec. 5.1, we found The uncertainty estimation proceeds in the same way as in the isovector channel up to minor modifications. Since error bands for the doubly-virtual TFFs in all kinematic configurations are not available in the literature, we estimated uncertainties by considering another TFF representation, namely the one based on Dyson-Schwinger equations (DSE) [52]. This yields ∆a η µ,ref = 2.10 × 10 −11 and ∆a η µ,ref = 4.18 × 10 −11 . The fact that individual results for η and η channels differ by 18 and 7 %, but the sum only by 3 % can be understood by comparing the mixing parameters against those obtained from the Canterbury parameterization. These coefficients enter Π (8/0)−η/η , asymp 1 quadratically, which leads to a reshuffling between a Since NLO results are not available for the η/η TFFs, we estimated the NLO OPE uncertainty by simply rescaling the one in the pion channel by the ratio of the reference outcomes. Due to the smallness of this uncertainty, this is expected to be sufficiently accurate.
For the range [m min , m max ] and the minimal allowed value Σ t for Σ match in the Monte Carlo simulation for P (x, y), we took the results from the isovector channel. We rescaled the π(1300) term below the matching surface by the ratio of reference results when adding this contribution to δ m a (8/0)−η/η µ . This is justified by the fact that the first excited pseudoscalars in the three flavor channels have similar masses, despite the large mass difference of the pseudo-Goldstone bosons.
All results are collected in Tab. 1 and our final estimate for the longitudinal shortdistance effects in a The relative contributions to the uncertainties are similar to the pion case illustrated in Fig. 4. A more precise description of η − η -mixing would of course help better separate the two isoscalar channels but would not play an important role in their sum leading to negligible shifts in the total contribution from longitudinal SDCs.

Sum over the flavor channels and comparison with literature
Combining the results from Secs. 5.1 and 5.2, obtained under the assumption that the ground-state pseudoscalar mesons dominate the low-energy region, our estimate for the total effect of the longitudinal SDCs on HLbL amounts to where we have combined the three uncertainties linearly since they originate from the same sources in all three channels. This result is remarkably close to what is expected based on flavor symmetry considerations. If the U (3) symmetry emerging in the combined chiral and large-N c limit is assumed, then ∆a  µ as a function of a lower limit on Q 2 3 in Eq. (9): our reference result and corresponding error band against the tower of excited pseudoscalars in the large-N c Regge model 1 of Refs. [32,33] and the curve from the MV model [25] evaluated using the up-todate dispersive pion TFF. At small non-vanishing Q 2 3,min , our reference curve is constant due to the finite Σ match , which for P (x, y) = 0 corresponds to m 2 = Q 2 3 = const., whereas the Regge model has a slope due to the absence of such a cutoff. The upper end of our error band shows a slope because of the inclusion of the π(1300) contribution in that region.
For this reason we do not expect that a more refined analysis of the subtler isosinglet contributions is going to change substantially our final results.
Refs. [32,33] have recently studied the possibility of saturating SDCs away from the chiral limit by including a tower of excited pseudoscalar states in the context of a Regge model matched to the pQCD quark loop. Their outcome is ∆a long µ = 13(6) × 10 −11 , which is well compatible with ours within errors. For the η -channel, the Regge model yields ∆a (8/0)−η µ = (6.5 ± 2.0) × 10 −11 , which is somewhat larger than our result but still compatible within errors. 16 This can partly be explained by the different value for C η used in Refs. [32,33], namely C η = 0.239, which results from imposing that Eq. (24) holds exactly. Fig. 5 shows ∆a (3) µ as a function of a lower cutoff on Q 2 3 in our approach as well as the large-N c Regge model 1 of Refs. [32,33]. In order to obtain this plot, we calculated the integral in Eq. (9) as a function of a lower limit on Q 2 3 (which depends on Σ, r and φ) for both the full a (3) µ as well as the pion pole contribution (cf. Eq. (40)). The Regge model result lies within our error band for all Q 2 3,min . Our estimate of longitudinal short-distance effects as well as the one in Refs. [32,33] are smaller than the shift obtained in Ref. [25], ∆a long µ = 23.5 × 10 −11 , which even increases to about 38 × 10 −11 if up-to-date TFF input is used [33]. These large values are due to two features of the MV model: the fact that the singly-virtual TFF is set to a constant over the whole integration region and not only in the OPE regime, and the fact that in the symmetric asymptotic limit the parametric momentum dependence is correct but its coefficient is too large. Both of these features can be clearly seen in Fig. 5 and are responsible for the discrepancies in the slope at small Q 2 3,min and the values at large Q 2 3,min , respectively. Refs. [40,41] have studied how the inclusion of an infinite tower of axial-vector mesons could help satisfy the OPE SDCs, focusing for this purpose on the relevant TFFs in the context of holographic QCD models. According to Ref. [40], the tower of axial-vector mesons contributes (29-41) × 10 −11 to a HLbL µ of which (57)(58) % are attributed to a long µ . Using instead holographic QCD input only for the momentum dependence of the TFF and fixing its normalization from experiment reduces the estimate of the contribution to a HLbL µ from the tower of axials to 22(5) × 10 −11 . Ref. [41] finds 14 × 10 −11 for the effect of axials on a long µ . Thus, the results of these studies appear to be at the high end of our range in Eq. (62). However, we stress that comparing these numbers with our result is not properly justified. Indeed, while in these models the parametric Σ-dependences implied by pQCD and the OPE in the respective limits are correctly reproduced, the coefficients thereof are typically too small. In addition, the lightest multiplet of axials significantly altersΠ 1 at small photon virtualities, which implies that in our approach its contribution should be included in the low-energy representation. This aspect will be discussed in the next section, also to show how information on additional states in the 1 GeV region can be incorporated in our analysis.

Including ground-state axial mesons at low energies
Here we adopt a model-dependent approach to illustrate the application of our procedure to the case of the inclusion in the low-energy region of the lightest of the axial-vector mesons, for which no dispersive treatment in the BTT formalism is available yet. According to the holographic QCD models in Refs. [40,41] and using the notation of Ref. [40], the contribution toΠ 1 of an axial meson of mass M A in the flavor channel a can be written as Π (a), axial 1 where A(Q 2 1 , Q 2 2 ) is the axial TFF. Among the various scenarios analyzed in Ref. [40], the hard-wall model by Hirn and Sanz (HW2) [81], which was also studied in Ref. [41] with different parameters, reproduces best the measured mass, the measured equivalent two-photon decay width and the singly virtual momentum dependence measured by L3 for the lightest multiplet [82,83]. Furthermore, it yields asymptotic axial TFFs whose momentum dependence is consistent with the behavior derived in Ref. [29]. The infinite tower of axials has the correct momentum scaling in the asymmetric asymptotic regime dictated by the OPE constraint, but the coefficient is 38 % too small [40].
We focused on the isovector channel, which is sufficient for our illustrative purposes, and thus on the inclusion of the a 1 meson. Based on the HW2 model, we obtained a a 1 µ,HW2 = 3.3 × 10 −11 for the a 1 contribution to a long µ . The rest of the tower of isovector axial mesons in this model yields ∆a (3), A, HW2 µ = 0.8 × 10 −11 , implying that in this framework about 80 % of the total effect comes from the lightest state. 17 By matching the interpolants in Eq. (36) to the contributions from the pion and the holographic a 1 with the reference set of assumptions in Sec. 5.1, we obtained 18 This result is more than twice as large as the resummed tower in HW2, ∆a (3), A, HW2 µ , but the significance of this discrepancy could only be assessed by a more sophisticated analysis including uncertainties, which is beyond the scope of this work. However, in the holographic model the infinite tower of axials does not fully saturate the pQCD nor the OPE constraints, which suggests that additional degrees of freedom besides axials should be included in a more realistic model.
We then considered the choice of parameters made in Ref. [41] and referred to as HW2(UV-fit) in Ref. [40]. This model is constructed to obey the OPE constraint exactly, but fails to describe low-energy physics like the ρ-meson mass, the pion TFF and the axial TFFs measured by L3. The longitudinal contribution from a 1 in this case amounts to a a 1 µ, HW2(UV-fit) = 3.4 × 10 −11 and the tower of states increases the value by ∆a Neglecting issues related to intrinsic model dependence in the low-energy input, our method based on interpolants that by construction satisfy all constraints indicates that the effects of longitudinal SDCs are relatively small compared to the dominant low-energy contributions, and what is crucial in order to achieve higher precision is to gain control over the latter. We stress that a reliable prediction with a robust uncertainty estimate of the effects of axial meson exchanges would require model-independent input information.

Conclusions
In this paper we have introduced a novel approach to incorporate longitudinal SDCs into the calculation of the HLbL contribution to the muon g − 2. At variance with the previous estimates based on hadronic models, we have constructed general functions interpolating between low-, mixed-and high-energy regions, without resorting to specify which and how hadronic intermediate states are responsible for saturating the constraints. Furthermore, our method allows us also to study in detail the role played by parameters and assumptions in a transparent and numerically efficient way.
Our main premise is that an accurate low-energy representation of the longitudinal functionΠ 1 (Q 2 1 , Q 2 2 , Q 2 3 ) entering the HLbL integral can be obtained by taking into account only intermediate states that are under good theoretical and numerical control. For the π 0 , due to the location of its pole, the form of this low-energy representation can be straightforwardly extended even to large Q 2 1 and Q 2 2 as long as Q 2 3 stays small. Using available input for the π 0 -pole term, we find that the shift due to longitudinal SDCs on the isovector part of a HLbL µ is in the range (2.6 ± 1.5) × 10 −11 . By including in the analysis also the isoscalar components, which the η-and η -poles are assumed to dominate at low energies, we obtained that longitudinal SDCs increase a HLbL µ by (9.1 ± 5.0) × 10 −11 in total. The quoted ranges encompass uncertainties in the low-energy input, perturbative 18 We had to shift C 2 π by 1.7 % in order to account for the a1 at small Q 2 3 and asymptotic Q 2 1 ∼ Q 2 2 .
corrections and fitting errors at asymptotic momenta, parametric variations of the functional form of the interpolants and of the matching surface, at which these functions are matched to the low-energy input, with the latter dominating the total uncertainty. Thus, according to our analysis, infinite towers of states heavier than 1 GeV, albeit crucial for the saturation of SDCs, give a relatively small contribution to a HLbL µ and this effect can be estimated with sufficient precision using our method. Conversely, states with masses around 1 GeV contributing significantly to the low-energy region play a decisive role also in a precision determination of short-distance effects.
Our result for the effects of longitudinal SDCs on a HLbL µ agrees with recent model estimates [32,33], fulfills the accuracy goal set by the forthcoming experimental results and is significantly smaller than the earlier model result of Ref. [25], especially when up-to-date TFF input is used. Furthermore, neglecting issues concerning intrinsic model dependence and the fact that holographic QCD calculations in Refs. [40,41] do not completely saturate the SDCs, we find in agreement with these studies that the infinite tower of axials has a relatively small impact on the longitudinal part of a HLbL µ if the lightest multiplet is treated explicitly as a low-energy contribution.
It will be straightforward to incorporate in our approach model-independent information on further intermediate states as well as higher-order corrections to asymptotic expressions once these become available. Furthermore, our method can be generalized to the case of transversal SDCs. Therefore, it paves the way for a combination of all available low-and high-energy information on HLbL into one model-independent, accurate numerical estimate of this contribution to the muon g − 2.
Let us now examine under which conditions the pole inΠ ) occurs for Σ pole < Σ match /2. This relation is independent of the low-energy input but depends on the pseudoscalar mass inΠ (a), asymp 1 . In our reference interpolant we set P (x, y) = 0 in Σ match and m 2 = 0.5 GeV 2 . For these choices the zero inΠ (a), asymp 1 is at sufficiently low Σ for all flavor channels and all (r, φ). ForΠ to 0.094 GeV 2 without violating the requirement Σ pole < Σ match /2. In all cases this does not place a serious limitation on the values for m we consider in the uncertainty estimation in Sec. 5.1.4. We have also checked that in the Monte Carlo sampling of the polynomial P (x, y) in Sec. 5.1.4, the requirement does not lead to a further constraint beyond Σ match (r, φ) > Σ t . Only inΠ (8/0)−η , asymp 1 , there is a very small region around (r, φ) = (1, π) where values of Σ match (r, φ) slightly larger than Σ t have to be excluded. In this region, however, the pole at (r, φ) = (1, π) ensures that Σ match (r, φ) is sufficiently large for all sampled parameter sets also without imposing that additional condition. SinceΠ  (36) for N → ∞ and the requirement Σ pole < Σ match /2 is trivially fulfilled due to Σ pole = 0. The convergence of interpolant 3 given in Eq. (37) also easily follows from that of interpolant 1, because the logarithmic term can be written as a Taylor series at finite (Σ match ) −1 so that the parameter b 1 in Eq. (37) is redundant as N → ∞.