Longitudinal short-distance constraints for the hadronic light-by-light contribution to (g − 2)μ with large-Nc Regge models

While the low-energy part of the hadronic light-by-light (HLbL) tensor can be constrained from data using dispersion relations, for a full evaluation of its contribution to the anomalous magnetic moment of the muon (g − 2)μ also mixed- and high-energy regions need to be estimated. Both can be addressed within the operator product expansion (OPE), either for configurations where all photon virtualities become large or one of them remains finite. Imposing such short-distance constraints (SDCs) on the HLbL tensor is thus a major aspect of a model-independent approach towards HLbL scattering. Here, we focus on longitudinal SDCs, which concern the amplitudes containing the pseudoscalar-pole contributions from π0, η, η′. Since these conditions cannot be fulfilled by a finite number of pseudoscalar poles, we consider a tower of excited pseudoscalars, constraining their masses and transition form factors from Regge theory, the OPE, and phenomenology. Implementing a matching of the resulting expressions for the HLbL tensor onto the perturbative QCD quark loop, we are able to further constrain our calculation and significantly reduce its model dependence. We find that especially for the π0 the corresponding increase of the HLbL contribution is much smaller than previous prescriptions in the literature would imply. Overall, we estimate that longitudinal SDCs increase the HLbL contribution by ΔaμLSDC=136\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \varDelta {a}_{\mu}^{\mathrm{LSDC}}=13(6) $$\end{document}× 10-11. This number does not include the contribution from the charm quark, for which we find aμc−quark\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {a}_{\mu}^{c- quark} $$\end{document} = 3(1) × 10−11.


Introduction
Current Standard Model (SM) evaluations of the anomalous magnetic moment of the muon, a µ = (g − 2) µ /2, differ from the value measured at the Brookhaven National Laboratory [1] a exp µ = 116 592 089(63) × 10 −11 , (1.1) by around 3.5 σ. In the near future, the new Fermilab E989 experiment [2] will be able to reduce the experimental uncertainty by a factor 4, and the E34 experiment at J-PARC [3] will provide an important cross check, see ref. [4] for a comparison of the experimental methods. Therefore, the theoretical calculation of a µ needs to be improved accordingly. The uncertainty of the SM prediction mainly stems from hadronic contributions, such as hadronic vacuum polarization (HVP), see figure 1 (a), and HLbL scattering, see figure 1 (b). Since the HVP contribution can be systematically calculated with a data-driven dispersive approach [5][6][7][8][9], lattice QCD [10][11][12][13][14][15][16], and potentially be accessed independently by the proposed MUonE experiment [17,18], which aims to measure the space-like finestructure constant α(t) in elastic electron-muon scattering, the HLbL contribution may end up dominating the theoretical error. 1 Figure 2: Pseudoscalar-pole contribution to (g − 2) µ . The cyan dots indicate the TFF of the pseudoscalar meson.
Apart from lattice QCD [27][28][29], recent data-driven approaches towards HLbL scattering are again rooted in dispersion theory, either for the HLbL tensor [30][31][32][33][34][35], the Pauli form factor [36], or in terms of sum rules [37][38][39][40][41]. In all approaches, the most important HLbL contributions are the π 0 -pole and other pseudoscalar-meson-pole contributions, see figure 2. The strength of these pseudoscalar poles is determined by the transition form factors (TFFs), which in turn can be reconstructed from dispersion theory [42][43][44][45][46][47], leading to [46,47] a π 0 -pole µ = 62.6 +3.0 −2.5 × 10 −11 , (1.2) in agreement with determinations from lattice QCD [48], Canterbury approximants (CA) [49], and Dyson-Schwinger equations (DSE) [50]. Since the central value (1.2) is close to earlier model-based calculations, e.g., within lowest-meson-dominance+vector (LMD+V) models [51], the second-most important aspect of the dispersive approach apart from rigorous uncertainty estimates is the clear definition of the pseudoscalar intermediate states in terms of physical, on-shell form factors, in contrast to earlier notions of a pion-exchange contribution, see ref. [52], which involve the model-dependent concept of an off-shell pion. This becomes particularly important when combined with other intermediate states, ensuring that the pseudoscalar poles are consistent with, for instance, the dispersive definition of two-pion intermediate states [34,35], which in turn are determined by the corresponding on-shell quantities, in this case the helicity amplitudes for γ * γ * → ππ [53][54][55][56][57][58]. However, in contrast to HVP there is no closed formula that resums all possible intermediate states (in terms of the cross section for e + e − → hadrons [59,60]), in such a way that the consideration of exclusive channels will break down eventually, irrespective of the complications when extending the dispersive formalism to higher-multiplicity intermediate states. Therefore, to control the regions in the (g − 2) µ integral where either two or all three independent photon virtualities become large, additional constraints are required. In close analogy to HVP, where perturbative QCD (pQCD) becomes applicable in the high-energy tail of the dispersive integral, such constraints arise from the OPE and pQCD. In the regime where all three virtualities are large it was shown recently [61] that the pQCD quark loop indeed arises as the first term in a controlled OPE, with the next order suppressed by small quark masses and condensates. For the case in which one virtuality remains small, the leading OPE constraint was derived in ref. [62], by reducing the HLbL tensor in this limit to the triangle anomaly and its known non-renormalization theorems [63][64][65][66]. 2 The latter constraint decomposes into longitudinal and transversal contributions. As noted in ref. [62], the longitudinal part is intimately related to the pseudoscalar poles, but cannot be saturated by π 0 , η, η alone, nor by any finite number of poles. As a remedy it was suggested to drop the momentum dependence of the TFF at the vertex to which the external photon is attached, see figure 2, which leads to a substantial increase of the pseudoscalar-pole HLbL contribution. Based on an LMD+V model for the π 0 and vectormeson-dominance (VMD) models for η, η from ref. [51], this increase was found to be 13.5 × 10 −11 for the π 0 and 5 × 10 −11 each for η and η . This shift has been included, in one way or another, in subsequent estimates of the total HLbL value [52,67]. In fact, as we will show below, with modern input for the TFFs the corresponding increase would become even larger. While there is no doubt that the SDC is important-it is, in fact, one of the few constraints on the mixed-energy regions in which one photon virtuality remains small-modifying the expression for the pseudoscalar poles in this way is not compatible with the dispersive description of the four-point HLbL tensor [31][32][33][34][35] and spoils consistency with other intermediate states in the same framework.
In this work we address the question from a different standpoint: already in ref. [62] it was observed that while a finite number of poles cannot saturate the SDC, an infinite tower of them potentially can-dropping the TFF at the external vertex has in fact been described in ref. [62] as a model for the resummation of the tower of pseudoscalar states. Here we present explicit constructions that implement the Melnikov-Vainshtein (MV) constraint in terms of an infinite tower of excited pseudoscalars, constraining their properties from Regge theory, all available SDCs, and phenomenology wherever possible [68]. Given all these constraints the resulting models for the HLbL tensor prove remarkably rigid, without altering the low-energy properties. Since phenomenological information on excited pseudoscalars, especially their TFFs, is scarce, we do not attempt to construct TFF representations that apply for arbitrary kinematics, but concentrate on minimal models that cover the space-like region needed for (g − 2) µ and at the same time are able to fulfill all SDCs. Systematic uncertainties are estimated by comparing two such representations, either based on a truncated or untruncated Regge sum for the TFFs themselves, as well as the available phenomenological constraints, see appendices D and E. Moreover, our model is only needed for the low-energy part of the (g − 2) µ integral: above the energy where the matching occurs, we calculate the integral with the quark loop. This strategy also ensures that the estimate of the asymptotic region still applies in the chiral limit, in which the excited pseudoscalar states decouple, see section 5.3, while at low energies phenomenological input is needed either way to account for the effect of quark-mass corrections. All this leads to a more reliable estimate for the impact of the OPE constraints on the total HLbL contribution. To this end, we first review the expression for the pQCD quark loop and the known OPE constraints on the HLbL tensor in sections 2 and 3, adapting the conventions to the language suitable for the decomposition of the HLbL tensor from refs. [33,35], in which the expressions for both the pseudoscalar poles and the pQCD quark loop become remarkably simple. Next, we present in section 4 the explicit construction of large-N cthe discussion of the various SDCs in sections 2 and 3.
inspired Regge models 3 implementing the OPE constraints and derive the consequences for HLbL scattering and (g − 2) µ . In section 5, we match the resulting expressions for the HLbL tensor to the pQCD loop to obtain a first estimate of the scale where the description of the HLbL tensor in terms of hadronic intermediate states and its asymptotic properties should meet. A more detailed comparison to the results obtained with the MV model is provided in section 6, before we summarize our main results and discuss future developments in section 7. Technical details and alternative evaluations that are used to estimate the systematic uncertainty are collected in the appendices.

Pseudoscalar poles
The pseudoscalar poles only appear inΠ 1 (and by crossing symmetry inΠ 2,3 ) with F P γ * γ * (q 2 1 , q 2 2 ) the doubly-virtual TFF, F P γγ * (q 2 ) = F P γ * γ * (q 2 , 0) the singly-virtual TFF, and P = π 0 , η, η (see ref. [33] for the detailed derivation). The TFFs are normalized to the two-photon decays according to They are defined by the matrix element In the chiral limit, the non-singlet normalizations are determined by the Adler-Bell-Jackiw anomaly [73][74][75] with Gell-Mann matrices λ a , λ 0 = 2/3 1, 20) and decay constants defined through F a P : which is in general a 3 × 3 matrix. Ignoring for simplicity any possible mixing between the π 0 and the other two states, this takes the form: where, after the second equality sign, we have already introduced the standard two-angle mixing scheme between η and η . For the pion a = 3, the corresponding low-energy theorem is very close to phenomenology, while for η, η chiral corrections and mixing effects need to be taken into account. In particular, we stress that due to the renormalization of the singlet current F 0 P is not actually an observable quantity, and the corresponding α s corrections [76,77] need to be considered when relating the normalization, asymptotic constraints, and ηη mixing parameters [78][79][80][81][82]. In the present work, we will take the η, η normalizations from experiment, so that the singlet corrections become most relevant when comparing the asymptotic constraints and η-η mixing parameters in different schemes. As described in section 4.2, we studied the impact of different such determinations on the numerics, with the result that the corresponding variations were numerically irrelevant in view of the accuracy anticipated for the pseudoscalar TFF models discussed in the following sections.
In addition to the normalizations, the pseudoscalar TFFs are subject to the (leading) asymptotic constraint [83][84][85] which for the asymptotic wave function φ a P (x) = 6x(1−x), and again ignoring α s corrections for the singlet component, produces the limits 4 (2.25) Figure 3: Quark-loop contribution to HLbL scattering.
In view of (2.19), multiplying these limits by F P γγ and summing over P one obtains an expression which depends neither on decay constants nor on mixing angles. Moreover, the block form of the matrix (2.22) leads to two separate combinations with such a property: and similarly for the asymptotic limit of the singly-virtual TFF. Beyond the singlet α s corrections that describe the scale dependence of F 0 P [78][79][80], there are genuine pQCD corrections to the TFFs suppressed by α s at scales related to the photon virtualities [87]. The impact of such next-to-leading-order pQCD corrections was studied in ref. [47] in the context of the pion TFF, with the result that even for the ground-state pion the effect is small and safely covered by the uncertainty estimated from the onset of the asymptotic region.

The perturbative QCD quark loop
The quark-loop contribution to HLbL scattering is shown in figure 3, indicating the different permutations that need to be considered. Compact expressions for the BTT scalar functions can be obtained as follows: we use a Feynman parameterization for the loop integrals and project the result onto the scalar basis functionsΠ i [33,35]. We find all necessary BTT functions Π i in the limit q 4 → 0 by taking this limit in the appropriate order, so that the Tarrach poles drop out. Then we combine the functions Π i according to (2.7) to obtain the scalar functionsΠ i . Due to the limit q 4 → 0, one integral can be carried out and we are left with a two-dimensional Feynman-parameter integral. The result for the integrands contains spurious kinematic singularities, but the residues of these poles vanish when the Feynman integrals are carried out. Therefore, we can subtract these poles and obtain a representation that is manifestly free of kinematic singularitieŝ where and In principle, it is also possible to extract the results by projecting onto the singly-on-shell basisΠ i [35]. However, it turns out that this method is less straightforward, because different spurious kinematic singularities appear, which have to be subtracted again and make the calculation more complicated. As a cross check of (2.28) we have evaluated light-quark loops for q = u, d, s with (constituent) quark mass m q , including a factor N c q=u,d,s Q 4 q = 2/3 as well as the lepton loops. The latter agree well with the known analytic expressions [88], while apart from the electron loop the results are well reproduced from the heavy-mass expansion [89]. Throughout, for the matching to our Regge models in section 5, we use the pQCD quark loop with m q = 0, given that even in configurations where chiral corrections for the light quarks q = u, d, s can be controlled within pQCD, they only enter at subleading orders.
As a first application, we consider the contribution from the charm quark. Assuming that this contribution is fully perturbative, with mass m c = 1.27(2) GeV [90], the quark loop evaluates to a c-loop µ = 3.1(1) × 10 −11 . In analogy to the light quarks, one would expect the most important non-perturbative effect to be related to the pole contribution from the lowest-lying cc resonance, the η c (1S) with mass m ηc(1S) = 2.9839(5) GeV and twophoton width Γ(η c (1S) → γγ) = 5.0(4) keV [90]. Using a VMD-type form factor with scale set by the J/Ψ (as suggested by a significant branching fraction BR(J/Ψ → η c (1S)γ) = 1.7(4)% [90]), this leads to the estimate a ηc(1S) µ = 0.8 × 10 −11 (this estimate agrees with the LMD result a ηc(1S) µ = 0.9(1) × 10 −11 from [91]). Given the relatively low scale set by m c one may also expect α s corrections in a similar ballpark. Altogether, we estimate while the b-quark contribution is already suppressed to the level of 10 −13 and the t-quark loop to 10 −15 .

