A data-driven approach to $\pi^{0}, \eta$ and $\eta^{\prime}$ single and double Dalitz decays

The dilepton invariant mass spectra and integrated branching ratios of the single and double Dalitz decays $\mathcal{P}\to\ell^{+}\ell^{-}\gamma$ and $\mathcal{P}\to\ell^{+}\ell^{-}\ell^{+}\ell^{-}$ ($\mathcal{P}=\pi^{0}, \eta, \eta^{\prime}$; $\ell=e$ or $\mu$) are predicted by means of a data-driven approach based on the use of rational approximants applied to $\pi^{0}, \eta$ and $\eta^{\prime}$ transition form factor experimental data in the space-like region.


Introduction
Anomalous decays of the neutral pseudoscalar mesons P (P = π 0 , η, η ) are driven through the chiral anomaly of QCD. Of historical importance is the process P → γγ, which apart from being the main decay channel of the π 0 and the η, its experimental discovery confirmed, for the first time, the existence of anomalies. In this case, the two final state photons are real and the transition form factor (TFF) encoding the effects of the strong interactions of the decaying meson is predicted to be a mere constant, the value of the axial anomaly in the chiral and large-N c limits of QCD, F π 0 γγ (0) = 1/(4π 2 F π ) in the case of neutral pion, where F π 92 MeV is the pion decay constant. On the contrary, if one of the two photons is virtual, the corresponding TFF is no longer a constant but a function of the transferred momentum to the virtual photon F Pγγ * (q 2 ), whereas when both photons are virtual the TFF depends on both photon virtualities and it is represented by F Pγ * γ * (q 2 1 , q 2 2 ). A single Dalitz decay occurs through the single-virtual TFF after the conversion of the virtual photon into a lepton pair, while double Dalitz decays proceed with the TFF of double virtuality involving two dilepton pairs in the final state. Dalitz decays are attractive processes to improve our knowledge of the not yet exactly known TFFs of the Pγ ( * ) γ * vertices. This is the main motivation of this work together with predicting the invariant mass spectra and the branching ratios (BR) of the decays P → + − γ and P → + − + − with P = π 0 , η, η and = e or µ.
On the theory side, the effort is focused on encoding the QCD dynamical effects in the anomalous Pγ ( * ) γ * vertices through the corresponding TFF functions. The exact momentum dependence of these TFFs over the whole energy region is not known, we only possess theoretical predictions from chiral perturbation theory (ChPT) and perturbative QCD (pQCD), thus constraining the low-and space-like large-momentum transfer regions, respectively. The TFF at zero-momentum transfer can be inferred either from the measured two-photon partial width, 1) or the prediction from the axial anomaly in the chiral and large-N c limits of QCD, as mentioned before, while the asymptotic behavior of the TFF at Q 2 ≡ −q 2 → ∞ should exhibit the right falloff as 1/Q 2 [12] 1 . Furthermore, the operator product expansion (OPE) predicts the behaviour of the double-virtual TFF in the limit Q 2 1 = Q 2 2 ≡ Q 2 → ∞ to be the same as for the single one, that is, 1/Q 2 [19] 2 . For the intermediate-momentum transfer region, the most common parameterization of the TFF, widely used by experimental analyses, is provided by the Vector Meson Dominance model (VMD). The dispersive representation of the TFF in terms of q 2 , where q 2 is the photon virtuality in the time-like momentum region, can be written as where s 0 is the threshold for the physical intermediate states imposed by unitarity and ρ(s) = ImF Pγγ * (s)/π is the associated spectral function. To approximate this intermediateenergy part of the spectral function one usually employs one or more single-particle states.
As an illustration, the contribution to the spectral function of a narrow-width resonance of mass M eff reduces to ρ(s) ∝ δ(s − M 2 eff ), which yields where F Pγγ (0) serves as a normalization constant and Λ(= M eff ) is a real parameter which fixes the position of the resonance pole on the real axis. However, this simple and successful single-pole approximation given in eq. (1.3) breaks down for q 2 = Λ 2 . One may cure this limitation by taking into account resonant finite-width effects as proposed by Landsberg in ref. [20] when considering the transitions P → + − γ in a VMD framework. According to this model, these transitions occur through the exchange of the lowest-lying ρ, ω and φ vector resonances and their contributions to the TFF are written as where F Pγγ * (q 2 ) = F Pγγ * (q 2 )/F Pγγ (0) is defined as the normalized TFF, g V Pγ and g V γ are the V Pγ and V γ couplings, respectively, M V the vector masses, and Γ V (q 2 ) the energydependent widths. Despite the notorious success of VMD in describing lots of phenomena at low and intermediate q 2 , particularly useful for the decays we consider in this work, this model can be seen as a first step in a systematic approximation. Padé approximants are used to go beyond VMD in a simple and model-independent manner also incorporating information from higher energies, allowing an improved determination of the low-energy constants relative to other methods [21]. For this reason, we make use in our study of the works in refs. [17,22], where all current measurements of the space-like TFFs γ * γ → P [23][24][25][26][27][28], produced in the reactions e + e − → e + e − P, have been accommodated in nice agreement with experimental 1 Perturbative QCD predicts lim Q 2 →∞ Q 2 F π 0 γγ * (Q 2 ) = 2Fπ. Alternative values to this result exist, see for instance refs. [13,14], though they seem to be disfavored, as pointed out in refs. [15,16]. For the η and η , see the asymptotic values obtained in refs. [17,18]. 2 The OPE predicts for the case of the pion lim Q 2 →∞ Q 2 F π 0 γ * γ * (Q 2 , Q 2 ) = 2Fπ/3. data using these rational approximants. We benefit from these parameterizations valid in the space-like region to predict the transitions Pγ ( * ) γ * used in the time-like region for the Dalitz decays we are interested in, with the primary aim of achieving accurate results for these decays. Different parameterizations existing in the literature are based on resonance chiral theory [29,30] and dispersive techniques [31,32], among others [33][34][35][36][37][38][39][40][41][42]. The paper is organised as follows. In section 2, we introduce our description for the π 0 , η and η transition form factors using the mathematical method of Padé approximants. Sections 3 and 4 are devoted to the analysis of single and double Dalitz decays, respectively, and predictions for the several invariant mass spectra and branching ratios are given. Finally, in section 5, we present our conclusions.

