Radiative corrections to the $\eta^{(\prime)}$ Dalitz decays

We provide the complete set of radiative corrections to the Dalitz decays $\eta^{(\prime)}\to \ell^+\ell^-\gamma$ beyond the soft-photon approximation, i.e. over the whole range of the Dalitz plot and with no restrictions on the energy of a radiative photon. The corrections inevitably depend on the $\eta^{(\prime)}\to\gamma^*\gamma^{(*)}$ transition form factors. For the singly virtual transition form factor appearing e.g. in the bremsstrahlung correction, recent dispersive calculations are used. For the one-photon-irreducible contribution at the one-loop level (for the doubly virtual form factor), we use a vector-meson-dominance-inspired model while taking into account the $\eta$-$\eta^\prime$ mixing.


I. INTRODUCTION AND SUMMARY
Looking for effects of new physics is certainly one of the major contemporary goals of particle physics. The discovery and quantification of new phenomena at this frontier is a very complicated task: it goes along with having our present knowledge under control. Due to the fact that experimental devices and techniques are getting more and more precise, theorists should provide sufficiently low uncertainties together with their predictions. Only in this way an eventual discrepancy can be clearly and correctly revealed. New physics has therefore no meaning without the 'old one' being fully explored. Unfortunately, the low-energy sector of strong interactions remains a significant challenge. The correct incorporation of radiative corrections in the QED sector might help to extract information about the strong sector for the chosen processes. It is the subject of this paper to study the next-to-leading-order (NLO) radiative corrections to the Dalitz decays η (′) → ℓ + ℓ − γ.
Unlike in the neutral pion case, the η (′) Dalitz decays do not belong to those with the highest branching ratio, since due to the higher η (′) rest masses hadronic decay channels are open. Nevertheless, studying these decays provides a way to access the electromagnetic transition form factors and consequently information about the structure of related mesons. The form factors in turn represent a valuable input for precision predictions of some other quantities like the anomalous magnetic moment (g − 2) of a muon. Moreover, these decays are used as normalization channels in rare decay (η (′) → ℓ + ℓ − ) searches.
In [1], motivated by contemporary experimental needs of the NA62 experiment [2], we revisited the classical paper of Mikaelian and Smith (M&S) [3]. We concentrated on the detailed recalculation and completion of the full set of NLO QED radiative corrections to the Dalitz decay of a neutral pion, i.e. to the process π 0 → e + e − γ. Doing so, we have avoided simplifications connected to neglecting higher orders in the final-state lepton mass and thus retaining generality for a future use -in particular for Dalitz decays including muon pairs. The onephoton-irreducible (1γIR) contribution at the one-loop level, which turned out to be indeed non-negligible in view of the two-fold differential NLO decay width, was also included. Finally, we have discussed contributions, which were numerically irrelevant for the pion case but were expected to became important e.g. for η decays. Here, muon loops as a part of the virtual radiative corrections need to be taken into account. Most importantly, we have also touched the fact that the deviations due to the slope of the η transition form factor cannot be overlooked and provided expressions covering a related correction.
It is the aim of the present paper to discuss in detail not only the case of η decays, but also the η ′ decays. In the former case, we could conveniently and extensively draw from the previous work [1], which governed the neutral pion Dalitz decay, since it was already written in a sufficiently general way. Consequently, we would proceed along very similar lines and after properly treating the η-η ′ mixing and generalizing the one-photon-irreducible contribution beyond the effective approach, we could immediately provide the relevant tables and plots. What really brings the current topic to a different level of difficulty is a desire to tackle the radiative corrections for the η ′ decays. The resulting framework is, of course, directly applicable for the η case and one can check that the results are (backward) compatible with the previous, although somewhat generalized, approach used in [1]. Needless to say, the same framework introduced here applies for the pion case as well and it provides some mi-nor corrections compared to the numerical results given in [1]. This stems from the fact that therein we have intentionally neglected the form-factor dependence of the bremsstrahlung correction. Let us only mention that the numerical results obtained for the pion decay using the new framework are indeed compatible with the formfactor slope correction suggested at the end of Section V of [1]. There is thus no particular need to use the presented framework for the pion case: one gains a correction to the correction at the level of 1 % and it would be rather an overkill. For the Dalitz decay of a neutral pion, the approach shown in [1] is sufficient.
Let us briefly discuss the subtleties and difficulties which one encounters and needs to deal with when facing the Dalitz decays of η (′) mesons and associated NLO radiative corrections and which are mainly driven by the properties of the η ′ meson. The main differences compared to the pion case stem from the following facts. First, it is the higher rest mass, which in the case of η is above the production of a muon pair and in the case of η ′ even above the lowest-lying resonances ρ and ω, the former of which is a broad resonance in ππ scattering. This is connected to the fact that the form-factor slope parameter is not negligible as it was in the pion case: the form factor cannot be scaled out anymore and its particular model is required to be taken into account. We then need to distinguish between two separate cases. Similarly to the leading-order (LO) decay width, in the case of the bremsstrahlung correction the singly virtual transition form factor appears. The calculation of this contribution includes integration over angles and energies of the bremsstrahlung photon. For these integrals to be well-defined in order to obtain reasonable results, including the width of the lowest-lying vector-meson resonances becomes necessary. Due to the fact that such a calculation will be unavoidably sensitive to the width of the broad ρ resonance, we have decided to incorporate the recent dispersive calculations [4,5]. In the case of the 1γIR correction, one needs to take into account the doubly virtual transition form factor. Here we do not expect any substantial dependence of the result on the vectormeson decay widths and we use a simple model, which incorporates the strange-flavor content of η (′) mesons and the η-η ′ mixing.
Naïve radiative corrections for the η → e + e − γ process were already published [6] soon after the work [3]: compared to [3], the numerical results presented in [6] correspond to the case in which only the numerical value of the physical mass of the decaying pseudoscalar was changed. Other work related to this paper is [7], where the twophoton exchange contributions to the cross sections of e + e − → η (′) γ processes were calculated. In the current work, we provide a complete systematic study of the NLO radiative corrections to the differential decay widths of the four processes under consideration: η → e + e − γ, η → µ + µ − γ, η ′ → e + e − γ and η ′ → µ + µ − γ.
As it was in the pion case, radiative corrections for these processes are crucial in order to extract relevant information from the data. This goes together with the fact that currently an ambitious experimental η (′) program aiming for an accuracy never reached before is running, for instance, at experiments BES-III [8], A2 [9] or GlueX [10]. Note that in [1] and also throughout the present work we study fully inclusive radiative corrections, i.e. no momentum or angular cuts on the additional bremsstrahlung photon(s) are applied. Consequently we are free of the collinear divergences (sensitivity to the smallness of the electron mass) that would otherwise appear; see, e.g., the corresponding discussion in [11].
The main message of the present work is the completion of the list of the NLO corrections in the QED sector and improving the previous approach [6]. Compared to [6], which relates to the case of the η → e + e − γ decay, we took into account muon loops and hadronic corrections as a part of the vacuum polarization contribution, 1γIR contribution, higher-order final-state lepton mass correction and form-factor effects. Moreover, we treat three additional processes including η ′ decays: η → µ + µ − γ, η ′ → e + e − γ and η ′ → µ + µ − γ. All the formulae necessary for the calculation of the considered correction are listed in the present paper or, whenever a repetition should occur, the reader is referred to the previous work [1]. Let us mention that in contrary to the pion case, the eventual NLO Monte Carlo event generator would not be able to profit from the real-time code calculating the desired correction at a given kinematical point. This is due to the much higher CPU time consumption caused by higher complexity of the involved integrals. On the other hand, sufficiently dense tables of values suitable for interpolation are submitted together with this text in a form of ancillary files.
Our paper is organized as follows. We recapitulate first some basic facts about the LO differential decay width calculation in Section II. Then we proceed to the review of the NLO radiative corrections in the QED sector in Sections III, IV and V. In particular, in Section III we discuss the virtual corrections, in Section IV we introduce the one-photon irreducible contribution and in Section V we describe the bremsstrahlung correction calculation. Some technical details together with extensive results concerning the bremsstrahlung and 1γIR contributions to the NLO correction have been moved to appendices. The building blocks for the 1γIR matrix element in terms of scalar form factors can be found in Appendix A. In Appendix B we present explicit results regarding the bremsstrahlung matrix element squared. This is related to Appendix C, where new basic integrals are listed which are necessary to append to the previously used basis presented in [1] due to a necessary generalization. Appendix D shows the partial fraction decompositions used to simplify the bremsstrahlung matrix element squared. In Appendix E we briefly describe how a simple vector-meson-dominance (VMD)-inspired model for the η and η ′ electromagnetic transition form factors is derived and provide its phenomenological test. Finally, we make use of the doubly virtual transition form fac-P p q k Figure 1. Leading-order diagram of the decay P → ℓ + ℓ − γ in the QED expansion.
tor from Appendix E and show in Appendix F a simple example of the approach discussed in Section IV on the case of the P → ℓ + ℓ − decays.