OPE for the asymptotic region
The first term in the OPE for the kinematic configuration in which all three momenta are large coincides with the pQCD quark loop. This has long been suspected in the literature, including ref. [62], but was only demonstrated recently in ref. [61], by working out the next order in the expansion. While at leading order all quark masses can simply be put to zero, this is no longer true at subleading orders. In fact, it is the presence of quark masses and condensates that numerically suppresses the next-to-leading order corrections. In the limit q 2 1 = q 2 2 = q 2 3 ≡ q 2 the expressions for the pQCD quark loop simplify tô where the Clausen function is defined as This result again includes the factor 2/3 due to N c and quark charges, after summing over q = u, d, s.

OPE for the mixed regions
The OPE constraint derived in ref. [62] applies to the case where one virtuality remains smaller than the others, Q 2 3 Q 2 1 ∼ Q 2 2 , also referred to as the mixed regions. This constraint traces back to non-renormalization theorems for the V V A correlator [63,64], which had been used before in the context of the electroweak contributions to (g −2) µ [65,66]. Explicit pQCD calculations at two-and three-loop order exist [92,93], but the main argument in ref. [62] was that the non-renormalization theorems allow one to address the regions in which both perturbative and non-perturbative aspects might be important. We first review this derivation, while casting the results in a form suitable for the BTT decomposition of the HLbL tensor.
The central object is the OPE for two electromagnetic currents: We consider large momentaq = (q 1 − q 2 )/2 flowing through the currents and expand the operator product into a series of local operators. For |q| Λ QCD , the coefficients can be calculated in perturbation theory. At leading order in α s , only two-quark operators are generated, hence the matching can be easily obtained by inserting the operator (3.3) into external quark states and expanding the diagrams for large momentaq: with 0123 = +1. Introducing the scalar density S(x), the axial vector j µ 5 (x), and the energy-momentum tensor θ µν (x) with flavors weighted by the squared electric charges, with the quark-mass matrix M = diag(m u , m d , m s ) as well as the derivative ∂ − = − → ∂ − ← − ∂ , we read off the matching for the OPE: The first term reproduces the expansion given in ref. [62], but differs in sign just because of different conventions (they use 0123 = −1).
Applying the OPE to the HLbL tensor in the limit Q 2 1 ∼ Q 2 2 Q 2 3 , Q 2 4 we then find at leading order where the correlator W µνρ is defined as Introducing the vector and axial-vector currents where {λ a , λ b } = 4/3δ ab + 2d abc λ c , we also define the correlator Performing the flavor decompositions with C a defined in (2.19) and D a = 1 2 Tr(Qλ a ), we write the correlator (3.9) as (3.14) The Lorentz decomposition of the V V A correlator is chosen as [64] with the following Lorentz structures: In the massless limit, one finds at one loop the contribution of the axial anomaly [73,74] w (a) (3.17) For the non-singlet contributions a = 3, 8, this result is modified neither by higher-order perturbative [94] nor non-perturbative contributions [95], while the singlet contribution is affected by the gluonic U (1) anomaly.
In the chiral limit, the factor C 2 a in (3.13) arises naturally due to the flavor decomposition. One factor of C a stems from (3.12), the second factor emerges as follows. We consider singlet and octet parts of the axial current and define W bc0 µνρ (q 1 , q 2 ) =: which implies The coefficients can be simplified to hence both singlet and octet components lead to another factor C a . In pQCD and in the chiral limit, the following non-renormalization theorems were derived in [64] for the non-singlet part of the axial current: The transversal functions in these relations are subject to non-perturbative corrections.
In the following, we will use the OPE constraints as they arise at leading order and in the chiral limit. Both the anomaly constraint (3.17) and the non-renormalization theorems (3.21) receive quark-mass corrections [96].

Projection onto BTT
In this section, we derive the asymptotic constraints that the leading-order expression (3.8) of the OPE imposes on the scalar BTT functionsΠ i (2.7) entering the master formula for a µ . One might be tempted to simply project the OPE expression (3.8) onto the BTT scalar functions. However, there are several problems with such an approach. First of all, the leading-order expression of the OPE is not manifestly gauge invariant: the contraction with (q 1 − q 2 ) µ vanishes, but the one with (q 1 + q 2 ) µ does not. Due to q 1 = −q 2 + O(1) this does ensure gauge invariance at O(1/q), while for the subleading orders relations with the matrix elements of the energy-momentum tensor are needed to restore gauge invariance. At leading order, gauge invariance could be restored by applying a gauge projector which does not alter the O(q −1 ) terms of the OPE expression. The subsequent projection onto BTT and extraction of the scalar functionsΠ i could then be performed immediately. However, this procedure is not uniquely defined: the BTT structures themselves become degenerate depending on the order of the expansion for largeq. This implies that the leading-order OPE only constrains certain linear combinations of scalar functionsΠ i . In an assignment of these constraints to individual scalar functions ambiguities are introduced. For the longitudinal amplitudes the linear combination of scalar functions that is uniquely constrained happens to coincide withΠ 1-3 , but for the transversal amplitudes the situation becomes more complicated. We proceed as follows in order to determine this ambiguity explicitly and to work out the exact form of the OPE constraint at the level of BTT functions. First, we remember that the HLbL tensor is linear in the external momentum q 4 . Due to the relation following from gauge invariance, it is enough to consider the derivative with respect to q ρ 4 and then take the limit q 4 → 0, as required for (g − 2) µ kinematics. The BTT functionsΠ i in this limit are unambiguously defined, hence they have their own proper expansion in 1/q, which we would like to constrain using the OPE. The derivatives of the Lorentz structureŝ T i multiplying the functionsΠ i in the tensor decomposition (2.5) however contain several terms with different scaling for largeq. For instance, for the tensor structureT µνλσ 5 , one finds where the first term scales as O(q 0 ), the second and third terms are of O(q), and the last term is of O(q 2 ). This illustrates that a certain coefficient of the expansion in 1/q of a scalar functionΠ i can contribute to different orders in the expansion in 1/q of the full HLbL tensor, i.e., to different orders of the OPE. Vice versa, in order to determine the leading-order OPE constraint on the BTT functions, we e.g. have to consider terms up to and including O(q −3 ) inΠ 5 . We now write the scalar functionsΠ i as a generic expansion in 1/q and sum up the scalar functions times (derivatives of) tensor structures. Collecting in the resulting tensor terms according to the scaling withq and requiring equality with the leading-order OPE limit (3.8) determines the expansion coefficients of the 19 BTT functions relevant in the (g − 2) µ limit aŝ (3.26) and the remaining ambiguities are parameterized by functions c (n) i behaving as c (n) i ∼ 1/q n , which are subject to certain crossing-symmetry relations following from (2.8) and (2.10). Note that the small dimensional quantity that makes the expansion parameter dimensionless can be any of the small scales, e.g., the small momentum or matrix elements of the operators in (3.7). Due to the scaling of the tensor structures, the neglected terms in (3.25) affect the HLbL tensor first at O(1/q 2 ) and therefore cannot interfere with the leading-order OPE. This result specifies the configuration Q 2 1 ∼ Q 2 2 ≡ −q 2 = −q 2 Q 2 3 . The related limits for small q 2 1 or q 2 2 follow directly from crossing symmetry. Since the longitudinal amplitude w L only contributes toΠ 1-3 , we will refer to these scalar functions as the longitudinal ones, and accordingly to the remainingΠ i as the transversal contribution. The non-trivial constraint on the non-singlet part of the latter emerges from the corresponding limit of (3.21) (3.27) but in contrast to the anomaly condition (3.17), which is exact in the chiral limit, this relation does receive non-perturbative corrections. As noted above, the projection (3.25) shows that only the OPE constraints on the longitudinal amplitudes are free from ambiguities, whereas all those on the transversal ones are affected by them. The presence of these ambiguities is not a problem per se: it simply means that at leading order the OPE only constrains certain linear combinations of BTT functions. We also note that the ambiguities would be moved to higher orders if the next terms in the OPE were included.
For asymptotic values of q 2 3 , the OPE constraints can be compared with the pQCD quark loop evaluated in the chiral limit and for q 2 These expressions perfectly agree with (3.25) if we use the non-renormalization theorem in the form (3.27) and set the ambiguities c (n) i to zero, which demonstrates that the OPE constraint and the pQCD quark loop coincide in the appropriate kinematic limit [97] (neither chiral effects nor α s corrections related to the gluon anomaly in the singlet channel matter in this limit). We stress that one could impose the OPE constraints on the transversal functions without having to deal with these ambiguities by first building linear combinations of the BTT functions that are free from them. In principle, one could even use the freedom in the projection at a given order to simplify expressions, e.g., at leading order one could choose the c (n) i in such a way that the only non-vanishing contribution arises inΠ 50 =Π 51 = 2/q 2 w + T +w − T (q 2 3 , 0, q 2 3 )f (q 2 ). However, such a simplification would no longer hold at subleading orders, therefore, we keep here the general form (3.25) that shows directly how the OPE limit corresponds to the pQCD quark loop (3.28) evaluated in the same kinematics. In the following we will focus on the OPE constraint on the longitudinal contribution, which can be unambiguously assigned to the BTT functionsΠ 1−3 already at leading order.

Relation to pseudoscalar poles
Separating the longitudinal OPE constraint into flavor components one findŝ (3.29) (3.30) which due to (2.26) matches precisely onto (2.16) when the meson masses and, crucially, the momentum dependence of the singly-virtual form factor are neglected. This is the basic premise of the model suggested in ref. [62]. We stress that for the non-singlet component these relations are exact in the chiral limit, see section 5.3 for an extensive discussion of this point. For non-vanishing quark masses and, in the case of the singlet, due to the gluon anomaly they do receive corrections. However, at low energies, where such corrections matter most, we always use the full dispersive result that automatically corresponds to physical quantities, while (3.29) and (3.30) are only implemented for asymptotic values of the virtualities q 2 i . The OPE constraint becomes potentially valuable in the context of the mixed-energy regions, where both a description in terms of hadronic intermediate states and pQCD have limited applicability. In practice, the constraint is rigorous once all momenta are large compared to Λ QCD to ensure that quark-mass corrections can be neglected. The Regge approach in the next section is our proposal for an explicit implementation of the OPE constraint, following a remark made in ref. [62]: while the 1/q 2 3 behavior in (3.29) and (3.30) cannot be obtained with any finite number of pseudoscalar poles, an infinite sum over excited states can produce the required asymptotics.

Regge models for the pseudoscalar-pole contribution
Assuming confinement, in the large-N c limit of QCD [98] the spectrum of the theory in any sector (set of quantum numbers) reduces to an infinite number of narrow resonances. One should not expect the spectral functions in this limit to be close to those of QCD with N c = 3 locally: a series of δ-functions does not look like the continuum observed in nature for any spectral function. On the other hand, one expects the large-N c limit to provide a good approximation to QCD on average, and in particular to reproduce to a reasonable accuracy its global properties such as asymptotic limits. There is a vast literature on the subject that shows that these theoretical considerations can be used with good success to build large-N c models that simultaneously satisfy low-and high-energy constraints [99][100][101][102][103][104].
The aim of the present section is to construct a large-N c Regge model in the pseudoscalar and vector-meson sectors of QCD that allows us to satisfy the SDCs discussed above via an infinite tower of pseudoscalar-pole contributions. The logic we follow in the construction of the model is very simple: we seek minimal models, in terms of algebraic form and number of free parameters, that are able to satisfy all known constraints, both of experimental as well as of theoretical nature, i.e., phenomenological constraints wherever available and all known high-and low-energy limits. Accordingly, we construct these large-N c Regge models with the application to HLbL scattering in mind and thus work with physical quark masses. We will comment on the chiral limit and the potential role of axial-vector resonances in section 5.3.

Large-N c Regge model for the pion transition form factor
The pion TFF describes the transition of a pion into two photons. VMD, LMD, and LMD+V models for the pseudoscalar TFFs are widely used [51,105], cf. figures 4 and 8. In this work, we use an untruncated large-N c model for the TFF, in which the pion couples to the photons through a tower of isovector, I G = 1 + , and a tower of isoscalar, I G = 0 − , vector mesons, J P C = 1 −− , e.g., the ρ and ω, respectively. Here, a tower of ρ (ω) mesons means an infinite sum over radially-excited ρ (ω) mesons. The contributions from a φ instead of an ω are subdominant, see appendix A, and thus will be neglected for the pion. The standard large-N c ansatz for the pion TFF (see refs. [106,107]) reads:

