Radiative corrections to double-Dalitz decays revisited

We revise the radiative corrections to double-Dalitz decays ($P\to\ell\bar{\ell}\ell'\bar{\ell}'$), completing the full next-to-leading order calculation in QED as compared to existing calculations. As a result, we find mild differences with respect to previous studies, that might be relevant for extracting information about the mesons transition form factors. The latter play an important role in determining the hadronic light-by-light contribution to the anomalous magnetic moment of the muon. Finally, we outline an ongoing extension of this work to the $e^+e^-\to e^+e^- P$ processes.


Introduction
The transition form factors (TFFs) of pseudoscalar mesons (P = π 0 , η, η ), commonly denoted as F Pγ * γ * (q 2 1 , q 2 2 ), describe the interaction of the former with two-possibly virtual-photons i d 4 x e iq 1 ·x 0| T { j µ (x) j ν (0)} |P = ε µνρσ q 1ρ q 2σ F Pγ * γ * (q 2 1 , q 2 2 ). (1.1) They have been an object of interest in the past and present years. From one point of view, such an exclusive process with a single-and a priori the simplest-QCD bound state, represents in principle the cleaner one that can be described in pQCD (see Refs. [1,2] for pioneering works).
As such, it serves as a probe of our understanding of QCD dynamics and the pseudoscalar meson structure itself. From a different point of view, and more recently, these have regained interest in connection with their role in the anomalous magnetic moment of the muon, (g−2) µ [3,4,5], which at the moment shows an interesting 3σ deviation [3] and for which new-more precise results-are expected soon from the new experiment at Fermilab [6] and also, in the future, from J-PARC [7].
A key ingredient here is the description of these form factors at the low energies relevant for (g − 2) µ [5,8]. While a wealth of experimental data is available in the single-virtual region [i.e. for F Pγ * γ * (q 2 , 0)] from e + e − → e + e − P collisions and P → γ + − decays, it is only in the past year that the first doubly-virtual measurement from e + e − colliders for the η became available [9]. Still, it would be desirable to gather complementary measurements at lower energies in order to test the most compelling frameworks devised for describing such TFFs in the region relevant to (g − 2) µ [5,10], and also to compare against the recent lattice results [11]. Experimentally, such a possibility can be brought by the double-Dalitz decays, P → γ + − + − , especially for the η and η cases due to their larger phase space [12,13,14,15,16,17,18]. 1 Notoriously, the required statistics for such measurements would be at hand at the proposed REDTOP facility [21]. Still, an accurate extraction of the TFF would not be possible without a proper implementation of QED radiative corrections. The first calculation of such kind for the double Dalitz decays appeared in [22], that included all diagrams necessary to cancel IR divergencies, but did not represent a full NLO calculation, with some diagrams/topologies missing. In our work [23], we revised such calculation, performing a full NLO evaluation in the soft-photon approximation. As a result, we find corrections to the previous estimate that are likely relevant for the extraction of the doublyvirtual TFF.

