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 μ 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 μ HLbL . Taking as input the known ground-state pseudoscalar pole contributions, we obtain that the longitudinal SDCs increase a μ HLbL by ( 9.1 ± 5.0 ) × 10 - 11 , where the isovector channel is responsible for ( 2.6 ± 1.5 ) × 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 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 (i.e. 0.14 ppm) [2], make it even more crucial and timely to further scrutinize and improve control over theory predictions.

Together with the hadronic vacuum polarization contribution, the hadronic light-by-light (HLbL) is the major source of theoretical uncertainty in the Standard Model [3-5]. In the last years, significant efforts have been devoted to improve the determination of and reduce model dependence by using analytic approaches based on dispersion relations [6-17] as well as lattice QCD [18-24]. In particular, the dispersive framework for the HLbL tensor in Refs. [6-11] has enabled accurate data-driven determinations with controlled error estimates of the contributions from one- and two-pion intermediate states.

In this framework 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 for asymptotic values of (subsets of) the photon virtualities. Since one of the photons corresponds to the static electromagnetic source in the definition of , one asymptotic regime is realized when the remaining three space-like photon virtualities are comparable and much larger than , and the other when two space-like photon virtualities are much larger than both the third and . 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 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-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--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 .

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 (), whose contribution is under firm theoretical control thanks to a dispersive evaluation [12-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 Sect. 2 we review the relevant constraints on HLbL and the assumptions made in their derivations. Section 3 describes our interpolation between the OPE and pQCD asymptotic constraints while in Sect. 4 we discuss its smooth connection with the low-energy region. In Sect. 5 we present our numerical analysis with particular emphasis on the error estimation. Conclusions are drawn in Sect. 6. Appendix A is devoted to the analysis of the convergence properties of our interpolants.

Longitudinal short-distance constraints on HLbL

Master formula for 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 is governed by the fourth-rank vacuum polarization tensor for fully off-shell photon-photon scattering in pure QCD,with momenta assigned as . In this expression, the electromagnetic current for the light quark triplet is given byBy 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 are free of kinematic singularities.

The HLbL contribution to can be derived from the tensor by using projection operator methods [5, 45, 46] in the static limit . After performing a Wick rotation to Euclidean momenta, angular averages [47, 48] lead to the master formula [11]where denote the magnitudes of the Euclidean loop four-momenta, , and is the cosine of the angle between these vectors. The scalar functions are linear combinations of the previous for . The analytic expressions of the integration kernels are given in Ref. [11].

Parameterizing the three-dimensional integration domain by the coordinates [49]which are related to the non-vanishing photon virtualities bywill prove very useful in the following discussion about asymptotic constraints on HLbL. The master formula in Eq. (4) then takes the formIn terms of the coordinates, the integration domain amounts to a cone with tip at the origin. In terms of , a given point in this cone is specified by the distance to the tip of the point’s projection on the symmetry axis (), and by the polar coordinates r and on the plane containing the point and orthogonal to the symmetry axis, normalized such that corresponds to the surface of the cone.

In the master formula, a special role is played by the scalar functions , which fulfillwhere the crossing operator 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 and also / poles yield sizable effects. Furthermore, distinctively, the OPE SDCs on do not involve other scalar functions [33].

The functional form of in specific kinematic regimes is constrained according to analytic QCD results, which we will fully exploit to estimate the impact of (longitudinal) SDCs onwhere the shift in the variable in corresponds to the crossing operation on . Thus for our analysis, we only need to study one BTT scalar function in the kinematics.

For the purpose of later discussion, we stress here that a pole term in due to a single-particle intermediate state of mass M yielding the denominator leads to a hierarchy among contributions in the space-like momentum region relevant for . For small values of larger masses get suppressed, while for 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 , 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 .

The lightest state contributing to is . The unitarity relation for a single pseudoscalar intermediate state yieldswhere 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 elementwith . If we set , then amounts to the pseudoscalar pole contribution . In the case, this has been evaluated within a few percent accuracy via a data-driven dispersive approach [12-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 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 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 [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. . 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 corresponds to and while the symmetric configuration occurs in a neighborhood of (see Fig. 1).

Fig. 1

The circle represents the boundary of the integration domain for a fixed value of . The angles , and correspond to , and , respectively. The colored regions denote where SDCs on hold at large . The blue domains yield contributions to 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

In the following we will review the relevant constraints on at large and describe in detail our method to provide general families of interpolants for between low- and high-energy regions in the integral.

The asymmetric asymptotic region: OPE constraints

For large Euclidean values of , one can expand the time-ordered product of two electromagnetic currents, which defines the tensorinto a series of local operators. At leading order in , by matching single-quark matrix elements and omitting the unit operator, which does not contribute to the connected HLbL tensor in Eq. (1), one obtains [59]where the axial current is defined by with charge matrix given in Eq. (2). The ellipsis denotes sub-leading terms suppressed by powers of . This result implies that, at leading order in the OPE and at leading order in , the HLbL tensor can be expressed in terms of the correlator of two vector currents with an axial current,for . This three-point function, which also appears in the calculation of fermion loop electroweak contributions to [60, 61], 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. [62]). The corresponding longitudinal scalar function determines the asymptotic behavior of in the asymmetric region and is fixed by the axial Adler–Bell–Jackiw anomaly up to chiral corrections and the gluon anomaly. Neglecting these effects, which will be discussed below, this translates into the following constraint [33] for the singlet and octet flavor components of , defined by the decomposition of the axial current,where in terms of the charge matrix and Gell–Mann matrices . In particular,Since Eq. (16) relies on a perturbative calculation of the VVA correlator, it holds in the kinematic limit . For the non-singlet components (), since perturbative [63] as well as non-perturbative [64, 65] corrections are absent in the chiral limit, the hierarchy between and can be dropped. In contrast, the singlet channel () is affected by the gluon anomaly, even in the chiral limit. This does not modify Eq. (16) for [33], but the extrapolation to small is only valid if in addition the large- limit is considered, where the anomaly vanishes. Furthermore, in the crossed kinematics ( and ), the leading-order OPE contributions to vanish.

Let us now compare against the pseudoscalar pole contributions, focusing on the pion pole first. In the chiral limit and using the fact that at leading order in [66, 67], one findsThis expression has a pole at since . The location of this pole as well as its residue agree with in Eq. (16) (cf. Ref. [25]), which is consistent with the pion being the only massless isovector state in the chiral limit.

For finite quark masses, the pole in is shifted from to , which lies outside the integration domain for . The closest point in the integration region for fixed asymptotic (see Fig. 1) is at , whereHere few percent chiral corrections to [68-72] have been neglected. This is still close to the actual pole, which leads to the enhancement by . Since no other contribution receives the same enhancement, the last expression is expected to provide an excellent approximation to the true 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 contributionThis is also consistent with the OPE result in Eq. (16) for , where chiral corrections are sub-leading. Thus, Eq. (20) is exact for , relies on the assumption of pion dominance at and has the correct chiral limit Eq. (16) for all . We extend it to the / channels and writeHere but cannot be directly identified with 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 . This assumption implies that can be read off from the pole contributionsOne can show that in the chiral limit and neglecting the gluon anomaly [33]At this point we note that, besides the corrections to the TFFs and OPE coefficient discussed in Sect. 2.3, which affect all ground-state pseudoscalars in the same way, the gluon anomaly induces a running of the flavor singlet decay constant [73-75]. 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 [76, 77].

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 implyup to the above-mentioned anomaly-induced scale- dependence, which leads to a violation of this equality (cf. Sect. 5.2).

For , the additional -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 Sect. 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 Sects. 5.3 and 5.4, we will compare the outcome of our analysis against these estimates of the effects of longitudinal SDCs.

corrections to the OPE

The derivation of Eq. (14) has been performed at leading order in . Since no other operator of dimension 3 can appear in that OPE, corrections only affect the OPE coefficient of the axial-vector current. In Refs. [78-80], this coefficient has been calculated to next-to-leading order (NLO). Including this contribution in Eq. (14) leads toIt follows that the NLO version of Eq. (21) readsThe two-current operator product 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. In fact, the symmetric asymptotic pion TFF has been calculated to NLO in Refs. [14, 81],7which is consistent with Eq. (25). The fact that the corrections agree between the HLbL tensor in the asymmetric asymptotic limit and the symmetric asymptotic pion TFF guarantees that the pion pole saturates at in the chiral limit also beyond leading order in .

A comment on the renormalization scale dependence of the terms in Eq.(25) is in order here. The non-singlet components of the axial current are conserved (up to quark mass corrections) and thus their anomalous dimensions vanish. This is not true for the singlet component due to the gluon anomaly [82, 83], but we neglect this effect here because it starts at . Therefore, since the perturbatively expanded dimensionless part d of the Wilson coefficient is scale () independent,and the terms () can be resummed using as input the -function and the one-loop result with evaluated at the scale (see also [84]).

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 -kinematics. At leading order in this OPE and at leading order in [33],In the symmetric limit, neglecting terms that are suppressed by powers of ,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 decays like , (towers of) hadronic contributions beyond ground-state pseudoscalar poles have to be responsible for the behavior shown by Eq. (30). 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 satisfyingSince Eq. (24) is violated, we definewhere the parameter is chosen such that Eq. (31) holds.

Interpolating between asymptotic constraints

We approximate the true 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 at small determined by single-particle intermediate states. Here we work at leading order in . Perturbative corrections will be discussed in our numerical analysis in Sect. 5.

The relevant constraints on at large are given by Eq. (21) for and Eq. (30) for . Both expressions as well as the vanishing result of the leading-order OPE contribution in the crossed kinematics are compatible withif . Thus, Eq. (33) interpolates between symmetric and asymmetric asymptotic limits. According to Sect. 2.4, parameterizes the anomaly corrections to the singlet VVA correlator and the resulting shift in with respect to . Since a term proportional to and independent of does not change the leading behavior at and thus does not spoil compatibility with the OPE constraint, we subtract from Eq. (33) in the case of /.

Obviously, the choice made in Eq. (33) 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 and we demand it to be finite and analytic for all .8 Therefore it can be approximated by a Taylor series in and truncated after order M,wherefor integer j, due to the pQCD constraint and crossing symmetry.

Up to now, we have applied the quark-loop result only at . However, the fact that Eq. (29) holds also in a neighborhood of this point can be used to fix the coefficients . To this end, we fitted Eq. (34) at fixed asymptotic with to Eq. (29) for .9 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 in the pion channel readand are all at most , as expected since in Eq. (34) they parameterize relative corrections. This holds true also for the channels, where the numerical values are different. In Sect. 5.1.2 we will discuss uncertainties due to the chosen fitting range, the number of parameters in the fit and corrections to the asymptotic constraints.

Interpolating between low and high energies

The next step is to smoothly connect our representation of for given by Eq. (34) to an accurate low-energy description. We achieve this by adding suitable terms to 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 at a suitably defined surface . In Sect. 4.2 we will discuss how is related to the mass scale at which intermediate states beyond the ones explicitly considered start to affect .

Interpolation functions and matching procedure

For we consider the following two interpolation functionswhose leading terms at asymptotic are given in Eq. (34), whereas below the matching surface we set . In Appendix A we will show that these functions converge to the true in the limit when matched to exact low-energy input using the convergence property of a Taylor series. The two different forms given in Eq. (37) will be used to estimate the sensitivity of our numerical results on the specific choice of interpolation between low and high energies.10

The coefficients are fixed from the requirement that the have the same value and the same -derivatives as the low-energy representation if evaluated at for each . 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 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 . Since it is desirable to match at least one derivative in order to have differentiable at the matching point, we will use 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 , where M is some non-perturbative mass scale. In fact, the Regge model considered in Refs. [32, 33] leads to interpolants containing terms like for , where 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 interpolantand use again .

The matching surface

The remaining crucial ingredient in our procedure is the function , which determines the value of at which the matching is performed for given r and . Choosing it too low leads to important modifications of at low energies with consequent overestimation of .11 Conversely, choosing too high assumes the low-energy input to dominate beyond what is expected according to mass and phase-space considerations and thus leads to underestimate .

For small values of , the , , poles are assumed to dominate independently of , due to the pole at (see Sects. 2.1 and 2.2). This implies that no matching is needed in this regime, i.e. . The most general function that is analytic for all except for a (first-order) pole at can be written aswhere determines the matching scale at and P is a polynomial with two arguments and no constant term. The transformation property of under crossing specified in Eq. (8) restricts P to contain only even powers of .

The parameter sets the absolute mass scale of and should thus be related to the masses of the states affecting beyond the ones explicitly included, namely , , here. In the following, we will assume that contributions to in the 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.12 This is realized for example in the large- limit of pure QCD: since the short-distance expressions for in both symmetric and asymmetric limits scale like , these can indeed be saturated by single-meson exchanges (see Ref. [85]). Non-resonant contributions from multi-hadron intermediate states (like , 2K, , ,...) are sub-leading for large and thus cannot contribute to the SDCs. Since scalar mesons have no impact on , the lightest states beyond the ground-state pseudoscalars that are the most relevant at small (see Sect. 2.1) are the axial mesons like and the tensor mesons like , with masses in the 1–2 GeV region, whose effects on can presently be estimated only using hadronic models.

For , corresponds to . Since a state of mass M ceases to be suppressed by the denominator compared to lighter states when approaches , should be chosen well below . At the same time, it should not be taken too small, because we do not expect any large contribution to at . We thus regard as a good starting point for our analysis. In Sect. 5.1.4 we will discuss a range of choices for this parameter as well as the effects of the polynomialwhich we have estimated by means of a Monte Carlo sampling over the coefficients .

Numerical results and error analysis

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 in this channel is the , 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 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 . In Sect. 5.2 we will discuss the inclusion of and in Sect. 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 singly- and doubly-virtual TFFs [13, 14], while the corresponding correction is included in the uncertainty. As reference interpolating function, we used with (see Eq. (37)), which turned out to yield results that are central in the range spanned by the interpolants 1, 2 and 3 and (cf. Eq. (49) below). For the asymptotic function, we included information from pQCD away from in the way explained in Sect. 3, while corrections contribute to the uncertainty. For the matching surface we used Eq. (39) with and . The resulting function, which we call , is shown in Fig. 2 for 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 due to the longitudinal SDCs in the isovector channel iswhere comes from using in the master formula Eq. (9), and is given in Eq. (12) according to Refs. [13, 14]. In Sect. 5.1.5 we will argue that our final result does not strongly depend on the choice of the reference set of parameters.

Fig. 2

The pion pole contribution and associated uncertainty from Refs. [13, 14] vs. the reference interpolant and its error band which includes all sources of uncertainty considered in the present analysis (see discussion in Sects. 5.1.1–5.1.4 below)

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. (41), we obtained an asymmetric error band around the reference result with boundary valueswhich correspond to the asymmetric error for the dispersive TFF. Given the smallness of these uncertainties, the (negative) correlation between them and the uncertainties of 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 , the Canterbury approximant with of Ref. [51] and the Dyson–Schwinger TFF from Ref. [52]. We obtainedwhich 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 (see Eq. (34)), which are related toWe start by discussing the fit to the quark-loop result. In Sect. 3, we chose , which leads to a strongly improved fit quality compared to . Considering a larger value of M gives an estimate of the errors made by approximating the pQCD result by a polynomial at at fixed asymptotic and by extrapolating to the regime , which is unknown except for the OPE constraint. Choosing shifts the result for by only indicating that the truncation at is sufficient. We also studied the effects of a substantial reduction of the radius, namely from down to , where no logarithm of ratios of squared momenta is larger than 1. We found that this leads to a small shift . We did not consider 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 giveswith respect to the reference contribution of longitudinal SDCs to the pion pole input in Eq. (41).

the choices made in the fit to the quark-loop result that lead to Eq. (36), namely the degree M of the polynomial and the radius of the fitting domain;

corrections to the OPE constraint as given by Eq. (26);

corrections to the quark loop.

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. (26))Here the Heaviside step function ensures that the perturbative correction only affects a region around , whose size can be varied via the free parameter A. By setting , this region does not extend into the domain. The choice of the renormalization scale is the same as in Ref. [14] and agrees with our discussion in Sect. 2.3 up to a factor of 2 since in the relevant regime. In our numerical analysis, for the running of we used the three-flavor one-loop beta function and matched to .

According to the discussion in Sect. 2, the OPE constraint in the chiral limit is saturated by the pion pole at 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. (37) by and matched the correspondingly modified to the pion-pole contribution with TFFs including effects [14, 81].13 Using this interpolant and this pion-pole result, our outcome for is larger than the reference result Eq. (41) byFor , the domain where the correction applies extends down to , but nevertheless the shift of turns out to be and thus still negligible. The smallness of these shifts can be understood from the large values of 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 corrections. We assumed in analogy with Eq. (45),and as in the leading-order quark loop fit, we set . Using this expression in Eq. (37) for the matching to the pion-pole with leading-order dispersive TFF, we obtained a shift of compared to the reference result. Even when inflating this uncertainty by a factor of 2,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. (37) and (38) 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 arewhere given by Eq. (41) 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. Settingall values listed above are within the range .

Choice of

The function in Eq. (39) contains the mass parameter m and the polynomial P(x, y), which has been set equal to zero so far. We have argued in Sect. 4.2 that m should be chosen considerably smaller than the mass M of the lightest resonances contributing to in addition to the ground-state pseudoscalar mesons. For this reason, for the reference interpolant we set . Here we discuss the effects of alternative choices for this parameter within a range between and .

Since according to Sect. 4.2 a conservative choice for the upper end of the range is , we set . In order to determine an appropriate value for , one has to estimate isovector contributions beyond the -pole. Following our argument in Sect. 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 for the following reasons. Models for tensor mesons around 1 GeV give similar or smaller contributions to [16, 86, 87]. For ground-state axials, recent studies based on different approximations and hadronic models yield quite different numerical results, see e.g. Refs. [40, 41, 86, 88], leading to large uncertainties. If future model-independent analyses show that axial-meson exchanges are responsible for significant effects in 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 below . Neglecting issues related to model dependence, in Sect. 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 with two photons can be described by a TFF, which determines the contribution to as in Eq. (10). In our analysis we used as input the TFF derived in Ref. [89] in the framework of Resonance Chiral Theory (RT) [90]. We fixed the free parameters in Eq. (69) of Ref. [89] by requiring that (i) the -TFF satisfies the Brodsky–Lepage condition, i.e. and (ii) the two-real-photon limit of the excited pion TFF is in the range , argued for in Ref. [33] based on experimental results [91, 92]. 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 can be illustrated by means of Fig. 3. Given a value of , at large enough , the range of interpolants (orange band) spanned by safely includes the green band representing the sum of the contributions from and , including errors on the latter due to the range for . Since this is not the case at small , the contribution to from this region is underestimated in our approach. To gauge this effect, we calculated the integral in Eq. (9) withusing the maximal 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 at a fixed . Repeating this calculation for different values of yields the function and by inverting this, we determined by fixing to values well below the accuracy goal set by forthcoming experimental results. For we obtained and for , .

Fig. 3

The figure displays the dispersive pion pole contribution, the reference interpolant and the (orange) band corresponding to the various choices of the parameter m. The blue line indicates the value of the matching surface for . The green band shows the sum of the - and -pole contributions, where the latter has been calculated using input from RT and phenomenology, including errors

Numerically, leads to the shift and the two values and yield and , respectively. If we add to the latter numbers, we obtain the conservative estimatesIn the following we will use for the main results and keep as an alternative, even more conservative uncertainty.

We also considered a different parameterization for the TFF, namely the one given by the Regge model in Refs. [32, 33]. Using the empirical instead of the Regge-model value of 1.36 GeV used in these references and following the same procedure discussed above, we obtained and . This leads towhere we again added to the uncertainties in the upward direction. For our final result we use the more conservative uncertainty estimates given in Eq. (52).

In order to study the effects of the polynomial in , Eq. (40), we set and sampled the free parameters according to a standard normal distribution. Since the pion pole gives an excellent approximation of for very small at any , we only allowed for parameters giving for all , where is defined as the smallest value of such thatholds for some . With calculated using RT, we obtained . This condition ensures that there are no large contributions from our interpolation at points where RT predicts a very small excited pion contribution. From this we calculated a distribution of results for , 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 errors for a Gaussian. This givesWe 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. (54).

Estimate of the effects of longitudinal SDCs in the isovector channel

The -column of Table 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 getSince we do not regard the reference parameterization as the central value, we symmetrized the uncertainty to finally obtain the rangeUsing instead in Eq. (52), the final result would be . 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.

Table 1

The effects on 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 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

	π^0	η	η^′
Δa_μ,ref×10^11	2.56	2.58	3.91
δ_TFFΔa_μ×10^11	-0.13+0.06	0.47	0.30
δ_pQCDfitΔa_μ×10^11	0.09	0.08	0.14
δ_NLOOPEΔa_μ×10^11	-0.00+0.01	-0.00+0.01	-0.00+0.02
δ_NLOpQCDΔa_μ×10^11	0.36	0.36	0.55
δ_intΔa_μ×10^11	0.62	-0.65+0.61	-0.84+0.74
δ_m,1Δa_μ×10^11	-1.20+1.31	-1.17+1.27	-1.60+1.68
δ_P(x,y)Δa_μ×10^11	-0.32+0.39	-0.33+0.31	-0.43+0.32
δ_totΔa_μ×10^11	-1.44+1.55	-1.50+1.56	1.97
Δa_μ×10^11	2.6±1.5	2.6±1.5	3.9±2.0

Figure 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 and the interpolant is performed. The uncertainties , and 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 in the interpolants in Eqs. (37) and (38). The uncertainties related to the perturbative corrections are considerably smaller. While we do not expect that calculations of 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.

Fig. 4

Relative contributions to the total uncertainty in the isovector channel. For asymmetric errors the mean of the squared errors is used

We have also checked that our results are robust against the choice of different reference sets of parameters. For example, if we set equal to the previous boundary values for the uncertainty in the reference configuration, namely and and choose as the range for the error estimation as before, we getwhere 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.14 We determined the parameters encoding -mixing as explained in Sects. 2.2 and 2.4 and obtainedwhich shows that is indeed sizable.

Following the same procedure for the construction of the reference interpolant as in Sect. 5.1, we foundThe 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 and . The fact that individual results for and channels differ by 18 and 8%, but the sum only by 3% can be understood by comparing the mixing parametersagainst those obtained from the Canterbury parameterization. These coefficients enter quadratically, which leads to a reshuffling between and . Due to Eq. (24) this effect drops out in the sum up to the anomaly correction affecting the OPE regime. As TFF contribution to the uncertainty on we took the absolute value of the differences between the results from the Canterbury and DSE TFFs.

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 and the minimal allowed value for in the Monte Carlo simulation for P(x, y), we took the results from the isovector channel. We rescaled the term below the matching surface by the ratio of reference results when adding this contribution to . 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 Table 1 and our final estimate for the longitudinal short-distance effects in readsThe 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 Sects. 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 towhere 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- limit is assumed, then . Using our isovector uncertainty and adding linearly a standard 30% U(3) breaking effect, we obtainFor this reason we do not expect that a more refined analysis of the subtler isosinglet contributions is going to change substantially our final results.

References [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 , which is well compatible with ours within errors. For the -channel, the Regge model yields , which is somewhat larger than our result but still compatible within errors.15 This can partly be explained by the different value for used in Refs. [32, 33], namely , which results from imposing that Eq. (24) holds exactly.

Figure 5 shows as a function of a lower cutoff on in our approach as well as the large- 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 (which depends on , r and ) for both the full as well as the pion pole contribution (cf. Eq. (41)). The Regge model result lies within our error band for all .

Fig. 5

as a function of a lower limit on in Eq. (9): our reference result and corresponding error band against the tower of excited pseudoscalars in the large- Regge model 1 of Refs. [32, 33] and the curve from the MV model [25] evaluated using the up-to-date dispersive pion TFF. At small non-vanishing , our reference curve is constant due to the finite , which for corresponds to , 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 contribution in that region ss

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], , which even increases to about 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 and the values at large , respectively.

References [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 to of which (57–58)% are attributed to . 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 from the tower of axials to . Reference [41] finds for the effect of axials on . Thus, the results of these studies appear to be at the high end of our range in Eq. (63). 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 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 of an axial meson of mass in the flavor channel a can be written aswhere is the axial TFF.

Among the various scenarios analyzed in Ref. [40], the hard-wall model by Hirn and Sanz (HW2) [93], 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 [94, 95]. 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 meson. Based on the HW2 model, we obtained for the contribution to . The rest of the tower of isovector axial mesons in this model yields , implying that in this framework about 80% of the total effect comes from the lightest state.16

By matching the interpolant in Eq. (37) to the contributions from the pion and the holographic with the reference set of assumptions in Sect. 5.1, we obtained17This result is more than twice as large as the resummed tower in 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 in this case amounts to and the tower of states increases the value by . Our reference interpolant leads to , which is again larger than the model result. However, also in HW2(UV-fit) the pQCD constraint is not fully fulfilled by the tower of axials.

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 . 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 entering the HLbL integral can be obtained by taking into account only intermediate states that are under good theoretical and numerical control. For the , due to the location of its pole, the form of this low-energy representation can be straightforwardly extended even to large and as long as stays small. Using available input for the -pole term, we find that the shift due to longitudinal SDCs on the isovector part of is in the range . 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 by in total. The quoted ranges encompass uncertainties in the low-energy input, perturbative 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 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 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 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 .