2)
F Vρ and F Vω , represented by blue dots in figure 4, are the current-vector-meson couplings and G πVρVω , the cyan dot in figure 4, is the coupling of two vector mesons to the neutral pion. We stress that the couplings in (4.1) are Q 2 independent as required by the large-N c approach (combined with analyticity): the contribution of a given intermediate state is fixed by its imaginary part, which for narrow resonances is a δ-function, which freezes any Q 2 dependence. Indeed the latter could be interpreted as coming from the continuum between resonances, which is suppressed in the large-N c limit. Another potential source of Q 2dependent corrections to (4.1) is related to subtraction terms: while (4.1) follows from an unsubtracted double-spectral representation for the TFF, introducing subtractions would produce single-propagator terms and, eventually, a polynomial. However, for a δ-function subtractions are not necessary, and even before taking the large-N c limit the pQCD behavior of the TFF implies that an unsubtracted representation holds. As argued in ref. [47], the advantage of using an unsubtracted dispersion relation, in favor of a subtracted variant that could suppress some of the high-energy input, is precisely that it allows one to manifestly incorporate the correct pQCD asymptotics. The large-N c ansatz (4.1) corresponds to this scenario. In the following, the vector-meson spectra are assumed to obey a radial Regge model, see figure 5: where σ ρ and σ ω are the slope parameters of the Regge trajectories, n ρ and n ω are radial excitation numbers, and the ground-state masses are M ρ = M ρ(770) = 775.26 (25) MeV and M ω = M ω(782) = 782.65 (12) MeV [90]. The green dotted lines with slope σ 2 = 1.11 GeV 2 are based on the lattice calculation of ref. [110]. The orange dot-dashed lines with slope σ 2 = 1.87 GeV 2 are based on the ρ → 2π decay [106,111].
satisfied after resumming the series of vector resonances. To this end, the coupling constants must be arranged in such a way that the Q −4 behavior of the individual terms becomes a Q −2 behavior after resummation. That this is possible was shown in refs. [106,107]. In addition, the pion TFF has been measured quite well in the singly-virtual case [115][116][117][118], and our model for the TFF would have to describe the data. For the doubly-virtual case a recent dispersive analysis has shown that data for related processes and theoretical arguments constrain the behavior of the TFF in that kinematical region [44,46,47]-a constraint we will also take into account.
Imposing all these constraints on the model (4.1) by adjusting its free parameters is technically cumbersome, especially if we consider that we must still add a third sum over the tower of pseudoscalar mesons, J P C = 0 −+ , cf. figure 2, with which we aim to change the large-Q 2 behavior of the whole HLbL tensor. In particular, we will implement the SDCs onΠ 1 introduced in (3.1) and (3.29): 5 SDC for the mixed region [62]: lim SDC for the asymptotic region: lim Here, F π(n)γ * γ * is the TFF of the n-th radially-excited pion and a radial Regge model is assumed for the pion masses starting from the first excitation, see figure 6: where M π = 134.9770(5) MeV is the π 0 mass [90]. 6 Given the complexity of implementing all these constraints simultaneously in terms of the general couplings of the Regge model, we therefore adopted a different approach: 1. we allow the ground-state pion to couple only to the ground-state ρ and ω mesons, and the n-th pion excitation to couple only to the n-th ρ and ω excitations; 2. we subsume the effect of the vector-meson excitations that we have just eliminated into a Q 2 i dependence of the numerator multiplying the resonance propagators; 3. the latter Q 2 i dependence will be parameterized in simple terms with as few free parameters as necessary to satisfy the constraints listed above. 5 Note that while for the MV SDC, which is derived based on the V V A triangle, the flavor decomposition into pion, η, and η is unambiguously given by C 2 a , see (3.29), the decomposition presented here for the SDC in the asymptotic region (3.1) is not unique. We choose to adopt the same separation as for the MV constraint. 6 Note that the ground-state is treated separately because for the Goldstone bosons a strong non-linearity of the Regge trajectory is expected [119].
The first step is motivated by the fact that non-diagonal couplings are suppressed by the reduced overlap of radial wave functions with different numbers of nodes [106]. For the same reason we are only considering the leading S-wave vector-meson trajectories and neglecting the D-wave daughter trajectories. In appendix B we will consider an alternative model that already for the pQCD limit of the TFF itself, cf. (4.6), uses the Regge resummation from ref. [106], but for the main text we restrict the presentation to the most economical form sufficient to fulfill all constraints simultaneously. This strategy leads us to and The second mass scale M diag is determined by fitting the experimental data, it parameterizes the doubly-virtual behavior of the TFF. With this parameterization, the three conditions for the TFF of the ground-state pion from (4.4)-(4.6) can be expressed as follows: anomaly: symmetric pQCD limit: (4.14) Since the mass scales M ρ , M ω , M +, 0 , and Λ as well as π 2 F π are of about the same order, all coupling constants that appear in the constraint equations (4.12)-(4.14) multiplied by ratios of these mass scales are expected to be of O(1). M −, 0 , on the other hand, is much smaller, so that the coupling constants multiplied by it (c BL and c B ) are expected to be of O(Λ 2 /M 2 −, 0 ) ∼ 100, otherwise their role in the equations would become irrelevant. Of course these three conditions are not sufficient to determine all five model parameters. Two more constraints follow from resumming the contributions of all excited pseudoscalars to the HLbL amplitude, see (4.7) and (4.8). Details on the evaluation of infinite sums over rational functions can be found in appendix C. The MV SDC for the mixed region translates into:  Table 1: "Natural" model parameters of the large-N c Regge models for the pion, η, and η TFFs. Note that here we rescaled the η ( ) parameters c A , c B , and c BL with a factor of where c diag from (4.14) has been used and The second SDC concerns the limit Q 2 i = Q 2 → ∞, for all i = 1, 2, 3. It also involves c A and c B , but now in a different combination together with c diag : with In appendix F.1, this system of equations is solved analytically. Here, we discuss numerical values for all parameters, also summarized in table 1, based on the following choice of Regge slopes [108]: 7 Furthermore, we use F π = 92.28 MeV, Λ = 1 GeV, and other input from the PDG [90]. 7 From figure 5 of ref. [108] we extractedMπ = 766 MeV.   [46,47]. The data are from CELLO [115], CLEO [116], BaBar [117], and Belle [118].
The constant c diag is independent of all the others and is directly determined by (4.14): Once this is fixed, equations (4.15) and (4.17) determine c A and c B . Since the second equation is quadratic it has two solutions, but one can be readily discarded because the two-photon couplings of the excited pions become unreasonably large, and so do the values of the constants c A and c B . The physical solution gives: Having determined c A and c B , (4.12) determines c anom to the value: and finally (4.13) fixes the remaining parameter: 47]. To find the best fit, we minimize the estimated variance: where j max is the length of the data set and p is the number of fit parameters. Here, f is our model and f data is the dispersive TFF. 8 The sum is over The singly-virtual TFF of the ground-state pion is shown in figure 7. The large-N c Regge model presented above is labeled as "Model 1." In appendix B, an alternative TFF model, to which we refer as "Model 2," is introduced, based on a Regge resummation for the TFF itself. Both models give a reasonable description of the experimental data, while Model 1 shows better agreement with the dispersive TFF in the intermediate-Q region. In appendix F.2, both models are shown also in the doubly-virtual region and further compared to the dispersive TFF [44,46,47], a prediction from lattice QCD [48], and a result from DSE [50]. We stress that neither model should be evaluated for other than purely space-like virtualities, both are constructed in such a way as to provide an efficient implementation of all constraints relevant for the space-like region, but do not properly incorporate the analytic structure required to continue to time-like virtualities. In addition to our fits to the dispersive π 0 TFF, we also checked that the π 0 contribution to (g − 2) µ is reproduced correctly Finally, as detailed in ref. [47], effective-field-theory constraints on the pseudoscalar-pole contributions [120,121] are automatically encoded in the TFF phenomenology, for the leading constraint in its normalization, for the subleading one in the momentum dependence.

Large-N c Regge model for the η and η transition form factors
Analogously to the pion case, our large-N c Regge model for the η and η TFFs shall satisfy the following five low-and high-energy constraints, cf. BL limit [84,85]: lim symmetric pQCD limit [114]: lim SDC for the mixed region [62]: lim SDC for the asymptotic region: lim Here, as compared to (2.25), we now use the notation Furthermore, we introduced: as follow by separating the η and η contributions to (2.26) according to (2.22). These coefficients fulfill C 2 η + C 2 η = C 2 0 + C 2 8 = 1/12. Switching to the η and η we face the problem that, since these are I = 0 mesons, they couple to isovector-isovector and isoscalar-isoscalar vector mesons, so to same-mass vector mesons only (ignoring φ-ω mixing), see figure 8. Taking the limit M ω(n) = M ρ(n) = M V (n) in our parameterization of the pion TFF (4.10), we obtain a significant simplification: Table 2: Pseudoscalar-vector-vector couplings derived in appendix A: C P V 1 V 2 with P = π, η, η and V i = ρ, ω, φ as defined in (A.9) and (A.10).
Since two free parameters dropped out, this parameterization cannot satisfy all relevant low-and high-energy constraints.
Fortunately, via vector-meson mixing in the isoscalar sector, there is a possible contribution of a mixed φ-ω term to the TFFs of the η ( ) , which would be absent in the case of ideal mixing. The φ-ω coupling to η ( ) will certainly be small when compared to the same-mass vector-meson couplings, see table 2, but since it contributes where the others cannot, it is important to retain.
In summary, our large-N c Regge model for the η ( ) TFFs reads: where the two parts parameterize the same-mass and mixed vector-meson contributions, respectively: . (4.36) Here we again use the short-hand notations from (4.11), with the modification that All meson spectra are assumed to follow a radial Regge model. For the φ meson, we use the analog of (4.3). For the η and η mesons, we distinguish: ) (GeV -1 ) Figure 9: Comparison to the doubly-virtual η TFF data from BaBar [122]. The large-N c Regge model, "Model 1" (4.34), is indicated by the pink bands with the dashed lines. Our alternative TFF model, "Model 2" (B.1), is indicated by the solid cyan lines. and with the ground-state masses M η = 547.862 (17) MeV and M η = 957.78 (6) MeV [90].
The normalization coefficient is defined as: where C P V 1 V 2 are the pseudoscalar-vector-vector couplings derived in appendix A, and C η φω is the parameter that measures the deviation from ideal mixing. By construction, each vector-meson pair contributes (up to normalization) exactly C P V 1 V 2 to F η ( ) γγ . To simplify the parameterization, equation (4.36) only contains terms which are unique to the φ-ω contribution, and the n-dependence has been removed from the numerator of (4.35). In this way, (4.36) is used to satisfy the BL limit for the ground-state η ( ) TFF as well as the two SDCs on the HLbL tensor.
The constraint equations following from (4.27) and (4.28) read: symmetric pQCD limit: where the same equations hold for the η with the obvious replacements (including C 8 → C 0 ). The MV SDC for the HLbL tensor in the mixed region translates to: where c diag has already been inserted. The SDC for the HLbL tensor in the asymptotic region becomes more complicated due to the presence of additional mass scales: In appendix G.1, the above system of equations is solved analytically. In the following, we discuss numerical values for all input parameters. The couplings C η ( ) V 1 V 2 are collected in table 2. They are calculated based on (A.9) with the phenomenological η-η mixing parameters [123,124]: 45) and the φ-ω mixing angle θ V = 36.4 • [90]. The parameters C 2 η ( ) , which describe our choice for the splitting of the SDCs on the HLbL tensor into η and η contributions, evaluate to: as follows from (4.32) with the η-η mixing parameters in (4.45). The decay constants F η ( ) , on the other hand, are not deduced from the η-η mixing parameters, but fit to experimental data for the singly-virtual η ( ) TFFs: with an estimated variance of χ 2 ∼ 1.1 and χ 2 ∼ 0.9, respectively. The errors are increased in order to cover the Padé approximant predictions from refs. [49] and [81] for η and η , respectively. The large error on F η may be partly due to the fact that it is not clear when the asymptotic BL limit sets in, accordingly, we will keep the full range in the error analysis. M diag is fit to the recent BaBar data for the doubly-virtual η TFF [122], see figure 9. The resulting value M diag = 898 MeV (with χ 2 ∼ 1.6) is used for both the η and η large-N c Regge model. Furthermore, we use the Regge slopes collected in (4.19) as well as [108]: 9 The final model parameters (c A , c B , c diag , c BL ) are summarized in table 1, where we rescaled the numerical values with C η φω /N ∼ 0.0129 and C η φω /N ∼ −0.0017, respectively, to show that all parameters are of "natural" size.
The singly-virtual TFFs of the ground-state η and η are shown in figures 10 and 11. The large-N c Regge model presented above is labeled as "Model 1." The alternative TFF model, introduced in appendix B, is referred to as "Model 2." The model error of the large-N c Regge TFFs is propagated from the errors of the input parameters σ P , σ V , F P γγ , F P , as well as the η-η mixing parameters, see (4.19), (4.26), (4.48), (4.45), and (4.47). While Model 2, for which we do not provide an error estimate, runs outside the error band of the Model 1 η (η) TFF for some low (intermediate) values of Q 2 , both models give a good description of the experimental data and, thus, come out close to the results from CAs [49] and fits within resonance chiral symmetric theory (RCST) [125]. 10 In appendices G.2 and G.3, both models are further compared to CA, RCST, and DSE [50], and the decomposition of Model 1 into 2ρ, 2φ, 2ω, and φω contributions is illustrated. We stress again that neither model should be evaluated for other than purely space-like virtualities.

