Pseudoscalar pole light-by-light contributions to the muon $(g-2)$ in Resonance Chiral Theory

We have studied the $P\to\gamma^\star\gamma^\star$ transition form-factors ($P=\pi^0,\eta,\eta'$) within a chiral invariant framework that allows us to relate the three form-factors and evaluate the corresponding contributions to the muon anomalous magnetic moment $a_\mu$, through pseudoscalar pole contributions. We use a chiral invariant Lagrangian to describe the interactions between the pseudo-Goldstones from the spontaneous chiral symmetry breaking and the massive meson resonances. We will consider just the lightest vector and pseudoscalar resonance multiplets. Photon interactions and flavor breaking effects are accounted for in this covariant framework. This article studies the most general corrections of order $m_P^2$ within this setting. Requiring short-distance constraints fixes most of the parameters entering the form-factors, consistent with previous determinations. The remaining ones are obtained from a fit of these form-factors to experimental measurements in the space-like ($q^2\le0$) region of photon momenta. The combination of data, chiral symmetry relations between form-factors and high-energy constraints allows us to determine with improved precision the on-shell $P$-pole contribution to the Hadronic Light-by-Light scattering of the muon anomalous magnetic moment: we obtain $a_{\mu}^{P,HLbL}=(8.47\pm 0.16)\cdot10^{-10}$ for our best fit. This result was obtained excluding BaBar $\pi^0$ data, which our analysis finds in conflict with the remaining experimental inputs. This study also allows us to determine the parameters describing the $\eta-\eta'$ system in the two-mixing angle scheme and their correlations. Finally, a preliminary rough estimate of the impact of loop corrections ($1/N_C$) and higher vector multiplets (asym) enlarges the uncertainty up to $a_\mu^{P,HLbL} = (8.47\pm 0.16_{\rm sta}\pm0.09_{1/N_C}{}^{+0.5}_{-0.0}{}_{\rm asym})\cdot 10^{-10}$.


Introduction
The electron anomalous magnetic moment, a e = (g e − 2)/2, is the most precisely measured [1] and predicted [2] observable in nature. There is, however, a greater interest in the muon anomalous magnetic moment, a µ = (g µ − 2)/2, since heavy physics beyond the Standard Model (SM) would have an effect of order (m µ /m e ) 2 ∼ 4 × 10 4 times bigger in a µ than in a e [3,4,5,6]. Similarly, a τ should be ∼ 280 times more sensitive to heavy new physics than a µ . However, the significantly lower mean lifetime of the τ lepton makes particularly difficult to measure this property (see e.g. refs. [7,8]), which by now is still consistent with zero [1].
The most accurate measurement from Brookhaven [9] of a µ seems at odds with the SM prediction [10,11,12,13,14,5,15,16] with a 3.5σ discrepancy [1] 1 and, both, FNAL muon g-2 and the J-PARC E34 collaborations have announced new experiments that will reduce the current Brookhaven error by, at least, a factor 4 [20,21]. Therefore, it becomes necessary to make a more accurate prediction for this observable with reduced uncertainty so as to be able to confront the forthcoming results from both collaborations with a SM result of comparable accuracy.
The main source of uncertainty comes from the hadronic contributions, so these are where our activity should focus on (see in ref. [22] an early discussion of this issue). All such contributions can be divided into two, the Hadronic Vacuum Polarization (HVP) and the Hadronic Light-by-Light Scattering (HLbL). The HVP can be completely data-driven through dispersion relations [23,24] 2 , while the latter cannot be obtained completely in this way yet. However, there have been remarkable advances in determining the HLbL part in a model independent way, by means of Lattice QCD [26,27,28] and dispersion relations [29,30]. Therefore, it is the HLbL contribution to the muon anomalous magnetic moment, a HLbL µ , which calls for a dedicated theory effort [31].
We will focus on the leading contribution to the HLbL [5], which is given by the pseudoscalar exchange 3 , a P,HLbL µ . To evaluate such contribution, a necessary ingredient is the P → γ ⋆ γ ⋆ transition form-factor (TFF), F P γ ⋆ γ ⋆ , which cannot be computed analytically in the underlying theory of strong interactions. However, as for the HLbL part, there have also been remarkable advances in determining the TFF in a model independent way by means of Lattice QCD [34] and dispersion relations [35]. We calculate it by means of the extension of Chiral Perturbation Theory (χPT) [36] that incorporates the lightest resonance multiplets in a chiral invariant framework [37], called Resonance Chiral Theory (RχT). We work within the large-N C limit and assume a 1/N C expansion, in such a way that we have a spontaneous chiral symmetry breaking pattern U (3) L × U (3) R /U (3) L+R , which gives place to a nonet of chiral (pseudo) Goldstones. We consider a RχT Lagrangian, L RχT , that includes U (3) symmetric operators and terms that introduce quark mass corrections in the even and oddintrinsic parity sectors [38,39,40,41]. Previous works [41] and [42] analyzed U (3) symmetric TFFs, given in the chiral limit. The novelty of the present approach is that the TFFs are studied beyond the massless pseudo-Goldstone limit, accounting for its leading order corrections in powers of the pseudo-Goldstone bosons squared masses, m 2 P , which explicitly break U (3) flavor symmetry. We will see that all but eight parameters (including among them the four η − η ′ mixing parameters) are fixed by short distance constraints. These eight unknown couplings will be determined through a fit to the experimental data in the space-like region of photon four-momenta q 2 ≤ 0 (note that the pseudoscalar pole contribution to a µ can be written in terms of just the space-like TFFs [43]). Since the data for the π 0 transition form-factor given by BaBar collaboration [44] seems to be at odds with its Brodsky-Lepage high-energy limit [45] and the η and η ′ TFFs, we will discuss its validity and consistency with other data. We will perform various alternative analyses and take as our reference fit the one without this set of data.
Our approach intends to introduce the U (3) breaking through a chiral invariant Lagrangian, where quark mass effects enter in a covariant way. This gives slightly simpler expressions for these form-factors and the short-distance constraints for the Lagrangian parameters, contrary to the strategy followed in ref. [46], which allows a completely general U (3) breaking pattern. Despite giving an accurate description of data, is a little less straightforward to employ and does not rely on chiral symmetry (where quark mass correction must enter through a more restricted pattern). Likewise, ref. [46] also used data from both time-and space-like regions, whereas the present article will only rely on space-like data. We will work in the large-N C limit and hadron loop effects will be neglected, so time-like observables (e.g., partial widths) will be excluded from our fits, as they require a dedicated analysis of the 1/N C corrections. In addition, photon radiative corrections to P → γ (⋆) γ (⋆) decays also play a non negligible role [47].
The paper is organized as follows: In section 2, we recall the relevant pieces of the Resonance Chiral Lagrangian needed for our study. Particularly, we explain the flavor-breaking corrections that we include for the first time in this kind of analysis (see, however, ref. [48] for a related discussion within vectormeson dominance, constituent-quark loops, the QCD-inspired interpolation by Brodsky-Lepage, and Chiral Perturbation Theory). Then, in section 3 we collect our results for the transition form-factors including the contribution from intermediate vector and pseudoscalar resonances. In section 4 we discuss the shortdistance constraints on the Resonance Chiral Lagrangian parameters that are obtained by demanding the QCD short-distance behaviour to the VVP Green's function and the pseudo-Goldstone form-factors (P -TFF). Our fits to data are presented in section 5, where we also quote our results for the parameters describing the η − η ′ mixing in the double angle scheme, including their correlations. We also compare our approach to other recent articles on the subject in section 5.2. Finally, the branching ratio predictions for P → γ (⋆) γ (⋆) processes are compared to their respective experimental measurements in section 5.3. The corresponding P pole contributions to a µ are obtained in section 6, with a careful statistical treatment and highlighting the influence of π 0 transition form-factor BaBar data. Finally, our conclusions are summarized in section 7. The three appendices collect, respectively, the Wess-Zumino-Witten Lagrangian, the relevant formulae used to evaluate the P pole contribution to a µ , and the correlation matrices for the two alternative fits considered in our analysis.

Relevant operators for the TFF
In modeling the TFF we make use of RχT, an extension of χPT that also includes the lightest resonance multiplets [37]. The Wess-Zumino-Witten action [49] describes the local P γγ interaction vertex, whereas a dressed photon description is needed to assure the ultraviolet convergence of the pseudo-Goldstone exchange contribution to a HLbL µ [5]. Within our approach, this is done by considering the vector resonances exchange between the pseudo-Goldstone and the final state photons. To include such interactions one needs to consider the oddintrinsic parity sector with operators including a pseudo-Goldstone and, either two vector resonances or one vector resonance and one external photon.
The complete odd-intrinsic parity basis involving pseudo-Goldstones and either two vector resonances or a vector resonance and a photon was given in [41]. We will rely, however, on the operators given in ref. [40], since these conform a complete basis for describing vertices involving only one pseudo-Goldstone. Also, this basis is simpler for the problem at hand (optimized for the study of this type of processes) and, as proven in [50], equivalent to the complete basis of Ref. [41] for a single pseudo-Goldstone field. In order to compare and employ results from [41], pseudoscalar resonances P ′ are also included in our description of the TFF. Ref. [42] found that the effect of the latter on the form-factors was equivalent to adding a second vector multiplet. Such P ′ interactions are taken from ref. [41] within the chiral limit, which give their first corrections to the TFF at O(m 2 P ) and will be addressed in subsection 2.2.
The Lagrangian describing the interaction between the lightest multiplet of pseudoscalars (the chiral pseudo-Goldstones) and massive meson resonances, R, can be organized according to the number of heavy fields: We will not consider operators with four or more resonance fields since they cannot contribute to the P → γ ⋆ γ ⋆ TFFs at the tree-level -considered in this analysis-. We will describe the spin-1 resonance fields in the antisymmetric tensor representation [37]. In this formalism, as a general feature, the simplest operators in the even-intrinsic parity sector contain an O(p 2 ) chiral tensor (with two derivatives or equivalent scales) [37] in addition to the resonance field, while in the odd-intrinsic parity sector the light-pseudoscalar tensor accompanying the resonance fields usually starts at O(p 4 ) [40,41]. Furthermore, the RχT Lagrangian will be divided in even and odd-intrinsic parity sectors, the latter containing the anomalous Wess-Zumino-Witten Lagrangian, L W ZW . We will present now the RχT operators that provide the relevant interactions between photons, γ, pseudo-Goldstones, φ a , and vector resonances, V : • Operators without resonance fields: The operators in L non−R only contain pseudo-Goldstone fields φ a and start in the even-intrinsic parity sector at O(p 2 ) (given by Gasser and χP T Lagrangian) and at O(p 4 ) for the anomalous oddintrinsic parity sector (where the lowest order contribution is provided by L W ZW ). However, as we want to study the most general possible U (3) breaking we will also consider operators with the structure of the oddintrinsic parity O(p 6 ) χPT Lagrangian [51]. Thus our non-resonant part of the RχT Lagrangian is given by where A stands for the trace in flavour space of A, the operators O W j can be found in table 1 and the well-known -though lengthy-form of L W ZW is given in App. A [49,51]. Although this non-resonant Lagrangian was considered in the chiral and large-N C limit VVP Green's function analysis [41], an appropriate description of physical P → γ ⋆ γ ⋆ processes requires further pseudo-Goldstone bilinear terms not shown above [52,103,59], which dress the φ a wave-functions and induces the η − η ′ mixing. This details are discussed in the later section 2.2. In order to explore the U (3) breaking in the η − η ′ system in more generality we have allowed here the presence of the double-trace operator C W 8 . Although it is subleading in 1/N C and it should be dropped from our analysis, 1/N C corrections play an important role in the description of the η − η ′ mixing so one could argue that this operator might be numerically relevant. We will show, nonetheless, that after demanding that the TFFs follow the highenergy QCD behaviour this subleading coupling C W 8 vanishes. We want to emphasize that L non−R in eq. (2) is not the low-energy χPT Lagrangian; it belongs to RχT Lagrangian which describes the meson interactions in the range of high and intermediate energies and the C W j and C W, χPT j couplings must not be confused.
• Operators with one vector resonance field: At tree-level, the P -TFF will have contributions from V −γ and V −γ −φ a vertices, with the vector multiplet V µν = 8 a=0 1 √ 2 λ a V a µν described in the antisymmetric tensor formalism [37] 4 : A quark mass correction to the V − γ transitions is allowed in L even V to account for the U (3) breaking in this vertex, where the only single-trace operator at O(m 2 P ) is given by the λ V term, given in ref. [38] by O V 6 with coupling λ V 6 = λ V / √ 2. Table 2 provides the full list of O j V JP operators [40], containing those (j = 1, 2, 3, 5, 6) that contribute to the V −γ−φ a vertex. Notice that the c 3 term is the only one that explicitly breaks the U (3) symmetry, via the tensor χ − , proportional to the quark masses. Nonetheless, the other c j couplings will also enter in the U (3) breaking contributions to the TFFs once the external pseudo-Goldstones are set on-shell and their wave function renormalizations and mixings are taken into account.
• Operators with two vector resonance fields: Operators with two vector fields give quark mass corrections to the resonance mass term (although we have preferred to keep it as a perturbation 4 The operator ∆L = i [37] is not included here in L even V since it does not contribute to the studied processes at tree-level. The same applies to other Lagrangian operators from refs. [37,41,38] not quoted here.  [40], containing those (j = 1, 2, 3, 5, 6) that contribute to the V − γ − φ a vertex.
instead of including it in L Kin V ) and contribute to the V − V − φ a vertices: The quark mass splitting induces a U (3) mass breaking in the vector nonet, in fair agreement with the phenomenology [39]. We note that in ref. [39] the e V m Lagrangian is loosely written as e R m RRχ + for a generic resonance. Table 3 provides the full list of O j V V P operators [40], containing those (j = 1, 2, 3) that contribute to the V − V − φ a vertex. Only the d 2 operator breaks explicitly U (3).
The chiral building blocks are defined as [51] where the φ a fields provide the U (3) nonet of the chiral pseudo-Goldstone fields (the η ′ becomes a chiral Goldstone under the chiral and large-N C limits), the scalar-pseudoscalar densities are set to the quark masses, s+ip =diag(m u , m d , m s ) (isospin symmetry will be assumed from now on, so m u = m d ≡ m q = m s ), and F is the pion decay constant (F π = 92.2 MeV [1]) in the chiral limit. The left and right sources are respectively ℓ µ = v µ − a µ = eQA µ + ... and .. and F µν R = F µν L | ℓ→r = eQF µν + ... are their associated field strengths (with the photon field strength tensor F µν = ∂ µ A ν − ∂ ν A µ ), and the dots stand for gauge boson fields not contributing to the P → γ ⋆ γ ⋆ TFF. For the Levi-Civita tensor we use the standard convention ǫ 0123 = 1.
The main contribution to the TFF is provided by the previous operators, as the pseudoscalar resonances P ′ can only enter via the P ′ −φ a mixing, suppressed by the ratio m 2 P /M 2 P ′ of the light pseudo-Goldstone and heavy pseudoscalar resonance masses. However, the P ′ multiplet is crucial to recover the right behaviour prescribed by the OPE for the VVP Green's function [41,50]. In order to use the short-distance constraints therein, for consistency, one needs to also include the pseudoscalar resonance operators in [41] that contribute to the TFF [37,41]: Our final results will crucially rely on the asymptotic behaviour imposed by QCD on the TFFs at high-energies [45].

U(3) breaking in the Φ and V nonets
A more accurate analysis of the P -TFF can be done by including corrections up to O(m 2 P ), with m P referring to the mass of the pseudo-Goldstones. In order to do this, one needs to take into account all possible corrections to interaction vertices, resonance masses and field renormalizations.
For the vector resonance nonet we will assume an ideal mixing, such that (V µν 11 , V µν 22 , V µν 33 ) = (ρ 0 µν + ω µν )/ √ 2, (−ρ 0 µν + ω µν )/ √ 2, φ µν , as prescribed by the large-N C limit. We will include only the lowest order terms in m q/s that break the U (3) symmetry. The quark mass corrections to the vector masses stemming from the operator of L even V V do not modify this mixing. In the large-N C and isospin limits, this Lagrangian yields the vector resonance mass pattern [39], in fair agreement with the phenomenology, with ∆ 2 2Kπ = 2m 2 K − m 2 π . The K * mass is provided to complete the information on the U (3) splitting of the vector nonet, although it will not be relevant for this work, obviously. Similar relations can be obtained for the pseudoscalar resonances (see ref. [39]), however these will not be needed. These effects would only correct the amplitudes at O(m 4 P ) and are beyond the scope of this work.
In the same way, the λ V term in L even V Lagrangian will introduce an U (3) breaking in the V − γ transitions, such that we will have to make the replacements: To compute the transition P → γγ one has to establish the pseudo-Goldstone mass eigenstates and, in particular, the η − η ′ mixing. For this we will use a two-angle mixing scheme [52,53] consistent with the large N C limit of QCD [54], where the octet and singlet bare fields φ a are related to the physical fields through 1 √ 2 φ a λ a is the pseudo-Goldstone matrix in terms of the bare pseudoscalar fields φ a . These fields are related to the physical ones through a rescaling of the form φ 3 = (Φ 11 − Φ 22 )/ √ 2 = C π π 0 , with C π = F/F π in the large-N C limit [55,56,57]. In the most general case, when loops and other corrections are taken into account C 2 π is just the wave-function renormalization Z π .
Within the two-angle mixing scheme, the mixing constants are parametrized in the most general form through in terms of two couplings f 8/0 and two mixing angles θ 8/0 [58]. The physical fields can be identified through the mixing parameters as a linear combination of the 'bare' SU (3) octet (bare η 8 = φ 8 ) and singlet (bare η 0 = φ 0 ) fields in Φ. In the large N C limit, one has η q = η and η s = η ′ , where the former is an isosinglet combination uū + dd and the latter is a pure ss state. These parameters have been computed in U (3) large-N C χPT at LO [52] (where θ 8 = θ 0 ≈ −20 • and f 8 = f 0 = F ), next-to-leading order (NLO) [52] and next-to-next-to-leading order (NNLO) [59]. These calculations point out an important issue in the treatment of the η and η ′ mesons with the large-N C and/or chiral limit: 1/N C and quark mass corrections to the U (3) pseudo-Goldstone properties happen to be of a similar order of magnitude and compete with each other. For instance, the lowest contribution to the η ′ mass is dominated by the O(N −1 C m 0 q ) contribution from the U (1) A chiral anomaly whereas the O(N 0 C m 1 q ) contributions mostly determine the π 0 mass. It is illustrative to observe the two limits separately: in the large-N C limit, the physical states are given by the ideal mixing angles θ 8/0 = − arcsin 2/3 ≈ −55 • and couplings f 8/0 = F , which imply C q = C ′ s = 1 and C ′ q = C s = 0, all up to quark mass corrections; in the chiral limit, Therefore, in the present work, we will always keep the full mixing coefficients C ( ′ ) q/s , not expanded; on the other hand, the large-N C limit will be assumed everywhere else in our space-like TFF analysis and resonance widths, meson loops or multi-trace RχT operators will be considered negligible for our study. Obviously, these subleading effects must be properly taken into account in the analysis of time-like observables (vector and pseudoscalar branching ratios, e + e − → P γ production, etc.). We are, however, modelling the large-N C limit of RχT with the infinite tower of resonant states truncated at the first multiplet. Still, we expect that this approximation will affect very mildly our prediction for a P,HLbL µ . The reason for this is that the integration kernel of the these contributions is completely dominated by the [0.1, 1] GeV 2 region for both photon virtualities q 2 1,2 [43], which suppresses the contributions from higher resonance excitations. Nonetheless we will use the [0, 100] GeV 2 region for both photon momenta in the form-factor when computing a P,HLbL µ . According to our results, ∼ 85% of the whole contribution comes from the [0.1, 1] GeV 2 region for both photon virtualities. Furthermore, if the upper limit in both squared momenta is set to 400 GeV 2 , the prediction for a P,HLbL µ changes by less than 0.005% with respect to our reference value. In addition to this, we will see in section 5 that current TFF data (extending up to roughly 35 GeV 2 ) do not show any hint of a sizable contribution from excited resonances and that their effect can be captured by a small shift of the lightest vector resonance parameters. It is illustrative to compare the present results with more involved studies including higher resonance poles [46,43,60].

Transition form-factors in RχT
The transition amplitude P (p) → γ ⋆ (q 1 , ǫ 1 )γ ⋆ (q 2 , ǫ 2 ) with an on-shell light pseudoscalar and two photons with virtualities q 2 1 and q 2 2 (polarizations ǫ ⋆ 1 (q 1 ) and ǫ ⋆ 2 (q 2 )) can be parametrized in terms of a scalar function, the form-factor where Bose symmetry implies This amplitude can be divided into three kinds depending on the number of intermediate vector resonances contributing to each diagram. First we have contributions involving no resonances, as that shown in Figure 1, which can be computed by means of the Wess-Zumino-Witten functional [49] (of chiral order p 4 ) and the O(p 6 ) χPT Lagrangian [51]. This gives where the constants C P0 , C P7 and C P8 depend on the pseudo-Goldstone boson considered and are given in table 4.  where the constants C V are and the constants C d 1R and C m 2R are given in table 5.
where the values for the constants C d 2R and C m 2R can be read from table 6. Table 6: Values of C m 2R (up) and C m 2R (bottom) for the couplings between different mesons.
Adding up the various π 0 → γ ⋆ γ ⋆ contributions we obtain the π 0 -TFF, where D R (s) = M 2 R − s. Note that the pion wave-function renormalization changes the global F −1 factor to F −1 π in the large-N C limit [55,56,57].
For simplicity, in addition to the couplings c 3 and d 2 , we will use the following combinations of constants: For the η meson we find the TFF, The result for the η ′ can be obtained from the previous one by substituting C q → C ′ q , C s → −C ′ s and m η → m η ′ . Notice that by taking the chiral limit one recovers the expressions given previously in refs. [41,42].
We would like to remark that the value of pseudoscalar decay constant in the chiral limit, F , is not required in our RχT description. In the π 0 -TFF it combines into the global F −1 π factor in eq. (16), with F π = 92.2 MeV [1]. In the same way, the η and η ′ form-factors only depend on the combinations C q /F , C s /F , C ′ q /F and C ′ s /F which, in turn, only depend on the mixing angles θ 8 , θ 0 and the value of f 8 and f 0 (not on their ratios f 8 /F and f 0 /F ).
Since we want to make use of the V V P Green's function short-distance constraints [41], the analysis is incomplete if the pseudoscalar resonance multiplet, P ′ , -included in ref. [41]-is not considered. The d m operator introduces a mixing between P ′ and the pseudo-Goldstone states P proportional to the quark masses. Thus, the contributions from the heavy pseudoscalar nonet to the P -TFFs begin at O(m 2 P ) and, at that order, take the form in terms of the P ′⋆ (p) → γ ⋆ (q 1 )γ ⋆ (q 2 ) transition amplitude in the chiral and large-N C limits, provided here for the isospin I = 1 resonance by where C P7 = O(m 2 P ) and M P ′ is the P ′ mass in the chiral and large-N C limits. U (3) breaking effects in the P ′ multiplet would contribute to our P -TFFs at higher orders in m 2 P . This result for the off-shell P ′ ⋆ γ ⋆ γ ⋆ transition is in agreement with ref. [41] and the U (3) limit of eq. (19) (given by m π = m K in C P7 ) coincides with U (3) symmetric contribution to the P -TFFs in ref. [42]. Comparing these expressions for the P ′ contributions with those in eqs. (16) and (18), it becomes evident that these contributions can be included merely by a redefinition of three parameters, namely, Note that the RχT couplings c 1,2,5,6 and d 1,3 do not need to be corrected. However, in order to take the P ′ contributions into account, we will also need to shift our auxiliary combinations

Short Distance Constraints
We will constrain the RχT parameters by imposing on the F P γ ⋆ γ ⋆ TFFs (P = π 0 , η, η ′ ) the high energy behaviour prescribed by QCD [45,61]: with Q 2 = −q 2 . These conditions are applied order by order in the m 2 P expansion, separately: first in the massless limit, at O(m 0 P ), and then at O(m 2 P ). We emphasize that in this work we are not matching the precise QCD predictions for the coefficients of the 1/Q 2 expansion [45,61], just requiring the vanishing of the TFF at high momentum transfer. We will see that these results are found 5 One does not need to take into account the corrections from λ V and e V m in eq. (20) for these redefinitions, as we are only interested in the P ′ contributions to the P -TFFs at O(m 2 P ).
to be compatible with previous results [41,50] from the analysis of the VVP Green's function and other observables. The constraints from O(m 4 P ) or higher are not taken into account, as we did not considered RχT operators contributing at that order.
The high-energy analysis of the TFF leads to the relations: • Additional η-TFF constraints, O(m 2 P ): Notice that the implicit κ P V 3 cancels out in the last equation so it should be actually read as c 3 = c 1235 /8. This, in combination with the constraint (24b), yields c 1 = c 2 − c 5 = 0.
• No additional η ′ -TFF constraints: up to O(m 2 P ), the analysis of the η ′ -TFF casts the same relations provided by the π 0 and η form-factors.
We note that we did not need to perform the m 2 P expansion of the mixing coefficients C ( ′ ) q/s . One obtains the previous set of consistent relations for any value of C ( ′ ) q/s : our conditions are actually requiring a good high-energy behaviour to the bare φ a → γ ⋆ γ ⋆ TFFs, which is inherited by the physical π 0 , η and η ′ TFFs after taking into account the field renormalizations and mixings provided by C π and C ( ′ ) q/s . This allows us to by-pass the cumbersome problem of the large-N C limit in the η − η ′ mixing [52], where 1/N C and quark mass corrections are found to be of a similar numerical order [52,53,59].
We will supplement these TFF constraints with the additional conditions derived from the short-distance analysis of the VVP Green's function Π(p 2 , q 2 , r 2 ) [41]. By matching the leading terms of the QCD OPE at p 2 , q 2 , r 2 → ∞ within the chiral and large-N C limits, ref. [41] obtains 6 , In the last line, we added the VVP constraints for the RχT couplings C W j in L non−R . It is remarkable that these relations are compatible with our TFF relations above, being the relations for c 125 , c 1256 , C W 8 and C W 22 identical. Combining these VVP constraints [41] with the previous ones from the high-energy TFF analysis up to O(m 2 P ), we conclude the constraints The relations obtained in this way fix all the parameters in the form-factors except for d ⋆ 123 , d ⋆ 2 , M V , e V m and the four η − η ′ mixing parameters. Particularly, among the operators violating U (3) explicitly, those with coefficients c 3 , C W 7 ⋆ , C W 8 and λ V vanish, according to eq. (28). Although there is a relation for d 123 from the OPE study of the VVP Green's function, d ⋆ 123 is not well known since the κ P V V coupling remains unconstrained [41]. The unrestricted parameters are fitted in this work using F P γ ⋆ γ ⋆ (q 2 , 0) experimental data for π 0 , η and η ′ . The results are given in the next section. There is one more constraint for κ P V 3 from the combination of the TFF and VVP relations for d 3 in the chiral and large-N C limits [41], which, however will not be considered in detail in this 6 In this article we use the notation from ref. [40]. It relates to the Lagrangian [41] through work: In principle, this, together with the previous constraints, leads to a negative vector mass splitting parameter e V m , in agreement with previous phenomenology [39] and our later fits. This prediction e V m = −2π 2 F 2 /(N C M 2 P ′ ) is, however, one order of magnitude smaller for F ∼ F π and M P ′ ∼ 1.3 GeV than the values required by the phenomenology. Nonetheless, one should be very cautious with the theoretical determinations of the P ′ couplings: under the lightest P ′ and V resonance multiplet approximation it is possible to obtain a good high-energy behaviour for the π 0 -TFF and VVP Green's function in the chiral limit but at the price of a P ′ → γ ⋆ γ TFF (20) which does not vanish a large momentum transfer. For this reason we have decided to leave κ P V 3 -or equivalently e V m -as a phenomenological parameter in our later fits, postponing further studies on the structure of the P ′ Lagrangian for a future work.
After using the SD constraints, the π 0 -TFF results simplified into the form Also the η-TFF is simplified to The constrained η ′ -TFF is given by this expression with the same replacements indicated before after eq. (18): In the chiral and large-N C limits one has m P → 0 and we recover the previous result from ref. [43]: and The result for the η ′ TFF is analogous with C q → C q ′ and C s → −C s ′ [42] 7 . We find again at O(m 2 P ) the 7 The values of C q/s will depend on whether we take first the chiral or the large-N C issue pointed out in previous studies of the TFF with only one vector multiplet at O(m 0 P ) [43]: it is not possible to match the 1/q 2 coefficient prescribed by QCD at deep virtual momentum Q 2 → ∞ at the same time for both kinematical configurations F P γ ⋆ γ ⋆ (q 2 , 0) ≈ 2F/q 2 [45] and F P γ ⋆ γ ⋆ (q 2 , q 2 ) ≈ 2F/(3q 2 ) [61]. We will see that, even though the first limit is not imposed, our TFFs fitted to data approximately recover the asymptotic behaviour 2F/q 2 . On the other hand, for the doubly off-shell TFFs we always find F P γ ⋆ γ ⋆ (q 2 , q 2 ) . This deficiency of our description might lead to a slight underestimation of our result for a P,HLbL µ that will be later discussed.
We would like to remark that the P γ ⋆ γ ⋆ and V V P short-distance relations provided here are fully compatible with the constraints from τ → P − γν τ [62] and τ → (P V ) − ν τ [63]. The high-energy conditions from τ → (KKπ) − ν τ [64] and τ → ηπ − π 0 ν τ [65] are also compatible with the relations in the present article provided F V = √ 3F [50]. Nevertheless, the outcomes in this article will not rely on this last relation, even though it will be useful for comparison. We will determine these 8 parameters by means of a fit to the available space-like experimental data. Notice, however, that our fit will be completely insensitive to F V , since it only appears in the form-factors multiplying d ⋆ 2 and d ⋆ 123 . Therefore, for convenience and to ease the comparison we will show the results for the combinations limit, i.e, on the hierarchy of the associated scales. Nonetheless, in the U (3) symmetric limit mq, N −1 C → 0 (also for the a HLbL µ integration kernels), the sum of the η and η ′ contributions to the anomalous magnetic moment is proportional , regardless of the order in which the limits are taken. 8 Instead of F 2 V d ⋆ 123 and F 2 V d ⋆ 2 we could have chosen κ P V V and d 2 for our set of independent fit parameters, related through eqs. (21c) and (28). However, we have preferred to use the former ones for sake of clarity.
Since we are mainly focused on the application of the transition Form Factor in the pseudoscalar exchange contribution to the a HLbL µ , we will only perform fits to data with negative squared photon momenta, q 2 ≤ 0. There are two reasons for this, the first one is that being a NLO effect in the 1/N C expansion, the width of the resonances is needed to obtain a finite result in the q 2 > 0 region of momenta, and to be consistent, we would need to consider also all NLO terms in the 1/N C expansion; the other one is that the observables in the time-like region may get large contributions from disregarded photon radiative corrections [47].
The Γ(P → γγ) decay widths give relevant information of the form-factor at very low energies and help reduce the error in the η − η ′ mixing parameters. The values for the Γ(P → γγ) decays are taken from the Particle Data Group (PDG) [1]. The relation between the form-factor and the on-shell photons decay width is given by We use CLEO [66], CELLO [67], LEP [68], BaBar [44,69] and Belle [70] data of the form-factors along with the decay width to real photons. Since the PDG takes into account LEP data at q 2 = 0 for the η ′ → γγ transition, we have removed this bin to avoid double counting. We fit our form-factors for π 0 , η and η ′ to all these data, simultaneously.
In order to stabilize the fit, we have added a contribution to the χ 2 given by the values of the parameters θ 8/0 and f 8/0 from previous determinations [52,53,58], namely This gives a very small contribution ∆χ 2 ≈ 1.5. These inputs could be understood as a prior distribution for the mixing parameters given by [52,53,58]. If one, however, takes the numbers given in ref.   Table 7: Values for the fitted parameters with (first column) and without (second column) BaBar data on the π 0 form-factor. We fit the parameters MV , e V m and the mixing couplings together with P1 and P2, notd2 andd123; these are reconstructed from the fitted values for P1 and P2 and eqs. (36a) and (36b) and shown here for illustration. The third column shows the parameters fitted excluding BaBar π 0 data and fixing MV and e V m (marked with †) according to refs. [39]. The last row collects the χ 2 per degree of freedom (dof). a P,HLbL µ still compatible with the result from our best fit (see section 7). If one uses F instead of F π in eqs. (35c) and (35d), the values for the parameters are completely compatible with those from our best fit. With this input we also get a result compatible with our best fit a P,HLbL µ (see secs. 6 and 7). If, furthermore, we remove the stabilizing conditions, eqs. (35a)-(35d), we end up with values for the mixing parameters that are incompatible with previous determinations [52,53,58,59]: θ 8 = (−27 ± 4) • , θ 0 = (−5 ± 8) • , f 8 = (220 ± 40) MeV and f 0 = (110 ± 30) MeV (though this leads to a P,HLbL µ = (8.81 ± 0.16) · 10 −10 , compatible with our best fit in secs. 6 and 7). The reason for this is that, as can be seen in table 8 (and tables 14 and 15), there is a large correlation between P 1 and θ 8 and between P 2 and θ 0 and f 0 (in addition to the strong correlations between the four mixing parameters) 10 . Therefore, the fit to the TFF data alone is not really sensitive to the four mixing parameters but to combinations of them and P 1,2 . The stabilizing conditions in eq. (35) lessen the effect of such high correlations avoiding, thus, a non-physical region for the mixing parameters.
We now focus on the fit to TFF data and eqs.   fit with f 8/0 , θ 8/0 , M V , e V m andd 2 andd 123 shows a correlation between the latter two close to one. Hence, instead of consideringd 2 andd 123 , we will fit the parameters P 1 and P 2 provided by the linear combinations with α 2 = −0.0272, α 123 = −0.233, σ d2 = 10.49 · 10 −3 , σ d123 = 81.28 · 10 −3 and r = 0.995. A fair estimate of the constants σ d2 , σ d123 and r that minimize these fit correlations can be extracted from the errors and correlation between d 2 andd 123 , respectively, in the preliminary fit. This transformation greatly reduces the correlation between both fit parameters, improving the efficiency and reliability of the numerical fit. We note that these constants σ d2 , σ d123 and r will change if one considers a different data set or fixes some parameters, e.g., M V and e V m . In this case, the fit to all the available space-like TFF data casts the mean values and errors given in table 7 (first column of results;d 2,123 were reconstructed from P 1,2 by means of eqs. (36a) and (36b)). Further details on the correlations for this 'fit 1'have been relegated to table 14 in App. C.
However, this first fit points out that the π 0 BaBar data [44] is in conflict with η and η ′ -TFF data, related through chiral symmetry, and with Belle data (see e.g. refs. [71,72,73,74,75,76,77,78,79]). Likewise, it is at odds with the asymptotic π 0 form-factor behaviour F π 0 γ ⋆ γ ⋆ ∼ O(Q −2 ) expected in highenergy QCD [45]. Indeed, this set of π 0 data is exclusively responsible for the important deviation in χ 2 (150.) from the number of degrees of freedom of this 'fit 1' (101). Once the BaBar π 0 data set is removed from the analysis, the ratio χ 2 /dof drops and our theoretical expressions provide a good description for all the π 0 , η and η ′ experimental data. Thus, we will take this 'fit 2' as our best fit. The corresponding mean values and marginal standard deviations are given in table 7 (second column of results) and their correlations can be found in table 8. In order to reduce the correlations in'fit 2' (π 0 BaBar data excluded), we have used transformations (36a) and (36b) with the same values for σ d2 , σ d123 and r as 'fit 1' (π 0 BaBar data included).
The values of M V obtained from these fits seem rather at odds with what one expects from the measured mass of the lowest-lying resonances: M ρ = M ω ≈ 807 MeV and M φ ≈ 1126 MeV. In principle, this is not a problem for our analysis, as we are considering a purely space-like analysis within the lightest resonance multiplet and large-N C approximation. Thus, the parameters employed here are not truly those of the full theory with the large-N C infinite tower of resonances. Thus, the lightest vector couplings become slightly shifted to compensate for the heavier resonances missing in our RχT description. This could mean that resonances from heavier multiplets (e.g., ρ ′ ) might be giving a contribution to the experimental data of the form-factor which is not fully negligible in our approach. Therefore, we decided to make yet another fit fixing the mass and the U (3) splitting parameter, excluding the π 0 data from BaBar, in order to study how this affects our results. This 'fit 3' gives the values shown in the third column of table 7, and the correlation among parameters in table 15 in app. C. As done in the previous fits, thed 2 andd 123 parameters were transformed in the same way, however, for this fit the values of the decorrelating transformation constants read α 2 = −2.95 · 10 −2 , α 123 = −0.254, σ d2 = 6.88 · 10 −3 , σ d123 = 5.04 · 10 −2 and r = 0.992. The comparison of the second and third columns of results in table 7 shows only small changes in the fitted parameters and the fit quality, which supports the values obtained for M V and e V m in the first two fits of table 7. Based on this discussion, we consider the second column of results in table 7 as our reference fit.
In fig. 4, we show the comparison between our best fit (excluding the BaBar π 0 TFF) and the fit to all the data (Q 2 = −q 2 ). The π 0 -TFF from the fit to all data is represented by the darker red band with dashed borders and the fit after removing BaBar π 0 data is given by the clearer green band with dotted borders. These bands provide the TFF 1σ uncertainty stemming from the corresponding fit, taking into account correlations. These results are shown together with the experimental data from BaBar [44,69] (full diamonds), Belle [70] (squares), CLEO [66] (empty diamonds) and CELLO [67] (triangles). The outcomes from both fits are very similar for the π 0 TFF and essentially identical for the η and η ′ form-factors. We have also plotted the QCD asymptotic 1/q 2 coefficient fig. 4 (horizontal line), estimated by 2F × C P0 with the central values from our best fit and F ≈ F π for illustration [45] 11 . Since the experimental TFFs (barring BaBar π 0 data) are compatible with this asymptotic behaviour, our theoretical form-factors also approximately agree with it, even though the matching of the asymptotic 1/q 2 coefficient from QCD was not among our high-energy constraints.
The comparison of our best fit with that disregarding the π 0 data from 11 A detailed analysis of the singlet contribution [80] tends to reduce (increase) slightly the value for the asymptotic 1/q 2 coefficient in the η (η ′ ) TFFs (see also the discussion in Ref. [81]), in closer agreement with the trend shown by data. Figure 4: Comparison between the π 0 (top), η (middle) and η ′ (bottom) TFFs from the fit to all data (darker red band with dashed borders) and the fit after removing BaBar π 0 data (clearer green band with dotted borders). See the text for details. Figure 5: Comparison between the π 0 (top), η (middle) and η ′ (bottom) TFFs from the fit with MV and e V m fixed (darker red band with dashed borders) and with these parameters included in the fit (clearer green band with dotted borders). BaBar π 0 data is not fitted in both cases. The remaining notation is the same as in fig. 4.  Table 9: Correlation matrix for the η − η ′ mixing parameters, mean values and marginal standard deviations according to our best fit analysis. These results are compared to the corresponding values from ref. [46].
BaBar but fixing M V and e V m is given in fig. 5. The TFF from the fit with M V and e V m fixed is given by the darker red band with dashed borders. It is compared to our best fit, where M V and e V m are fitted together with the other six parameters, provided by the clearer green band with dotted borders. The remaining notation is the same as in the previous plots in fig. 4. The results for the η and η ′ TFFs happen to be compatible -though less precise for our best fit-. On the other hand, there is some small discrepancy in the π 0 -TFF: the 1σ bands do not overlap and the asymptotic value for the π 0 form-factor is slightly smaller in the case with fixed M V and e V m , yielding a poorer χ 2 /dof. Notice that the worsening of the fixed-(M V , e V m ) fit is due to the π 0 -TFF.
The departure of all three values ford 123 in table 7 from the prediction d 123 = 1/24 ≈ 4 · 10 −2 [50] shows the impact of pseudoscalar resonances in this particular coupling. Our fitted values ofd 2 exhibit a similar deviation with respect to those obtained in ref. [82], where d 2 ∈ [4, 7] · 10 −2 (see also ref. [83]), pointing again to the relevance of pseudoscalar resonance contributions to V V P vertices.
We also evaluated the mixing parameters C (′) q and C (′) s according to our best fit (BaBar π 0 data excluded). Since our fit determines f 8 and f 0 rather than their respective ratios f 8 /F and f 0 /F with the chiral limit decay constant F , we define the quantities C q/s × (F π /F ). These quantities only depend now on the angles θ 8/0 and the ratios f 8/0 /F π and can be easily translated into the actual mixing coefficients C ( ′ ) q/s once the value of the ratio F π /F is provided. They are given in table 9, where the mean values, the marginal standard deviations and the correlation matrix are quoted. The errors have been propagated from the fit parameters taking into account correlations by means of a Monte-Carlo.

Comparison with other recent TFF determinations
Our description of the Transition Form Factor provides fairly simple expressions that can be implemented in the computation of P γγ observables, despite the lack of NLO contributions in the 1/N C expansion (barring those in the η − η ′ mixing). This derives from the also simple short distance constraints (28), contrary to those in ref. [46], where much more involved ultraviolet restrictions are reported. These are the result of a more general U (3) breaking pattern, not introduced in the chirally covariant form proposed in this article, and the inclusion of several vector multiplets. Their fits are done to observables in both q 2 regimes, time-like and space-like, without considering photon radiative corrections [47] which may become important for some observables.
Czyż et al. [46] This work This model is analogous to ours in the sense that all the parameters in our model can be expressed in terms of those in ref. [46], when restricted to the first multiplet given by i = 1 12 . In our approach the Lagrangian couplings are quark mass independent and the structure of m q/s corrections is dictated by the operators and the kinematics ((q 1 + q 2 ) 2 = p 2 = m 2 P ). On the other hand, the effective parameters in [46] encoded an important part of the quark mass corrections. Note, however, that the U (3) splitting in ref. [46] is not fully general and assumes some restrictions: their vector-photon-pseudoscalar (h Vi ) and vector-vector-pseudoscalar (σ Vi ) couplings are the same for the three pseudo-Goldstones. The equivalences are given in table 10. In order to do the comparison, we have translated our resonance Lagrangian in the antisymmetric tensor V αβ formalism into the Proca four-vectorV ρ realization considered in [46] through the transformation V µν i → −(∂ µV ν i − ∂ νV µ i )/M Vi for each vector V i = ρ, ω, φ, ρ ′ ... [85,86]. Once our L RχT is expressed in terms of Proca fields, the operators are set with their particles on-shell and identified with the corresponding notation in [46].
In table 11 we give the numerical values of the parameters in table 10 for both descriptions ('fit 2' results are used in the comparison). It can be seen that, in spite of the differences in both descriptions, a very good agreement is found among the numerical values for table 11. It is worth mentioning that, despite the very different results for θ 8 , θ 0 , f 8 and f 0 , the mixing parameters C (′) q are in good agreement but C s and C ′ s show some discrepancy (more significant for the latter), as it can be seen in table 9. However, if the relations of eqs. (35) are not imposed in our fits, the values obtained for the mixing parameters in this way get closer to those in ref. [46].
Czyż et al. [46] This work f V1 0.2020 ± 0.0008 0.198 ± 0.003 F ω1 0.88 ± 0.01 1 F φ1 0.783 ± 0.005 0.72 ± 0.04 h V1 0.0377 ± 0.0008 0.0326 ± 0.0005 H ω1 1.02 ± 0.03 1 A π0 On the other hand, descriptions of the TFF such as those in ref. [60] by means of Padé approximants provide a neat and simple approach, which can also incorporate asymptotic QCD information [45,61]. However, the lack of a Lagrangian in such method makes it complicated to combine information from the π 0 and η ( ′ ) TFFs. As a result, one needs to perform a separate Padé analysis for each quantity, needing a set of similar statistical quality data to describe a closely related process with a comparable accuracy.

B(P → γ (⋆) γ (⋆) ) predictions vs. data
As a quality check of the fitted parameters and the RχT description, we give our estimates for different branching ratios related to P → γ (⋆) γ (⋆) processes, uncertainties are naively obtained by independently varying the RχT parameters within their 1σ ranges without considering correlations among them 13 . These are provided by our best fit in the table 12 (all parameters are floated and BaBar π 0 data is excluded). The width of the resonances is needed to obtain finite results. A momentum dependent width is employed for the ρ meson, following ref. [87]; for the ω and φ, we use the total widths from PDG [1] neglecting any off-shellness dependence, as they are rather narrow. This is accomplished by changing the denominator of the corresponding resonance propagator to D R (s) = (M phys R ) 2 − s − iM R Γ R (s), with the physical masses M phys R taken from the PDG [1]. Except for the resonances masses in the propagators D R (s), all parameters correspond to those from our best fit. For the computation of our predicted branching ratios we have divided our RχT results for the partial widths by the experimental total decay width (no theoretical prediction is considered here for the latter).
The second and third column of table 12 agree at less than two standard deviations in all cases (theory undertainties are not discussed in these estimates), which corroborates our best fit results. We note that no experimental limit is known for the decays in the last three lines of table 12. We emphasize that our space-like analysis is missing several features which may be crucial in time-like observables, such as subleading 1/N C effects as, e.g., resonance widths and the impact of higher resonance multiplets. Thus, the RχT branching ratios in table 12 should be taken with a grain of salt. Nonetheless, the fair agreement with data is remarkable, giving support to the hypotheses and approximations assumed in this work and the final results for the anomalous magnetic moment in the next section.

Pole prediction with one vector resonance multiplet
For the evaluation of the pseudoscalar pole contribution to the HLbL, a P,HLbL µ , we used the expressions for the loop integrals given in ref. [43] 14 . Dispersion relations [29,30] show that the HLbL is determined by the various physical absorptive channels of the V V V V Green's function of four electromagnetic currents. The lightest absortive cut is given by the meson pole topologies γ ⋆ γ ⋆ → P → γ ⋆ γ ⋆ . Additional intermediate states (P P cuts, multiparticle states or with heavy resonances, etc.) will add further corrections to a HLbL µ . These were effectively included before in the so-called meson-exchange contributions (with off-shell mesons) [5]. In a dispersive framework, however, the analyticity structure of the amplitudes is encoded in their poles and cuts alone, in such a way that their residues and imaginary parts (discontinuities) are uniquely related to on-shell form factors and scattering amplitudes [29,30].
Notice that the O(m 2 π ) correction to the form-factor from considering a nonvanishing pseudo-Goldstone mass given in ref. [42] is automatically included in our analysis. Since no U (3)-splitting parameter was considered in the analysis of the form-factor in ref. [42], such correction was underestimated. This contribution accounts for a relative variation of the form-factor of ∆ ∼ −2.5 · 10 −2 , 15 while in ref. [42] this correction is reported to be ∆ ∼ 5.9 · 10 −3 . However, as stated there, additional corrections come as further suppressed powers of m 2 P .
The total pseudo-Goldstone pole contribution is estimated by means of a MonteCarlo run with 5000 events, which randomly generates the eight fit parameters with a normal distribution according to their mean values, errors and correlations between parameters. It integrates at the same time all three contributions from π 0 , η and η ′ exchanges, thus taking into account the correlations between the π 0 , η and η ′ TFF contributions. Our best fit ('fit 2' -excluding BaBar π 0 data-) leads to our final result a P,HLbL µ = (8.47 ± 0.16) · 10 −10 , which perfectly agrees with previous results [41,42,60,43,89,90,99,91,92,93,94,95,96,97,98,26,27] 16 and with a reduced uncertainty. In general, all evaluations are in agreement, as can be seen from table 13. Noticeable exceptions are (we consider the last result from every group as the reference one) those using Melnikov-Vainshtein short-distance constraints [99], obtaining ∼ 13.5 · 10 −10 and the results obtained within the non-local chiral quark model by Dorokhov et al., ∼ 5.85 · 10 −10 [95]. Although the lattice result of ref. [26] (5.35 ± 1.35) · 10 −10 (see also refs. [28,34,100,101]), might seem at odds with other determinations, one still has to take into account potentially large finitevolume systematics, finite-lattice-spacing corrections and the limited statistics for the leading disconnected contribution, all of which are being refined and could bring in reasonable agreement this result with the remaining predictions. From these numbers and comparing ours with those previously quoted, it seems that further orders in m 2 P might be neglected in our analyses. For comparison we also provide the pole contribution for the anomalous magnetic moment stemming from the other two types of fits done in this article: Including BaBar π 0 data 'fit 1': a P,HLbL µ = (8.58 ± 0.16) · 10 −10 .
Fixing M V and e V m 'fit 3': a P,HLbL µ = (8.50 ± 0.13) · 10 −10 , The values given by BaBar for the π 0 -TFF at high energies are much larger than the form-factor derived from the η and η ′ TFFs through chiral symmetry. This leads to a slightly higher value for a P,HLbL µ in comparison to our best fit, with BaBar π 0 data excluded, which shows a very good compatibility between the π 0 , η and η ′ data. Both results are still compatible within errors. On the other hand, the fit where the vector masses are fixed gives a 1σ confidence interval fully included in the 1σ interval of our best determination (37), although with a slightly smaller error, as expected. Given the abovementioned tension between QCD-driven predictions and BaBar π 0 transition form-factor data, we consider the result (37) excluding these data from the fits as our reference value, a P,HLbL µ = (8.47 ± 0.16) · 10 −10 .
Since the separate contributions of the P mesons are interesting in their own (i.e., some papers only consider the π 0 contribution, which is highly restricted by chiral symmetry, contrary to the η (′) cases; data quality varies from channel to channel, etc.), we quote their values in the following. The pole contributions from each separate pseudo-Goldstone exchange are computed for our best fit via a 5000 event MonteCarlo in the same way described before, giving the following values for our best fit: Notice that despite the slightly higher (absolute) uncertainty in the π 0 contribution 18 , the η has an uncertainty reduced by a factor of four compared with previous determinations [42]. Also, the η ′ contribution has a mildly reduced uncertainty. We note that the sum in quadrature (uncorrelated) of the errors for the individual contributions in the preceding equations yields an error estimate (∼ 0.13 · 10 −10 ), which is a bit smaller than the one in eq. (37) because of the 17 Since a Gaussian distribution accounting for correlation is used to compute a P,HLbL µ , there are some points of the generated phase space that give negative squared resonance masses, leading to spurious divergences. Therefore, these points were dropped in order to obtain a finite result. 18 Our result for this contribution is in agreement at the 1. correlations, accounted in the latter: performing the simultaneous integral of the three contributions leads to eq. (37). As a check of the size of these corrections, we have additionally computed the P pole contributions using our form-factor in the chiral and large-N C limits. We have taken the central values of the vector mass M V and mixing parameters from our best fit, though keeping the physical pseudo-Goldstone masses in the integration kernels. As a result, we find (F/F π ) 2 a P,HLbL µ = 8.27 · 10 −10 , where the change is essentially given by the a η,HLbL µ contribution. For F ≈ F π , the small change in the mean value of the pseudo-Goldstone pole contribution in the chiral and large-N C limits (∼ −2.5%, up to corrections in F/F π ) suggests that NNLO corrections, suppressed by further powers of m 2 P , must be negligible. An effect on this ballpark could already be inferred by comparing a P,HLbL µ in the U (3) symmetric analysis of ref. [42] with our present result. To check the compatibility between the present fit and that done in ref. [42] we compare the values obtained for the relevant combination F 2 V d 3 in the TFF: F 2 V d 3 /(3F 2 π ) = (−116.2 ± 1.8) · 10 −3 for this work and F 2 V d 3 /(3F 2 π ) = (−110.7 ± 8.3) · 10 −3 for that in ref. [42], which show a good agreement between both determinations.
The study of the so-called off-shell pole contribution to a HLbL µ will not be discussed in this article. The P -TFF with an off-shell pseudo-Goldstone (p 2 = m 2 P ) is by construction an ill-defined quantity (arbitrary off-shell contributions can be obtained through pseudoscalar field redefinitions). One should actually rather analyze the Green's function of four electromagnetic currents T {J µ EM (0)J ν EM (x)J α EM (y)J β EM (z)} , which is free of these ambiguities. However, its study within the RχT framework is a more cumbersome problem than just accounting for the pseudo-Goldstone tree-level exchanges provided by the TFFs, and has been postponed for a future work.

Further error analysis
As said previously, only the leading order terms in 1/N C have been considered in the computation of a P,HLbL µ . However, one must include a non-zero width in the vector resonance propagators in order to get finite values for the branching fractions of subsection 5.3. In particular, the ρ meson width plays the most important role. Intermediate ππ and KK loops account for the main NLO corrections in 1/N C to the ρ propagator [87]: where being σ P (q 2 ) = 1 − 4m 2 P q 2 . Note that A P (q 2 ) is real for q 2 < 4m 2 P . Thus, the ρ propagator provided by eq. (40) is real in the whole space-like region q 2 < 0.
We perform the NLO replacement in the ρ propagators given in eq. (40) of the functions g Vi (q 2 ) in app. B, which enter in the a P,HLbL µ integral representation [43]. The ω and φ propagators in the g Vi (q 2 ) are left unchanged. Likewise, we keep f (q 2 ) = 0, as it is found at LO in 1/N C after imposing the short distance constraints. This leads to a decreasing in the theoretical prediction, a P,HLbL µ | LO+NLO − a P,HLbL µ | LO = −0.09 · 10 −10 , essentially dominated by the contribution in the virtuality range [0.1, 1] GeV 2 , as before. This is, nonetheless, just one of the possible NLO corrections in 1/N C to the anomalous magnetic moment. One-loop modifications to the π 0 V V ′ vertex can be, e.g., equally important in the space-like domain and may lead to a positive contribution to a P,HLbL µ . Thus, we take the absolute value of this shift as a crude estimate of the 1/N C effects: Other source of unaccounted error originates in the lightest meson dominance assumption of the present work. As shown in [43], a TFF with only one vector multiplet fails to reproduce at the same time the correct asymptotic behaviour q 2 F P γ ⋆ γ ⋆ (0, q 2 ) = −2F [45] and q 2 F P γ ⋆ γ ⋆ (q 2 , q 2 ) = −2F/3 [61] for q 2 → −∞. In our one-multiplet study, we were able to reproduce F P γ ⋆ γ ⋆ (0, q 2 ) ≈ −2F/q 2 at large momentum transfer but the doubly off-shell form-factor behaved like . However, both QCD limits can be correctly recovered by considering a second multiplet of vector resonances. In the chiral limit, the π 0 -TFF including ρ and ρ ′ takes then the form form-factor follows exactly the QCD asymptotic behaviour from refs. [45,61], the two-multiplet TFF only depends on two free resonance couplings, which here have been chosen to be c (the superindex shows the multiplet to which it couples to). One of the short distance conditions is that If one assumes that the lightest multiplet dominates in this high-energy relation, then one would expect that each multiplet contribution cancels on its own, yielding c = −N C M 2 ρ /(64π 2 ) from eq. (28) in sec. 4. We have taken the chiral limit π 0 -TFF and the inputs M ρ = 0.77 GeV, M ρ ′ = 1.45 GeV and F ≃ F π in the a π 0 ,HLbL µ integral representation [43], while keeping the physical pion mass in the integration kernels. Hence, we observe that the second vector multiplet increases a π 0 ,HLbL µ by an amount ∼ 0.5 · 10 −10 . This increasing is not unexpected: since our single-resonance F π 0 γ ⋆ γ ⋆ (q 2 , q 2 ) vanishes too fast at high energies, it provides an underrated prediction in comparison with the asymptotic QCD TFF and, consequently, an underrated anomalous magnetic moment prediction. We therefore expect that all modifications of the latter TFF that arrange its asymptotic behaviour will push a P,HLbL µ upwards. We thus have the rough estimate, (∆a P,HLbL We assume the correction to the π 0 contribution as the dominant effect and neglect here the corrections in a η,HLbL µ and a η ′ ,HLbL µ [43]. A more rigourous and detailed analysis will be performed in a future work [88], where we will consider all the short distance constraints beyond the chiral limit and extract the free parameters from a fit to the experimental data. This issue has also been discussed in [43]: the difference between the VMD result (reproduced here by (32) and failing to fulfill the OPE constraint [61]) and the LMD+V rational approximation (with a good short-distance behaviour) was used to estimate an uncertainty ∆a π 0 ,LbL µ = ±1.0 · 10 −10 . This error was essentially dominated by the uncertainty in one of the LMD+V parameters, namely h 2 [43].

Conclusions
We have given a more accurate description of the form-factor by including terms up to order m 2 P , for the first time in a chiral invariant Lagrangian approach with the quark mass corrections introduced covariantly. In addition to chiral symmetry, the high-energy conditions from the TFF (up to O(m 2 P )) and the VVP Green's funtion (at O(m 0 P )) were crucial to fix the various unknown Lagrangian parameters. These short distance constraints were found to be consistent with previous determinations in the chiral and large-N C limits. We have also been able to fix various V JP couplings for the first time: c 1 = c 2 − c 5 = c 3 = 0.
By fitting to the different sets of data for the TFF of the pseudo-Goldstones we were able to confirm that BaBar data for the π 0 -TFF is not compatible with measurements of the η (′) transition form-factors. However, in order to stabilize the fit, some prior distribution for the mixing parameters were provided. Otherwise, the fit leads to an η − η ′ mixing in strong conflict with current phenomenology -though still yielding a fair a P,HLbL µ determination-. After assuming previous phenomenological determinations [52,53,58] as inputs, we obtain as a by-product of our analysis, the correlated η − η ′ mixing parameters, given in table 9, which may be useful for future analyses involving these mesons. The contribution from the four mixing inputs f pheno 8/0 and θ pheno 8/0 in eqs. (35a)-(35d) to the χ 2 turns out to be very small (∆χ 2 ∼ 1.5) and yielded values very similar to previous determinations of such parameters, correlated and with a reduced marginal error. These results depended very mildly on the mixing inputs.
Comparing our results with other recent determinations, we find that the TFFs presented in this article turn into rather simple expressions after demanding the correct high energy QCD behaviour while, at the same time, they contain the corresponding O(m 2 P ) corrections and implement chiral symmetry. We also have obtained values for the parameters that are compatible with those obtained in ref. [46] despite the lack of heavier copies of the resonances in our approach. This shows that, even though these heavier multiplets are not taken into account, our result should be compatible with the a P,HLbL µ considering these resonances.
We have determined the pseudo-Goldstone pole contribution to the a HLbL µ with an improved precision, a P,HLbL µ = (8.47 ± 0.16) · 10 −10 , with respect to previous works by using a more accurate description of the Resonance Chiral Lagrangian: we have considered corrections up to O(m 2 P ) in the form-factor within the lightest V and P resonance multiplet approximation. We have performed some checks which suggest that contributions suppressed by higher powers of m 2 P are negligible at the present level of precision.
Future works will be directed towards an improved estimate of the impact from higher resonance multiplets and 1/N C corrections, assumed negligible in most of this work. Nonetheless, a rough estimate of possible further uncertainties was performed in subsection 6.2, yielding a P,HLbL µ = ( 8.47 ± 0.16 sta ± 0.09 1/NC +0.5 −0 asym ) · 10 −10 , where the first error (sta) comes from the fit, the second one (1/N C ) from the estimate of NLO effects in 1/N C and the last one (asym) stems from the incorrect F P γ ⋆ γ ⋆ (q 2 , q 2 ) asymptotic behaviour and the impact of higher resonance multiplets. Likewise, a more throrough HLbL determination should go beyond the meson-pole contribution, demanding a future study of the V V V V four-point Green's function.
B The P form-factor in the evaluation of a HLbL µ For the evaluation of a HLbL µ the form-factor is conveniently written as follows [43] F P 0 γ ⋆ γ ⋆ (q 2 1 , q 2 2 ) = where, for P 0 = π 0 the f and g functions are given by where M V ′ i = M ω δ ρVi + M ρ δ ωVi , being δ the Kronecker delta. Notice that since in the π 0 decays, no photon comes from a φ meson at the large N C limit (as can be seen in the values of C d 1R , C m 1R , C d 2R and C m 2R in tables 5 and 6) the sum in the π 0 -TFF must only contain the ρ and ω resonances.
For the η-TFF one gets the following f and g functions where f ω (q 2 ) = 9f ρ (q 2 ), and Note that the f ρ,ω,φ and h ρ,ω,φ in these equations do not refer to the constants employed in ref. [46] and discussed in Tables 10 and 11.
Just as explained for the case of the η ′ −T F F below eq. (18), one can obtain the expressions for the η ′ in the form given by eq. (50) by replacing C q → C ′ q , C s → −C ′ s and m η → m ′ η including ∆ 2 ηπ → ∆ 2 η ′ π and ∆ 2 2Kπη → ∆ 2 2Kπη ′ . It is worth noticing that, in agreement with refs. [41,42], by imposing the short distance constraints, one obtains f (q 2 ) = 0 , for the three pseudo-Goldstone bosons, which greatly simplifies the numerical computation of a P,HLbL µ .

C Additional information on the fit correlations
We collect in this appendix the correlations between fitted parameters in the two main alternative fits discussed in the paper: i) 'fit 1', in which all sets of data are analyzed ; ii) 'fit 3', where BaBar π 0 TFF data is excluded from the fits with M V and e V M fixed to the values from [39]. These can be read in Tables  14 and 15, respectively.  Table 14: Correlation of fitted parameters including the π 0 form-factor BaBar data, 'fit 1'.  Table 15: Correlation of fitted parameters with fixed M V and e V m excluding the π 0 form-factor BaBar data, 'fit 3'.