II. LEADING ORDER
In what follows we will stick to the notation used in [1]. Let us briefly recapitulate it for completeness. Throughout the text we denote the four-momenta of the decaying pseudoscalar meson (of mass M P ), lepton (mass m ℓ ), antilepton and photon by P , p, q and k, respectively. For a parent meson we have in mind η or η ′ . Traditionally, we introduce kinematical variables x and y defined as where x is a normalized lepton-antilepton pair invariant mass squared and y has the meaning of the rescaled cosine of the angle between the directions of the outgoing photon and antilepton in the ℓ + ℓ − center-of-mass system (CMS). If we introduce a lepton-specific constant ν ℓ = 2m ℓ /M P and associated CMS lepton speed we can write the limits on x and y as Consequently, these depend on the final-state lepton mass. The leading-order diagram of the decay P → ℓ + ℓ − γ is shown in Fig. 1. The shaded blob corresponds to the singly virtual electromagnetic transition form factor F ((p + q) 2 ), which is closely related to the doubly virtual transition form factor by (4) The two-fold differential decay rate reads Above, we have used the LO expression for the most prominent (electromagnetic) decay rate of a neutral pseudoscalar meson P : Integrating (5) over y, we find the one-fold differential decay width Moving beyond the leading order, it is convenient to introduce the NLO correction δ to the LO differential decay width, which allows us to write schematically dΓ = (1 + δ + . . . ) dΓ LO . In particular, in the case of the twofold differential decay width we define and in the one-fold differential case we have Concluding the definitions closely related to the previous work [1], such a correction can be divided into three parts emphasizing its origin Here, δ virt stands for the virtual radiative corrections, δ 1γIR for the one-photon-irreducible contribution, which is treated separately from δ virt in our approach due to the reasons of historical development, and δ BS for the bremsstrahlung. As a trivial consequence of previous equations, having knowledge of δ(x, y) allows for obtaining δ(x) using the following prescription: x .
(11) In the following sections we discuss the individual contributions one by one. We mainly point out the differences in comparison to the pion case.

III. VIRTUAL RADIATIVE CORRECTIONS
What we call the virtual radiative corrections is obtained from the interference terms of the LO diagram shown in Fig. 1 and the diagrams in Fig. 2a and 2b. In the pion case, the vacuum polarization was dominated by the electron loop. It turns out though that in the high-invariant-mass region of the photon propagator the hadronic effects become significant, which should be taken into account for the η (′) decays. Thus, in general, we shall deal with the photon self-energy in the form where Π L and Π H stand for the leptonic and hadronic parts, respectively. We discuss these contributions separately in the following. At NLO, the contribution of the vacuum polarization (12) to the virtual radiative correction then reads Note that δ virt Π (x, y) does not depend on y and consistently with (11) we have δ virt Π (x) = δ virt Π (x, y). Regarding the lepton part, all the necessary formulae connected with the contributions of lepton loops to the vacuum polarization are stated in Section III of [1] and hold also in the current cases. However, we shall point out some interesting details. It was already discussed in [1] that in the case of the η (′) -meson decays not only the electron loop but also the muon loop should be taken into account as a part of the vacuum polarization contribution. The contribution of the leptonic vacuum polarization insertion -i.e. the contribution of the lepton loops to the photon propagator -to the correction δ virt can be expressed as (cf. the formulae in Section III of [1]) Here, M 2 P x is the invariant mass of the final-state lepton pair and ℓ ′ stands for the leptons circulating in the loop. Consistently with definition (2) we have β ℓ ′ ≡ β(M 2 P x; m 2 ℓ ′ ) and |β ℓ ′ | = |β 2 ℓ ′ | . Let us note that β ℓ ′ no longer stands for the CMS speed of the final-state leptons as was the case for β ℓ in (2). Since β depends on x, which for the given process of a decay to a lepton pair of flavor ℓ satisfies the limit x ∈ [ν 2 ℓ , 1], we can run into different kinematical regimes, which are explicitly covered by formula (14); for this purpose the Heaviside step function θ was used. Whenever β 2 ℓ ′ < 0, i.e. β ℓ ′ itself becomes imaginary, no on-shell lepton-antilepton pair can be created in the loop and the diagram in Fig. 2a lacks the imaginary part. This happens if the invariant mass of the final-state lepton pair is under the production threshold of the two loop leptons of flavor ℓ ′ : ν 2 ℓ ≤ x < ν 2 ℓ ′ . In turn, this condition can be, of course, only realized in the case of the muon loop contribution to the NLO decay widths of the π 0 → e + e − γ decay and in the process η (′) → e + e − γ, where the kinematical condition ν 2 e ≤ x < ν 2 µ is met within a part of the kinematically allowed region.
The hadronic contribution to the photon self-energy can be expressed via a dispersive integral [12] Here, σ H is the total cross section of the e + e − annihilation into hadrons and is related to the ratio R(s) as σ H (s) = (4πα 2 /3s)R(s). In order to obtain a smooth curve from the data [13], we fit the experimental values in a similar way as it was done in [12]. Concerning the upper bound of the integral in (15), in practice it is sufficient to use the scale Λ cut ≃ 4 GeV. It is set in such a way that all the important resonances are covered while the higher-energy region does not contribute significantly for s < ∼ M η ′ . The resulting contribution to the radiative corrections is plotted as one of the constituent curves in Fig. 4.
Instead of (13), one might take into account the whole geometric series of vacuum polarization insertions and sum it in order to find Eq. (13) is then recovered by taking the linear expansion of (16), since |Π(s)| ≪ 1. Note that in (16) one should use a full form of the vacuum polarization including its imaginary part, similarly as it was e.g. defined in (20) of [1] for the lepton case: with γ ℓ ≡ (1 − β ℓ )/(1 + β ℓ ). But since |Π(s)| ≪ 1, also | Im Π(s)| ≪ 1 + Re Π(s) and 2 Re Π(s) ≫ |Π(s)| 2 . One could thus safely use only the real part of Π(s), as it was defined in (13), for the purpose of the numerical evaluation of (16). This is equivalent to effectively using Im Π(s) → 0. On the other hand, the numerical difference between the definitions (13) and (16) is ∆δ virt Π (x) ≃ 0.25 % at s ≃ M 2 ω . Therefore we take into account the precise formula (16) together with (15) and (17).
For completeness, let us revise at this point what we mean by the virtual correction. In agreement with Eq. (16) in [1], we use where Π(s) contains not only single electron and muon loops, but also the whole hadronic contribution to the photon self-energy. For the form factors F 1 and F 2 , which arise from the vertex correction diagram in Fig. 2b, the reader is referred to seek the expressions in [1].