Comparison of two-photon couplings
Apart from the Regge slopes for the trajectories of pion, η, η , as well the vector mesons, phenomenological input for the excited-pseudoscalar contributions could in principle be provided by their TFFs. Even though the normalizations are poorly known, it is still important to verify that the two-photon couplings implied by our Regge models compare  [49]. The dark blue dotted curve is the RCST result from ref. [125]. The data are from L3 [127], CELLO [115], CLEO [116], and BaBar [126].
reasonably to the available phenomenological constraints. For the first excited state in the pion trajectory, there is a limit see appendix D, which is indeed satisfied by F π(1300)γγ = 0.050 GeV −1 from our Regge model. Nothing is known about the two-photon coupling of the π(1800) and even heavier excited pions. The situation is more involved in the η ( ) sector. As alluded to in the caption of figure 6, the spectroscopy of excited η ( ) states is contentious, especially regarding the role of the states below 1500 MeV. Table 3 collects two possible assignments of states to Regge trajectories. The first interpretation, favored by ref. [90], considers the η(1295) the lowest η excitation and differentiates between η(1405) and η(1475) states. The latter is considered as the first η excitation, while the η(1405) is described as a glueball candidate. In contrast, ref. [113] argues that there is only a single state below 1500 MeV, the η(1440), which should be interpreted as the first η excitation. The X(1835) is identified as suitable candidate for the first η excitation, although its quantum numbers are not yet established. In both cases, the η(1760) emerges as the second η excitation.
The available constraints on the two-photon couplings of η ( ) are collected in table 4, see appendix D for details. The results from the second column are valid under the assumption Table 4: Constraints on the two-photon couplings of excited η ( ) states, as collected in appendix D. The column labeled "direct" includes constraints that follow directly from branching fractions listed in the PDG, while the column labeled "indirect" lists results obtained when assuming that the channels from the last column are dominant.
that the branching fractions listed in the PDG are accurate, while the third column assumes, in addition, dominance by some decay channels (given in the last column). For the η(1295) only an indirect limit is available, in the first assignment the two-photon coupling of the η(1295) comes out slightly larger. Note, however, that the very existence of the η(1295) is called into question in ref. [113], with the fact that in contrast to the η(1475) this resonance has not been seen in the γγ reaction as one of the arguments. The η(1405) is discarded in either assignment of Regge trajectories. However, in the second assignment the η(1440) would be interpreted as a single state instead of η(1405) and η(1475), see also ref. [129], in such a way that for the comparison the two-photon couplings of both states should be considered. Remarkably, the measured value for the η(1475), F η(1475)γγ = 0.041(6) GeV −1 , agrees perfectly with F η(1440)γγ = 0.035 GeV −1 from Assignment 2. In Assignment 1, where the η(1475) is considered the first η excitation, there is still reasonable agreement. Next, the experimental result for the η(1760) nicely confirms the two-photon coupling implied by both assignments, since a tiny correction beyond the dominant η ππ channel would suffice to bring the numbers into complete agreement. Finally, the two-photon coupling of the X(1835) in Assignment 2 fulfills the direct limit but not the one assuming dominance of η ππ, which may indicate that in case this assignment is correct, other channels besides η ππ may play a role (as indeed suggested by other decay channels listed in the PDG). Moreover, the significance of the two-resonance fit from ref. [130] used to obtain the much stricter limit is only quoted at 2.8 σ. Taken together with the fact that not even the quantum numbers of the X(1835) are firmly established, it thus seems difficult to draw meaningful conclusions on its two-photon coupling at this point. Altogether, we conclude that the two-photon couplings implied by our large-N c Regge models are well compatible with the phenomenological constraints. In particular, for the cases where measurements and not just limits exist, the η(1475) and the η(1760), the resulting couplings are close to the ones that our large-N c Regge models would imply. The same is true for our alternative TFF model (B.1), see figure 19, whose couplings are similar to the ones of the large-N c Regge models. We stress that the detailed comparison depends on the assignment of observed states to Regge trajectories, but in both variants considered there is reasonable agreement with the two-photon phenomenology of excited η ( ) states.

Excited-pseudoscalar contributions to (g − 2) µ
The ground-state pseudoscalar-pole contributions to (g − 2) µ , calculated based on our TFF models, are given in (4.25), (4.49), and (4.50). The uncertainty on the predictions from Model 1 is the propagated error from the input parameters σ P , σ V , F P γγ , F P , and the η-η mixing parameters. In all cases we observe good agreement with the literature [46,47,49,128], which demonstrates that in addition to fulfilling the various SDCs, our Regge models capture the properties of the TFFs most relevant for the g − 2 integral.
In the following, we derive the contribution to (g −2) µ originating from radially-excited pseudoscalar mesons. The large-N c Regge models introduced in the preceding sections and the alternative model discussed in appendix B are constructed in such a way as to describe not only the ground-state pseudoscalar TFFs, but also the TFFs of excited pseudoscalar mesons. Phenomenological input on these excited states enters mainly in terms of their masses as contained in the Regge parameters, while the infinite sum restores the correct asymptotic properties of the HLbL tensor, which cannot be achieved with a finite number of pseudoscalar-pole contributions. Moreover, for some of the excited states limits on their two-photon couplings are available, see appendix D as well as the discussion in the previous subsection, which shows that the couplings implied by our Regge model are consistent with the available constraints from phenomenology.
With the large-N c Regge model, we can calculate the pseudoscalar-meson tower exactly, i.e., we can perform the infinite sum over pseudoscalar-pole diagrams with excited pseudoscalars. For Model 2, we sum over the lowest n = 100 radially-excited pseudoscalars (P = π, η, η ) and then fit a saturation curve, in order to extrapolate to infinity. Here, we defined: with the infinite-summation result denoted as: The saturation curve procedure is illustrated in figure 12, where for reasons of clarity only every other data point is plotted above n = 4. The fits start from n max = 1 and describe the data perfectly. The dotted lines indicate the extrapolated values for ∆a P -poles µ and illustrate the good convergence of the summation already at n max = 100. This procedure has been verified with the large-N c Regge model, for which the sum is already saturated at n max = 100.
For the full pseudoscalar-pole contributions to (g − 2) µ , we obtain:  (4.57) With the alternative assignment of η ( ) excitations in the radial Regge trajectories [113], see dot-dashed purple lines in figure 6, we obtain: indicating that the net effect is remarkably insensitive to the assignment of the η, η Regge trajectories. Expressing C 2 η ( ) through the experimental F η ( ) γγ and F η ( ) , see left-hand side of (2.26), instead of the η-η mixing parameters, leads to a decrease of a η ( ) -poles µ that is well within the uncertainty quoted in (4.57). Our final result for the sum of pion, η, and η states is:   The simplest and most instructive matching to the massless pQCD quark loop proceeds at the level of the (g − 2) µ integral. The asymptotic pQCD region where all Q i are large can be captured by imposing the condition that all Q i be larger than Q min . To be able to add the mixed regions, where one virtuality is smaller than Q min , in the quark-loop integration, one needs to dampen the contribution in the additional integration region, since it is already partly covered by the ground-state pseudoscalar poles. To this end, we introduce the suppression factor Q 2 /(Q 2 + Λ 2 ) for the virtuality Q < Q min , while retaining the cut that at least two Q i ≥ Q min . In this way, the limit Λ → ∞ reproduces a common lower cutoff on all Q i . The results are shown in figure 13 for the total (Π 1-12 ) as well as for longitudinal (Π 1-2 ) and transversal (Π 3-12 ) contributions separately.
The Regge models in the preceding section predict a ratio ∆a η,η µ /∆a π 0 µ near the expectation (C 2 0 + C 2 8 )/C 2 3 = 3. Similarly, we obtain ∆a η µ /∆a η µ close to 2, as suggested by the scaling with C 2 η /C 2 η (4.32). To first approximation, the implementation of the various asymptotic constraints on the HLbL tensor thus reproduces the simple scaling that Q min (GeV) a μ ⨯ 10 11 Figure 13: Contribution of the pQCD quark loop (with vanishing quark mass) to a µ from the region Q 1,2 ≥ Q min with the contribution from Q 3 below Q min damped by Q 2 3 /(Q 2 3 +Λ 2 ) (plus crossed). The total contribution fromΠ 1-12 is shown in black, together with the partial ones fromΠ 1-2 (red) andΠ 3-12 (blue). The pQCD contribution with common lower cutoff in all Q i is reproduced in the limit Λ → ∞.
originates from the weight factors (2.20) appearing in the V V A triangle. For the mixed regions this behavior is exact due to (4.7) and (4.29), as long as the low-energy properties of the HLbL tensor are not disturbed, while for the asymptotic region it is a consequence of the flavor decomposition chosen in (4.8) and (4.30). The fact that the results from the summation of excited pseudoscalars confirm these expectations indicates that the pQCD quark loop dictates, if not the overall size of the effect, at least its decomposition in the various isospin channels. This is an encouraging sign that the model dependence which is intrinsic in the approach we are following here is mitigated by the QCD constraints. To understand even better the extent to which this mitigation occurs we analyze here in detail the matching between the Regge models and the quark loop integral, after introducing appropriate cutoffs.
For ∆a PS-poles µ ∼ 13 × 10 −11 figure 13 suggests scales Λ and Q min around 1.4 GeV. In addition, the pQCD quark loop would predict an additional increase from the transversal amplitudes around 4 × 10 −11 , but for these scales the interplay with axial-vector resonance contributions needs to be studied in more detail. In the following, we will instead focus on the comparison of our Regge model and the pQCD quark loop in the longitudinal amplitudes.

Matching of short-distance contributions
Beyond the matching at the level of (g − 2) µ , it could be instructive to also compare the specific contributions to the BTT functions in the various kinematic domains. However, once the respective scaling with the virtualities is factored out, we find that the coefficient converges relatively slowly to its asymptotic value. We conclude that it is rather the convolution with the kernel functions T i that becomes important to assess the relevant scales of the SDCs for the HLbL contribution. This is illustrated in figure 14, which shows various contributions to a µ as a function of a lower cutoff on all three virtualities Q i , as well as in figure 15, which shows the opposite case of an upper cutoff on all three virtualities Q i . The ground-state pseudoscalars are saturated by 90% for Q max = 1.5 GeV, while for the excited pseudoscalars only about 25% of the total contribution comes from this energy region. By construction, their contribution asymptotically matches onto the one from the pQCD quark loop, and figure 14 shows how fast that asymptotic limit is reached after convolution with the (g − 2) µ integral kernels: at 1.5 GeV it is saturated by 70%, or about 80% if the tail of the ground-state pseudoscalars is included.
The mixed region is more difficult to illustrate, especially for the corresponding OPE constraint, because in addition to the hierarchy Q 2 3 Q 2 1 ∼ Q 2 2 the small virtuality still needs to be large compared to Λ QCD , otherwise chiral corrections will become important. For that reason, the low-energy part of the integration region was suppressed by the second cutoff Λ in figure 13. To obtain some measure of the size of the mixed-region contribution, figure 16 shows the remainder if for a given cutoff Q cut both the regions where all Q i ≤ Q cut and all Q i ≥ Q cut are subtracted. For the ground-state pseudoscalars at Q cut = 1.5 GeV, this produces the remaining 10% beyond the low-energy region, while the asymptotic region Q i ≥ Q cut is already largely negligible. For the sum of excited pseudoscalars, it is instructive to further scrutinize the decomposition at this scale into low-energy (25%), mixed (40%), and asymptotic (35%) regions. As concerns the contribution from the lowest excitations, in Model 1, the low-energy region is entirely saturated by the sum of the first three excitations, the mixed region by 80%, but for the asymptotic part of the integral the higher excitations make up about 50%. This pattern suggests the interpretation that indeed the lowest excitations are most important for the low-energy and mixed regions, while the infinite tower of resonances restores the correct asymptotic behavior. In fact, we find that the numerical impact of the integration regions where the OPE constraint strictly applies, i.e., where both Q 3 Q 1,2 and Q i Λ QCD , is already very small, so that in practice its main effect lies in constraining the TFF Regge models.
Altogether, this discussion indicates that at some point around 1.5 GeV the description of the HLbL tensor in terms of hadronic intermediate states should be matched onto the one from pQCD. In particular, the implementation of the SDCs in terms of excited pseudoscalars gives an indication how big an impact the intermediate regime may have (in the longitudinal amplitudes): while from pQCD alone one may have guessed a contribution around 5×10 −11 from the asymptotic region, for a value of Q min chosen at 1.5 GeV, the excited pseudoscalars with masses in the same region will add a contribution of similar size, covering also the mixed regions of the (g − 2) µ integral.
To quantify the matching between the quark loop and the description in terms of hadronic states, one would need to define a concrete criterion for the matching scale. One way to define an optimal scale could be to consider the difference between Regge model and quark loop as a function of Q min in combination with the uncertainties of each description for a particular cutoff. For the Regge model, we can estimate this uncertainty as before, but for the quark loop one would need to know the α s corrections and/or higher orders in the OPE, which when compared to the leading-order quark loop would already entail information about the scale where pQCD becomes an efficient description of the HLbL tensor. Absent such calculations, we may obtain a first estimate by comparison to similar  Figure 16: Mixed-region contribution to a µ , defined as the full integral minus the contributions from the low-energy (all Q i ≤ Q cut ) and high-energy (all Q i ≥ Q cut ) regions.
pQCD uncertainties in inclusive τ decays [131][132][133][134][135][136], given that we expect a matching scale not too far off the τ mass, which would suggest an uncertainty around 20%. Based on the combined uncertainties of the Regge model and the pQCD quark loop, we then find a preference for a matching scale around Q match = 1.7 GeV, leading to the decomposition  Note that, in the first line, we also subtracted the very small contribution from the groundstate pseudoscalars from the integration region Q i ≥ Q match to avoid any double counting. As expected, the comparison to (4.59) confirms that for the asymptotic part of the integral it does not matter whether a description based on hadronic intermediate states or pQCD is employed: this means that about one third (i.e., the second line) of (5.1) is a modelindependent part of the effect we have calculated. But how model dependent is the rest and can we adequately cover this model dependence with our uncertainty estimate? There are different uncertainties which need to be considered and we summarize all of them here: • 3.6 units coming from the uncertainties in the parameters of Model 1, as given in (5.1), obtained by stretching the uncertainties in the Regge slopes by a factor three; 11 • 1.7 units are obtained by varying the matching point by 0.5 GeV (the main effect comes from the lowest Q match , which we vary to as low as 1.2 GeV); • 3.8 estimated from the difference between Model 1 and 2, cf. (4.59). 11 We thereby aim to cover scenarios in which other hadronic states could be used to implement the SDCs, in which case the Regge slopes would differ; e.g., according to ref. [108], the Regge slopes of the axial-vector a1 and f1 trajectories are σ 2 a 1 = 1.36(49) GeV 2 and σ 2 f 1 = 1.27(64) GeV 2 .
All these uncertainties concern essentially the contribution below the matching point of Q match = 1.7 GeV, as estimated in the previous section. The outcome of our analysis for this part therefore reads: with a 65% uncertainty, which we consider as sufficiently conservative and, in addition, covers the systematic effects related to the asymptotic behavior of the excited-state TFFs as discussed in appendix E. Another observation corroborating this conclusion is that the contribution to the central value due to the first three pseudoscalar excitations (whose masses and, in part, two-photon couplings are constrained by phenomenology) amounts to and stress that the contribution of the higher excitations (n > 3), which has been calculated with our Regge model and is the most uncertain and model-dependent part of our calculation, amounts to only less than 10% of the total. We conclude that our final result (5.3) has a generously estimated uncertainty that we expect to cover the remaining model dependence.

