Pseudoscalar transition form factors: $(g-2)$ of the muon, pseudoscalar decays into lepton pairs, and the $\eta-\eta'$ mixing

We present our model-independent and data-driven method to describe pseudoscalar meson transition form factors in the space- and (low-energy) time-like regions. The method is general and conforms a toolkit applicable to any other form factor, of one and two variables, with the potential to include both high- and low-energy QCD constraints altogether. The method makes use of analyticity and unitary properties of form factors, it is simple, systematic and can be improved upon by including new data. In the present discussion, the method is used to show the impact of experimental data for precision calculations in the low-energy sector of the Standard Model. In particular, due to its relevance for New Physics searches, we have considered the hadronic light-by-light scattering contribution to the anomalous magnetic moment of the muon (the pseudoscalar exchange contribution), the pseudoscalar decays into lepton pairs, and the determination of the mixing parameters of the $\eta$ and $\eta'$ system. For all of them we provide the most updated results in a data-driven manner.


Introduction
The pseudoscalar meson transition form factors (TFF) describe the effect of the strong interaction on the γ * γ * − P transition (where P = π 0 , η, η · · · ) and is represented by a function F Pγ * γ * (q 2 1 , q 2 2 ) of the photon virtualities q 2 1 , and q 2 2 . From the experimental point of view, one can study the TFF from both space-and time-like energy regimes. The time-like TFF can be accessed from a single Dalitz decay P → l + l − γ process which contains an off-shell photon with the momentum transfer q 2 1 , a F Pγ * γ (q 2 1 , 0), covering the 4m 2 l < q 2 < m 2 P region. The space-like TFF can be accessed in e + e − colliders by the two-photon-fusion reaction e + e − → e + e − P, see Fig. 1. The common practice is to extract the TFF when one of the outgoing leptons is tagged and the other is not, that is, the single-tag method. The tagged lepton emits a highly off-shell photon with the momentum transfer q 2 1 ≡ −Q 2 and is detected, while the other, untagged, is scattered at a small angle and its momentum transfer q 2 2 is almost zero, i.e., F Pγ * γ (Q 2 ) ≡ F Pγ * γ * (−Q 2 , 0). Experimental information on its doubly virtuality can be obtained through the double-tag method which demands the tagging of both leptons in the final-state. Due to the decrease of the cross section for both virtual photons, and the difficulties of the detection of all the particles entering into the process, data for the doubly virtuality are not yet available for the lowest pseudoscalars.
Theoretically, the limits Q 2 = 0 and Q 2 → ∞ are well known in terms of the axial anomaly in the chiral-limit of QCD [2] and perturbative QCD [3], respectively. The TFF a e-mail: masjuan@kph.uni-mainz.de b e-mail: sanchezp@kph.uni-mainz.de is then calculated as a convolution of a perturbative hardscattering amplitude and a gauge-invariant meson distribution amplitude which incorporates the non-perturbative dynamics of QCD [3]. At that point, a model needs to be used either for the distribution amplitude or for the TFF itself [3]. The discrepancy among different approaches reflects the model-dependency of that procedure. On top, when the desired object should be parameterized for the whole space-like range, and available experimental data are going to be used, the asymptotic scale of QCD should be fixed and a contact with the low-energy description performed. Theoretical efforts have not yet succeeded on performing such an endeavor, even worse for ascribing a theoretical errors to the procedure. A different strategy might be, then, desirable.
Let us remark that form factors are not interesting by themselves as represent the knowledge of QCD in a nutshell, but also for their important role on precision calculations of low-energy Standard Model observables such as the anomalous magnetic moment of the muon, the pseudoscalar decays into lepton pairs, the radius of the proton, etc, where deviations between the theoretical predictions and the experimental measure point towards an indirect search of New Physics phenomena.
The experimental information on the TFFs together with the theoretical knowledge on their kinematic limits yield the opportunity for a nice synergy between experiment and theory in a simple, easy, systematic, and userfriendly method. This synergy is the purpose of our work.
In this respect, our proposal can be summarized taking the following considerations: • We want a method, not a model. Our proposal is actually a TOOLKIT.
• It is a method based on the analyticity and unitary of the TFF, important ingredients when heading towards errors at the 10% level.
• Should be a method as simple as possible, maximally transparent.
• If possible, the method should not use any assumption, only approaches. That means, a method improvable without ad hoc statements.
• We shall provide a systematic method in two different senses: easy to update whenever new experimental data or new theoretical calculations are available; capable of provide a purely theoretical error from the approaches performed.
• Finally, should be predictive and checkable.
The answer to this catalogue of wishes can be found within the Theory of Padé approximants. The connexion with the mathematical problem is given by the welldefined general rational Hermite interpolation problem. This problem corresponds with the situation where a function should be approximated but previous information about it is scarce and spread over certain information on a given set of points together with a set of derivatives. Our goal is to make a contact with this problem from our needs and provide a model-independent and data-driven parameterization of the TFFs useful for the problems we have at the precision low-energy frontier of the Standard Model.

