Contribution of neutral pseudoscalar mesons to $a_\mu^{HLbL}$ within a Schwinger-Dyson equations approach to QCD

A continuum approach to Quantum Chromodynamics (QCD), based upon Schwinger-Dyson (SD) and Bethe-Salpeter (BS) equations, is employed to provide a tightly constrained prediction for the $\gamma^{*} \gamma^{*} \rightarrow \{ \pi^0, \eta, \eta', \eta_c, \eta_b \}$ transition form factors (TFFs) and their corresponding pole contribution to the hadronic light-by-light (HLbL) piece of the anomalous magnetic moment of the muon ($a_\mu$). This work relies on a practical and well-tested quark-photon vertex Ansatz approach to evaluate the TFFs for arbitrary space-like photon virtualities, in the impulse approximation. The numerical results are parametrized meticulously, ensuring a reliable evaluation of the HLbL contributions to $a_\mu$. We obtain: $a_{\mu}^{\pi^0-\textrm{pole}} = (6.14 \pm 0.21) \times 10^{-10}$, $a_{\mu}^{\eta-\textrm{pole}} = (1.47 \pm 0.19) \times 10^{-10}$, $a_{\mu}^{\eta'-\textrm{pole}} = (1.36 \pm 0.08) \times 10^{-10}$, yielding a total value of $a_{\mu}^{\pi^0+\eta+\eta'-\textrm{pole}} = (8.97 \pm 0.48) \times 10^{-10}$, compatible with contemporary determinations. Notably, we find that $a_{\mu}^{\eta_c+\eta_b-\textrm{pole}} \approx a_{\mu}^{\eta_c-\textrm{pole}} = (0.09 \pm 0.01) \times 10^{-10}$, which might not be negligible once the percent precision in the computation of the light pseudoscalars is reached.


I. INTRODUCTION
More than half a century after the advent of the Standard Model (SM) of particle physics, it has successfully withstood a continuous barrage of innumerable experimental tests. Many of us are keenly interested in high precision measurements of quantities which can be theoretically best calculated in order to zoom into the very limits of this model, hunting for the possible discrepancies. Measurement and calculation of the muon anomalous magnetic moment, a µ = (g µ − 2)/2, provide precisely such battleground [1][2][3]. The most recently reported value by the Brookhaven National Laboratory (BNL), 116592091(63) × 10 −11 [4] shows a persistent 3.5 standard deviations away from the SM prediction 116591823(43)×10 −11 [5]. Well deserved attention is currently being paid to this anomaly 1 due to the ongoing rigorous experimental endeavours to pin it down with increasing precision. The dedicated FNAL experiment will reach a fourfold improvement of the current statistical error within about two years from now [10]. Later on, J-PARC also plans to achieve a comparable accuracy [11]. If this deviation does not wither away, it would be highly desirable to reduce the SM calculational uncertainty as much as possible to be able to associate the discrepancy with possible new physics. What stand on the way are the hadronic contributions which are hard to tame and severely restrain our efforts to make predictions with the desired exactitude.
The SM prediction includes quantum electrodynamics (QED) corrections up to five loops [12][13][14], two-loop (and leading-log three-loop) electroweak ones [15,16] and hadronic contributions, the latter saturating the error of the SM precision quoted above. These are divided into hadronic vacuum polarization and hadronic light-by-light (HLbL) contributions. While the former could be related to data already in 1961 [17], a similar data-driven extraction is not yet possible for the HLbL piece, although a dedicated effort [18][19][20][21][22][23][24][25][26] has made remarkable advances towards reaching this goal in the near future.
In this paper, we compute the HLbL contributions coming from the light neutral pseudoscalar transition to two photons (and the first ever estimation for η c and η b ) to the anomalous magnetic moment of the muon. We fol-low a novel Schwinger-Dyson and Bethe-Salpeter equations (SDEs, BSEs) approach to compute γγ * → M [58][59][60] transition form factor for arbitrarily large space-like momentum for the first time, in a unique framework with a direct connection to quantum chromodynamics (QCD). Such an approach is known to unify those form factors with their corresponding valence quark distribution amplitudes [61][62][63], charged pion and kaon form factors [64][65][66], their parton distribution functions [67][68][69] and a wide range of other hadronic properties (masses, decay constants, etc.) [70][71][72]. We extend the SDE-BSE treatment of Refs. [58][59][60] to account for arbitrary space-like virtualities of both photons.
Different plausible parametrizations of the numerical data are discussed. In particular, the flaws and strengths of Vector-Meson and Lowest-Meson Dominance (VMD, LMD) parametrizations [42] as well as Canterbury Approximants (CAs) [73] are analyzed. In the context of a HLbL µ , the latter were presented in Ref. [55] and also in a recent, but different, SDE-BSE approach [74]. As explained therein and below, we find the CAs parametrization more adequate.
This article is organized as follows. In Sec. II, we introduce the basics of our extension of Refs. [58][59][60] to the doubly-off shell (DoS) case of the TFFs. The proposed parametrizations of the numerical data are presented in Sec. III, together with the implications of the low and high energy behavior of the TFFs and the corresponding constraints. This framework is then applied to π 0 , η and η cases. This section ends with the corresponding description of the heavy η c and η b mesons. In Sec. IV we discuss our results for a HLbL µ . Based on this analysis, we present our conclusions in Sec. V.

II. SDE-BSE APPROACH
The transition γ * γ * → M is described by a single form factor. In the impulse approximation [58], where Q 1 , Q 2 are the momenta of the two photons and (Q 1 + Q 2 ) 2 = P 2 = −m 2 M (m M is the mass of the pseudoscalar). The kinematic arrangement is q 1 = q + Q 1 , q 2 = q − Q 2 ; with q being the integration variable 4 . In addition, e M = e c M is a charge factor associated with the valence quarks of the given meson (e is the charge 4 For simplicity in the notation, we have defined q ≡ of the positron) 5 . The other symbols carry their usual meanings: is the propagator of the f -flavoured quark. It is determined from its SDE (namely, the gap equation): where K(q, p) is the kernel of the gap equation and {r, s, t, u} are color indices (not displayed when obvious). Quark propagator, and every other Green function involved in its SDE, are renormalized at the resolution scale of ζ = 2 GeV := ζ 2 .
• Γ M (p; P ) is the Bethe-Salpeter amplitude of the pseudoscalar meson M , obtained from its BSE: where P is the total momentum of the bound state and χ M (q; 1]. No physical observable depends on η, the definition of the relative momentum. Here M(q, p; P ) is the renormalized, fully amputated, two particle irreducible, quark-antiquark scattering kernel. It is related to K via the Axialvector Ward-Takahashi identity [75,76].
In conjunction with Eq. (2), we employ the so-called Rainbow-Ladder truncation (RL), which is known to accurately describe the pseudoscalar mesons [58-60, 64, 69]. This entails: where k = p−q, D 0 µν (k) is the tree-level gluon propagator in the Landau gauge and G(k 2 ) is an effective dressing function. We employ the well-known Qin-Chang interaction [78], compatible with our modern understanding of the gluon propagator [79][80][81][82]: its dressing function saturates in the infrared and monotonically decreases as the momentum increases and recovers the perturbative QCD running coupling in the ultraviolet. With the interaction strength (ωD = m 3 G ) fixed, all physical observables are practically insensitive to variations of ω ∈ (0.4, 0.6) GeV [66,78]. Moreover, sensible variations of m G do not alter the output observables significantly either, a fact that is illustrated through the computation of the pion mass and decay constant with two parameter sets (all these details are taken into account in our final result). Here on, we shall employ isospin symmetry m u = m d := m l . Typically, m G and the current quark masses are fixed such that ground-state masses and decay constants are properly reproduced [74,78]. Thus, for our first set of parameters (RL-I), we employ {m π , m K , f π , f K } as benchmarks; for the second set (RL-II), we use {m π , m ηc , m η b , f π } instead. Precise input parameters, computed masses and decay constants are given in Table I.

A. Quark-photon vertex
In principle, the QPV can be obtained from its inhomogeneous BS equation. This process automatically incorporates vector meson poles [77,83] in the vertex and guarantees the preservation of the Abelian anomaly [84][85][86]. However, it also limits the domain in which a direct evaluation of the TFFs is possible [74,77,87,88]. Thus we follow an alternative route. We introduce a reliable QPV Ansatz based upon gauge covariance properties and multiplication renormalizability of the massless fermion propagator. In conjunction with this approach, we incorporate the non-Abelian anomaly at the level of the BSE. This approach sets apart our work (and our previous ones [58][59][60]) from the recent SDE approach of [74,77].
A kindred version of the vertex Ansatz we employ was first introduced in [64], for the calculation of the pion elastic form factor (EFF) and subsequently adapted in [58] for the TFFs. The QPV is expressed completely via the functions which characterize the dressed quark propagator (q = k f − k i ,s = 1 − s): where Longitudinal pieces alone do not recover the Abelian anomaly, since it turns out impossible to simultaneously conserve the vector and axial-vector currents associated with Eq. (2). Thus a momentum redistribution factor is introduced: where M f E = {p|p 2 = M 2 f (p 2 ), p 2 >0} and M f (p 2 ) is the quark's mass function. As s is exponentially suppressed, it does not affect the large Q 2 behavior of the TFFs. To account for Q 2 2 = 0, the simplest symmetrization corresponds to the replacement Q 2 1 → Q 2 1 + Q 2 2 , which clearly recovers all the limits 6 . The way the specific values of s f are established is addressed in Sections II C and IV.
The following pattern is observed: Given the proposed form in Eq. (6), the QPV is determined by the quark propagator dressing functions. Consequently, the TFFs are fully expressed in terms of quark propagators and BS amplitudes, obtained in the RL truncation. We then employ perturbation theory integral representations (PTIRs) for those objects, as previously done in the calculation of the pion distribution amplitude [61] and its EFF [64]. The particular representations were presented in [61] and discussed herein in the appendix. PTIRs allow us to write Eq. (2) in terms of objects which have q-quadratic forms in the denominator. Thus, after introducing Feynman parametrization and a suitable change of variables, the 4-momentum integrals can be evaluated analytically. Subsequently, integrations over the Feynman parameters and the spectral density are performed numerically. The complete calculations basically require a series of perturbation-theory-like integrals. This expedites the computation considerably and allows a direct evaluation of the TFFs in the whole domain of space-like momenta. We thus managed, for the very first time, to compute the γ * γ * → All lowest-lying neutral pseudoscalars TFFs in this domain. In using the proposed QPV Ansatz, we overcame the inconvenience stemming from solving its BS equation and it expedited the computation of the TFFs. Despite lacking explicit non-analytic structures, associated with vector meson poles in the time-like region, the effects in the space-like region are appropriately reproduced (later discussed in connection with the charge radius). Thus we expect our approach to be valid by essentially maintaining the key quantitative details. Our Ansatz follows from using the gauge technique [89]. Thus it satisfies the longitudinal Ward-Green-Takahashi identity (WGTI) [90][91][92], is free of kinematic singularities, reduces to the bare vertex in the free-field limit, and has the same Poincaré transformation properties as the bare vertex.

B. The η − η case
Our dealing with the η − η mesons is now discussed. First, we use a flavour basis to rewrite the BS amplitudes as follows: where we keep using the isospin symmetric limit, such that l = u, d. The RL kernel by itself does not produce any mixing between the pure ll and ss states. Thus, the Bethe-Salpeter kernel is improved by including the non-Abelian anomaly kernel (see Refs. [60,93]): with χ l = M l (0) and θ A controlling the relative strength between the γ 5 and γ 5 γ · P terms; It models a dependence on U (3) flavour-symmetry breaking arising from the dressed-quark lines which complete a 'U-turn' in the hairpin diagram (see Fig. 1 in Ref. [60]). The strength of the anomaly is controlled by Here ω ξ and D ξ provide a momentum dependence for the anomaly kernel, as a generalization to that introduced in Ref. [93]. The set of Dirac covariants which describe Eq. (10) can be inferred from the axial-vector WGTI [93]; we keep those which dominate. The rest of the pieces are not determined by the WGTI, but they can be driven by phenomenology [60]. Since the RL truncation does not produce any mixing by itself, it is natural to require more input to describe the η − η system. In particular, given the anomaly kernel of Eqs. (10)-(11), we fix D ξ , ω ξ and cos 2 θ ξ to provide a fair description of m η,η and f l η,η 7 . More weight is given to the masses, which are better constrained empirically. Input and output values, together with the RL counterpart, are listed in Table I. As a reference, if a single mixing angle scheme (and a pair of ideal decay constants) is assumed, our results yield f l ≈ 1.08f π , f s ≈ 1.49f π and φ ηη = 42.8 • . Moreover, note that D ξ = 0 turns off the non-Abelian anomaly and produces an ideal mixing with pure ll and ss states, which implies m η = m π = 0.135 GeV and m η = m ss = 0.698 GeV (f l := f π = 0.093 GeV, f s := f ss = 0.134 GeV).

C. Kinematical limits
We now turn our attention to the transition form factors, which we define in the standard way (see e.g. Ref. [95]), focusing on large Q 2 behavior to start with. It is well known that above a certain large scaleQ 2 0 >Λ 2 QCD , 7 Thus, in addition to the usal setting of the free RL parameters [74,78], 3 more are introduced to obtain 6 new observables.
the TFFs take the form [87,95,96]: where Q 2 >Q 2 0 and φ q M (x; Q 2 ) is the q-flavour valence quark distribution amplitude of meson M . For notational convenience, and in order to match experimental normalization, the TFFs have been rescaled as . This normalization will be employed from this point onward. In the asymptotic domain, Q 2 → ∞, where the conformal limit (CL) is valid, one arrives at: from which the corresponding limits of the SoS and equally off-shell (EoS) TFFs are obtained: For the pion, c 2 π = e 2 (4 − 1)/9, thus arriving at the well-known limit F ∞ M = 2f π [95,96]. To account for the flavour structure of the η − η systems, Eq. (15) is modified as where c l = 5/9, c s = √ 2/9, c 8 = 1/ √ 3, c 0 = 2/3; and f 8,0 η,η are the decay constants in the octet-singlet basis [97]. Owing to the non-Abelian anomaly, the singlet decay constant, f 0 η,η , and thus f l,s η,η , exhibits scale dependence [98]. Moreover, although the RL gives the correct power laws as perturbative QCD, it fails to produce the correct anomalous dimensions. This is readily solved by a proper evolution of the BS wave function. The way QCD evolution is implemented in our calculations is detailed in Refs. [58][59][60]; this process entails that the correct limits of Eqs. (15)- (16) are numerically reproduced.
In the opposite kinematical limit of Q 2 , the Abelian anomaly dictates the strength of F M (0, 0) for the Goldstone modes. In the chiral limit, this entails: where the index '0' denotes chiral limit value. The nonmasslessness of the π 0 and η mesons produces slight deviations from the above result [5,99]. Supplemented by our value f 0 π = 0.092 GeV, Eq. (18) can be employed to fix s 0 = 1.91 s l in the QPV, Eq. (6). I: RL parameters (left and central panels) are fixed to produce the ground-state masses and decay constants. We restrain ourselves to ω = 0.5 GeV, the midpoint of the domain of insensitivity [78]. The η − η values follow from RL-I parameters plus the anomaly kernel (Eq.(10)) inputs given in the right panel. Experimental (Exp.) PDG values are taken from [5], ( * ) lattice QCD results from [94] and ( † ) phenomenological reference numbers from [60]; here mπ and mK correspond to the average of the neutral and charged mesons. The mass units are in GeV.

RL-I
with α em = e 2 /(4π), the electromagnetic coupling constant. From the computed masses and decay constants in Table I, one can readily infer the corresponding decay widths for the η c and η b mesons [59]: This yields to the values: such that s c = 0.78 and s b = 0.23 in order to hold Eq. (19) true 8 . Current algebra is adapted to obtain the analogous for the η − η case [60]: Thus, from the values of Table I, one gets: predictions which are commensurate with empirical determinations, respectively [5]: 0.516(22) keV, 4.35 (36) keV. The results of (23) fix s s = 0.48 and demand a reduction of s l = 1.91 → 1.21, for the η − η case. This happens for two reasons: since we give more weight to 8 Experimentally, Γ γγ ηc = 5.0(4) keV. Nothing is gained for the TFF if sc is fixed to reproduce that value, [59], and the corresponding contribution to aµ would be contained within our final error estimate. the correct description of masses, our best set of parameters in Table I underestimates the value of f l η , as compared to phenomenology. Secondly, a key difference of η − η , with respect to π 0 , η c , η b , is the presence of the non-Abelian anomaly, which conceivably generates corrections to Eq. (2) at infrared momenta. This issue will be addressed elsewhere. Nonetheless, we can estimate the potential impact of our model inputs by following the criteria explained in Section IV.
A comparison of our results for F M (0, 0) to the corresponding measurements is given in Table II. Basically, the error bars are obtained by varying the strength of the transverse terms in the QPV (s f ), as well as the computed masses. This process does not alter any of the conclusions presented in Refs. [58][59][60], nor the agreement those results have with the empirical data; instead, it allows us to provide an error estimate in a quantity that could be sensitive to small variations of the inputs, such as a µ . Notice that, while the F (0, 0) values of π 0 , η and η c show an accurate match with the empirically inferred results, the η is underestimated and produces a larger error bar. A similar pattern has been observed for the charge radii (r M ), which is essentially the slope at Q 2 = 0. While the π and η c charge radii are obtained with desired accuracy [58,59] (this also occurs with the pion EFF [64,66]), and the corresponding comparison with the experiment is within the 1.5% level [5], the η − η system suffers from a larger uncertainty [60]. This is attributed mostly to the presence of the non-Abelian anomaly and the failure of Eq. (2) to incorporate beyond RL effects. Firstly, it is worth mentioning that although our TFF lacks (dynamical) poles in the timelike region, the space-like behavior at low-Q 2 is compatible with VMD, F M (Q 2 , 0) ∼ (Q 2 + m 2 V ) −1 , as can be read from Figure 1. Second, our vertex Ansatz is further validated by the neat agreement with Ref. [83], in which the connection between the QPV and the pion charge radius is clearly established. While our analysis yields r π = 0.675 (9), the most complete result in [83] gives r π = 0.678 fm. Finally, our π 0 , η c , η b , predictions (obtained with QPV Ansatz and PTIRs) have been proven entirely compatible with those approaches that solve the vertex BSE instead and perform a direct calculation [87,88].
Given the set of SDE-BSE inputs and results, in the next section we discuss the parametrizations for the obtained numerical data.

III. PARAMETRIZATIONS AND CONSTRAINTS
Regardless of the approach one takes to compute TFFs, it is highly convenient to look for certain types of theory-driven parametrizations for those form factors, such that the corresponding integrals of a HLbL µ can be computed with relative ease and with a minimum error following standard methods [42,52].
Besides accurately fitting the numerical data, we look for parametrizations that reproduce the low and high Q 2 constraints to the fullest extent possible. We now discuss VMD, LMD and CAs parametrizations.

A. LMD parametrizations: flaws and strengths
For considerable time, VMD and LMD type of parametrizations have been quite popular. Among other attractive aspects, they allow us to rewrite the a HLbL µ related integrals in such a way that there is no dependence on Q 1 · Q 2 [2,42].
Nevertheless, VMD and LMD fail in reproducing the large Q 2 limits, yielding an incorrect power law in both the SoS (Eq. (15)) or EoS (Eq. (16)) cases. Extensions of LMD that include one or more additional vector mesons (LMD+V or LMD+V+V') can potentially fulfill such requirements. Due to a higher number of parameters, such attempts can provide a more reliable fit in a larger domain of momenta.
In general terms, with x = Q 2 1 and y = Q 2 2 , one can define the LMD+V N −1 parametrizations as follows: However, for any finite y 0 = 0, xF M (x, y 0 ) diverges linearly with c N −1,1 + y 0 c N −2,2 as x → ∞. Consequently, the asymptotic limits cannot be recovered for arbitrary y 0 and the accuracy of the fit is compromised as x increases, which makes this parametrization unsatisfactory.

B. Canterbury Approximants
A more convenient approach to the problem at hand is through the so called Canterbury Approximants [73], which have been recently employed to evaluate a HLbL µ [55,74]. The latter reference follows another SDE treatment to evaluate the pole contributions of π 0 , η, η to a HLbL µ . We explore this alternative to parametrize our numerical solutions and calculate the respective contributions to a HLbL µ . Consider a function f (x, y) symmetric in its variables and with a known series expansion CAs are defined as rational functions constructed out of such polynomials P N (x, y) and R M (x, y): whose coefficients a ij , b ij fulfill the mathematical rules explained in detail in Ref. [55]. We shall employ a certain C 1 2 (x, y) such that the TFFs can be written as: P (x, y) = a 00 + a 10 (x + y) + a 01 (xy) , (29) R(x, y) = 1 + b 10 (x + y) + b 01 (xy) + b 11 (x + y)(xy) The large number of parameters can be reduced straightforwardly: • a 00 = F M (0, 0), low energy constraint.
• a 10 = b 20 (1 + δ BL )F ∞ M , fully asymmetric limit. It has been seen that the pion TFF, xF π (x, 0), marginally exceeds its asymptotic limit in the domain x>20 GeV 2 [58] (subsequently recovering it as x continues to grow). Thus, we have included a parameter δ BL to improve the quality of the fit for our given set of numerical data. This is by no means an implication that the Brodsky-Lepage limit of Eq. (15) is violated; it is rather a numerical artifact to obtain a better interpolation.
The parametrization of Eq. (31) cannot be recast in any way so that the Q 1 · Q 2 dependence in the a HLbL µ integrals disappears, but alternative methods can be implemented [52,55]. Unlike LMD+V N −1 parametrization, it also has a well defined limit when one of the variables is finite (but not zero) and the other tends to infinity. Since the large Q 2 domain is well under control, it enhances its reliability even far beyond the domain that contributes the most to a HLbL µ , and that of the available data set.

C. LMD: ηc and η b
As explained in Ref. [59], the η c,b TFFs lie below their corresponding asymptotic limits even at very large values of momentum transfer 9 . Imposing any asymptotic constraint is useless and potentially harmful for the accuracy of the fit. Moreover, those form factors are harder due to the larger masses of the η c,b mesons. In fact, the curvature of η b TFF is only very pronounced above a couple of hundred GeV 2 [88].
Therefore, a simple LMD-like form can be employed: where we find M V0 := 3.097 GeV = m J/ψ , c 00 = 6.5613 GeV 3 and c 10 = 0.0611 GeV for η c ; M V0 := 9.460 GeV= m Υ , c 00 = 30.8424 GeV 3 and c 10 = 0.0426 for η b . Here the flaws of the large-Q 2 behavior of the LMD parametrizations are irrelevant: those appear far beyond the domain of integration. Notably, the LMD representation of the η b TFF, for x, y<20 GeV 2 , reproduces the numerical result within 1% error. In the next section we present our numerical results for a µ . 9 In a less noticeable way, this also happens for the η . Thus, we found convenient to redefine F ∞ η → x 0 F (x 0 , 0), where x 0 = 70 GeV 2 . We note that the symmetric form factor is not affected by this. Particularly, it is recovered exactly, to our numerical precision, for x = y = 10 GeV 2 .

IV. RESULTS
We display the γγ * → π 0 , η, η TFFs in Fig. 1 and their respective comparisons with available low-energy experimental data [5,[99][100][101][102]; a keen agreement is exhibited. The η c result is plotted in Fig. 2 and the analogous for η b in Fig. 3. It is seen that the η c prediction matches the experimental data [103] and the corresponding for η b is in accordance with the non-relativistic QCD (nrQCD) approach [104]. In the domain of interest, the LMD representations of η c and η b TFFs accurately reproduce the numerical SDE calculations. The corresponding results for all the pseudoscalars till much higher values of the probing photon momentum can be consulted in Refs. [58][59][60]. Notably, the DoS extension we present in this work ensures F π (x, 0) and F π (x, x) converge to their well-defined asymptotic limits, as can be observed in Fig. 4. The CAs faithfully accommodate this numerical behavior. Additionally, the charge radius is reproduced to 1.5 % accuracy.
As we have discussed, all the pieces in our SDE-BSE treatment, in particular the QPV Ansatz, ensure an accurate description of the π 0 , η c and η b mesons; but the case of η − η is not completely satisfactory. This occurs mostly due to the presence of the non-Abelian anomaly which, in principle, could introduce infrared corrections to the impulse approximation [60]. To account for the influence of the model for the QPV and other assumptions, firstly we vary the strength of the transverse terms in the QPV such that: 1) we reproduce (as much as possible) the empirical values of F η,η (0, 0) and 2) a rather large uncertainty is included in the more sensitive domain, around Q 2 ∼ 0.4 GeV 2 . From the computed decay constants, one gets a value of F η (0, 0) which is about 10% smaller than the empirical one, thus producing a broader band for the η meson. In the case of π 0 and η c , without the presence of the non-Abelian anomaly, the goal of this minimal variation is to produce an error band in the vicinity of Q 2 = 0, such that the uncertainty associated with F π,ηc (0, 0) is comparable in size to that reported in PDG [5]. The η b TFF has not been measured yet. To be on the safe side, we include error bars on the charge radius by resorting to the nrQCD result. This produces a 5% error around F η b (0, 0). Any additional but reasonable change in the Bethe-Salpeter kernel parameters has a sufficiently small impact [66,74] and we find it to be contained within those bands. Furthermore, small variations of the meson masses also have negligible effects on the TFFs (in fact, one can take m π = 0 with impunity), but they exhibit moderate to large impact on a µ . Thus we allow ourselves to vary m π ∼ 0.135 − 0.140 GeV, m η ∼ 0.548 − 0.560 GeV and m η ∼ 0.956 − 0.960 GeV. and CLEO [101] collaborations (we have also included L3 data [102] for the η ). Additionally, we display the most recent x = 0 values from PDG [5,99]. The mass units are in GeV.
Putting all together, we obtain: Our SDE prediction of a π 0 −pole µ is compatible with other reported values [42,52,55,56,105,106]. For example, a π 0 −pole µ = (5.81 ± 0.09 ± 0.09 +0.5 −0 ) · 10 −10 according to Ref. [56] (resonance chiral lagrangians), while Hoferichter et al. obtain a π 0 −pole µ = (6.26 +0. 30 −0.25 ) · 10 −10 (dispersive evaluation). Ref. [55] reports a π 0 −pole µ = (6.36 ± 0.26) · 10 −10 (CAs) and a recent SDE evaluation [74] obtains Assuming a two-angle mixing scheme in the flavour basis [107], and a chiral approach [52]: one can write the η − η TFFs in terms of F π (x, y). This simplification yields a η+η −pole µ = (2.86 ± 0.42) · 10 −10 , which is consistent with our result albeit with a larger error. The SDE result from Ref. [74] follows a symmetrypreserving RL approach to compute F l (x, y) and F s (x, y), such that the physical states η − η are obtained from there after assuming a two-angle mixing scheme. It is shown that whether one takes the chiral approach or not, the sum of the η − η contributions to a µ remains the same (although the individual contributions are different). This is due to cancellations that occur because of the structure of the mixing matrix. However, in our SDE-BSE treatment, the mixing between the l and s flavours is produced directly due to the presence of the non-Abelian anomaly kernel in the BSE, Eq. (10). It is the anomaly kernel that produces the mixing; no particular mixing scheme is assumed. For our data sets, limited to the range x, y ≤ 10 GeV 2 , we show the CAs parameters in Table III. Regarding the heavy mesons, although the value of η b is 3 orders of magnitude smaller, the η c is commensurate with the current experimental and theoretical error bars. Moreover, our obtained value is fully compatible with that reported in Ref. [108], a ηc−pole µ = (0.08)·10 −10 . Thus, this contribution might not be omitted when the theoretical calculations reach a higher level of precision. Also, it could serve as an estimate for potential non-perturbative corrections to the charm loop 10 .

V. CONCLUSIONS
It is highly timely to revisit the computation of a HLbL µ on the eve of the FNAL (and hopefully J-PARC) improved measurements. We calculate the dominant piece of this observable (coming mainly from the π 0 pole and secondarily from the η and η poles). For the first time, the sub-leading contributions of η c and η b poles were obtained. As a result of our analysis, we find: Our findings for the light pseudoscalars are compatible with previous determinations and have a comparable uncertainty. While the η b result is negligible, the magnitude of a ηc−pole µ (confirmed in Ref. [108]) is sizable as compared with the contemporary error bars; thus, it could promote more theoretical calculations on the topic.
Earlier and recent SDE works [46,74,109] have shown this continuum approach as a promising tool in understanding the QCD contributions to a µ . This is clearly supported by the consistency with our predictions and those from [74]. Moreover, the present work heavily relies on our earlier studies, [58][59][60], where we compute the pseudoscalar transition form-factors: γγ * → {π 0 , η, η , η c , η b }, all lowest-lying neutral pseudoscalars. Such calculations are based upon a systematic and unified treatment of QCD's SDEs. Several efforts have followed this approach to compute a plethora of hadron properties, with the resulting predictions invariably being in agreement with or confirmed by experimental data and lattice QCD simulations (see Refs. [72,110] for recent reviews).
Our previous research [58][59][60] and the resulting current work not only explain the existing data accurately but are also quantitatively predictive for the ones to be measured in modern facilities. Thus, we believe this work is useful in the collective effort to reduce the error of the SM prediction of a µ so as to maximally benefit from the forthcoming improved measurements and hopefully find indirect evidence for new physics in the future. thanks P. Masjuan and P. Sánchez-Puertas for useful conversations on this topic. This research was also partly supported by Coordinación de la Investigación Científica (CIC) of the University of Michoacan, CONACyT-Mexico and SEP-Cinvestav, through Grant nos. 4.10, CB2014-22117, CB-250628, and 142 (2018), respectively.

Appendix A: Quark propagator and BS amplitudes
It is convenient to express the quark propagator in terms of complex conjugate poles (ccp). Omitting flavor indices, it can be expressed as where z j , m j are obtained from a best fit to the numerical solutions, ensuring Im(m j ) = 0∀j, a feature consistent with confinement [61]. We find that j m = 2 is adequate to provide an accurate interpolation. The BS amplitude of a neutral pseudoscalar is written in terms of four covariants, namely: Each scalar function, F k = F ( k; P ), is split in two parts and can be parametrized in terms of PTIRs in the following way: F(k; P ) = F i (k; P ) + F u (k; P ) , (A3a) with ∆ Λ (s) = Λ 2 ∆ Λ (s), k 2 z = k 2 + zk·P , a − F = 1 − a F . The indices 'i' and 'u' denote the connection with the infrared and ultraviolet behaviors of the BS amplitude and, the spectral density: The interpolating parameters are obtained through fitting to the Chebyshev moments: with n = 0, 2, where U n is an order-n Chebyshev polynomial of the second kind. F 4 (k; P ) is small and has no impact, hence it is omitted in all cases. For similar reasons, F 3 (k; P ) might be omitted for η and η b as well. The forms given in Eqs.(A1)-(A3) have been proven undoubtedly useful; their specific interpolating values are presented in Refs. [58][59][60].