Chiral limit and role of axial-vector mesons
One may ask whether the implementation of the longitudinal SDCs adopted here would work in the chiral limit: in this limit, excited pseudoscalars have a vanishing coupling to the axial current and therefore would not be able to contribute to the fulfillment of the OPE constraint. 12 However, there are known cases in which the chiral and the large-N c limits do not commute, most notably in the context of baryon chiral perturbation theory. For instance, if one first takes N c → ∞ and then m q → 0, the entire tower of excited baryons contributes to the first non-analytic term in the quark-mass expansion of the nucleon mass, while in the opposite order only nucleon intermediate states appear [137,138], and similar subtleties arise elsewhere due to mass splittings of order 1/N c [139]. Further subtleties in the order of the chiral limit have been pointed out before even for the V V A anomaly itself: the discontinuity of the fermion triangle loop function in the axial-vector virtuality vanishes with the fermion mass, but in a dispersion relation the mass dependence is canceled and produces the anomaly that survives in the chiral limit [140,141]. While these examples show that care is required when exchanging the limits, at least the implementation of the large-N c Regge models described here does not allow any such subtleties to occur and is meant to be used only away from the chiral limit. If the excited pseudoscalar poles were to decouple in the chiral limit, an alternative solution could be provided by axial-vector intermediate states, which do contribute in the chiral and large-N c limits. For these mesons, however, only model-dependent calculations are available in the literature so far, 13 whereas a calculation based on a dispersive framework is still lacking. Such a framework would allow one to express the contribution to HLbL scattering in the most general way in terms of all TFFs of the axial vectors (which admit three), as is the case for the pseudoscalars (which admit only one). Another significant difference is that while for pseudoscalars the sum rules that guarantee the absence of ambiguities in the evaluation of the HLbL contribution [35] are automatically satisfied, this is not the case for axial-vector mesons. We believe that at present it is fair to say that even the ground-state contributions of the latter are poorly understood.
Besides these theoretical reasons, there are also phenomenological ones that favor a discussion in terms of pseudoscalars: while for the most relevant excited pseudoscalar resonances, those in the energy range between 1-2 GeV, there is at least some information on the phenomenology relevant for HLbL, the situation is even worse for the known axial-vector resonances in the same mass range. This is related to the fact that for axial vectors a decay into two real photons is forbidden by the Landau-Yang theorem [145,146]: hadronic channels such as three pions are dominant with respect to suppressed decays to virtual photons, which have been observed only for two ground-state axial-vector resonances [147,148].
If a viable implementation of the longitudinal SDCs in terms of axial-vector resonances were possible, it would have to look quite different from ours in terms of pseudoscalars excitations. Besides the fact that different TFFs contribute, we observe that the axial-vector contribution toΠ 1-3 does not resemble the pseudoscalar-pole contribution (2.16), in fact, both in a Lagrangian-based approach and in dispersion theory the pole in q 2 3 − m 2 A cancels in the longitudinal BTT amplitudes. Based on what is known about the axial-vector TFFs, we cannot preclude the possibility that a finite number of axial vectors could be used to construct such a solution. If that were possible while being consistent with phenomenology and the SDCs on the axial-vector TFFs, this would be an appealing solution, but the necessary theoretical framework for carrying out such an analysis is not yet available. For the moment we took the pragmatic point of view that we implement the longitudinal SDCs in terms of the hadronic intermediate states that we can control best, both theoretically as well as phenomenologically. Having adopted this strategy, we need to address the question of whether our estimate of the systematic uncertainty is large enough to cover the possibility of implementing the SDCs in terms of other hadronic intermediate states.
We believe that we can answer positively this question under the reasonable assumption that even in an alternative scenario the matching to pQCD will occur in the range we have considered. In this case the contribution from the quark loop will remain unchanged and all we need to discuss is the excited-pseudoscalar-pole contribution estimated as 8.7(5.5) × 10 −11 . About one unit out of nine comes from excited states with n > 3: if these were not needed to satisfy the SDCs, this contribution would have to be dropped, a possibility amply covered by our uncertainties. The bulk of the contribution comes from excited states with n ≤ 3, and as we have discussed in section 4.3 and in appendix D, the estimate of their two-photon couplings we have obtained by requiring that the longitudinal SDCs be satisfied is compatible with what is know from phenomenology. Our uncertainty estimate covers the present phenomenological uncertainties on these couplings and could be reduced if the phenomenological information on them were improved. In the end, even if the SDCs were to be implemented using axial-vector states, the first few pseudoscalar excitations would need to be included regardless, it is just that the pseudoscalar TFFs would not be constrained by the HLbL SDCs.
In conclusion, we believe that the uncertainty related to the nature of the hadronic states used in the implementation of the SDCs should be covered by the error assigned in (5.2), an expectation that has been supported more recently by models in holographic QCD [149,150], in which the SDCs are implemented by summation of an infinite tower of axial-vector resonances instead. Beyond the model context, assessing the role of axial-vector resonances in fulfilling the SDCs, especially the transversal ones, will first of all require an improved understanding of their ground-state contributions.

Comparison to the Melnikov-Vainshtein model
In this section, we compare our implementation of the longitudinal SDCs to the one from ref. [62], which is based on the observation that the modification of (2.16) ensures that both the normalization and the mixed-region OPE constraint (3.30) are fulfilled. 14 Since the form (2.16) of the pseudoscalar poles is a direct consequence of the dispersion relation for the HLbL tensor, which we suggested in refs. [31][32][33][34][35] for the case of general four-point kinematics, this modification is not compatible with the description of other intermediate states in the same framework. However, the replacement (6.1) could be justified by a dispersion relation for the HLbL amplitudes directly in the kinematic limit relevant for (g − 2) µ , i.e., for q 4 = 0. In this limit it is not possible to work with a dispersion relation in the Mandelstam variables s, t, and u at fixed q 2 i , because they cease to be independent: s = q 2 3 , t = q 2 2 , and u = q 2 1 . This means that when writing dispersion relations in the q 2 i for g − 2 kinematics, there is no clear separation between the singularities of the HLbL amplitude and those generated by hadronic intermediate states directly coupling to individual electromagnetic currents, e.g., two-pion states as in figure 17: in this framework the pseudoscalar TFFs can no longer be treated as external input quantities. The same holds for higher intermediate states, so that in general the factorization of form factors and scattering amplitudes of the intermediate state in question would cease to apply.
Nevertheless, we observe that, in principle, both forms of dispersion relations are perfectly legitimate-the transition from the dispersion relation for the four-point function [31][32][33][34][35] to a dispersion relation in the photon virtualities in the (g −2) µ kinematic limit amounts to a relabeling of contributions from different principal cuts. This is illustrated by writing primary secondary q 4 = 0 q 2 q 3 q 1 Figure 17: Dispersion relation for the HLbL tensor in the (g − 2) µ kinematic limit with singularities from primary and secondary channels (2π state and pseudoscalar pole). the pseudoscalar pole (2.16) in our framework aŝ where the first term reproduces the pole in the alternative dispersive framework for (g − 2) µ kinematics, while the second term does not contain a pole at q 2 3 = M 2 P . More precisely, the second term is the contribution due to intermediate states X in a cut through the TFF, with the discontinuity determined by the sub-processes γ * (q 3 ) → X and a pseudoscalar-pole contribution to γ * (q 1 )γ * (q 2 ) → γX, as illustrated in figure 17. This piece is present even in the alternative dispersive framework, which demonstrates that changing the dispersive framework simply amounts to a reshuffling of contributions between different principal cuts. Due to the mixed-region SDC, for large q 2 1 ∼ q 2 2 the second piece in (6.2) has to cancel against the contribution from the infinite tower of higher intermediate states up to chiral corrections. Since this is a key point in ref. [62], and the basis for the construction of the MV model, it is worthwhile discussing in detail how this cancellation has to work. For simplicity we concentrate on the pion contribution only (isospin-triplet component) and include all other contributions other than the pion pole toΠ 3 1 into a single function G: From the requirement that (3.29) provide the leading term for asymptotic q 2 1 ∼ q 2 2 , but that it be exact (in the chiral limit) as far as the q 2 3 dependence is concerned, it follows that the function G must have the following asymptotic behavior: We stress that (6.4) is exact in the chiral limit, a property inherited from (3.29), which is a remarkable and interesting result. The MV model (6.1) consists of taking the pion pole as the only contribution toΠ 3 1 . This effectively amounts to promoting (6.4) to an equation valid for any value of q 2 1 and q 2 2 : While this provides a simple way to exactly satisfy (6.4), there is no physical justification in support of such a very strong assumption. It is therefore not surprising that this leads to uncontrolled numerical effects, which we have been able to quantify here. In addition, we note that away from the chiral limit the residue in (6.2) contains F P γγ * (M 2 P ) instead of F P γγ , which at least for η ( ) entails significant chiral corrections.
Numerically, ref. [62] concluded an increase of 13.5 × 10 −11 for the pion and 5 × 10 −11 each for η ( ) , based on the modification in (6.1) and the TFFs from ref. [51] (LMD+V for the pion and VMD for η ( ) ): However, we note that with modern input for the TFFs this number would increase substantially: for the pion, our Model 1 implies an increase of 16.2 × 10 −11 , which would increase to 17.3 × 10 −11 if one used the dispersive TFF instead. Here, the change to the original MV number mainly reflects the differences between the LMD+V model from ref. [51] and the dispersive result for the π 0 TFF [46,47]. For η ( ) , the differences are more severe because the incorrect asymptotic behavior of the VMD form factors in the pQCD limit suppresses the impact of taking the singly-virtual form factor to a constant. We find 10.0 × 10 −11 and 12.1 × 10 −11 for η and η , respectively, which in total produces an increase of 38 × 10 −11 beyond the pseudoscalar ground-state contributions, nearly three times the result given in (4.59). Apart from the overall size, another key difference in our implementation concerns the hierarchy ∆a π-poles µ < ∆a η-poles µ < ∆a η -poles µ found with the excited pseudoscalars, see (4.57), while in ref. [62] the largest effect was found for the pion. The fact that ∆a η-poles µ comes out much smaller than ∆a η -poles µ can be partly explained by the two-photon couplings, F η(n)γγ < F η (n)γγ , and also through the scaling of the excited state TFFs in the BL limit, see figure 19.
This observation also has consequences for the matching to the quark loop. While in our case the scaling of the flavor components follows essentially the expectation from the weights C 2 a , this is not the case for the MV model, and therefore it is less clear how the matching to pQCD should proceed. In fact, as shown in figure 14, despite not being implemented explicitly, the MV model also comes close to the pQCD asymptotics: the main difference to our Regge model occurs in the low-energy region below 1 GeV. This matching onto pQCD asymptotics is coincidental, however: by construction, the model saturates the MV constraint also in the limit in which all virtualities are large and therefore exceeds the proper pQCD limit by a factor of 3/2. Since the asymptotic value is approached rather slowly, the resulting curve happens to be close to the pQCD quark loop for the range of Q min displayed in figure 14. Figures 14-16 also illustrate the origin of the difference between our implementations. For the reference scale of 1.5 GeV, the low-, mixed-, and high-energy regions contribute 75%, 20%, and 5%, respectively, which demonstrates that indeed the approximations of the MV model manifest themselves primarily in the low-energy region, where the dispersive framework provides the best constraints and the contribution of higher states only leads to a moderate uncertainty.