IV. ONE-PHOTON-IRREDUCIBLE VIRTUAL RADIATIVE CORRECTION
Considering the QED expansion, the NLO diagrams contributing to the one-loop one-photon-irreducible (1γIR) correction to the P → ℓ + ℓ − γ process are shown in Figs. 2c and 2d. We treat this contribution separately from the virtual correction to emphasize the fact that it was not included in the original approach [3]. Therein, it was considered to be negligible already for the pion case. This statement had been corrected in Ref. [14] many years before the debate about this issue was finally closed. Note also that in the pion case, the 1γIR contribution was studied in [15] in connection with the bremsstrahlung correction to the π 0 → e + e − rare decay.
In the 1γIR contribution one cannot factorize out the electromagnetic transition form factor. This correction becomes therefore unavoidably model-dependent already at the two-fold differential level and it is necessary to choose a particular model to evaluate the correction numerically. Compared to the previous calculations of the LO diagram and NLO virtual radiative corrections, a doubly virtual transition form factor F P γ * γ * (l 2 , (P − l) 2 ) is needed to be used. This form factor enters the loop and the integration over the unconstrained momentum l is then performed. Note that P stands for both the decaying pseudoscalar as well as for its four-momentum.
Let us now proceed further and consider that during the calculation of the 1γIR loop diagrams the following structure inevitably appears: By construction, the arguments of the form factor coincide with the photon propagators in the loop; l denotes the loop momentum as is usual. In what follows we consider a family of large-N c motivated analytic resonancesaturation models, which are in detail discussed in [16] (N c denotes the number of colors). Here F P γ * γ * is a rational function with vector-meson poles. We then realize that due to a use of the partial fraction decomposition one can perform -within the loop integrals appearing during the evaluation of the diagrams -the following substitution: (20) Note that from now on we will suppress the explicit writing of the '+iǫ' part. Above, In this way we come from the matrix element M 1γIR to α i M h 1γIR [h i ]; note that the normalization constant F π 0 γ * γ * (0, 0) = −N c /(12π 2 F π ) from decomposition (20) is already included in M h 1γIR [h i ] and that we have used a shorthand notation h i ≡ h(c 1,i , c 2,i , M 2 1,i , M 2 2,i ). In order to get results for the whole family of models it is necessary to analytically integrate over the loop momentum just once; this is the main advantage of this approach. At the end one can choose the particular model of the form factor by setting the parameters in the final matrix element appropriately. We can find the above used constants c 1 and c 2 from projecting on the product of the normalized form factor and the photon propagators: for instance in the case M V1 = M V2 we have This little trick is highly convenient when it is necessary to create a universal code for calculating radiative corrections within different models. Let us also note that substitution (21) does obviously not conserve term by term the desired property of the doubly virtual form factors -the symmetry in their arguments. This ansatz though works generally for the rational models mentioned above, which might have a rather complicated structure. The symmetry in question is then always restored in the final result after including all the pieces M h 1γIR [h i ] . Moreover, we realize that if we define we can immediately write This trivial observation simplifies the loop integration even further. Instead of (20), we can then write During the calculation of the amplitude it is only necessary to perform the following substitution: The desired final amplitude for the particular model is then obtained by writing a suitable combination (in spirit of (25)) of such amplitudes which are calculated using substitution (26). One only needs to insert the correct masses and coefficients into this combination, which goes along the lines of (24), (22) and (21). Lets discuss the previous procedure on a particular case and consider for a while the eta-meson decays: P = η. The simplest physically relevant model we can imagine is based on the VMD scenario, i.e. by an ansatz which assumes that the form factor is saturated by the lowest-lying multiplet of vector mesons. It has the following form: Above, φ is the η-η ′ mixing angle and f ℓ together with f s are the associated decay constants in the quark-flavor basis of the quark currents [17,18]; for further details concerning derivation of this model see Appendix E. It is then clear after counting of the loop-momenta powers that such a form factor guarantees the UV convergence of the loop integrals. The matrix element for such a form factor can be schematically written as Following the subsequent decomposition (24) and using linearity, one then finds As the only building block, one needs to calculate M h 1γIR g(M 2 1 , M 2 2 ) obtained in terms of substitution (26) in the original matrix element. For the purpose of the 1γIR contribution to the correction δ(x, y) at the one loop level, the same decomposition will work. Indeed, an interference of M 1γIR with the LO matrix element needs to be considered and thus linearity is preserved. Therefore it is necessary to take into account a linear combination of δ h 1γIR defined according to the prescription (26) in [1]: The explicit expressions for the building-block form factors A and T are shown in Appendix A. Let us have a look at the ratio of the electromagnetic form factors in (30). Clearly, irrespective of the model that we use for the doubly virtual form factor, it is convenient that the normalization of the form factor in such a model equals to F (0) appearing in the LO expression; see also (4). This leads to the fact that the correction will be independent of any overall normalization effects. In the next section we introduce a spectral representation of a normalized form factor; see (41). Following what we have just assumed, we can write and using the VMD-based scenario for instance for η as in (27), we find A simple example of the presented approach applied on the P → ℓ + ℓ − decays is included as Appendix F. There is also a complementary technique how to cover a whole set of form factors under consideration. Instead of putting the particular form factor into our diagrams, we can use the local Wess-Zumino-Witten (WZW) term. In other words we trade the form factor for the constant given by the chiral anomaly. It is then clear from simple considerations, that the contributions from Fig. 2c need counter terms to compensate UV divergences. The convergent part of such a counter term carries an undetermined constant χ (r) (µ) renormalized at scale µ, which can effectively mimic the high-energy behavior of the would-be complete form factor. Using a proper matching procedure, it is possible to acquire a numerical value of this constant χ (r) (µ) for a given form-factor model. Up to mass corrections we can use an approximate formula to estimate this effective parameter (see Eq. (46) in [16]). The question is, if this procedure can be used also for a box diagram in Fig. 2d which is already convergent for the local WZW form factor. It turns out that the corrections are of order m 2 ℓ /M 2 V and M 2 P /M 2 V . Hence, for the pion case this assumption works well. On the other hand, for η (′) it does not and one would need to introduce additional effective parameters in a consistent way. This is though not a trivial task and is beyond the scope of this work. Let us also mention that χ (r) (µ) enters the corrections being multiplied by ν 2 ℓ and its effect is thus negligible for the decays with electrons in the final state.
In the result section (Section VI) we will use the simple VMD-inspired model for the doubly virtual transition form factor to estimate the importance of the 1γIR contributions; for model details see Appendix E. For completeness, let us then present the correction δ 1γIR (x, y) within the model under consideration: For the η case, we have κ η ℓ ≡ 5 sin φ fs . The η ′ case then corresponds to the simultaneous interchange {cos φ → sin φ, sin φ → − cos φ}, of course followed by specifying the right M P within the building blocks δ h 1γIR . For the resonance masses we put M ℓ ≡ M ρ/ω (the average physical mass of ρ and ω mesons) for the light sector and M s ≡ M φ (the physical mass of the φ meson) for the strange sector. However, we would like to stress again that the framework presented in the present section is capable of dealing also with more sophisticated form-factor models. In general, the model dependence of part of the radiative corrections is, of course, a nuisance. But for the doubly virtual transition form factors of η and η ′ there is at present no alternative. In the next section the singly virtual transition form factors are required. Here the model dependence can be reduced by the use of dispersive methods [4,5].