Transition Form Factors
The usefulness of Padé approximants (PAs) as fitting functions for different form factors have been extensively illustrated [17,18,21,22,44,45]. The reader is referred to these references for details on the method though here we cover some important aspects for consistency. The PAs to a given function are ratios of two polynomials (with degree L and M , respectively), 3 constructed such that the Taylor expansion around the origin exactly coincides with that of f (q 2 ) up to the highest possible order, i.e., f (q 2 )−P L M (q 2 ) = O(q 2 ) L+M +1 . We would like to point out that the previous VMD ansatz for the form factor (1.3) can be viewed as the first element in a sequence of PAs which can be constructed in a systematic way. By considering higher-order terms in the sequence, one may be able to describe the experimental space-like data with an increasing level of accuracy. The important difference with respect to the traditional VMD approach is that, as a Padé sequence, the approximation is well-defined and can be systematically improved upon. Although polynomial fitting is more common, in general, rational approximants are able to approximate the original function in a much broader range in momentum than a polynomial [43].
Yet the success of PAs as fitting functions for space-like TFFs, some important remarks are in order. First, there is no a priori mathematical proof ensuring the convergence of a Padé sequence to the unknown TFF function, though a pattern of convergence may be inferred from the data analysis a posteriori 4 . For instance, the excellent performance of PAs in ref. [18] (see figures 2 and 7 there) seems to indicate that the convergence of the η TFF normalization and low-energy constants is assured (see also ref. [45] for the η case). Second, unlike the space-like TFF data analyses [17,22], one should not expect to reproduce the time-like TFF data since a Padé approximant contains only isolated poles and cannot reproduce a time-like cut 5 . However, if this right-hand cut is approximated by one or more single-particle states in the form of one or several narrow-width resonances, as stated before, then the Padé method may be still used up to the first resonance pole, indeed, up to neighbourhoods of the pole. The size of the region which is affected by the presence of the pole, a disk of radius ε, is not known but, as we will see later, may be deduced, thus fixing the range of application of the PAs for the time-like data. Third and last, the poles found in the PAs fitting the TFFs can not be directly associated to physical resonance poles in the second Riemann sheet of the complex plane. These, in turn, may be obtained following the prescriptions of refs. [46,47], which is beyond the scope of the present work.
We would like to emphasize that the use of PAs as fitting functions for some set of experimental data can be viewed as an effective mathematical method which intrinsically contains relevant physical information of the function represented by the data set. In this work, we benefit from the findings of refs. [17,22], where the π 0 , η and η TFFs were fixed in the space-like region from the analysis of the intermediate process γγ * → P by several experimental collaborations, to predict the time-like region of the same TFFs needed for the description of the reaction P → γγ * and therefore for the single and double Dalitz decays studied here. The extrapolated version of the TFFs, from the space-like region to the time-like one, used in this analysis are discussed case by case in the following.
2.1 π 0 → γγ * Given the small phase-space available in this transition, 4m 2 ≤ q 2 ≤ M 2 π , the π 0 TFF can be expressed in terms of its Taylor expansion, where F π 0 γγ (0) is fixed from eq. (1.1) and the values of the low-energy constants (LECs), slope and curvature, b π and c π , respectively, are borrowed from eqs. (12,13) in ref. [22] 6 , b π = −3.24(12) stat (19) sys × 10 −2 , c π = 1.06(9) stat (25) sys × 10 −3 , where the statistical error is the result of a weighted average of several fits using different types of PAs to the same joint set of π 0 TFF space-like data and the systematic error is attributed to the model dependence of the PA's method. In this way, the values obtained for the LECs can be considered as model independent. It is worth mentioning that the systematic errors ascribed to the LECs are quite conservative, in the sense that they are obtained from a comparison of the constants predicted by several well-established phenomenological models for the TFF and their counterparts extracted using various types of PAs from fits to pseudo-data sets generated by the different models. For each LEC, the systematic error is chosen to be the largest difference among these comparisons, making the whole approach reliable and with a stamp of model independence [22]. The Taylor expansion in eq. (2.2) can be safely used for the description of the π 0 TFF in the time-like region within the range of available phase-space since the first pole seems to appear, for all types of PAs considered, inside the region of ρ-dominance [22], thus well beyond the phase-space end point. Finally, we would like to remark that this expansion of the single-virtual π 0 TFF will be used for predicting both π 0 → e + e − γ and π 0 → e + e − e + e − decays, the latter by means of a factorisation of the double-virtual π 0 TFF in terms of a product of the single-virtual one (see subsection 2.4 for details).