Summary and outlook
In this work, we studied short-distance constraints (SDCs) for the hadronic light-by-light (HLbL) contribution to (g − 2) µ . We concentrated on the longitudinal constraints that are intimately related to pseudoscalar-pole contributions. Since the HLbL tensor can only be constrained from data in the low-energy region, but not in the mixed-and high-energy regions, SDCs are important for a model-independent approach towards HLbL scattering. In sections 2 and 3, the Lorentz decomposition of the HLbL tensor from refs. [33,35] was used to formulate the known expressions for the perturbative QCD (pQCD) quark loop and the operator-product-expansion (OPE) constraints on the HLbL tensor, respectively. The OPE constraint in the symmetric region with Q 2 1 = Q 2 2 = Q 2 3 ≡ Q 2 is given in (3.1), the Melnikov-Vainshtein constraint [62] for the mixed region with Q 2 3 Q 2 1 ∼ Q 2 2 in (3.25). Both are implemented including the singlet component, for which in addition to chiral corrections also perturbative corrections arise.
Subsequently, we focused on the longitudinal SDCs, related to the pseudoscalar-pole diagrams (2.16) by means of (3.30). While a finite number of poles cannot saturate the SDCs, an infinite tower of them can [62]. To that end, we have constructed two models for the transition form factors (TFFs) of ground-state and radially-excited pseudoscalar mesons: our large-N c Regge model for pion, η, and η is described in sections 4.1 and 4.2, and an alternative model using the Regge resummation from ref. [106] is introduced in appendix B to estimate the systematic uncertainty (see also appendix E). While applicable only in the space-like region as relevant for (g − 2) µ , both models satisfy all relevant low-and highenergy constraints for the TFFs-the chiral anomaly (normalization), the Brodsky-Lepage limit, and the symmetric pQCD limit, see (4.12)-(4.14) and (4.26)-(4.28)-give a good description of the experimental data, and reproduce the established results for the groundstate contributions to (g − 2) µ . In addition, with an infinite tower of excited pseudoscalars, they restore the correct asymptotic Q 2 -behavior of the HLbL tensor in the mixed-and high-energy regions, see (4.7) and (4.8), as well as (4.29) and (4.30).
Thus, it has been shown that the SDCs on the HLbL tensor, and in particular the MV constraint, can indeed be satisfied with an infinite sum over excited pseudoscalar-pole diagrams, while maintaining the correct low-energy behavior. Our result (4.59) ∆a PS-poles µ = ∆a π-poles µ + ∆a η-poles µ + ∆a η -poles µ = 12.6(4.1) × 10 −11 , (7.1) derived from the large-N c Regge models alone, is significantly smaller than the original estimate ∆a µ | MV = 23.5×10 −11 from ref. [62], which was obtained by removing the momentum dependence of the TFF at the external photon vertex. In fact, with modern input for the pseudoscalar TFFs this effect would increase further to ∆a µ | MV ∼ 38 × 10 −11 , demonstrating the dangers of ad-hoc modifications of the low-energy properties of the HLbL tensor. Indeed, we observe that by far the main part of the difference to our implementation originates from the low-energy part of the g − 2 integral. Furthermore, in contrast to ref. [62], we find ∆a π-poles µ < ∆a η-poles µ < ∆a η -poles µ . Accordingly, the flavor decomposition into excited π 0 , η, η states follows roughly the expectation from the coefficients determining the SDCs, motivating a matching of our hadronic imple-mentation onto a description in terms of the pQCD quark loop. This matching, illustrated in figures 13-16, shows that, as expected, the ground-state pseudoscalars are relevant only at low energies, but about half the excited-state contribution comes from the integration region of Q i ≥ 1 GeV, while the other half could be interpreted as an estimate of the mixed regions. Since, by construction, the excited-state contribution asymptotically matches onto the one from pQCD, we then replaced the hadronic formulation in favor of the quark loop in the asymptotic part of the integral, at a matching scale of Q match = 1.7 GeV obtained from our best estimates of the uncertainties in the Regge models and pQCD corrections. Due to the assumed pQCD uncertainties and variation of the matching scale, as well as the inflated errors for the Regge slopes in Model 1, see discussion between (5.1) and (5. slightly increases with respect to (7.1), the advantage being that the asymptotic part of the result is manifestly independent of the nature of the hadronic states in terms of which the correct asymptotic behavior was restored. In this way, our final result mainly relies on the Regge models for an estimate of potential contributions for which the asymptotic constraints do not yet apply, and, while data is scarce, this is the energy region where at least some phenomenological guidance for the excited pseudoscalar states is available. In particular, this strategy ensures that since the excited pseudoscalars decouple in the chiral limit, see section 5.3, our implementation of the asymptotic part of the integral remains valid for vanishing quark masses, while for the low-energy phenomenology chiral corrections are essential.
In the future, the matching to pQCD could be improved if explicit calculations of pQCD corrections became available, a first step in this direction was already taken in ref. [61]. Moreover, the phenomenological analysis would profit from further experimental information on the two-photon physics of hadronic resonances in the 1-2 GeV region, which holds true not only for the longitudinal amplitudes but in general. In fact, to address the transversal amplitudes, the effects of axial-vector resonances need to be understood in the context of dispersion relations, especially given that their masses are much closer to the typical matching scale found for the longitudinal SDCs in this paper. Work along these lines is in progress.

A Anomalous Pseudoscalar-Vector-Vector Coupling
For the pion TFF, we only considered the coupling of the pion to a ρω pair, see figure 4, and neglected the contribution given by a ρφ pair. In the following, we motivate why the ρφ pair can be neglected for the pion, and derive the relative strength of 2ρ, 2ω, 2φ, and φω contributions to the TFFs of the η and η , see figure 8.
In ref. [152], we find the Lagrangians for the anomalous pseudoscalar-vector-vector coupling, with g VVΦ = 3g 2 /32π 2 F 2 π , and the electromagnetic photon-vector interaction, with g V γ the individual coupling strengths. Φ stands for the neutral ground-state pseudoscalar mesons, denoted by π 0 , η (8) and η (0) : and V µ stands for the neutral ground-state vector mesons, denoted by ρ µ , φ (8) µ , and φ (0) µ : Note that the latter Lagrangian (A.2) for the neutral vector mesons is given in the ideal mixing situation: ρ ∼ 1/ √ 2 uū − dd , ω ∼ 1/ √ 2 uū + dd and φ ∼ −ss. In general, we use a φ-ω mixing: 90] which, however, almost corresponds to the ideal case (θ ideal . For the η-η mixing, we use the short-hand notation: where the mixing matrix in the standard two-angle mixing scheme is given by: with the mixing parameters introduced in (2.22). The coupling strengths of the pseudoscalar meson to two-photon interactions in the VMD picture, see figures 4 and 8, can be reconstructed from the above Lagrangians, taking into account the φ-ω and η-η mixings: A similar approach is chosen in ref. [153], where the contributions to the singly-virtual TFFs are analyzed through the combination of pseudoscalar-photon-vector and photon-vector interactions. The dependence of the electromagnetic photon-vector interactions (A.2) on the vector-meson masses are canceled out by the vector-meson propagators. Our final couplings read, for P = η, η : with I η = 1, I η = 2, and T IJ given in (A.7). Similarly, we define for the pion:  Table 5: Parameters of the alternative model for pion, η, and η TFFs.
Note that in (A.9) and (A.10) we divided all couplings by a common factor: . This is allowed because the relative strength of the different couplings does not change. The large-N c Regge model for the η ( ) TFFs is then constructed such that each vector-meson pair contributes exactly C P V 1 V 2 /N to F η ( ) γγ , where N is the normalization (4.40). Numerical values for C P V 1 V 2 can be found in table 2. One can clearly see that C π ρω C π ρφ , which is why we neglected the ρφ contribution to the pion TFF. Furthermore, one can see that the ground-state η ( ) TFFs are dominated by the 2ρ, while the contribution from φ-ω mixing is small. This is also illustrated in figures 30 and 33, where we show the 2ρ, 2ω, 2φ, and φω contributions to the singly-virtual and doubly-virtual η ( ) TFFs.

B Alternative model for pion, η, and η transition form factors
In this appendix, we present an alternative model for the pseudoscalar TFFs, which will help us to study the systematic uncertainty of our g − 2 result. This alternative model uses the Regge resummation from ref. [106] to satisfy the pQCD limit of the TFF, cf. (2.25). It reads: where P = π, η, η and the introduced mass spectra again follow a radial Regge ansatz: 2) The first term in (B.1), proportional to c 1 , corresponds to a variant of a large-N c Regge TFF model with equal mass spectra for all vector mesons. The second term in (B.1), proportional to c 2 , has an additional factor of x(1 − x) in the numerator, originating from the asymptotic wave function, cf. (2.24). This term is crucial for the model to satisfy the BL limit (2.25): lim and the fit F η , F η from (4.47). The exponential functions in the numerator, e −(Q 2 shall support the VMD in the region of small momentum transfers, and suppress the x(1 − x) correction which is only needed in the asymptotic region. Therefore, the real-photon limit (2.17) is proportional to c 1 : with ψ (1) the trigamma function and whereas the symmetric pQCD limit (2.25), to leading order in Q 2 , is proportional to c 2 : (B.7) In this way, c 1 and c 2 are fixed and the TFF model reproduces the chiral anomaly, the BL limit, and the symmetric pQCD limit exactly. Evaluating the SDCs of the HLbL tensor, cf. (3.1) and (3.29), we obtain for the mixed region: with Li 2 the dilogarithm, and for the asymptotic region: Thus, both the mixed region and the asymptotic region acquire the correct Q 2 behavior, as is discussed in detail in appendix C. The model parameters M V i , σ V i , and Λ are determined as follows, see table 5: • For the pion TFF we use M V 1 = 1 2 M ρ(770) + M ω(782) = 0.779 GeV; for the η and η TFFs we use for M V 1 the pole-mass parameters of a VMD ansatz fit to the CLEO data [116].
• For reasons of comparison, σ V 1 is chosen to reproduce the values of the two-photon couplings to the first excited pseudoscalars obtain with our large-N c Regge model: Alternatively, one could use the phenomenological constraints on the two-photon couplings listed in table 4; • σ V 2 is chosen to satisfy the MV SDC; • For the pion TFF Λ and M V 2 are adjusted to bring the model in line with the dispersive description of the π 0 TFF [44,46,47]; 15 for the η and η TFFs the same Λ as in the pion case is used, while M V 2 is fit to the available experimental data.
With the parameters in table 5, the MV SDC is satisfied to about ∼ 2 × 10 −3 relative accuracy or better, and the two-photon couplings of our large-N c Regge model are reproduced to about ∼ 3 × 10 −4 relative accuracy or better. The SDC for the asymptotic region, cf. (4.8) and (4.30), is not implemented in our alternative TFF model, however, even without further adjustments it is reproduced to 117 % for the pion, 124 % for the η, and 120 % for the η . Of course, one could also choose the model parameters differently and implement the SDC for the asymptotic region precisely and the MV limit approximately. Note that the parameters M V 1 and σ V 1 are close to the physical values for the masses of the lightest vector mesons and the slopes of their radial Regge trajectories, cf. figure 5. These physical values assure that the first term in (B.1) indeed resembles a large-N c Regge model.
In figure 18, the TFF presented in this appendix (Model 2) is compared to the large-N c Regge model (Model 1) from section 4 for the lowest radial excitations of pion, η, and   Figure 18: TFFs of the first n = 1, . . . 5 radially excited pion, η, and η states. Comparison of the large-N c Regge models from section 4, indicated by the solid curves, and our alternative TFF model (B.1), indicated by the dotted curves. The left panel shows the TFFs in the singly-virtual limit, the right panel shows the doubly-virtual region with η . A comparison to experimental data for the ground-state pseudoscalars is postponed to appendices F.2, G.2, and G.3. The two-photon couplings F P (n)γγ of the excited states come out in close agreement between both models, see also figure 19. For Model 2, we observe an enhancement of the excited-state TFFs in the low-Q region, especially for the doubly-virtual kinematics. This enhancement becomes weaker with increasing Λ, since it is an artefact of the interplay between the two terms in (B.1). Fitting both M V 2 and Λ to data for the ground-state η and η TFFs would lead to Λ < 1 GeV, and thus, exacerbate the enhancement of the excited-state TFFs at low Q. Therefore, we decided to use Λ = 1.318 GeV, as obtained  for the pion, also for η and η . Note that for Model 1 the derivatives of the TFFs in the limit of zero momentum transfer are not unique but depend on the direction, a consequence of the construction in terms of Q 2 − /Q 2 + in (4.10) as a minimal way to implement the different asymptotic limits. This can be seen when comparing the slopes of the singly-virtual and symmetric doubly-virtual TFFs in the left and right panels of figure 18. The modifications of Model 1 described in (E.3) and (E.4) reduce the direction-dependence of the derivative at the origin. However, the derivative of the TFFs is not needed for the evaluation of (g − 2) µ and the alternative implementation in Model 2 does not exhibit this issue: in figure 18, the slopes for Model 2 are always positive, but for Model 1 they change sign between the left and right panels. Accordingly, this will be another systematic effect estimated by the comparison of the two models. We stress again that neither model has the required good analytic properties to remain valid outside the space-like region relevant for (g − 2) µ , of which the zero-momentum-transfer limit of the derivatives in Model 1 is one particular manifestation.
In the right panel of figure 19, the BL limits of the excited-state TFFs are shown. For Model 1 this limit increases with the excitation number n until it reaches an asymptotic value, but for Model 2 it remains constant. Since the true asymptotic behavior for radiallyexcited pseudoscalar TFFs in the BL limit is unknown, the two models with different asymptotics will allow us to understand the systematic uncertainty of our prediction for the excited-state contributions to (g − 2) µ . The symmetric pQCD limit of the TFFs, on the other hand, is independent of the excitation number n for both models. The two-photon couplings, which enter dominantly into the (g − 2) µ integral, agree by default for n = 0 and n = 1, and also match perfectly for n > 2.

C Verifying short-distance constraints for the HLbL tensor
In this appendix, the mathematical formalism used to derive the behavior of the HLbL tensor in the mixed-energy region, cf. (4.15), (4.43), and (B.8), and the high-energy region, cf. (4.17), (4.44), and (B.9), is presented.

C.1 Polygamma functions and infinite sums over rational functions
The gamma function is defined on R * + as [154]: It can be analytically continued to a meromorphic function in the complex plane, with poles at non-positive integers. In order to deal with the infinite sums over pseudoscalar and vector-meson poles, we use the polygamma functions, which are defined on C as derivatives of the logarithm of the gamma function: They are meromorphic in the complex plane and admit the following series representation [154]: which is converging for any z ∈ C except negative integers. With this, we can express an infinite sum over rational functions. Let {f n } n∈N be a sequence of the form f n = p(n) q(n) where p(n) and q(n) are polynomials in n with deg(p(n)) < deg(q(n)). Let α k be the roots of the denominator q(n). If all the roots are simple, the fraction f n can be written as (partial fraction decomposition): where m = deg(q(n)). In general, if one or more roots α k have multiplicity m k ≥ 2, the formula becomes: wherem ≤ deg(q(n)) is the number of distinct roots. It follows that: provided that the series based on the sequence {f n } n∈N is converging. This can be used to compute the infinite sum over pseudoscalar poles, ∞ n=0Π P (n)-pole 1 , within our large-N c Regge model for the TFFs.

C.2 Euler-Maclaurin summation formula
A key ingredient in the discussion of the SDCs for the HLbL tensor is the Euler-Maclaurin summation formula, which describes the difference between an integral and a related sum, see for instance ref. [155, chapter 8]. Notably, it can be used to derive asymptotic expansions.
Let a < b and m > 0 be integers, and f be a function whose derivatives f (2m) (x) are absolutely integrable over the interval (a, b). Then the Euler-Maclaurin formula reads: where the remainder is given by: where x is the greatest integer smaller or equal to x, B s are the Bernoulli numbers, and B s (x) are the Bernoulli polynomials. 16 Using (C.11), we can find a bound for the remainder: In particular, if f (2m) (x) does not change sign in the considered interval, the remainder is bounded by 2 − 2 1−2m times the first neglected term in (C.7). The Euler-Maclaurin formula can be used to derive the asymptotic expansion of the polygamma functions (C.2) at large z ∈ R. To illustrate, consider the trigamma function ψ (1) (z) = ∞ k=0 1 (z+k) 2 . Inserting f (x) = 1 (z+x) 2 , a = 0, b = ∞, and m = 1 into (C.7) leads to the asymptotic expansion: where the notation of the remainder has been slightly modified compared to (C.7) in order to highlight the additional z-dependence. The derivatives f (2m) (x) = (2m+1)! (x+z) 2m+2 do not 16 The Bernoulli numbers can be generated through: and the Bernoulli polynomials can be constructed according to: For 0 ≤ x ≤ 1, they satisfy [155]: 11) change sign and the first neglected term in (C.13) is given by B 2 2! (0 + 2 z 3 ) = 1 6z 3 . This implies that |R 1 (∞, z)| ≤ (2 − 2 1−2 ) 1 6z 3 = 1 4z 3 . In the next subsection, we will be interested in the remainder generated by truncating the asymptotic expansion in (C.13) after the first term: (C.14) It follows from (C.13) that: where the last inequality holds when z ≥ 1. For a general n ∈ N, a similar procedure leads to [154]: The asymptotic expansion (C.16) in combination with (C.3) is sufficient to derive (4.15), (4.17), (4.43), and (4.44), and thereby fix the parameters of the large-N c Regge models to satisfy the required SDCs on the HLbL tensor.