V. BREMSSTRAHLUNG
In this section we briefly build on [1] -which among others discussed the bremsstrahlung correction calculation for the pion case -and show the differences which come into play due to the fact that we are now interested in η (′) mesons. We also show some techniques how to deal with the obstacles which arise. Concerning the notation, we will restrict ourselves to the one used in previous works.
The diagrams which contribute to the Dalitz decay bremsstrahlung are shown in Fig. 2e. Their contribution is (among others) important to cancel IR divergences stemming from the virtual corrections depicted in Fig. 2b. The corresponding invariant matrix element (including cross terms) can be written in the form where Here, we use l and k for the photons and p and q for the electron and positron four-momenta, respectively, and F ((p + q) 2 ) ≡ F P γ * γ * (0, (p + q) 2 ). Note that we use the shorthand notation for the product of the Levi-Civita tensor and four-momenta in which ε (k)... = ε µ... k µ . Inasmuch as an additional photon comes into play it is convenient to introduce a new kinematical variable which stands for the normalized invariant mass squared of the two photons It has the similar meaning as x in the case of the leptonantilepton pair. The contribution of the bremsstrahlung to the next-to-leading-order two-fold differential decay width can be written as The above used operator J is defined for an arbitrary invariant f (k, l) of the momenta k and l as follows: In the case of the pion Dalitz decay, the value of the slope parameter of the form factor is small: a π ≃ 0.03. Consequently, the form factor F ((l+p+q) 2 ) which enters (34) can be conveniently expanded in the following way: Therefore, F ((l + p + q) 2 ) can be approximated by F (M 2 P x) for the process π 0 → e + e − γ. This squared leads to the direct cancellation with |F (M 2 P x)| 2 appearing in the leading-order expression. The bremsstrahlung contribution to the radiative corrections then becomes effectively independent of the particular model of the pion transition form factor. However, in the case of an η meson, the slope parameter is a η ≃ 0.5, which is definitely not negligible anymore. One would need to include higher-order corrections in expansion (40), the convergence would be slower and things would become in general more complicated since additional terms would need to be treated; see also the discussion at the end of Section V of [1]. The real obstacles though appear with the η ′ -meson case. Due to the fact that a η ′ ≃ 1.4, expansion (40) is not applicable at all. One thus needs to calculate with the full form factor.
Although in such a case the situation is somewhat different compared to the one in which the form factor cancels out, in general it is possible to use a similar framework as in the pion case -at least in the sense of treating the kinematical integrals. Accordingly, one needs to rewrite the bremsstrahlung correction in terms of integrals which are known from the pion case. These need to be somewhat generalized due to the presence of poles in the form factor; for results see Appendix C. This becomes more important in the η case compared to the pion decay and needs to be taken into account explicitly. For the η ′ case, this procedure then becomes absolutely crucial since in the hadron spectrum the mass of the η ′ meson lies above the masses of the lightest vector-meson resonances ρ and ω and close to the φ-meson mass. The subsequent (numerical) integrations over all the relevant kinematical configurations of the bremsstrahlung photon become significantly nontrivial due to the running over these poles, which are regulated by incorporating physical widths of the resonances. The narrow resonances like ω and φ are somewhat straightforward to include. However, the width of the broad ρ resonance is sensitive to the π-π scattering. This can be also taken into account by the use of recent dispersive approaches [4,5]. To this extent, it is convenient to use the Källén-Lehmann spectral representation of the Feynman propagator, which allows for the use of a common spectral density function for all the mentioned resonances. The facts stated above make the bremsstrahlung contribution -especially in the case of the η ′ meson -significantly sensitive to the form-factor model. To minimize the model dependence, whenever possible, data and general principles of QFT, analyticity and unitarity, are used that give rise to a dispersive framework.
In the Källén-Lehmann spectral representation, the form factor has the following form: Concerning the integration range, we restrict ourselves to (4m 2 π , Λ 2 ). This is sufficient for our purpose and covers all the important physics. In the lower bound, m π is the mass of a charged pion π ± and 4m 2 π then constitutes the squared mass of the lowest hadronic state that can couple to a photon. The upper bound is governed by the cutoff Λ ≃ 1.05 GeV chosen in such a way just to cover the peak and width of the φ resonance. The spectral function has in the chosen energy range two main contributions distinguished by isospin, The narrow resonances ω and φ contribute to the isospin-zero part with the narrow-resonance spectral function (43) Note that after a full integration over s one indeed gets (up to a sign) the resonance propagator: Finally, let us mention, that a one-narrow-resonance VMD model for the form factor can be written in terms of (41) with A(s) = A V (s). The isospin-one part is governed by the broad ρ meson. In order to include the important effect of ππ scattering, we use the dispersive approach of [4,5]. There the spectral function has the form where F π = 92.2 MeV is the pion decay constant, Ω(s) is the Omnès function, a dispersive tool incorporating pion rescattering, and P (s) and R(s) = (1 + α V s) are polynomials related to the η → ππγ reaction amplitude and the pion vector form factor F V (s) = R(s)Ω(s), respectively. In the case of η ′ , we use P η ′ (s) = (1 + α η ′ s + β η ′ s 2 ) with α η ′ = 0.99(4) GeV −2 and β η ′ = −0.55(4) GeV −4 [5]. For the η spectral function, we take (based on [4]) P η (s) = (1 + α η s)R(s) with α η = 1.32(13) GeV −2 . Note that in view of (45) and taking into account only terms linear in s, P η (s) = 1 + (α η + α V )s. The rest of the parameters from (42) and (45) are given in Table I. For numerical reasons, the spectral function for the broad ρ resonance may be fitted, since it is much faster to numerically integrate the analytical expression compared to the dispersive data interpolation. To this extent, it is  Table I. Values for w are taken from [4] and values of κ were calculated using the following prescription therein: The numerical value of α η V was estimated from the fit in the upper panel of Fig. 1 in [4] and the value of α η ′ V is based on values tabulated in [5].
necessary to model the resonance peak behavior. The following function copies the dispersive shape of A 1 satisfactorily: (46) The energy dependence of the width of the ρ mesonassuming the main contribution comes from the 2π-decay -can be expressed as The advantage of the latter approximation is the fact, that when we insert it into (46) and subsequently into (41), the form factor F (q 2 ) can be evaluated analytically. This speeds up the numerics even further. The values of the fit are shown in Table II. For the η ′ case, the formula  Let us now show, how the matrix element squared looks like in the spectral representation. For the sake of writing down its structure, we recast (34) as where we have used the shorthand notation E = (l + p + q) 2 + iǫ and F = (k + p + q) 2 + iǫ and defined I E ≡ I(k, l) and I F ≡ I(l, k). After inserting the representation of the form factor (41) we get (49) An attentive reader might have noticed that in the previous formula the explicit iǫ is not necessary. At the same time, it might be sensible to explicitly write the infinitesimal imaginary part of the propagator in the following expressions. To avoid any potential confusion further on, let us define e ≡ Re E and, similarly, f ≡ Re F . After the operator J is applied on the matrix element squared and summed over all the spins and polarizations |M BS | 2 ≡ sp., pol. |M BS | 2 , we can actually use the sym-metry k ↔ l ⇔ E ↔ F term by term, which results into In order to perform the integrations of the J operator on the respective terms, we need to do a few fractionproduct decompositions; the procedure is described in detail in Appendix D. Taking also into account that Following the definition (41), we made use of the fact that With the previous simplifications (43), (46) and (47), the above integral can be evaluated analytically which speeds up the time-demanding numerics.
For the sake of integrating out the bremsstrahlung photon it is convenient to define the following three rescaled parts of the matrix element squared: One gets the remaining building blocks of (52) by performing the limit s → 0, i.e. taking Tr E (0) and Tr EF (0). Note that it does not matter if one uses 1/E * or simply 1/E for the Feynman propagators in the definitions of Tr E (0) and Tr EF (0) since these give real contributions. The name Tr is chosen to be in agreement with [3]. Note that the original work of Mikaelian and Smith (M&S) -which in some sense corresponds to A(s) = 0 since the form factor could have been factorized out -is connected with the previous definitions through where For numerical reasons, the integration of the resulting expression is performed term by term. We can name these terms for further convenience. Consequently, we write where Terms t 1a and t 1b are the only ones present in the M&S case. Terms t 1a and t 2 are dominant. By construction (cf. (D6)), terms t 3a and t 3b belong together since separately they develop peaks which are exactly compensated only in the sum. Let us mention that for the purpose of (59), in t 3a it is sufficient to take only the real part of I(V 0 ) since Tr EF (0) gives only a real-part contribution. For the correction δ BS (x, y) it holds Together with (D6) this suggests that we need to perform five subsequent integrations: two nontrivial analytical ones are implicitly hidden in the J operator and at least two need to be inevitably performed numerically. Let us briefly look at the explicit integral structurefor instance on the case of t 4 -and how it contributes to the bremsstrahlung correction. In terms of definitions (53) and (56) we have The term Re{I(V 0 − s)} can be evaluated analytically by taking the fit (46) of the spectral function instead of its numerical form (45) given by the term containing the Omnès function. This procedure significantly speeds up the numerics. The evaluation of J{Tr EF (s)} is performed analytically and the strategy is similar to the one used in [1]. Our goal is to rewrite the expressions under consideration into basic integrals, which are then listed in Appendix D of [1] and Appendix C of the presented work. For this purpose, symmetries and properties of operator J are used together with the reduction procedure described in detail in Section V of [1] and summarized in Appendix B therein. Due to the presence of the effective mass s in the denominators, this procedure needs to be slightly modified since now also E − s appears instead of simple E. As a trivial example, we choose the following reduction: Note also that the IR-divergent terms (those which diverge after integrating over x γ ) need to be treated separately. The current case is discussed in Appendix B. For further details and an introduction we refer an interested reader to [1].