η → γγ *
In order to describe the time-like region of the η TFF from the space-like data analysis in ref. [17], we will employ two PAs, the P 5 1 (q 2 ) and the P 2 2 (q 2 ). These are the highestorder PAs one can achieve when confronted with the joint sets of space-like experimental data. The sequence P L 1 (q 2 ) is used when the TFF is believed to be dominated by a single resonance, while the P N N (q 2 ) one is appropriate for the case the TFF fulfils the asymptotic behaviour 7 . A Taylor expansion equivalent to eq. (2.2), with b η = 0.60(6) stat (3) sys and c η = 0.37(10) stat (7) sys for the slope and curvature parameters, respectively, is better not to be used in this case because of the larger phase-space available, 4m 2 ≤ q 2 ≤ M 2 η . From the analysis in ref. [17], we also obtained that the fitted poles for the P L 1 (q 2 ) sequence are seen in the range (0.71, 0.77) GeV, beyond the phase-space end point, thus making again our approach applicable and the predictions reliable.
Our predictions for the modulus square of the normalized time-like η TFF, F ηγγ * (q 2 ), as a function of the invariant dilepton mass, √ s ≡ m , are shown in figure 1, together with the experimental data points from the A2 Collaboration on the decay η → e + e − γ [4] (black circles) and the NA60 experiment on η → µ + µ − γ [7] 8 (green squares). The predictions from the P 5 1 (q 2 ) (red solid line) and P 2 2 (q 2 ) (black long-dashed line) are almost identical and in nice agreement with the experimental data, whereas the Taylor expansion (blue dot-dashed line) is not so precise in the upper part of the spectrum. For this reason, we will use in our analysis both PAs indistinctly. The one-sigma error bands associated to P 5 1 (q 2 ) and P 2 2 (q 2 ) PAs are displayed in light-red and light-gray, respectively. These error bands are built from the uncertainty in the coefficients of the PAs 9 and the normalization factor extracted from the two-photon decay width.