C.3 Short-distance constraints for the alternative transition form factor model
For the alternative TFF model, introduced in appendix B, the situation is more complicated, because the summation over the pseudoscalar-pole diagrams involves three infinite sumsone additional sum over vector-meson towers per TFF (B.1)-and only two of them can be performed analytically. Therefore, one needs to use the Euler-Maclaurin formula to extract the asymptotic behavior.
We first consider the mixed-energy region Q 2 1 ≈ Q 2 2 Q 2 3 . The terms in our TFF model (B.1) with exponential weights in the numerator are suppressed and do not contribute to the MV limit: . (C.17) The integration over the Feynman parameter x is trivial. Since the f nij (y) in (C.17) do not contain any singularities in the space-like region and the integration domains are bounded, the convergence is uniform and the commutation of integrations and summations is justified. 17 Using (C.3), we can express the sums over i and j in (C.17) with trigamma functions: where we use the notations from (4.9), (4.38), (4.39), and (B.2), assuming for simplicity thatM P = M P . The remaining sum over n cannot be performed analytically. We can, however, rewrite the trigamma function as the first term in its asymptotic expansion and a remainder, see (C.14): Here, we defined: and included the remaining terms of (C. 19) in δH P MV . The sum over n in (C.20) can now be expressed in terms of polygamma functions, but the integral over the Feynman parameter y is difficult to perform analytically. Therefore, we expand in Q 2 and Q 2 3 before integrating over y. The assumption that the two operations commute will be checked a posteriori. We find (Q 2 Q 2 3 ):      3 , and the integration over y in (C.19). In all cases, we can see that when r → 1 (MV kinematics), the curves tend to the coefficients obtained by expanding in the virtualities first and then integrating over the Feynman parameter y. Note that the scales on the y-axis vary between the plots.
An appropriate choice of σ V 2 in a P MV therefore reproduces the MV limit. Let us now verify that expanding before integrating in (C.21) was justified and that δH P MV (Q 2 , Q 2 3 ) is subleading. The first issue can be addressed numerically. We use the coordinates defined in (2.14). The kinematics corresponding to the MV limit can be expressed in those coordinates by setting φ = π, then considering r close to 1, and finally taking the limit Σ → ∞. In figure 20a, we study the function Q 2 Q 2 3 H P MV (Q 2 , Q 2 3 ) in (C. 20) for different values of r using a numerical integration over y. One can see that with r getting closer to 1 the curves tend to a P MV and this justifies the commutation of expansion and integration in this case.
We are left to study the error made by considering only H P MV and not the remainder δH P MV , which can be decomposed into three terms: The notation can be understood as follows. From (C.18) to (C. 19), two trigamma functions, which stem from the doubly-virtual TFF F P γ * γ * (Q 2 , Q 2 ) (first index) and the singly-virtual TFF F P γγ * (Q 2 3 ) (second index), were expanded in a leading piece (0) and a remainder (1). In other words, δH P MV(1,0) combines the remainder of the trigamma function in Q 2 with the leading term in the expansion of the trigamma function in Q 2 3 : (C.25) Here, we used (C.12) to show that the term is bounded from above, as can be done for each term in (C.24). As before, it can be checked numerically that H P MV(0,1) indeed tends to the result obtained by expanding first and integrating second, see figure 20b. We proceed analogously for the two remaining terms: which is checked numerically in figure 20c and: (C.27) see figure 20d. The above considerations add up to: (C. 28) i.e., the error we make by keeping only the leading term in the expansion of the ψ (1) in (C.19) is subdominant. When considering the high-energy region Q 2 1 ≈ Q 2 2 ≈ Q 2 3 = Q 2 , the same technique can be applied, but the situation simplifies slightly, since there is only one large scale. The pQCD constraint on the HLbL tensor reads: Similarly to the previous case, since the sum over n cannot be performed analytically, we rewrite the polygamma function as the first term in its asymptotic expansion and a remainder. This leads to: Analogously to H P MV , the term H P pQCD is treated as follows: the sum over n is performed, the expression is expanded in Q 2 , and the integration over y is carried out. The commutation of the expansion and integration is checked numerically in figure 21. For the other terms, δH P pQCD(i,j) , we use (C.15) and then proceed analogously to the leading term, see figure 21 for the numerical checks.  Figure 21: Numerical checks for the validity of commuting the expansion in Q 2 and the integration over y in (C.30). In all cases, we can see that the curves tend to the coefficients obtained by expanding in the virtualities first and then integrating over the Feynman parameter y. Note that the scales on the y-axis vary between the plots.