VI. RESULTS
In the previous sections, we discussed the three main parts of radiative corrections as we referred to them in (10), i.e. the virtual corrections in Section III, the 1γIR correction in Section IV and the bremsstrahlung in Section V. In order to get the final result, one then simply sums over the partial results: Eqs. (18), (33) and (67). Let us now comment once again on these contributions in detail. The model-independent expression for the virtual corrections (18) includes the photon self-energy contribution (12), consisting of the exact (at NLO) leptonic (17) and the experimental-data-based hadronic (15) parts, and the vertex correction expressed in terms of the F 1 and F 2 form factors. The expressions for these form factors can be found in [1]; in particular see Eqs. (22) and (23) therein. The form factor F 1 then contains the IR-divergent piece which cancels with the corresponding term stemming from the bremsstrahlung contribution, as shown at the end of Appendix B. The bremsstrahlung contribution (67) itself consists of the fully model-independent ingredients (60), already present in the M&S case, and the pieces (61-63) that still contain some model dependence, though mitigated by our dispersive, data-driven input. All the necessary definitions are then presented in Appendices B and C. Finally, for the model-dependent 1γIR contribution (33) we used the VMD-inspired model discussed in Appendix E. Let us again stress at this point that we used a rather general approach applicable for a wide family of rational models. The final result within a particular model can then be found as an appropriate linear combination of the building blocks (30), with the form factors A and T defined in Appendix A. As another example of a suitable model we mention the one based on the Resonance Chiral Theory framework [19]. There are other -data-drivenmodels [20,21], which are then not directly applicable using the presented approach. However, we expect that different models would lead to compatible values in the numerical results.
The overall NLO radiative corrections to the two-fold differential decay widths and their comparison to the main constituents are shown, for the case y = 0, in Fig. 4. In light of a direct comparison of our results for the η → e + e − γ process to the previous work [6], we can mention a significant difference between both approaches. Taking into account all the discussed contributions, a table of values of corrections δ(x, y) to the process η → e + e − γ similar to the one provided in the original work [6] can be produced at the very same points; see Table III. The difference can then be seen even more explicitly when comparing directly the corresponding values in Table III with Table I in [6]. In the same manner, values for the decay η → µ + µ − γ are presented in Table IV. From the plots in the bottom panels in Fig. 4 we can see that for the case of an η ′ meson, such a table would need to be much denser to cover the wiggle behavior. For this reason, the table-like plots were created, which show the overall NLO radiative corrections for various values of y; see  Tables V and  VI. In Fig. 4 we can see that the 1γIR correction has an intriguing character. For the decays η (′) → e + e − γ it is negative for the whole range of values of x and enhances thus the effect of the M&S-like corrections which are also negative in a wide range of x. On the other hand, for the process η (′) → µ + µ − γ the situation is different due to higher masses involved. In this case, δ 1γIR is very close to zero, mostly slightly negative, but then it starts to be significantly positive. Its dominant behavior for large x contributes to the fact that δ(x, y) is positive over nearly the whole range of x.
The form-factor effects are significant and shown in Fig. 4 as large-spaced dotted lines. Out of these contributions, δ BS t2 is most prominent and responsible -together with the hadronic part of the vacuum polarization -for the wiggles in the case of η ′ decays. The terms δ BS t3a + δ BS t 3b and δ BS t4 are on the other hand -although being most complicated and time consuming to calculate -rather suppressed. Together with another numerically small contribution to the bremsstrahlung correction, δ BS t 1b , they can be even considered negligible in the case of the η (′) → µ + µ − γ. In the case of the decays into electrons, the considered correction is of order 0.1 % for small |y|, otherwise negligible, too. Except for the already mentioned form-factor-dependent term δ BS t2 , there is then only one significant contribution to the bremsstrahlung correction δ BS t1a , which comes with its separately treated IRdivergent counterpart δ BS D . In order to fully understand the effect of the amusingly looking wiggles showed in Fig. 4 and Fig. 5 and eventually appreciate the size of the calculated radiative corrections, we present the differential decay widths in Fig. 6. In these plots, we can see a direct comparison of the LO and NLO decay widths, together with their sum. In general, the total radiative corrections look rather negligible, with an exception of the process η ′ → e + e − γ. In this case, the corrections are significant and the NLO differential decay width is sizable compared to the LO width. The wiggles, which appear only in the η ′ case, then contribute in such a way that they change the shape of the resonance peaks and make them smaller. In particular, Fig. 6c shows that the height of the ω peak is considerably influenced by the radiative corrections. This might be interesting for the extraction of ω properties or of the ω-η ′ interplay. One might deduce such information from η ′ → ωγ → e + e − γ or from η ′ → ωγ → π + π − π 0 γ. It can be expected that the radiative corrections are different for these two decay branches. Ignoring such radiative corrections in the analyses of these decays might lead to contradictory conclusions.
Finally, let us shortly mention that the hadronic part of the vacuum polarization has, in contrary to the η case, a significant impact on the corrections for the η ′ decays, with up to a ∼5 % effect in the M 2 η ′ x ≃ M 2 ω region.
ACKNOWLEDGMENT T.H. would like to thank Uppsala University and University of Birmingham for their hospitality. We would also like to thank A. Kupść for providing us with the Omnès function data used in [5].
This work was supported by the Czech Science Foundation grant GAČR 18-17224S and by the ERC starting grant 336581 "KaonLepton".
is depicted as a dashed line. It would have directly corresponded to the M&S correction [6] if muon loops had been included to the vacuum polarization and higher orders in lepton masses (O(ν 4 ℓ )) had not been neglected in [6]. This contribution also includes the IR-divergent part of the bremsstrahlung. The hadronic contribution δ virt Π L to the virtual radiative corrections is shown as a small-spaced dotted line. The part of the bremsstrahlung correction which is dependent on the model of the form factor and is nonzero for any nonvanishing spectral function and which corresponds to the sum δ BS t 2 + δ BS t 3a + δ BS t 3b + δ BS t 4 is shown as a large-spaced dotted line. The one-photonirreducible contribution δ 1γIR is then shown as a dash-dot line. The divergent behavior of δ(x) near x = ν 2 ℓ has its origin in the electromagnetic form factor F1(x) and is connected to the Coulomb self-interaction of the dilepton at threshold.  Figure 5. The NLO radiative corrections δ(x, y) given in percent for η ′ → ℓ + ℓ − γ processes. Different lines correspond to different values of y. These start at y = 0 and are equidistantly separated by 0.1. The last value is y = 0.99 for the electrons and y = 0.95 for the muons due to kinematical restrictions. The corrections (integrably) diverge at x → ν 2 ℓ , which happens only for y = 0 since for y > 0, x = ν 2 ℓ is not kinematically allowed: indeed, for given x, |y| < β(x) and β(ν 2 ℓ ) = 0. These plots are complementary to Tabs. V and VI, in which the wiggle structure could not be precisely covered.  Figure 6. The two-fold differential decay widths dΓ(x, 0) at NLO (solid line) and its constituents for the η (′) → ℓ + ℓ − γ decays. The LO differential decay width for y = 0 is shown as a dotted line. The corresponding NLO contribution to the differential decay width is represented by a dashed line. Table III. The overall NLO correction δ(x, y) given in percent for a range of values of x and y (i.e. the Dalitz-plot corrections) for the process η → e + e − γ. It is sufficient to show the results only for positive values of y since the corrections are symmetric under y → −y. The larger values at the edge of the kinematical region (as x → 1) are naturally present due to the fact that the correction itself is defined as a ratio of the NLO and LO decay widths which both vanish for x → 1. The LO differential decay width, however, approaches zero much faster than the interference term of the LO graph with the 1γIR one. Let us also note that these larger values are beyond any practical significance in the final spectrum (see Fig. 6).  Table IV. The overall NLO correction δ(x, y) given in percent for the process η → µ + µ − γ. The dot-filled entries correspond to kinematically forbidden combinations of x and y. See also the caption of Table III   Appendix A: Matrix element of the one-photon-irreducible contribution The building-block matrix element for the one-photon-irreducible contribution can be expressed in terms of scalar form factors in the following (manifestly gauge-invariant) way (cf. (24) in [1]): The form factor P does not contribute to δ 1γIR , so we do not include it in this appendix. Note also that only the box diagram (see Fig. 2d) contributes to A(x, y). Let us also define the two-point (bubble) scalar one-loop integral as a reference point for the used notation: Above, we have introduced ε = 2 − d/2. Before providing the final expressions, let us define Using the dimensional regularization, the dimensional reduction scheme [22,23]) and Passarino-Veltman reduction [24], the explicit results for the two contributing form factors, A and T , in terms of scalar one-loop integrals read (we use again for simplicity m ≡ m ℓ , ν ≡ ν ℓ and M ≡ M P ) (A6) Appendix B: Bremsstrahlung matrix element squared In this appendix we provide explicit formulae for the building blocks of the bremsstrahlung invariant matrix element squared, as described in Eqs. (52-63). If we generalize (54) to we only need to provide its explicit form in order to cover Tr E 2 and Tr E (s) defined in (54) and (55). Indeed, one can obtain Tr E (s) by means of multiplying Tr E 2 (s) by the appropriate propagator: Tr E (s) = (E − s)Tr E 2 (s). For completeness, Tr E 2 = Tr E 2 (0). In the following expressions, as denoted below, the symmetries of the J operator were already taken into account and the reduction procedure (in the sense of Appendix B of [1]) into the basic integrals was performed as well. Before we get to the desired expressions, we need to define a complete set of the Feynman denominators (while suppressing the '+iǫ' part): A = l · q , B = l · p , C = k · q , D = k · p , E = (p + q + l) 2 , F = (p + q + k) 2 . (B2) For simplicity, in what follows we use m ≡ m ℓ , ν ≡ ν ℓ and M ≡ M P . The term which includes squares of the bremsstrahlung diagrams reads the IR-divergent part of which, up to G 2 (s) and the y → −y part, is given as This is in agreement with Appendix A of [1] up to an additional factor of 1/2 visible in matching (58). Above, we have used the following definitions: The second building block which arises from the interference of bremsstrahlung diagrams under consideration can be expressed as where we have defined the function Tr X(E−s) (s; X,x, ξ) Needless to say, whenever (y → −y) appears above it applies for entire expressions. Let us note that after inserting the expressions (B3) and (B10) into formula (58), one indeed recovers expressions presented in Appendix A of [1]. In (59), which we can reformulate into only t 1a and t 2 contribute to the overall IR-divergent part, which can be rewritten as Tr div = 4(t div 1a + t div 2 ) = 4 1 2 Tr div E 2 + ds A(s) We already know that J{Tr div E 2 } = J{Tr D }+(y → −y); note that G(0) = 1. In addition, by definition, the IR-divergent part of Tr E (s) is given by the divergent part of (E − s)Tr E 2 (s). To this, only Tr div E 2 can contribute. Moreover, we realize that E − s = M 2 x − s + 2A + 2B and that only its M 2 x − s + iǫ part does not suppress the essentially divergent terms 1/A 2 , 1/B 2 and 1/AB appearing in Tr div E 2 . Consequently, J{Tr div E (s)} = (M 2 x − s + iǫ)G 2 (s)J{Tr div E 2 }. We thus find J Tr div = 2J Tr div which, using (B5) and (41), translates into J Tr div = 2J Tr div Due to the symmetries of the double integral and the fact that in the end we are nevertheless interested only in the real part, the term in the square brackets can be further recast into Taking the real part of the previous result, we finally get (cf. definition (41)) The form factor squared cancels after normalizing on the LO differential decay width and, subsequently, we get the same δ BS D (x, y) as in [1], Appendix A. This is the desired result, since δ BS D (x, y) contains the terms which exactly cancel the IR divergences stemming from the virtual corrections: δ BS D (x, y) + 2 Re {F 1 (x)} is then IR-safe. Let us recall that the integration of J{Tr D } over x γ needs to be performed analytically. Note that it was only kept as a part of J{Tr E 2 (s)} since it also contributes to the IR-convergent part of J{Tr E (s)}.
where we have used v 1 = (A + C) and v 2 = (B + D) + g(s). Note that (C16) includes an IR-divergent part J[1/A 2 ]/g(s). Let us also mention that even though it might seem like that after inspecting definition (B3) of J{Tr E 2 (s)}, it is not necessary to define other integrals, e.g.
since only Tr E 2 (0) and (E − s)Tr E 2 (s) appear in the final expression for δ BS (x, y). For completeness, let us take advantage of our notation and add that J[1/(C(E − s) 2 )] = Q(0, v 1 , v 2 ; −g(s)) . All the other terms can be found in Appendix D of [1]. At this point, let us mention that in [1], terms (D20) and (D21) were not listed properly in light of definition (D1). It was not stated explicitly that the second version of (D1) was used for evaluating all the terms. This (second) version, however, equals to the original definition only for a > 0, as mentioned in (D1). When a < 0, an additional overall minus sign appears. This was not emphasized and shall at least be clarified at this point. Taking thus into account (C18), which is universal for any a, we find Other expressions in Appendix D of [1] require no changes. Let us also mention that one should indeed get now (C19) by putting s → 0 in (C14).
Appendix D: Partial fraction decompositions for the bremsstrahlung matrix element squared In this appendix we show how to rewrite the bremsstrahlung matrix element squared (50) into a form which will allow us to perform the integrations of the J operator on the respective terms. To achieve this we need to use a few fraction-product decompositions. First, we take the simplest case: Note that E represents a Feynman denominator corresponding to the virtual photon and in this and following applications the difference between E and E * plays no role. We also write s instead of s − i(ǫ + ǫ ′ ) since s = 0 due to the positive limits of the integration. Similarly, In the denominators, ǫ and ǫ ′ represent positive and infinitesimally small independent numbers and the numerical result remains the same even when we assume for simplicity that ǫ ′ → ǫ. Now, we can use the fact that this term is multiplied by two spectral functions A(s) and A(s ′ ) and integrated symmetrically over s and s ′ . After we rename s ↔ s ′ in the second term on the right-hand side of (D2), we realize that due to the symmetric integration we obtain the complex conjugate of the first term. Since we are anyway interested in the real part only, we just get the factor of 2. In the following we use the knowledge of the fact that e + f = M 2 (D3) Note that due to the presence of the J operator, we can substitute 1/F * by 1/E * on the right-hand side of the previous equation. Above, ǫ ′ is actually artificial and can be safely sent to zero. This can be also seen due to some other reasons. First, on the left-hand side we could have already written f instead of F * . Another, more general, reasoning is similar to the following case. The last term which needs to be treated can be decomposed as Again, ǫ and ǫ ′ are infinitesimally small and independent. For a moment we can assume these to be finite though small. The left-hand side does obviously not depend numerically on the fact if ǫ is greater or smaller than ǫ ′ so the apparent dependence of this type needs to be canceled on the right-hand side. For this purpose both of the terms in the brackets are necessary. Without loss of generality we can assume that ǫ > ǫ ′ and put ǫ = ǫ ′ + δ with δ > 0 small. Since ǫ and ǫ ′ can be arbitrary for (D4) to be valid, we can set ǫ ≃ ǫ ′ . In the corresponding limit δ → 0 we can apply the Sochocki's formula to find where p. v. stands for the principal value. The part containing the delta function then trivially vanishes (due to the first bracket) together with the dependence on the sign of the difference ǫ − ǫ ′ . Again, using the fact that the integration is symmetric in s and s ′ and that we can on the right change F * → E * , we obtain the sum of a term and its complex conjugate, which results in a factor of 2 under the real-part operator; note that I E I * F has no (nonvanishing) imaginary part since it is related to the tree diagram.
Next we define the correlator η A V V for each basis state (again A ∈ {ℓ, s}): The factors Z η A = 0|j A |η A are related to the pion case Z π = B 0 F π by Z η A = Z π f A : we have introduced the ratio of the decay constants f A ≡ F A /F π . For the ηV V correlator we can then finally write The √ 2 factor comes from (E11) and the mixing factors from (E12). To avoid difficulties connected with calculation of the decay constant values and defining the final ηV V correlator all the way through (E14), we will use the VMD ansatz These values are used throughout the paper.