η → γγ *
The description of the whole time-like η TFF by means of PAs is cumbersome. The available phase space, 4m 2 ≤ q 2 ≤ M 2 η , includes now an energy region where poles associated to these PAs can emerge. The analysis of the η TFF space-like data performed in ref. [17] revealed the appearance of a pole in the range (0.83, 0.86) GeV for the cases of a P L 1 (q 2 ) sequence. Consequently, we cannot employ the method of PAs for describing the timelike TFF in the entire phase-space region and a complementary approach must be used. 7 In ref. [17], the fit to space-like data is done for Q 2 |F (Q 2 )| and not for the TFF itself. As a consequence, PAs satisfying the correct asymptotic limit, that is, lim Q 2 →∞ Q 2 F (Q 2 ) = const., are represented by the sequence P N N (q 2 ). 8 We thank S. Damjanovic from the NA60 experiment for providing us with the time-like TFF data points obtained from η → µ + µ − γ. 9 The coefficients of the PAs along with their errors and the correlation matrix can be obtained from the authors upon request. , as a function of the invariant dilepton mass, √ s ≡ m . The predictions coming from the P 5 1 (q 2 ) (red solid line) and P 2 2 (q 2 ) (black long-dashed line) PAs, and the Taylor expansion (blue dot-dashed line) are compared to the experimental data from η → e + e − γ [4] (black circles) and η → µ + µ − γ [7] (green squares). The one-sigma error bands associated to P 5 1 (q 2 ) (light-red) and P 2 2 (q 2 ) (light-gray) PAs, and the QED prediction (gray short-dashed line) are also displayed.
Then, we propose to match the description based on PAs to that given by eq. (1.4) at a certain energy point 10 . Given the mass and the width of the ρ meson, the first of the resonances included in the VMD description, the region of influence due to its presence may be defined using the half-width rule as M ρ ± Γ ρ /2 [48], thus deducing the value of the radius ε mentioned earlier. The particular energy point located at √ s 0.70 GeV, the lowest value of the former region for M ρ 775 MeV and Γ ρ 150 MeV, fixes the optimal matching point 11 . Fixed this value, a representation valid in the whole phase-space domain is that given by the PA below the matching point an eq. (1.4) above it. In order to match both descriptions of the form factor at the matching point we have to rescale the VMD result accordingly. In this manner, we keep track of the resonant behaviour in the upper part of the spectrum where PAs cannot be applied, while the low-energy region is predicted in a more systematic way as compared to VMD by PAs established uniquely from space-like 10 To proceed with the matching, we have considered an energy-dependent width for the ρ resonance, , with σ(q 2 ) = 1 − 4M 2 π /q 2 , and a constant width for the ω and φ narrow resonances. Input values for the masses and widths as well as for the rest of the couplings entering eq. (1.4) are taken from ref. [1]. 11 The region of influence attributed to the ω and φ poles is negligible since these are narrow resonances and are placed far from the matching point. are used. The one-sigma error bands associated to P 6 1 (q 2 ) (light-red) and P 1 1 (q 2 ) (light-gray) PAs, and the QED prediction (gray short-dashed line) are also displayed.
data. This will allow us to integrate the whole spectrum and predict the branching ratio of the several η Dalitz decays considered here.
Our predictions for the time-like η TFF together with the experimental data points from the BESIII Collaboration on the decay η → e + e − γ [8] (black circles) are displayed in figure 2. The results from the P 6 1 (q 2 ) (red solid line) and P 1 1 (q 2 ) (black long-dashed line) are shown up to the matching point. The corresponding error bands are in light-red and light-gray, respectively. From the matching point on our predictions are replaced by a rescaled VMD representation based on the three lowest-lying vector resonances. Our PAsbased predictions are again in fine agreement with experiment. A Taylor expansion with b η = 1.30(15) stat (7) sys and c η = 1.72(47) stat (34) sys [17] is also included for comparison. It is worth mentioning that an extrapolation of the P 6 1 (q 2 ) PA beyond the matching point nicely passes through the last experimental point. This can be understood in the following terms. The VMD description includes three resonances whose poles are located at different places, while the P L 1 (q 2 ) PAs include only one. Making use of the single pole approximation in eq. (1.3) and the Taylor expansion in eq. (2.2) for the case of the η , the slope parameter is identified as b η = m 2 η /Λ 2 . Using the b η value deduced from eq. (1.4) one gets Λ = M eff = 0.822(58) GeV, where the error is due to the half-width rule and can be utilized as a measure of the region of influence of the pole. The former value is very similar to the one obtained from the pole position of the P 6 1 (q 2 ) PA, located at √ s = 0.833 GeV. Therefore, the region of influence of this pole can be estimated to be in the interval (0.77, 0.89) GeV. It is for this reason that the last experimental point would be also in agreement with the P 6 1 (q 2 ) prediction [45]. This is not so for the P 1 1 (q 2 ) PA, thus showing that increasing the Padé order allows for a better description of the data. In any case, for the numerical analysis of the different decays involving the η we also keep both PAs for the sake of comparison.
, depends on both photon virtualities, q 1 and q 2 . Due to Bose symmetry, it must satisfy F Pγ * γ * (q 2 1 , q 2 2 ) = F Pγ * γ * (q 2 2 , q 2 1 ). Its normalization is obviously the same as the single virtual TFF, F Pγ * γ * (0, 0) = F Pγγ * (0), and can be extracted either from the two-photon partial width by means of eq. (1.1) or from the axial anomaly. It must also satisfy that when one of the photons is put on-shell the double-virtual TFF becomes the single-virtual one, i.e. lim q 2 In addition, the double-virtual TFF can fulfil the following asymptotic space-like constraints, Due to the lack of experimental information in the case of double-virtual TFFs, our initial ansatz will be to use the standard factorisation approach, which in terms of normalized form factors reads -52]. This double-virtual TFF description may or may not satisfy the high-energy constraints above. For instance, the PA P 0 1 (q 2 ) = a 0 /(1 − a 1 q 2 ), corresponding to the single pole approximation in eq. (1.3) motivated by VMD, would induce a 1/q 4 term in the double-virtual TFF, which violates the last of the asymptotic constraints mentioned before (OPE prediction) [19,[53][54][55]. For this reason, we also use for our study the lowest order bivariate approximant which consists in a generalization of the univariate PAs named Chisholm approximants (CAs) [43]. The analysis of the π 0 → e + e − decay is a recent example that illustrates the application of these CAs [56] (see also P. Masjuan's contribution in ref. [57]). In eq. (2.4), a 0,0 is identified as the normalization F Pγ * γ * (0, 0) and then fixed from eq. (1.1), b 1,0 is the slope of the single-virtual TFF obtained in refs. [17,22], that is, b π from eq. (2.3) for the pion and b η ( ) from eq. (5) in ref. [17] for the η and η , respectively, and b 1,1 would correspond to the double-virtual slope which may be extracted in the future as soon as experimental data for the double-virtual TFFs become available. For the numerical analysis, we consider, as a conservative estimate, to vary b 1,1 from the value respecting the OPE prediction, b 1,1 = 0, to b 1,1 = 2b 2 1,0 , far from the factorisation result b 1,1 = b 2 1,0 . In this manner, we can test the sensitivity of our predictions to the double-virtual slope. We also encourage experimental groups to perform double-virtual TFF measurements in order to fix this parameter. In this work, we employ both descriptions indistinctly, the factorisation ansatz and the bivariate approximant in eq. (2.4). See also ref. [58] for a recent approach to the double-virtual TFF of the η meson based on the standard factorisation approach.

Single Dalitz Decays
Single Dalitz decays involve the single virtual TFF as described in section 2. The amplitude of the decays reads while the corresponding differential decay rate is given by where the TFF appears in its helpful normalized version in order to avoid misunderstandings due to different conventions on the definition of F Pγγ * (0) existing in the literature. Notice that the experimental measurement of the partial width to two photons appears as a normalization in any case. For our numerical calculations we employ the PrimEx Collaboration result Γ π 0 →γγ = 7.82 (14) 3.1 π 0 → e + e − γ The decay π 0 → e + e − γ was suggested for the first time by Dalitz in 1951 [59]. The first BR prediction arose from a pure QED radiative correction calculation neglecting the momentum dependence of the TFF [60] (radiative corrections have recently been revisited in ref. [61]). By looking at figure 3 where we compare, as a function of the invariant mass of the dielectron pair, our description of the decay rate distribution (green solid curve) with the QED result (gray dashed curve), this seems to be a reasonable approximation since the main contribution to the decay rate comes from the very low-energy part of the spectrum where the effect of the TFF is almost negligible. In fact, there is an almost perfect overlap between the two curves and only really tiny differences appear on the second half of the spectrum. Numerical results for the BR are presented in table 1. Our prediction is in very good agreement with the experimental measurements. The main source of the error we have quoted arise from the uncertainty associated to the low energy constants, eq. (2.3), which include also the error associated to the measured decay width to two photons [64].
Our results are also in agreement with other theoretical predictions existing in literature, refs. [20,32,34,41], as well as with the QED estimates of refs. [32,62,63]. The dimuon mode in the final state is not kinematically allowed in this case.

η → + − γ ( = e, µ)
Qualitatively, this is the same process as the π 0 Dalitz decay, with the novelty that a dimuon pair in the final state is also allowed since the larger mass of the η increases the upper kinematic limit. A priori, these decays are expected to be more challenging for testing the momentum dependence of the TFF because the energy released in the process is now larger, expecting higher deviations from the QED estimates. This is precisely what our predictions reflect in figure 4 where the decay rate distribution of η → e + e − γ (blue solid curve) and η → µ + µ − γ (black solid curve) are compared to the QED estimates (gray dotted and dashed curves, respectively). As a matter of example we have employed the P L 1 (q 2 ) Padé type in the figure, here and hereafter. Diagonal PAs, P N N (q 2 ), would produce very similar description in accordance with the transition form factors descriptions given in sections 2.2 and 2.3. One interesting feature of both decays modes is that, in absolute terms, the effect of the TFF is much more sizable in the muonic case than in the electronic one. The reason is because the shape of the distribution of the latter shows, as occurred in π 0 → e + e − γ, a strong peak in the low-momentum transferred region of the spectrum, where the effect of the TFF is small, which provides the main contribution to the BR. Noticeably, the high energy part of the spectrum of both modes is overlapped as it may be since the only difference between them is the dilepton threshold production. Numerically, we see from table 2 that the BR involving muons in the final state has increased by 50% with respect to the QED prediction while the effect is much less considerable when dealing with electrons (∼ 3.5%) where, predictions with and without considering TFF effects are compatible within errors at the current level of accuracy. The source of the associated errors arise from the error bands shown in figure 1. In all, our predictions are in good agreement  Source with present experimental measurements. Comparing with other authors results, we agree with: the QED estimates of ref. [63], the predictions of ref. [41] and the values of ref. [33], while tiny differences with ref. [34] are noticed.

η → + − γ ( = e, µ)
The large mass of the η increases the upper kinematical limit by ∼ 410 MeV with respect to the case of the η. The TFF description is given in section 2.3, where the effect of the intermediate vector resonances ρ, ω and φ is included. As shown in figure 5, the distribution of the decay η → e + e − γ (blue solid curve) evidences again a marked peak at low-energies which, despite the contribution coming from the resonance region, dominates the decay as occurred in π 0 (η) → e + e − γ. On the contrary, the effect of the TFF on the decay η → µ + µ − γ (black solid curve) is larger than in η → µ + µ − γ, increasing the BR by a factor of about 2. This is so because both phase space considerations and the effect of passing through a q 2 region where resonances may be produced on-shell. Interestingly, the contribution due to the ρ resonance bends the distribution while the inclusion of the ω Source BR(η → e + e − γ) · 10 4 BR(η → µ + µ − γ) · 10 4 this work [P 6   resonance accounts for the sharp peak at around 0.8 GeV. Numerical results are presented in table 3, where the source of the error comes from the error bands associated to the TFF. From the theory side, our predictions are in accordance with those of ref. [41], while they are slighlty below respect to both the recent experimental measurement of η → e + e − γ [8] and the old measurement of η → µ + µ − γ [9], though in agreement within errors in both cases. To sum up, the pattern of both η and η single Dalitz decays is notably similar: the impact of the TFF on the muonic channel is much larger than in the electronic ones as discussed in section 3.2.

Double Dalitz Decays
Double Dalitz decays involve the TFF of double virtuality as described in section 2. They implicate four particles in the final state which makes the phase space integration much more tedious. In case of having two pairs of non-identical particles, that is η ( ) → e + e − µ + µ − , the required diagram is shown in figure 6 (left diagram) and the amplitude of the decay reads (4.1) The corresponding decay rate distribution can be reduced to 12 where, in this case, S = 2 in agreement with the expression given in ref. [42]. The TFF appears again normalized to unity at the origin. On the contrary, in case of having two pairs of identical particles in the final state, that is P → e + e − e + e − or η ( ) → µ + µ − µ + µ − , one should consider both the direct and exchange diagrams of figure 6 (left and right diagrams). Therefore, the total amplitude of the process reads where the appearance of the minus sign is due to the exchange of two indistinguishable fermions in the final state. Then, squaring the amplitude of eq. (4.3) we arrive at where not only appear the contributions from both the direct and exchange diagrams but an interference term. We notice that the contribution to the partial decay width coming from the first and the second term of eq. (4.4) is obviously the same, that is Γ dir = Γ exch , because the integration variables are nothing more than dummy indices. In this way, the contribution coming from the sum of the direct and exchange diagrams, Γ dir+ex , is obtained through the use of eq. particles in the final state, has been taken into account. Regarding the interference term, its computation is much more cumbersome. We have relegated to appendix A the detailed expression due to its length but it is worth to comment that we have obtained an expression in terms of five invariant masses which has required a Monte Carlo (MC) simulation to be integrated.
The only possible double Dalitz decay of the neutral pion is π 0 → e + e − e + e − , other possibilities are not kinematically allowed. In view of the results from π 0 →→ e + e − γ, one may expect that the overall effect of the TFF will be again small. In figure 7, we show our results for the different contributions to the decay rate distribution as a function of the invariant mass of one dielectron pair of the direct diagram. Concretely, we display the curve corresponding to the direct diagram (green solid line), the curve of the contribution of the exchange diagram expressed in terms of the former dielectron invariant mass of the direct diagram 13 (red dotted line), the interference term (blue dotted line) and finally the total distribution (black dotted line). We want to note that the contribution from both direct and exchange diagrams integrates obviously the same and that the interference is small and destructive. Our BR predictions are shown in table 4 from which we corroborate 13 The curve of the exchange diagram expressed in its own variables would look equal as the green solid 3.46 (19) [10] Table 4. Branching ratio predictions for π 0 → e + e − e + e − compared with experimental measurements.
that the effect of the TFF is small because the main contribution to the BR proceeds from the very low-momentum transferred region where we a peak emerges, as already occurred in π 0 → e + e − γ. Our results are well in accordance with current experimental measurements. The source of the associated error comes from the uncertainty on the low-energy parameters eq. (2.3). Notice that the sensitivity of this decay to the variations of the double virtual slope parameter, b 1,1 , is at the fifth decimal number. In this sense, the high level of accuracy demanded to infer its value is unthinkable at the current precision level. It is also interesting to compare with other authors results. We are in good agreement with the results given in refs. [41,42,50] for the direct and exchange contributions while the result of ref. [34] is about 5% lower than our predictions. Regarding the interference term we have, a perfect agreement with refs. [41,50] and a value about 30% higher than ref. [42] while ref. [34] did not consider this term. Comparing with previous QED estimates we agree with refs. [62,63] for the direct and exchange contributions. For the interference term the former did not consider it and the later gave a result 5 times larger than us.
The double Dalitz decays of the η meson, η → e + e − e + e − , η → µ + µ − µ + µ − , and η → e + e − µ + µ − , are now kinematically allowed. Let us first analyze the latter for simplicity. In this case, the two dilepton pairs are different and consequently there is no interference phenomenon. Hence, the distribution rate is just given by eq. (4.2) and shown in figure 8 in two different manners, one expressed in terms of the dielectron invariant mass and the other in the dimuon variable (blue and red solid lines respectively), where, of course, both curves integrate the same. Our predictions are shown in table 5, where the source of the associated errors comes from the error bands associated to the TFF for the case of the factorisation approach, and from the uncertainty on the single virtual slope for the description employing CAs.
From the experimental side, we respect the current upper limit, while from the theory side, because of the appearance of a dimuon pair in the final state, the effect of the TFF increases the BR about 50% for the same arguments as explained in η → + − γ. This decay, though much more sensitive than π 0 → e + e − e + e − to the double virtual slope, b 1,1 ,

Source
Double virtual TFF BR(η → e + e − µ + µ − ) · 10 6 This work Chisholm approximants  would require accurate measurements as well as demanding a very precise description of the TFF, in order to diminish its associated error, for deducing b 1,1 , far from the present situation. Comparing with other authors, we agree with the predictions of ref. [41], while we have found discrepancies with the value 5.83·10 −7 of ref. [34], with the prediction 2·10 −7 of ref. [33] and with the estimate 7.84 · 10 −7 of ref. [63]. In the later case, the reason seems to be a typographical fault of a factor of 2 missing as pointed out in both refs. [41,50]. In such case, it would reproduce the QED result of table 5 as it should be, because they did not consider the momentum dependence of the TFF. The decays involving two identical dilepton pairs in the final state, η → e + e − e + e − and η → µ + µ − µ + µ − , require to consider eq. (4.4). Their distributions are given in figure 9 (left and right panels respectively) as a function of one dilepton invariant mass of the direct diagram. We explicitly show the contribution from the direct diagram (green solid curve), the curve of the exchange diagram expressed in terms of the former dielectron (dimuon) invariant mass of the direct diagram (red dotted curve), the interference term (blue dotted  Figure 9. Different contributions to the decay distributions of η → e + e − e + e − (left) and η → µ + µ − µ + µ − (right), respectively. Direct diagram (green solid curve), exchange diagram (red dotted curve), interference term (blue dotted curve) and the total distribution (black dotted curve) are displayed with respect to one dilepton invariant mass of the direct diagram.

Source
TFF This work  Table 6. Branching ratio predictions for η → e + e − e + e − and η → µ + µ − µ + µ − confronted to current experimental status. curve) and the total decay rate distribution (black dotted line). The integrated BR results are shown in table 6 where the error comes again from the error bands of the TFF description as given in figure 1, for the factorisation approach, and from the uncertainty associated to the slope, b η , for the description using CAs.
Comparing with present experimental status, our prediction for η → e + e − e + e − is compatible at less than 1σ with the KLOE measurement [11] as well as with the recent measurement value of the WASA@COSY collaboration [6], while our estimate for η → µ + µ − µ + µ − respects the current upper bound of ref. [5]. We have found the same trend as in η → + − γ that is, while the overall effect of the TFF on the electronic mode is small, increasing the BR of η → e + e − e + e − by 6% respect to the QED estimate, the impact on the muonic channel, η → µ + µ − µ + µ − , becomes important increasing the BR by a factor ranging (1.6 − 1.7) respect to the QED calculation. As a consequence of that, the sizable sensitivity of η → µ + µ − µ + µ − to the TFF of double virtuality makes it a good candidate to improve our knowledge on it. Interestingly, a precise experimental measurement of this mode at the per cent level of precision leaves us in position to estimate the value of b 1,1 .
For that purpose, it is also required to diminish the associated uncertainty to the TFF.
Here enters the ability of the Padé method we use for accommodating new experimental data as soon as released by experimental groups. On the contrary, this same exercise for η → e + e − e + e − would demand accurate measurements at the per mil level to unveil this quantity, far from the present situation. Our predictions are in good agreement with the results of ref. [41] for the electronic mode, while a (10 − 15)% over the muonic prediction.
Comparing with ref. [34] (who did not considered the interference term) we are a 10% over for the electronic case while his result for η → µ + µ − µ + µ − is 60% smaller. We are also in accordance with the estimate of ref. [33] for η → e + e − e + e − . Regarding the pure QED calculation of ref. [63], we are in perfect agreement for the electronic channel while tiny differences are found in the muonic decay, probably caused by the updated values of our inputs values.
Regarding the double Dalitz decays of the η , we have the same three possible final states as for the η. However, in this case we have only adopted the factorisation approach ansatz for describing the double virtual TFF of the η . The reason is because the use of Chisholm approximants, which may respect the appropriate asymptotic behavior q −2 , would only apply at low energies, concretely up to the matching point where PAs are applicable, while beyond, we are somehow forced to employ the factorisation approximation, through a VMD description, which induces a q −4 term. So, there is no gain respecting the high-energy behavior in the low-energy region if we violate it at high energies. We compute first η → e + e − µ + µ − again through eq. (4.2). Noticeably, it follows the same trend as η → e + e − µ + µ − , with the difference that this case is sensitive to the resonance region as can be read off from figure 10. Once more, the low-momentum region basically dominates the distribution when working with the electronic variable (blue solid curve) while it is a smooth falling function of the dimuonic momentum with a small bump and a sharp peak accounting for the effect of the ρ and the ω, respectively (red solid curve). Indistinguishable, both curves integrate the same BR. Our predictions are presented in table 7 without, in this case, any experimental reference to compare with. The effect of the TFF increases by a factor of about 2 the BR respect to the QED estimate, which is much notorious than in η → e + e − µ + µ − . The decay spectra for η → e + e − e + e − and η → µ + µ − µ + µ − shown in figure 11 (left and right panels respectively) have been computed by taking eq.    in table 8 which also reflect the tendency that the effect of the TFF is sizable and larger than for the case of the η. In particular, the BR of η → e + e − e + e − and η → µ + µ − µ + µ − have increased by 20% and by a factor of 2, respectively. On the experimental side, we neither have an observation to compare with, while on the theory side we have only found the predictions given in ref. [41] with which we are in good agreement for the cases of having two identical dilepton pairs in the final state, while we are slightly below for η → µ + µ − e + e − .  Figure 11. Different contributions to the dielectron and the dimuon invariant masses distribution for η → e + e − e + e − (left) and η → µ + µ − µ + µ − (right), respectively. Direct diagram (green solid curve), exchange diagram (red dotted curve), interference term (blue dotted curve) and the total distribution (black dotted curve) are displayed with respect to one invariant mass of the direct diagram.

Conclusions
The single and double Dalitz decays P → + − γ and P → + − + − (P = π 0 , η, η ; = e or µ) have been analysed by means of a data-driven model-independent description of the transition P → γ ( * ) γ * . We have benefited from (our) previous findings on the space-like single-virtual TFF γγ * → P obtained through the use of Padé approximants to represent these transitions in the time-like energy region where they are applicable. We have shown that this extrapolation from the space-like to the time-like is supported by current experimental data η and η TFFs obtained from η ( ) → e + e − γ and η → µ + µ − γ decays. This nice behaviour proves that these TFFs are well approximated by meromorphic functions. Regarding the TFF of double virtuality, besides the standard factorisation approach, we have motivated the use of bivariate approximants, which would satisfy the high-energy constraints and whose coefficients may be determined as soon as experimental data become available. From the phenomenological point of view, we have found that the invariant mass distributions involving electrons in the final state show strong peaks at the very lowmomentum transfer region, which mainly dominate the contribution to the branching ratios, hence suppressing the effect of the TFFs. On the contrary, distributions implicating muons in the final state are much more homogeneously distributed and clearly manifest the neat effect of the TFF, which, in particular, is enhanced for the η decays due to phase space considerations. Our central final branching ratio predictions are summarised in table 9, where a combined weighted average of the results shown in the different tables have been considered and the uncertainties symmetrised. The values of n σ in the table account for the number of standard deviations the experimental measurements are from our predictions. All these predictions are seen to be in accordance with present experimental measurements, only η → µ + µ − γ appears slightly in tension. To end, we would like to encourage once more experimental groups to measure these TFFs.

A.1 Four-body decay width in invariant variables
The partial decay width of a particle P of mass M P decaying into four particles p 1 p 2 p 3 p 4 reads [1] Γ(P → p 1 p 2 p 3 p 4 ) = dΦ(p P ; q p 1 , q p 2 , q p 3 , q p 4 ) (2π) 4 2M P M(P → p 1 p 2 p 3 p 4 ) 2 , (A.1) where dΦ(p P ; q p 1 , q p 2 , q p 3 , q p 4 ) is the four-body phase-space element given by dΦ(p P ; q p 1 , q p 2 , q p 3 , q p 4 ) = δ 4 p P − Following refs. [66,67], the phase space is expressed in terms of independent invariant masses (instead of using three-momenta and angles) as dΦ(p P ; q p 1 , q p 2 , q p 3 , q p 4 ) = 1 8π is the basic four-particle kinematic function. As argued in ref. [66], the limits of integration of the remaining variables,