D Two-photon couplings of excited pseudoscalars
In this appendix we collect the phenomenological information that is available on the twophoton couplings of the excited pseudoscalars listed in the PDG [90], in comparison to the two-photon couplings of the first radially-excited pion, η, and η states as shown in figure 19 for both the large-N c Regge and the alternative TFF model.
Phenomenologically, the two-photon couplings of the excited pseudoscalars are unknown, but for many states some information on these couplings can be extracted either from direct limits on the two-photon channel or from measurements of particular branching fractions. For the excited pion states, the only available information concerns the π(1300). The blue bar in the left panel of figure 19 indicates values excluded by the limit F π(1300)γγ < 0.0544(71) GeV −1 . (D.1) Since at present there is no measurement of the π(1300) width and the two-photon branching ratio, the above bound is an estimate based on the available empirical information. The π(1300) decays predominantly into 3π, e.g., into ρπ: 18 Γ(γγ)Γ(ρπ) Γ total < 85 eV [159], Γ(π(ππ) S-wave ) Γ(ρπ) = 2.2(4) [157], (D.2) whereas the πf 0 (1300) and γγ decays are suppressed [160]. Assuming Similar constraints exist for several excited η, η states. As discussed in the main text, the assignment of Regge trajectories is not settled, so here we simply reproduce the listing according to the PDG, see section 4.3 for a discussion of the phenomenological implications. We stress that given that even the identification of states is contentious, the experimental limits should be treated with caution and mainly serve as guidance that our Regge models do not assume implausible values for the two-photon couplings. Assuming that the branching fraction into other channels can be neglected, 19 we would conclude Γ(γγ) < 80 eV and thus Using the total width Γ η(1405) = 48(4) MeV as measured in the KKπ channel, the two limits involving this channel imply Γ(γγ) < 1.73 keV, while assuming that KKπ and ηππ 18 Note that the role of the S-wave component is not settled: while Γ(π(ππ)S-wave/Γ(ρπ) = 2.12 from ref. [156] agrees with ref. [157], ref. [158] found a negligible S-wave component Γ(π(ππ)S-wave/Γ(ρπ) < 0.15. In the latter case the limit on the two-photon decay width F π(1300)γγ would become stricter by a factor √ 3.2 ∼ 1.8, in which case there would be some mild tension with F π(1300)γγ implied by our Regge models. 19 The limit from ref. [161] already includes the conversion factor Γ(ηπ + π − )/(Γ(ηπ + π − ) + Γ(ηπ 0 π 0 )) = 2/3, which emerges from the combination of isospin and symmetry factors in Γ(ηπ 0 π 0 )/Γ(ηπ + π − ) = 1/2. which, assuming dominance of η ππ and including the neutral channel by means of the relation Γ(η π 0 π 0 )/Γ(η π + π − ) = 1/2, translates into In case other channels do contribute, this number would have to be considered a lower limit.

E Systematic uncertainties and decay constants of excited pseudoscalars
The systematic errors quoted for ∆a µ in section 5.2 are based on comparing results from two different models for the pseudoscalar TFFs introduced in section 4 and appendix B, respectively. This has then been added to a conservatively estimated uncertainty coming directly from the parameters of our models. By construction, our models link the TFFs of the different pion, η, or η states such that a resummation is at all possible, but it is clear that the details of the TFFs so obtained may turn out not to be realistic, at least for the lowest lying excited pseudoscalars. Here, we explore this specific question in particular for the first excited pseudoscalars, on which some information from the phenomenology is indeed available. As we will show, if we adapt the parameters of our models presented in sections 4.1 and 4.2 to be in agreement with phenomenology, or theoretical expectations, we obtain shifts in our results which are well covered by the present error budget.
In section 4.3 and appendix D, we confirmed that our TFF models are well compatible with phenomenological constraints for the two-photon couplings of η(1295), η(1405), η(1475), η(1760), and X(1835). In the following, we will be interested in the leptonic decay constants, F P (n) , of the excited pseudoscalars, limiting our analysis to n ≤ 3 states. Since all pseudoscalar mesons, except for the Goldstone mode, decouple from the axial-vector current in the chiral limit, the corresponding decay constants, defined in (2.21), are suppressed. Note that contrary to the decay constants, the two-photon couplings are non-vanishing in the chiral limit. At low Q 2 it is the latter that are most relevant.
We start by considering the symmetric pQCD limit of the TFFs, which for the groundstate pseudoscalars is given in (4.6) and (4.28). The same relations also hold for the excited pseudoscalars, replacing only the ground-state decay constants, i.e., F P → F P (n) . To change the asymptotic limit of F P (n)γ * γ * (−Q 2 , −Q 2 ) for n = 1, 2, 3, we modify c diag by replacing (F.1) with and (G.1) with and analogously for η , while keeping all other model parameters the same, cf. table 1.
Varying only the TFFs of the lowest excitations n ≤ 3, it is ensured that the SDCs remain intact. Applying this modification for n = 1 decreases ∆a µ by (0.22+0.27+0.45)×10 −11 = 0.94×10 −11 , where the individual numbers are the pion, η, and η contributions, respectively. Note that we are removing the contribution from Q i ≥ Q match , as this region is described by the pQCD quark loop, cf. section 5.2.
The BL limit of the pseudoscalar TFFs is not known for general n. Presently, this limit is therefore not fine-tuned in our large-N c Regge models, as can be seen in the right panel of figure 19. Here, we want to assume a BL limit constant in n, lim Q 2 →∞ F P (n)γγ * = lim Q 2 →∞ F P γγ * , which is also found for the alternative TFF model introduced in appendix B. Replacing (F.3) with and (G.2) with while keeping all other model parameters the same, this is achieved without changing the two-photon couplings or other SDCs. Looking at the n = 1 case, ∆a µ decreases by (0.49 + 1.18 + 1.92) × 10 −11 = 3.59 × 10 −11 . This shift predominantly comes from a change of the low-Q 2 shape of the TFFs and not from the different asymptotic limits. Combining the changes in (E.1) and (E.3), as well as (E.2) and (E.4), the total decrease of ∆a µ amounts to (0.62 + 1.27 + 2.10) × 10 −11 = 3.99 × 10 −11 if only n = 1 is modified, 5.26 × 10 −11 if also n = 2 is modified, and 5.82 × 10 −11 if n = 1, 2, 3 are modified.
As a last point, we study the effect of non-diagonal couplings.
To be more precise, in our large-N c Regge models from section 4, we allowed the n-th pion, η, or η excitation to couple only to the n-th ρ, ω, and φ excitations, whereas now we allow the first excited pseudoscalars to couple to the ground-state vector mesons. For the π(1300), we modify , and generate the two-photon coupling through ρ(770) and ω(782). For the first excited η ( ) , we only express the dominant isovector-isovector part of the two-photon coupling through ρ(770): .

(E.6)
These modifications lead to a decrease of (0.53 + 0.54 + 0.72) × 10 −11 = 1.79 × 10 −11 . Such non-diagonal couplings are also present in the alternative TFF model introduced in appendix B. All these modifications affect in one way or another the low-Q 2 behavior of the excitedstate TFFs, which could be constrained more rigorously if data were available. At present, we observe that the corresponding changes, which tend to lower ∆a µ , are well covered by our final uncertainty estimate in (5.2). Since, on the other hand, the consideration of Model 2 suggests systematic effects in the opposite direction, we leave the central value as derived from Model 1, with uncertainties as assigned in (5.2).

F Pion transition form factors
In this appendix, we describe the large-N c Regge model for the pion TFF, given in (4.10), in more details. A comparison to experimental data and other parameterizations available from the literature is postponed to section F.2. Based on the constraint equations in (4.14), (4.15), and (4.17), the model parameters should be replaced by: where the parameters related to the asymptotic limits of the HLbL tensor simplify to: and the auxiliary functions Here, the definitions from (4.16) and (4.18) are used.
F.2 Comparison of data and literature: F π(n)γ * γ * In this appendix, we compare our large-N c Regge model, "Model 1" (4.10), and our alternative model, "Model 2" (B.1), for F π(n)γ * γ * to data and other parameterizations available from the literature. In figure 7, the singly-virtual π 0 TFF is shown for Q 2 ∈ [0, 35] GeV 2 . In figure 23, we focus on the low-Q region and include a comparison to the recent lattice QCD [48] and DSE [50] results. Our π 0 TFF models, for which we do not display error estimates, are in good agreement with the dispersive and lattice QCD TFFs, while we observe some deviation of our Model 1 from the DSE prediction. However, the error quoted for the DSE result in ref. [50], as pointed out therein, is only a rough estimate based on the variation of their one model parameter and does not account for the total truncation error. Therefore, we conclude that our π 0 TFF models also agree with the DSE prediction in the singly-virtual region.
In figure 24, we show the doubly-virtual π 0 TFF in the low-Q region. Both lattice QCD and DSE are able to give much more accurate predictions of the (pseudoscalar) TFFs for doubly-virtual than for singly-virtual kinematics, as is obvious by comparing the error bands in figures 23 and 24. In the symmetric region, Q 2 1 = Q 2 2 = Q 2 , starting from ∼ 1 GeV 2 , the DSE predict a slightly larger π 0 TFF than lattice QCD, see left panel in figure 24. Our models for the π 0 TFF run just between these DSE and lattice QCD predictions. Note, however, that the discrepancy in figure 24 is visually enhanced by showing . For doubly-virtual kinematics away from the symmetric limit, see for instance Q 2 1 = Q 2 and Q 2 2 = 2Q 2 in the right panel of figure 24, our π 0 TFF models are in closer agreement with the lattice QCD prediction.  Figure 23: Singly-virtual π 0 TFF in the low-Q region. The large-N c Regge model, "Model 1" (4.10), is indicated by the dashed pink curve. Our alternative TFF model, "Model 2" (B.1), is indicated by the solid cyan curve. The gray band with the dotted curve is the dispersive result from refs. [46,47]. The blue band with the long-dotted curve is the lattice QCD result from ref. [48]. The green band with the dot-dashed curve is the DSE result from ref. [50]. The data are from CELLO [115] and CLEO [116].  Figure 24: Doubly-virtual π 0 TFF in the symmetric region Q 2 1 = Q 2 2 = Q 2 (left) and in the region where Q 2 1 = Q 2 and Q 2 2 = 2Q 2 (right). Legend is the same as in figure 23.
In figure 25, we show the doubly-virtual π 0 TFF in the symmetric region for Q 2 ∈ [0, 40] GeV 2 . Both models, but in particular Model 1, are in perfect agreement with the dispersive description [46,47]. In figure 22, Model 1 and 2 are shown in the full space-like region for Q 2 1 , Q 2 2 < 10 GeV 2 . One can see that their main difference is in the regions where at least one of the photon virtualities is small.   In this appendix, we describe the large-N c Regge model for the η and η TFFs, introduced in section 4.2, in more details. A comparison to experimental data and other parameterizations available from the literature is postponed to appendices G.2 and G.3.
All expressions are given for the η, but hold as well for the η after obvious replacements (including C 8 → C 0 ). Based on the constraint equations in (4.42), (4.43), and (4.44) the model parameters should be replaced by: where the parameters related to the asymptotic limits of the HLbL tensor simplify to: Here, the definitions from (4.16) and (4.18), the auxiliary functions from (F.6), as well as  [49]. The green band with the dot-dashed curve is the DSE result from ref. [50]. The data are from CELLO [115] and CLEO [116].
are used. Note that for the η one has to instead choose: as the physical solution for the quadratic equation (4.44).
G.2 Comparison of data and literature: F η(n)γ * γ * In this appendix, we compare our large-N c Regge model, "Model 1" (4.34), and our alternative model, "Model 2" (B.1), for F η(n)γ * γ * to data and other parameterizations available from the literature. The error band shown for Model 1 is generated by propagating the errors of the input parameters σ P , σ V , F ηγγ , F η , F 8 , F 0 , θ 8 , θ 0 . In figure 10, the singly-virtual TFF of the ground-state η is shown for Q 2 ∈ [0, 40] GeV 2 . In figure 28, we focus on the low-Q region and include a comparison to the DSE result [50]. One can see that our models agree with the experimental data from CELLO [115] and CLEO [116], but tend to a larger η TFF than CA [49] and DSE. In addition, Model 2 is larger than Model 1 for Q 2 < 2.4 GeV 2 . This low-Q enhancement explains why a η-pole µ | Model 2 > a η-pole µ | Model 1 , see (4.49).  Figure 29: Doubly-virtual η TFF in the symmetric region Q 2 1 = Q 2 2 = Q 2 (left) and in the region where Q 2 1 = Q 2 and Q 2 2 = 2Q 2 (right). Legend is the same as in figure 28. In figure 29, the doubly-virtual η TFF is shown for two kinematic situations: symmetric momenta, and Q 2 1 = Q 2 and Q 2 2 = 2Q 2 . Considering Model 1, we observe a slight tension with the DSE prediction in the region of Q 2 ∈ [0.2, 0.8] GeV 2 . This tension should, however, not be taken too serious, because both the DSE and our error band are only based on the variation of input parameters and do not take into account all possible error sources.
In figure 26, Model 1 and 2 are shown in the full space-like region for Q 2 1 , Q 2 2 < 10 GeV 2 . One can see that their main difference lies, similarly as for the π 0 TFF, in the regions where at least one of the photon virtualities is small.  [49]. The green band with the dot-dashed curve is the DSE result from ref. [50]. The data are from L3 [127], CELLO [115], and CLEO [116].
In the left panel of figure 30, the ground-state η TFF is decomposed into the contributions from 2ρ, 2ω, 2φ, and φω vector mesons. As expected, the TFF is dominated by the isovector-isovector 2ρ contribution, followed by the isoscalar-isoscalar 2φ contribution. The φω contribution (4.36), which was needed to generate enough freedom in our large-N c Regge model to satisfy the BL limit of the TFF and the two SDCs on the HLbL tensor, is small.
In the right panel of figure 30, we show the TFF of the first (n = 1) radially-excited η state. In the doubly-virtual region, the relative strength of vector-meson pairs is comparable to what one finds for the ground-state η. The 2ρ contribution is now slightly smaller than the total TFF, and the φω contribution is now larger than the 2ω contribution. In contrast, the singly-virtual TFFs of the radially-excited η states will be dominated by the φω contribution. This enhancement is generated by the n-dependence in the numerator of (4.36) through terms proportional to M +, n . The two-photon couplings and BL limits of the excited-state η TFFs are shown in figure 19.
G.3 Comparison of data and literature: F η (n)γ * γ * In this appendix, we compare our large-N c Regge model, "Model 1" (4.34), and our alternative model, "Model 2" (B.1), for F η(n)γ * γ * to data and other parameterizations available from the literature. The error band shown for Model 1 is generated by propagating the errors of the input parameters σ P , σ V , F η γγ , F η , F 8 , F 0 , θ 8 , θ 0 .  Figure 32: Doubly-virtual η TFF in the symmetric region Q 2 1 = Q 2 2 = Q 2 (left) and in the region where Q 2 1 = Q 2 and Q 2 2 = 2Q 2 (right). Legend is the same as in figure 31.
In figure 11, the singly-virtual TFF of the ground-state η is shown for Q 2 ∈ [0, 40] GeV 2 . In figure 31, we focus on the low-Q region and include a comparison to the DSE result [50]. One can see that Model 1 agrees with the experimental data from L3 [127], CELLO [115], and CLEO [116], as well as the CA [49] and DSE results. Model 2 tends to a larger η TFF for Q 2 < 2 GeV 2 . This low-Q enhancement explains why a η -pole µ | Model 2 > a η -pole µ | Model 1 , see (4.50).
In figure 32, the doubly-virtual η TFF is shown for two kinematic situations: symmetric momenta, and Q 2 1 = Q 2 and Q 2 2 = 2Q 2 . Model 1 is in slight tension with the DSE prediction for Q 2 ∈ [0.2, 1.6] GeV 2 . This tension should, however, not be taken too serious. A comparison of our models with the CA result shows perfect agreement for symmetric momenta. For large photon virtualities, both models agree with each other and give a reasonably good description of the recent doubly-virtual η TFF data from BaBar [122], see figure 9.
In figure 27, Model 1 and 2 are shown in the full space-like region for Q 2 1 , Q 2 2 < 10 GeV 2 . One can see that their main difference is, similar as for the η TFF, in the regions where at least one of the photon virtualities is small.
In the left panel of figure 33, the ground-state η TFF is decomposed into the contributions from 2ρ, 2ω, 2φ, and φω vector mesons. As expected, the largest contribution to the TFF is coming from the isovector-isovector 2ρ mesons, followed by the isoscalar-isoscalar 2φ mesons. Unlike in the case of the η TFF, the 2φ mesons gives a positive contribution to the η TFF, just like the 2ρ, 2ω, and φω mesons. Thus, since the 2ρ contribution does not need to cancel out a negative 2φ contribution as it does in the η TFF, it appears to be smaller than the total η TFF. The φω contribution (4.36), generated through φ-ω mixing, is small.
In the right panel of figure 33, we show the TFF of the first (n = 1) radially-excited η state. In the doubly-virtual region, the relative strength of vector-meson pairs is similar to what one finds for the ground-state η . The φω contribution is now larger than the 2ω contribution, and the 2φ contribution at low Q. In contrast, the singly-virtual TFFs of the radially-excited η states will be dominated by the φω contribution, while all other contributions are of negligible size. This enhancement is generated by the n-dependence in the numerator of (4.36) through terms proportional to M +, n . The two-photon couplings and BL limits of the excited-state η TFFs are shown in figure 19. s = q 2 3 and t = q 2 2 , but by no means does it imply q 2 3 = M 2 π . Of course one is free at that point to separate the pure pole in q 2 3 (with only its residue in the numerator) from non-pole terms. Between the two different dispersive representations, a simple reshuffling takes place, see (6.2) and the whole discussion in section 6.
2. As discussed in section 6, the MV model is based on an unjustified extrapolation to low q 2 1,2 of the constraint at high q 2 1,2 . We have called this a "distortion" of the lowenergy behavior of the HLbL tensor in ref. [68], a description considered unjustified in ref. [151]. Figures 14 and 15 very clearly illustrate this distortion. An alternative solution to the SDCs based on a tower of axial-vector mesons in holographic QCD has been presented in two papers [149,150], which appeared after ref. [151]. These alternative solutions to the SDCs have a very similar behavior as the curves corresponding to our model in figures 14 and 15, and confirm that the MV model [62] leads to a low-q 2 behavior of the HLbL amplitude that cannot be explained in terms of any other physical states-in other words, a "distortion." 3. After (18) in ref. [151] it is stated that: the model in the present manuscript "violates the above equation and claims, effectively, that c ρ L ∼ 1 also in the chiral limit." This statement is incorrect: in our model we are not able to take the chiral limit simply because it is formulated in terms of effective parameters that are fit to data or theoretical constraints. It is not the point of the model to make any claim about the behavior in the chiral limit. The underlying philosophy is to fulfill the SDCs only for large q 2 3 , where the chiral limit becomes irrelevant, and not to rely on it at low q 2 3 , because it would be a bad approximation and we can instead use known phenomenological constraints. Such a strategy is best carried out with excited pseudoscalars because, unlike for axial-vector resonances, there are no ambiguities regarding their dispersive definition and because at least some information from phenomenology is available. Section 5.3 explicitly discusses the issue of the chiral limit in our model. 4. The conclusion of ref. [151] contains the following three statements: "there is no doubt that: (a) this region (Q 2 1,2 Q 2 3 ) provides the largest contribution to a HLbL µ ; (b) it allows for an exact non-perturbative analysis of the longitudinal structure function in the chiral limit and (c) it supplies strong evidence that corrections to the chiral limit are small." The first statement is plainly wrong, in particular for the model by MV, which actually receives most of its corrections from the region Q 2 1,2,3 < Q 2 match , as can be clearly seen from figure 14. It is precisely this observation that leads to the conclusion that the modifications in the low-q 2 region are unphysical. Point (b) is correct, by construction, but the question is whether the chiral limit is a useful approximation at low q 2 , which is the most important region for a µ . This is claim (c), which, unfortunately, is also not correct: the difference between the original MV model (6.1) and the first term in (6.2) is such a quark-mass correction. In the case of the pion the two expressions give contributions to a HLbL µ that indeed differ by a small amount (about 10%), but in the case of η and η the difference is much larger, about 100%, as anticipated in section 6.