Model relevancy check
Now we have to check if such a VMD-inspired model is phenomenologically successful and if it is suitable for calculating 1γIR contribution to the NLO virtual radiative corrections.

a. Decays including vector mesons
We may start with processes containing vector mesons. To investigate this we need to take into account few more quantities. We now introduce the overlap between a V meson (V ∈ {ω, ρ 0 , φ}) and the vector current j V µ , i.e.
Here we assume that ρ 0 and ω mesons contain only the light quarks (no hidden strangeness) and the φ meson contains only the strange quarks. This is a fairly good approximation to the real world [13]. For later convenience we also define the γ-V coupling strength With this at hand we can obtain F V from the V → e + e − decay processes. The direct calculation from the Lorentz invariant matrix element yields in the limit m 2 e ≪ M 2 V after averaging and summing over polarizations Above, κ V is the overlap of the electromagnetic current j em µ with the meson-related current j V µ , For a little more details concerning the ω case see [16]. Necessary numerical inputs and results are shown in Table VII.
782.65 (12) [13]. In the upper part of this table we see masses, decay widths and branching ratios of the V → e + e − processes for the vector mesons under consideration. In the lower part we then see the overlaps κV (coefficients in (E5)), decay widths of V → e + e − calculated from values provided in the upper part and finally the resulting coupling strengths FV evaluated according to (E20).
For the following, it is convenient to introduce the η (′) VV correlator (E21) Consequently, the η (′) Vγ * form factor within the VMD model can be written as The additional factors κ η (′) V , which come from the η-η ′ mixing, can be found in Tabs. VIII-IX. Using these form factors we can calculate the two-body decay widths containing one pseudoscalar meson (η (′) ), one vector meson (ω, ρ 0 or φ) and one photon. Depending on the masses of the particles involved, we use one of the following two prescriptions: The results together with a comparison with the experimental data can be found, for the η case, in Table VIII, and for the η ′ case in Table IX. We see that the VMD  [13] for the decays V → ηγ. The √ 2 factor in κ ηφ comes from (E11).
model and the experimental data are in agreement. The last class of the processes containing vector mesons which we will investigate is a counterpart of the previous decays, when the photon in the final state is virtual and turns into the lepton pair. The decay width can then be expressed as a form-factor-dependent integral over the dilepton invariant mass where λ denotes the Källén triangle function defined as After inserting numerical values we get which should be compared with Γ exp φ→ηe + e − = (0.49 ± 0.05) keV. For the above evaluation we have used B(φ → ηe + e − ) = (1.15 ± 0.10) × 10 −4 [13].