Full NLO analysis
At LO, the amplitude for double-Dalitz decays (P → + − + − ) reads where the phase-space variables {s i j , y i j , φ }, and also λ (i j) , have been defined in [23,22]. Such result is modified at higher orders in QED. Particularly, denoting the generic amplitude as the NLO contribution to the differential branching ratio stands for the interference of the LO and NLO amplitudes above. Representative diagrams contributing to the NLO amplitudes are shown in Fig. 1. Among these diagrams, the vertex, vacuum polarization, and five-point amplitudes, all of them containing IR divergencies, where accounted for already in [22], while the three-and four-point amplitudes where missing. The last types constitute, together with the revision of previous topologies, our main contribution. Further, in order to obtain a finite infrared (IR) result, bremsstrahlung diagrams need to be incorporated as well (see Fig. 1). Employing the soft-photon approximation, these can be expressed as where Q k and p k stand for the charge and momentum of the k-th particle respectively, and I(p i , p j ) are well-known integrals [24] that can be found, adapted to our case of interest, in [23] 2 . In the following, we revise each of the contributions as compared to Ref. [22]. Concerning the vertex corrections, these amount to shifting the γ µ Dirac matrices in Eq. 2.1 to where F 1 (q 2 ) and F 2 (q 2 ) are the well-known NLO contributions to the Dirac and Pauli form factors. The first one needs to be regularized, together with the self-energies, and contains IR divergencies to be compensated by bremsstrahlung graphs. Its inclusion amounts to a multiplicative factor with respect to the LO result and is in good agreement with Ref. [22]. Concerning the second one, while we agree on the expression for F 2 , we find differences with respect to the resulting matrix element squared. While its impact on the total branching ratio (BR) is negligible [23], it might affect (softly) the differential distributions. Regarding the vacuum polarization, we find good agreement. Concerning the three-and four-point functions, these have been computed for the first time. In comparison to the vertex functions, these involve an integral over the TFF that, for computational purposes, has been assumed that can be taken as constant or be expressed in a propagator-like form (modulo powers of momenta). Further, for the sake of simplicity, a simple factorized Padé approximant has been employed for numerical results. Analytic expressions for them have been provided, for = , in [23] in terms of three-and four-point tensor integrals-the full result has been made accessible as ancillary files on arXiv. The three-point topologies deserve further comments since, unlike for the four-and five-point ones, the amplitude itself is divergent for a constant TFF. In our work, two possibilities have been investigated: to use the form factor or to use chiral perturbation theory (χPT), where a counterterm that is common to P → + − decays is required to render the amplitude finite. For light pseudoscalars and leptons, the three-point topology is very similar in both approaches, while the performance of χPT (at this order) deteriorates for heavier pseudoscalars and leptons.
Finally, the five-point topology has been recomputed employing different techniques as compared to Ref. [22], that uses a tensor decomposition of the five-point loop integrals together with spinor amplitudes. In particular, in our final calculations, we compute directly the NLO piece, after summing over polarizations. This allows to decompose the integral in terms of the scalar five-point function and lower-point tensor ones. As a cross-check, we computed the NLO amplitude in terms of five-point tensor integrals that were decomposed into lower-point ones using the method in [25]. Such method overcomes some numerical instabilities encountered in Ref. [22] connected to vanishing Gram determinants, and resulted in identical numerical results as the first method. Since such an integral contains IR divergencies, we also checked that, once combined with bremsstrahlung graphs, the result was IR-finite. Since final expressions are lengthy, it was not possible to compare to the previous results in [22] that might be a source of the mild discrepancies that we found.

Numerical results
In the following, we provide results on the correction to the total decay width with a softphoton cut-off energy E c , corresponding to x 4 = 0.9985, in analogy to [22]. 3 For the calculation, we employed LoopTools [26,27] to evaluate the loop integrals that we checked provided the same results as our implementation of the methods in [25]. For the numerical integrals we employed Mathematica's NIntegrate method as well as the MonteCarlo routines in [28] for each of the contributions, 4 providing compatible results. As a cross-check, the total contribituon of the five-point amplitudes was compatible with 0 within integration errors, providing an additional cross-check of the correctness of the code [23]. The main results for the BR are given in Table 1, that displays the quantity δ (NLO) = BR NLO /BR LO − 1, where NLO superscripts refers to the results containing LO and NLO contributions to the BR. The second row, δ (NLO), is our main result, obtained for a non-constant TFF, which effect is outlined in the third row. The fourth includes only those contributions in [22] and a constant TFF, so that it can be readily compared to [22], fifth row.
3 E c = m P (1 − x 4 )/2 and x 4 = p 2 4 m −2 P . For relating to different x 4 values, we refer to [23]. 4 For details on results for each contribution separated, we refer to [23].  Table 1: Numerical results for the radiative corrections (second row). The third shows the effects of including a TFF on the LO result. The fourth row contains the same corrections included in Ref. [22] (in fifth row), that show the disagreement. The last row stands for the BR at NLO accuracy.  As advanced, we find some discrepancies-of the order of 1% compared to [22]. The source of this is unclear: a few possibilities have been outlined, but numerics in their MC generator could contribute too. Finally, we note that for heavier pseudoscalars the previously ommitted three-and four-point topologies represent an additional 1% effect-likely relevant for extracting information about the doubly-virtual TFF [17], yet a quantitative statement would require a MC simulation.

Conclusions and Outlook
In this work, we have revisited and completed the full NLO calculation for double-Dalitz pseudoscalar decays in the soft-photon approximation. In doing so, we find small numerical differences with respect to the topologies already accounted for in the existing calculation, the source of which could not be determined precisely. In addition, the new topologies have a small but likely non-negligible effect for the extraction of the doubly-virtual TFF. The results have been provided in a Mathematica notebook available at arXiv as an ancillary file. Also, this work is relevant for completing the NLO corrections for pseudoscalar production at e + e − colliders, which is currently a work in progress.