Why Padé Theory?
In front of other more sophisticated models and methods, it may seem curious to appeal now to an old tool such as Padé Theory to perform precise calculations within the low-energy frontier of the Standard Model. In this respect, two considerations must be taken into account. First, Padé Theory is an active field of research in applied mathematics. Few of the results considered here have been the subject of research in the last years only. Second, the corpus of Padé Theory defines precisely the problem we have at hand (as we already announced, the general rational Hermite interpolation problem) and provides with a solution in the form of convergence theorems and tools.
From a theoretical point of view, one should notice that for precise Standard Model calculations at low energies considered here, what is needed is not directly the TFF but an integral over the TFF with a particular weight. (The exception to this point discussed in this talk is the extraction of the η − η mixing parameters from the corresponding TFFs.) Another interesting assertion is the fact that parameterizations such as the Vector Meson Dominance, Lowest Meson Dominance (and extensions), models from holographic QCD [4] are already a certain kind of Padé approximant (the so-called Padé-Type approxiamnts [5] where the poles of the Padé are given in advanced). These parameterizations should take advantage of the theory of Padé approximants if a robust calculation is to be claimed. For example, the uncertainties due to the truncation of the PA sequence which are never accounted for when using these parameterizations, do not need to be small to be neglected [6] 1 . On top, it is proven that Padé Type approximants converge slower than PAs, specially when considering integrals up to infinity (cf. second Ref. in [7]). In this respect, heading towards 10% errors for the pseudoscalar contributions to the HLBL (the typical size of the foreseen experimental uncertainty in the new experiment at Fermi-Lab), all these considerations are a must.
The PAs feature an interesting relation with dispersion relations. The state-of-the-art dispersive formalism (DR) accounting for the HLBL (see Ref. [8] together with the discussion in Ref. [9]) has to face the lack of high-energy constraints coming from QCD (cf. as well [4]). In this respect, PAs can provide a natural tool to extend the DR results and interpolate up to the high-energy region, including the well-known pQCD behavior into the game as it should. Actually, the PA themselves can be used by the DR as a parameterization where to impose the dispersive constraints, obtain the PA coefficients, and then extrapolate far away in the space-like region. From another point of view, the complementarity between PA and DR may come from the interplay between TL and SL regions. In this case, DR after accounting for the TL data (a region not accessible by the PA built from the SL), can provide the LECs to be used for reconstructing the PA in the SL. In this situation, PA would represent an extension of the DR into the SL. In an interactive procedure, SL experimental data can be used, as well, to constrain the PA parameters for later on being used back in the DR and so on.
Last but not least, as a user-friendly and simple tool, PA can be used in the analysis of experimental results as a corpus of fitting functions. In this situation, instead of using a VMD to fit the experimental data, the highest PA would yield the most accurate theoretical result including systematic errors without the need to use a particular model.