b. Electromagnetic transition form factor
Using previously defined meson-specific factors, we can finally define the doubly virtual electromagnetic transi-tion form factor of an eta meson in the quark-flavor basis: (E28) For A(V) above we simply substitute A(ρ 0 ) = A(ω) = ℓ and A(φ) = s . In the VMD case (inserting ansatz (E16)) the form factor becomes .

(E29)
To get the η ′ form factor, it is only necessary to perform the following substitution: As a simple application, we can have a look at the two-photon decay of a pseudoscalar P ∈ {π 0 , η (′) }. The decay with of such a process can be expressed as follows: Of course, for a neutral pion we would have κ π 0 = 1 . In the η (′) case we can write which becomes (cf. (E29) for p 2 = q 2 = 0) cos φ f s ≃ 1.23 (10) .
Note that e.g. κ η ≃ 1 is consistent with experiment, although it differs significantly from a naïve WZW-based calculation [25] for which κ η = 1/ √ 3. The results are shown in Table X. We see that there is a very good agreement between the prediction of the VMD model and the experimental results.
Next, we should investigate the singly virtual transition form factor. One possibility is to compare the decay widths of the η (′) Dalitz decays. These can be written as the following integral:  A simple model we discuss here is however (in the singly virtual mode) not suitable for the Dalitz decays of η ′ due to the unregulated pole which needs to be integrated over. We thus present results only for the η case. These are shown in Table XI. The other way is to directly plot the transition form factor over the respective data. The measurements of the η transition form factor in time-like region are shown in Fig. 7. In the space-like region the results can be found in Fig. 8.
Finally, we can test the transition form factor in its whole power, in the doubly virtual mode. Using the VMD model (E29), we can calculate the decay width of one of the rare processes η → µ + µ − . Note that data for the rare processes η → e + e − and η ′ → ℓ + ℓ − are not available yet; for predictions for these decays within the model under discussion see Table XII in Appendix F. Experimentally, the branching ratio was found to be B(η → µ + µ − ) = (5.8±0.8)×10 −6 , which translates into Γ exp η→µ + µ − = (7.6 ± 1.1) × 10 −6 keV. Within the VMD model we find Γ VMD η→µ + µ − = (6.0±1.4)×10 −6 keV, which is in agreement with the experimental value. Note that we have calculated this value using the convenience of the approach shown in the following appendix. It is similar to the one we have used to calculate the 1γIR contribution to the radiative corrections in Section IV. Finally, let us also mention that for the decays under consideration, more advanced models were developed [31], which could also mimic ππ rescattering effects.
Appendix F: Form factors in P → ℓ + ℓ − decays In this appendix we apply the approach explained in Section IV to the P → ℓ + ℓ − decays. We would like to show how the building block for the matrix element looks like in this case and calculate the coefficients for the specific transition form-factor models.
On account of the Lorentz symmetry and parity conservation, the on-shell matrix element of the P → ℓ + ℓ − process can be written in terms of just one pseudoscalar form factor in the following form: Subsequently, the decay width reads Taking into account only the LO contribution in the QED expansion, we find for the pseudoscalar form factor 4 F P γ * γ * (p 2 , q 2 )λ(M 2 P , p 2 , q 2 ) p 2 q 2 (l 2 − m 2 ℓ ) .
(F3) Here, p = l−q 1 and q = l+q 2 , where q 1 and q 2 are the lepton momenta and λ is the triangle Källén function. For the rational resonance-saturation models, we will use in The normalized singly virtual η electromagnetic transition form factor squared within the VMD-inspired model compared to data: time-like region. The shaded area corresponds to the uncertainty of our fit for the VMD model. Data are taken from [26,27]. agreement with substitution (20) the following definition: h(c 1 , c 2 , M 2 V1 , M 2 V2 )λ(M 2 P , p 2 , q 2 ) (l 2 − m 2 ℓ ) .
After recalling (29) we know that the previous expressions might be obtained in terms of linear combinations of the building blocks P h P →ℓ + ℓ − g(M 2 V1 , M 2 V2 ) . Using the dimensional regularization, the dimensional reduction scheme [22,23] and Passarino-Veltman reduction [24], the explicit result of the necessary loop integration in terms of scalar one-loop integrals reads For completeness, we list the predictions for the branching ratios of the η (′) → ℓ + ℓ − decays in Table XII.
In what follows, we would like to provide some basic examples of the decomposition of the loop integrals containing various models for transition form factors in the case of the rare decay of a neutral pion. Lets start η → e + e − η → µ + µ − η ′ → e + e − η ′ → µ + µ −  Table XII. Branching ratios for the η (′) → ℓ + ℓ − decays within the VMD-inspired model. For instance, the value for the η ′ → µ + µ − process is in good agreement with the value B(η ′ → µ + µ − ) = 1.4(2) × 10 −7 calculated in [7].
with some definitions. The VMD ansatz for the electromagnetic transition form factor of a neutral pion takes a simple form .
(F9) The lowest-meson dominance (LMD) model [32], where also the lowest-lying pseudoscalar multiplet was taken into account, gives the following result: (F10) As the last example we introduce the two-hadron saturation (THS) model proposed in [16], which for the P V V correlator takes into account two meson multiplets in both vector and pseudoscalar channels: .

(F11)
In terms of decomposition (20) we can write for the amplitudes where the coefficients c model 1,i and c model 2,i have the following form: We can find the above-listed constants from projecting on the product of the normalized form factor and the photon propagators: for instance we have c THS 2,2 = lim p 2 ,q 2 →M 2 V 1 F THS π 0 γ * γ * (p 2 , q 2 ) F π 0 γ * γ * (0, 0) (p 2 − M 2 V1 )(q 2 − M 2 V1 ) p 2 q 2 .