Transition form factors from rational approximants
The toolkit we have introduced is a powerful method to deal with functions of one and two variables for which experimental data points are known and theoretical constraints available. The method is, of course, general. We developed it focusing on the need to properly introduce experimental data in a systematic manner in observables that manifest a discrepancy with the Standard Model calculation. For the HLBL and pseudoscalar decays into lepton pairs, the object are meson TFFs. In this section we collect our main results for π 0 , η and η . We proposed in Refs. [11,12] to use a sequence of PA P N M (Q 2 ) = P N (Q 2 ) P M (Q 2 ) [7], to fit the space-like data [1,[15][16][17]. Since PAs are constructed from the Taylor expansion of the Figure 1. e + e − → e + e − P space-like process. q 1 and q 2 represent the photon virtualities and the black blob in the interaction point represents the meson transition form factor F Pγ * γ * (q 2 1 , q 2 2 ). Table 1. π 0 , η, and η slope b P , curvature c P , and asymptotic limit (Q 2 → ∞) from Refs. [11][12][13] , from the fits we can obtain the derivatives of the F Pγ * γ (Q 2 ) at the origin of energies in a simple, systematic and model-independent way [11,12]. The idea behind is that the TFF can be expanded as: where a P 0 is related to the decay of the pseudoscalar into two photons through b P and c P in Eq. (1) are the slope and curvature of the TFF respectively. The low-energy parameters are fundamental quantities for constructing the PAs and any hadronic model to be used to evaluate hadronic contributions to the observables we are discussing.
Since the analytic properties of TFFs are not known in general although the time-like region at low energies show the well-known unitary ππ cut, the kind of PA sequence to be used is not determined in advance. We consider two different sequences and the comparison among them should reassess our results. The first one is a P L 1 (Q 2 ) sequence inspired by the success of the simple vector meson dominance ansatz [10], and the second one is a P N N (Q 2 ) sequence which satisfy the pQCD constrains Q 2 F Pγγ * (Q 2 ) ∼ constant. After combining both sequences' results, slope and curvature from Refs. [11][12][13][14] are shown in Table 1 is also shown. The results are shown in Fig. 2.  Figure 2. π 0 (upper panel), η (middle panel), and η (lower panel) TFFs. Green-dot-dashed lines show our best P L 1 (Q 2 ) fit, and black-solid lines show our best P N N (Q 2 ) fit. Black-dashed lines display the extrapolation of the P N N (Q 2 ) at Q 2 = 0 and Q 2 → ∞. Experimental data are from CELLO (red circles), CLEO (purple triangles), and BABAR (orange squares) Colls. [15]. The π 0 figure contains also data from BELLE (blue diamonds) [16]; and the η figure data from L3 (blue diamonds) [17]. Figures form Refs. [11,12].

Applications of the transition form factors 3.1 Hadronic light-by-light scattering contribution to the muon (g − 2)
The HLBL cannot be directly related to any measurable cross section and requires the knowledge of QCD contributions at all energy scales. Since this is not known yet, one needs to rely on hadronic models to compute it. Such models introduce systematic errors which are difficult to quantify. The large-N c limit of QCD [18] provides a very useful framework to approach this problem and has been used to perform what nowadays are the reference numbers for the HLBL [4]. This, however, has a shortcoming. Cal- culations carried out in the large-N c limit demand an infinite set of resonances. As such sum is not known in practice, one ends up truncating the spectral function in a resonance saturation scheme, the so-called Minimal Hadronic Approximation [19]. The resonance masses used in each calculation are then taken as the physical ones from PDG [20] instead of the corresponding (but unknown) masses in the large-N c limit. Both problems might lead to large systematic errors not included so far [5,11,21]. It was pointed out in Ref. [5] that, in the large-N c framework, the Minimal Hadronic Approximation can be understood from the mathematical theory of PAs to meromorphic functions. Obeying the PA rules, one can compute the desired quantities in a model-independent way, within the large-N c , and even be able to ascribe a systematic error to the approach [22].
Our proposal is, then, to use a sequence of PAs built up from the low-energy expansion obtained in the previous section. The form factor needed in this calculation is the one of double virtuality. For this reason we extend the PA method to the Canterbury approximants (CA) [23]. Now, the Taylor expansion of the doubly-virtual TFF reads: where a 1,1 corresponds to the slope of double virtuality and we have made use of the Bose symmetry to simplify our expression. Having experimental data in some energy region, even without high precision (30% or even 50% statistical error), one can attempt the TFF reconstruction in a systematical way by a sequence of doubly virtual approximants fitted to such data. CAs can accommodate the highenergy constraints from QCD as well [23]: Knowing the Taylor expansion of the F(Q 2 1 , Q 2 2 ), Eq. (2) would be unique: a 0 is determined from the Γ(π 0 → γγ) as before, b p is the slope of the single virtual TFF, and a 1,1 corresponds to the doubly-virtual slope.
Experimental data for F(Q 2 1 , Q 2 2 ) is not available yet and we cannot extract a 1,1 from them. The OPE tells us that lim Q 2 →∞ F(Q 2 , Q 2 ) ∼ Q −2 and implies a 1,1 = 2b 2 p . On the other hand, at low energies, the factorization approach F Pγ * γ * (Q 2 1 , Q 2 2 ) = F Pγ * γ (Q 2 1 , 0)×F Pγγ * (0, Q 2 2 ) works well as it is already seen in the construction of the approximant (2) where the parameter a 1,1 enters as a correction to the 2b 2 p . Taking the factorization into account, then a 1,1 = b 2 p . To be on a conservative side, we consider for the numerical analysis b 2 p ≤ a 1,1 ≤ 2b 2 p . This range should effectively take into account our ignorance on a 1,1 . As soon as experimental data will be available, this range will be immediately shrunk. Our results for the pseudoscalar pole contribution to the HLBL for π 0 , η and η are collected in Table 2. Let us remark that this is the most robust calculation of the η and η contribution to the HLBL up to now. The range (7.33 ÷ 9.94) × 10 −10 is as large as the pion loop contribution which is subleading in 1/N c . There is a large effort to estimate the pion loop at the 10% precision using DRs [8]. We notice, however, that this efforts will be in vain if the range we quote in Table 2, last column, cannot be substantially reduced.
On top, the TFF is considered to be off-shell. To match the large momentum behavior with short-distance constraints from QCD, calculable using the OPE, we consider the relations obtained in Ref. [4]. In practice this will amount to extend the CAs we are using to match the set of high-energy OPE constraints. These results are still preliminary and will be given elsewhere [25]. We anticipate, however, an increase of the results in Table 2, last column, by about 20%, again of the order of the pion loop contribution.
The HLBL is dominate for Q 2 ranging from 0 to 2 GeV 2 , in particular above around 0.5 GeV 2 . Therefore a good description of TFF in such region is very important. Such data are not yet available, but any model should reproduce the available one. That is why in [11,12], in contrast to other approaches, we did not used data directly but the low-energy parameters of the Taylor expansion for the TFF and reconstructed it via PAs. The LECs certainly know about all the data at all energies and as such incorporates all our experimental knowledge at once. This procedure implies a model-independent result together with a well-defined way to ascribe a systematic error. In other words, it is the first procedure that can be considered an approximation, in contrast to the assumptions considered in other approaches.
Let us emphasize the role of experimental data by using the LECs in Table 1 together with the π 0 → γγ to match the free parameters of the LMD+V model introduced in [26]. We would obtain a HLBL,π 0 µ = 7.5 × 10 −10 which contrasts with the original a HLBL,π 0 µ = 6.3 × 10 −10 [26]. The impact of the new experimental data is then clear, inducing a 20% effect. On top, since the LMD+V is not a PA but a well-educated model, it is difficult to ascribe a systematic error due to the large-N c approach. PAs in turn, already account for such corrections Flavour changing and conserving processes π 0 (q) which in Refs. [11,12] were found to be of the order of 5% − 10% less dramatic than the naive 30% from the N c counting, but still required to be taken into account.

pseudoscalar decays into lepton paris
Pseudoscalar decays into lepton pairs provide a unique environment for testing our knowledge of QCD. As such decays are driven by a loop process, encode, at once, low and high energies. For the π 0 decay, the process (neglecting electroweak corrections) proceeds in two steps as shown in Fig. 3. The loop does not diverge due to the presence of the pseudoscalar transition form factor on the π 0 → γ * γ * anomalous vertex [2], the F Pγ * γ * (k 2 , (q − k) 2 ) with k 2 , (q − k) 2 space-like photon virtualities. For η and η intermediate ππ cuts must be taken into account as well [37].
The branching ratio (BR) for this decay may be expressed in terms of the two photon decay width as where A(q 2 ) represents the loop amplitude (see Ref. [23] and references therein for details) which is unknown as far as the normalized TFF (F Pγγ (0, 0) = 1) is unspecified. The role of the doubly virtual TFF is actually rather important as the given integral, similarly to the HLBL case, is UV-divergent. Remarkably, for the π 0 case, it is possible to go further without a single clue on the TFF. Being the lightest hadron, it is not possible to find any additional intermediate hadronic state which may be on-shell, and contribute therefore to the imaginary part. Consequently, its imaginary part is solely given by the intermediate γγ state, the well-known unitary bound discussed by Drell [29], which gives Due to the presence of the photon propagators, the kernel of the integral is peaked at around the electron mass. Then, one can expand that kernel in terms of m e /m P but also m e /Λ, being Λ the cut off of the loop integral, or the hadronic scale driven by the TFF. [34] resummed such power corrections and found them negligible [35]. Then, using a Vector Meson Dominance for the TFF they found BR(π 0 → e + e − ) = (6.23 ± 0.09) × 10 −8 [35], 3σ off the KTeV result.
In [23], we investigated the role of the TFF of doubly virtuality using the CAs. Beyond considering the eventual effects of different New Physics scenarios (which are negligible), we consider for the numerical analysis that 0 ≤ a 1,1 ≤ 2b 2 P , and obtained the new Standard Model value. Extending our calculations to the η and η decays now including as well the µµ mode, the preliminary results are found in Table 3.
The error comes from P → γγ and b π together with the evaluation of the systematic error from our approximation [23], and the two main numbers from the ranging of a 1,1 . To shrink the window here provided, experimental data would then be very welcome. For the π 0 , this final number represents still a deviation of the measured BR between 2.6 − 1.4σ.
Forcing our approximant (2) to reproduce the KTeV result and then used for the π 0 contribution to the HLBL [4], we obtain a HLBL;π 0 µ = 2.9 × 10 −10 , a deviation of 2 − 3 × 10 −10 with respect to the standard result. Taking into account that the global theoretical SM error for the muon (g − 2) is 6 × 10 −10 [4], the role of the π 0 → e + e − is certainly remarkable, and never been considered so far. Similar effect is also found for the η → µ + µ − decay, indicating that the current precision of the SM error on the muon (g − 2) may be underestimated.

η − η mixing parameters
The physical η and η mesons are an admixture of the S U(3) Lagrangian eigenstates [27]. Deriving the parameters governing the mixing is a challenging task. Usually, these are determined through the use of η( ) → γγ decays as well as vector radiative decays into η( ) together with Γ(J/Ψ → η γ)/Γ(J/Ψ → ηγ) [27]. However, since pQCD predicts that the asymptotic limit of the TFF for the η( ) is essentially given in terms of these mixing parameters, we use our TFF parametrization to estimate the asymptotic limit and further constrain the mixing parameters with compatible results compared to standard (but more sophisticated) determinations [12].
In this section we reanalyze [14] the η − η mixing as we did in Ref. [12,13]. In these works, we took advantage of the flavor basis, where the η and η pseudoscalar decay constants defined in terms of the axial current, J a 5µ ≡ qγ µ γ 5 λ a 2 q, as 0| J a 5µ |P ≡ ip µ F a P (where a = q, s refers to light and strange quarks, resp.) can be expressed as  [32,33] η → µ + µ − (1.37 ÷ 1.49)(33) × 10 −7 - Table 3. Preliminary results for the range a 1,1 ∈ (b 2 P ÷ 2b 2 P ), i.e. (OPE ÷ Fact). The errors refers to the statistical error for BR(P → γγ), the error from b P and the systematic, combined in quadrature.
theory NLO predictions yield [38,39] where F π is the pion decay constant, and Λ 1 is an OZIrule-violating parameter expected to be small. Assuming Λ 1 to be neglected, Eq. (4) implies φ q = φ s ≡ φ, an approximation which has been successful in phenomenological applications. Large-N c ChPT predicts also that Here, phenomenological studies [12-14, 38, 39] do not support Λ 1 = 0 since they find F q > F π . Therefore, to be consistent we will consider the most general case φ q φ s and work in the so-called octet-singlet basis, where the decay constants are defined as In this basis, Large-N c χPT at NLO predicts [38,39] where F K is the kaon decay constant. To derive (6) and (7) we made use of the relation between the low-energy constant L 5 appearing in the Large-N c χPT Lagrangian at NLO [40] and the ratio F K /F π . F 0 is renormalization group dependent. This is connected to the J 0 5µ anomalous dimension implying [27,41] where N F is the number of active flavors. Solving this equation at order α s , In the singlet-octet basis, the different limiting behaviors of the TFF, F Pγγ ≡ F Pγ * γ (0), and P ∞ ≡ lim Q 2 →∞ Q 2 F Pγ * γ (Q 2 ) take the simple form are charge factors and δ = −0.17 [13] accounts for the F 0 running up to Q 2 → ∞ [41]. In addition, we have included the OZI-violating parameter Λ 3 , which has been neglected in our previous studies since it enters at the same level as Λ 1 .
The set of Eqs. (8-11) form a system of 4 equations with 5 unknowns (F (8,0) η ( ) , Λ 3 ). Then it could seem that, at least taking Λ 3 = 0, we may solve the system. However, as explained in [13], such a system is not linear independent as there is the relation which is free of mixing parameters. Indeed, Eq. (12) determines Λ 3 once its left-hand-side is (experimentally) known but, we still have to face the fact that our system is linear dependent. In order to overcome this problem, we notice that at NLO in large-N c χPT , Eqs. (6,7), provides a clean prediction both for F 8 and (θ 8 − θ 0 ) in terms of the well-known value for F K /F π [20]. Taking either F 8 or (θ 8 − θ 0 ) constraint, would add an additional equation to the previous system, which would provide a unique solution. Taking both, would lead to an overdetermined system, which in general has no solution. For this reason, we adopt a democratic procedure in which we perform a fit including both, F 8 and (θ 8 − θ 0 ) constraints, together with Eqs. (8)(9)(10)(11). In addition, we include the theoretical uncertainties from Large-N c χPT predictions, Eqs. (6,7) by noticing that F K /F π typically receives 5% corrections from the NNLO [40]. Consequently, we add this error in quadrature on top of the one from [20] for our fitting procedure.
The solutions are collected in Tab. 4. The value for χ 2 /DOF is excellent which indicates a good agreement with large-N c χPT but predicts non-negligible NNLO corrections accounted here as a 5%. Without this 5%, the χ 2 /DOF would grow up to 1.5. In addition we can use Flavour changing and conserving processes (7) (7) 1.0 Table 4. Predictions for the mixing parameters in the singlet-octet basis. F π = 92.2MeV and θ 8,0 in degrees.
In Fig. 4 we collect our results (orange crosses) [14] and compare them to different predictions in the literature [27,38,39,42] (red dots), together with our previous results [13] in blue-empty squares. We see that the main difference with respect to our previous work [13], where we did not use the η TFF asymptotic value and assumed φ q = φ s , appears in F 0 . This is to be expected as the inclusion of Λ 1 and Λ 3 affects the singlet part exclusively. In addition, we have reduced our errors thanks to the constraints from Large-N c χPT with respect to our previous work. Our prediction for Λ 3 may be compared with the one in [42], Λ 3 = −0.03 (2). Both of them point towards a small value for this parameter, though they differ in sign. We find that Λ 3 actually plays an important role not only in fulfilling the degeneracy condition, Eq. (12), but in the η(η ) → γγ decays as well. In addition, the Λ 1 term is rather important and affects specially the η results, where deviations of order 10% appear if this is omitted. Finally, we stress that the use of the RG equation for F 0 is fundamental, increases η ∞ and diminishes η ∞ , bringing in agreement experiment and theory. We encourage then the future analysis to account for it.

Conclusions and Outlook
Pseudoscalar meson transition form factors are a good laboratory to study the properties of pseudoscalar mesons. Their interest goes, however, much beyond the mesons themselves as they play a key role on precision calculations of Standard Model observables at low energies where hadronic contributions are the cornerstone of the error evaluation. We propose the method of Padé approximants as a toolkit to analyze them. The method is easy, systematic, user-friendly, and can be improved upon by including new data. Provides, as well, information about the underlying structure of the TFF and can be used to extrapolate experimental information to extract the low-energy parameters of the form factor together with their asymptotic limits. The former makes contact with the experimental measurement of the pseudoscalar decay into two photons as well as calculations of the slope and curvature, and the later provides insights into the perturbative QCD regime.
Once LEPs and asymptotic values are known, the TFFs allow us to study the η − η mixing in a framework consistent with chiral perturbation theory within the large-N c limit at NLO including OZI-rule-violating parameters.
The most relevant feature of the method here described is their excellent performance as an interpolation tool. As such, it is the most compelling method to provide an accurate parameterization for the TFF in the whole spacelike region. This is of primordial importance for the correct assessment, in a data-driven approach, of the pseudoscalar contribution to the HLBL [6,24]. But this does not stop here. Since the approximants can as well penetrate into the time-like region below the first resonance, precise experimental data can be easily incorporated. Beyond the knowledge the approximants bring with respect to the features of the TFFs at time like, with this extension pseudoscalar decays into lepton pairs can be studied. In this regards, our method provides the most accurate datadriven and model-independent result consistent not only with the well-known QCD features at high and low energies, but as well with a well with a mathematical theory. This excellent performance is here exemplified by calculating the pseudoscalar pole contributions to the HLBL and pseudoscalar decays into leptons pairs (for π 0 , η and η ) resulting in the most updated and data-driven results for these quantities up to now.   Figure 4. η − η mixing parameters in the octet-singlet basis from L [27], FKS [39], BDO [42], EF [38], EMS(14) [12], EMS(15) [13].