A theoretical analysis of the doubly radiative rare decays $\eta^{(\prime)}\to\pi^0\gamma\gamma$ and $\eta^\prime\to\eta\gamma\gamma$

The rare, doubly radiative decays $\eta^{(\prime)}\to\pi^0\gamma\gamma$ and $\eta^\prime\to\eta\gamma\gamma$ are analysed in terms of scalar and vector meson exchange contributions using the frameworks of the Linear Sigma Model and Vector Meson Dominance, respectively. While our predictions for the decay $\eta\to\pi^0\gamma\gamma$ show a satisfactory agreement with the experimental values reported by the A2 and Crystal Ball collaborations, thus supporting the validity of our approach, our estimates for the partner reaction $\eta^{\prime}\to\pi^0\gamma\gamma$ show some dispersion. These are addressed and discussed by performing a first analysis of the associated BESIII data. We also provide a first prediction for the non yet measured decay $\eta^{\prime}\to\eta\gamma\gamma$ that should be taken as a first indication of the possible value of the associated branching ratio. We hope our study to be of interest for present and future experimental analyses of these decays.


I. INTRODUCTION
Measurements of η and η decays have reached unprecedented precision placing new demands on the accuracy of the corresponding theoretical descriptions. Among them, the rare, doubly radiative decay η → π 0 γγ in particular, has attracted much interest due to both the large uncertainties in the experimental data and in the theoretical calculations as we will see below. Also the analysis of the electromagnetic rare decays η → π 0 γγ and η → ηγγ is interesting for several reasons. First, it completes existing calculations on the brother process η → π 0 γγ, which has been studied in many different frameworks, from the seminal works based on Vector Meson Dominance (VMD) [1,2] and Chiral Perturbation Theory ChPT [3] to the more modern treatments based on the unitarization of the chiral amplitudes [4,5] and [6]. At present, while there is only a crude estimate of the branching ratio of the η → π 0 γγ decay [7,8], there is neither a calculation nor a prediction of the branching ratio associated to the decay η → ηγγ. Second, the BESIII collaboration has recently reported the first measurement of the decay η → π 0 γγ [9] making hence the topic of timely interest. Third, the analysis of these decays could help to extract information on the properties of the lowest-lying scalar resonances, in particular of the isovector a 0 (980) from η → π 0 γγ and the isoscalars σ(500) and f 0 (980) from η → ηγγ, thus complementing other analyses based on V → P 0 P 0 γ decays with V = (ρ, ω, φ) and P 0 = (π 0 , η) [10], D and J/ψ decays, central production, etc. (see for brevity the note on scalar mesons in Ref. [11]).
For these reasons, the purpose of this work is to provide a first detailed evaluation of the invariant mass spectra and the integrated branching ratio of the η → π 0 γγ and η → ηγγ decays. As a check of our approach, we also (re)evaluate the η → π 0 γγ decay, for which a measurement of the branching ratio and invariant mass spectrum already exists [12,13].
On the theory side, the η → π 0 γγ process has been seen as a stringent test of the predictive power of ChPT. In this framework, the tree-level contributions at O(p 2 ) and O(p 4 ) vanish because the pseudoscalar mesons involved are neutral. The first non-vanishing contribution comes at O(p 4 ), either from loops involving kaons, largely suppressed due to the kaon masses, or from pion loops, again suppressed since they violate G-parity and are thus proportional to m u − m d [3]. Numerically, one obtains Γ (4) π = 0.84 × 10 −3 eV, Γ (4) K = 2.45 × 10 −3 eV and Γ (4) π+K = 3.89 × 10 −3 eV, respectively, for the π, K and π+K loop contributions at O(p 4 ) to the decay width, two orders of magnitude below the measured width Γ exp η→π 0 γγ = 0.33 ± 0.03 eV [11,13]. The first sizable contribution comes at O(p 6 ), but the coefficients involved are not well determined and one must resort to phenomenological models to fix them. In this sense, for instance, VMD has been used to determine these coefficients by expanding the vector meson propagators and retain the lowest term. This leads, with an equal contribution from ρ and ω, to a decay rate Γ (6) ρ+ω = 0.18 eV, two times smaller than the "all-order" estimate keeping the full vector meson propagator, Γ VMD = 0.31 eV, in reasonable agreement with older VMD estimates [1,2] and Refs. [20,21]. The contributions of the scalar a 0 (980) and tensor a 2 (1320) resonances to the O(p 6 ) chiral coefficients were calculated in the same manner but no "all-order" estimate was given in any case. In addition, contrary to the VMD contribution where the coupling constant appears squared, the signs of the a 0 and a 2 contributions are not unambiguously fixed [3]. At O(p 8 ), a new type of loop effects taking two vertices from the anomalous chiral lagrangian appear. Pion loops are no longer suppressed because there is not G-parity violation in these vertices and the kaonloop suppression doesn't occur. Numerically, these loops give Γ the "all-order" VMD estimate, implies Γ χ+VMD η→π 0 γγ = 0.42 eV. Adding the contributions at O(p 6 ) from the a 0 and a 2 resonances with sign ambiguities, which do not represent an "all-order" estimate of these effects, one can conclude conservatively that Γ χ+VMD+a 0 +a 2 η→π 0 γγ = 0.42 ± 0.20 eV [3]. The further inclusion of C-odd axial-vector resonances raises this value to 0.47 ± 0.20 eV [22] (see also Ref. [23]). Other determinations of the O(p 6 ) low-energy constants in the early and extended Nambu-Jona-Lasinio models lead to 0.11-0.35 eV [24], 0.58 ± 0.30 eV [25], and 0.27 +0. 18 −0.07 eV [26]. A different approach based on quark-box diagrams [27,28] yields values of 0.70 eV and 0.58-0.92 eV, respectively. Finally, in the most recent analyses, the η → π 0 γγ process has been considered within a chiral unitary approach for the meson-meson interaction, thus generating the a 0 resonance and fixing the sign ambiguity on its contribution. This also allows to calculate the loops with one vector meson exchange.
Altogether, they get 0.47 ± 0.10 eV and 0.33 ± 0.08 eV in their 2003 [4] and 2008 [5] works, the difference being due to sizable changes of the radiative decay widths of vector mesons used as input in their calculations, but in very good agreement with Γ exp η→π 0 γγ = 0.33 ± 0.03 eV [13] in any case. As stated before, there is only a rough estimate of η → π 0 γγ [7,8] and no theoretical analyses for η → ηγγ.
Then, in view of the previous discussion and since we are mostly interested in evaluating the η → π 0 γγ and η → ηγγ decays and reveal the possible effects of the lightest scalar resonances on them, our approach is defined as follows. First, we compute the dominant contribution to these processes, that is, the exchange of intermediate vector mesons through the decay chain P 0 → V γ → P 0 γγ. Second, we determine in ChPT the next important contribution, namely, the O(p 4 ) chiral loops. This determination takes into account the input of the pseudoscalar singlet in the chiral amplitudes. Therefore, the large-N c limit of ChPT has been taken and the η 0 is regarded as the ninth pseudo-Goldstone boson of the theory.
To simplify, we work in the isospin and thus only the kaon loops are considered. Later, the chiral-loop prediction will be substituted by a Linear Sigma Model (LσM) calculation where the effects of scalar mesons are included explicitly. In this way, we provide in the given model an "all-order" estimate of these scalar effects, fix the sign ambiguity and test the relevance of including the full scalar meson propagators instead of integrating them out. In Ref. [26], it was shown for the η → π 0 γγ case that the effects from scalar and tensor resonances can account for as much as 20% of the total decay width, depending on the choice of sign for the different contributions. However, on general grounds, one would expect the a 2 (1320) effects to be smaller than the a 0 (980) ones due to the heavier mass involved in the propagators. In the chiral unitary approach, where the a 0 is generated dynamically and the sign ambiguity is then fixed, its resonant effects are seen to enhance the decay width by 2% [4]. This small effect is a consequence of the available phase space, m η − m π 0 413 MeV, for the a 0 (980) to resonate 1 .
For this reason, we will consider only the scalar meson contributions to the processes under analysis and provide an "all-order" estimate of these scalar effects based on a calculation performed in the LσM model. In this way, we will be able, first, to fix the sign ambiguity and, second, to test the relevance of including the full scalar meson propagators, in a given model, instead of integrating them out. However, for the sake of completeness, we start considering the dominant chiral-loop contribution, that is, the contributions containing two This article is organized as follows. In section II, we review the calculation of the chiral loops for η → π 0 γγ and provide the first chiral-loop prediction for the η → π 0 γγ and η → ηγγ decays. In section III, we present the main contribution to these processes, that is, the exchange of an intermediate vector meson through the decay chain P 0 → V γ → P 0 γγ. Later, in section IV, the chiral-loop prediction will be substituted by a Linear Sigma Model (LσM) calculation where the effects of scalar meson resonances are taken into account explicitly. As a check of our approach, in section V we first calculate these two contributions for the decays that concern us, η ( ) → π 0 γγ and η → ηγγ, and then we perform fits to the η ( ) → π 0 γγ experimental data. Finally, our conclusions are presented in section VI.
Preliminary results of these processes have been presented in Refs. [30,31].

II. CHIRAL-LOOP PREDICTION
We start discussing the η → π 0 γγ case. As stated before, the contribution from kaon loops is dominant and the pion loops vanish in the isospin limit. The amplitude is written (cos ϕ P + √ 2 sin ϕ P ) with ϕ P the η-η mixing angle in the quark-flavour basis, resulting from the loop computation 2 . It is important to notice that in the seminal work of Ref. [3] this chiral-loop prediction was computed taking into account the η 8 contribution alone and the mixing angle was fixed to θ P = ϕ P − arctan √ 2 = arcsin(−1/3) −19.5 • . Now, the η 0 contribution is also considered (in the large-N c limit where the pseudoscalar singlet is the ninth pseudo-Goldstone boson) and the dependence on the mixing angle is made explicit.
For the η → π 0 γγ case, the associated amplitude is that of Eq.
Finally, for the η → ηγγ case, two amplitudes contribute, one through a loop of charged kaons, as in the former two cases, and the other through a loop of charged pions, which in this case is not suppressed by G-parity. Again, the corresponding amplitudes are that of Eq. (1), replacing s K → s π and m K + → m π + for the pion loop, with 2 Not to be confused with the four-pseudoscalar scattering amplitude calculated in ChPT at lowest order. and The latter amplitude coincides with that of η → ηπ + π − when computed in the large-N c ChPT at lowest order [32]. Needless to say, the former amplitudes for η → π 0 γγ and η → ηγγ constitute the first chiral-loop predictions of these two processes.

III. VMD PREDICTION
Next to the chiral-loop amplitudes, there are also "all-order" estimates of the corresponding exchange of intermediate vector bosons which are calculated in the framework of VMD.
The full VMD amplitude was seen to produce the dominant contribution to η → π 0 γγ [3], and the same happens, as we see below, for η → π 0 γγ and η → ηγγ. Now, we review the calculation for the η → π 0 γγ case, with some improvements with respect to Ref. [3], and then calculate for the first time the full VMD amplitudes of η → π 0 γγ and η → ηγγ. The calculation of these amplitudes includes the diagrams depicted in Fig. 1.
For η → π 0 γγ, the amplitude is written as where t, u = (P − q 2, P is the four-momentum of the decaying particle, 1,2 and q 1,2 are the polarization and fourmomentum vectors of the final photons, and . For η → π 0 γγ and η → ηγγ, the related amplitudes are Eq. (6) with the replacements g V ηγ g V π 0 γ → g V η γ g V π 0 γ and g V ηγ g V π 0 γ → g V η γ g V ηγ , respectively, and m 2 η → m 2 η . The corresponding couplings are where G = 3g 2 /(4π 2 f π ) and g is the vector-pseudoscalar-pseudoscalar coupling constant of VMD which can be fixed from various ρ and ω decay data, and with ϕ V accounting for the ω-φ mixing angle.
Notice that, when the OZI-rule is applied, that is ω = (uū + dd)/ √ 2 and φ = ss, the expressions of these couplings are reduce to In Ref. [3], the VMD prediction for η → π 0 γγ was calculated assuming equal ρ and ω contributions and without including the decay widths in the propagators. In this case, these approximations were valid since the phase space available prevents the vector mesons to resonate. However, for η → π 0 γγ, the phase space allowed permits these vectors to be onshell and the introduction of their decay widths is mandatory. For this reason, we include, for all the three cases, the decay widths in the vector meson propagators.

IV. LσM PREDICTION
An "all-order" estimate of the scalar meson exchange effects to the processes under study can be achieved in the LσM where the complementarity between this model and ChPT can be used to include the scalar meson poles at the same time as keeping the correct low-energy behavior expected from chiral symmetry. This procedure, that we briefly recapitulate in the following, was applied with success to the related V → P 0 P 0 γ decays [10].
In order to quantify the contribution coming from scalar resonance exchanges we employ an U (3)×U (3) LσM [33][34][35]. In this context, the processes η ( ) → π 0 γγ proceed through kaon loops and by exchanging the a 0 (980) in the s-channel and the κ in the t-and u-channels. The reaction η → ηγγ is a little bit more complex allowing both kaon and pion loops with the σ(600) and the f 0 (980) exchanged in the s-channel and the a 0 (980) in the u-and t-channels.
The full amplitudes in a LσM for the three processes are computed according to with the quantities L(z), s π,K and {a} defined in section II. The participating amplitudes in the LσM turns out to be s, t and u dependent and are expressed in terms of the pion and kaon decay constants f π and f K , the masses of pseudoscalar and scalar mesons involved in the process and the scalar and pseudoscalar mixing angles in the flavour basis, ϕ S and ϕ P , respectively. For our analysis, we follow the procedure outlined in Ref. [10] by one of us to retrieve a consistent full s-dependent amplitude. This mainly consists of replacing the t-and u-channel contributions by the result of subtracting from the chiral-loop amplitude Eqs. (3), (4) and (5) the infinite mass limit of the s-channel scalar contributions. We refer the interested reader to [10] for further details.
In doing so, one finally gets the scalar amplitudes where D S (s) are the S = σ, f 0 and a 0 propagators defined in Appendix A. While a Breit-Wigner propagator is used for the σ, for the f 0 and a 0 a complete one-loop is preferable [10,36]. The required couplings in Eqs. (14) and (15) are given by where ϕ S is the scalar mixing angle in the quark-flavour basis defined as with σ q ≡ 1 √ 2 (uū + dd) and σ s ≡ ss. The couplings g σηη and g f 0 ηη can be written in several different equivalent forms, but we have chosen the ones involving the a 0 mass and the pion decay constant for the sake of clarity. These amplitudes, now only s-dependent, are then incorporated in A LσM η ( ) →π 0 γγ and A LσM η →ηγγ to give the result shown in Eqs. (9), (10) and (11). We can anticipate that, taking into account the scalar meson effects in an explicit way does not provide, as we will numerically see in the following section, a noticeable improvement with respect to the chiral-loop prediction, except for the case of η → ηγγ where the σ contribution turns out to be considerable.

V. PREDICTIONS AND FITS TO EXPERIMENTAL DATA
In this section we present our results. These are obtained using the standard formula for the three body decay rate [11] with the squared amplitude given by MeV [42]) and ϕ P = 38.5(2.9) • [42]. Finally, for illustrative purposes, we also consider a fourth set of parameters that uses the same values of set 3 but with the OZI-rule applied, that is ϕ V = 0 • , and therefore neglecting the contribution of the φ-meson in η ( ) → π 0 γγ (cf. Eq. (8)).
Our results are collected in Table I, where the predictions of chiral loops, the LσM, which replaces the former when scalar meson poles are incorporated, VMD, and the total decay width and branching ratio for the three processes are included. The comparison with the corresponding experimental results, if available, is also displayed. The total decay width is the result of adding the LσM and VMD contributions coherently. The associated errors come from the uncertainties associated to the aforementioned parameters.
For η → π 0 γγ, our calculation agrees with the "all-order estimate" of Ref. [3] and the more involved analysis of Refs. [4,5], thus giving support to our approach as a starting point for the determination of the other two processes. Our predictions, although showing some dispersion, are rather stable against the chosen set of parameters and are seen well in agreement among each other within errors. A graphical account of the two-photon invariant mass distribution obtained with the parameters from set 1 is shown in Fig. 2 and is compared to experimental data. In this plot, the individual contributions to the spectrum, that is, VMD, LσM, and their interference, are also displayed. As seen, our prediction agrees very well with the 2008 Crystal Ball data [12] in Different to η → π 0 γγ, the predicted branching ratios for η → π 0 γγ given in Table I show a significant dispersion going from a branching ratio that appears to be a factor of two bigger than the BESIII experimental value (set 1) to a prediction that agrees well with the measured ones (set 4). We would like to note that the branching ratio moves down towards the experimental value as the numerical value of the coupling |g| and of the ω-φ mixing angle ϕ V diminish. To further investigate this, in Fig. 3 we provide the spectrum in terms of the two-photon invariant mass as obtained with set 1 (solid red curve), set 2 (dotted red curve), set 3 (dashed red curve) and set 4 (solid black curve). Notice that, while the shape and trend of the experimental data is followed by all our predictions, the Decay chiral loops LσM VMD Γ BR th BR exp [11] η → π 0 γγ (eV) set 1 1.24 × 10 − height of the different curves move downwards towards the measured spectrum as |g| and the value ω-φ mixing angle ϕ V diminish. The prediction with the least bias with respect to the experimental data is provided by set 4. This tells us that a scenario with a value of |g| in the lower part of the ballpark of determinations of this coupling, and with a small contribution of the φ-meson, seems to be favored by data. In summary, the η → π 0 γγ decay is rather sensitive to the product of couplings g ωη γ g ωπ 0 γ (cf. Eq. (7)) since the ω contribution prevails with ∼ 80% of the total VMD signal, while the ρ contributes with ∼ 5%. Notice The individual VMD (dotted black curve), LσM (dashed black curve), and their interference (dotdashed black curve), contributions are also shown. Data is taken from Ref. [13] (A2) and Ref. [12] (Crystal Ball). See main text for details. mass measurements, we next run fits allowing this parameter to float. We have considered two different scenarios and the resulting fits are gathered in Table II. The first scenario (named fit 1 and fit I in the table, respectively, for η → π 0 γγ and η → π 0 γγ), employs ϕ P = 40.4(6) • and ϕ V = 3.32(9) • , while the second one (fit 2 and fit II) takes ϕ V = 0 • and allows the η-η mixing angle to float. We would like to notice that while the value of |g| as extracted from a fit to the η → π 0 γγ data, |g| = 4.52(8) (fit 1), tends to lie on the upper side of the ballpark of determinations of this coupling, this is found to be on the lower side when it is fitted to the η → π 0 γγ reaction, |g| We have also tried a simultaneous analysis of both decay channels and the result of the fit tends to the individual fit results of the mode η → π 0 γγ with a χ 2 dof ∼ 2.5. This is not a surprise since we have more data for η → π 0 γγ than for η → π 0 γγ, but the BESIII data of the latter decay is much more precise thus favoring the minimization. No improvement is seen and hence we refrain to show the corresponding fit results.
As seen, the situation for the decay η → π 0 γγ in particular, is not free of controversy Decay Assumptions Fit results VMD Γ BR th BR exp [11] η → π 0 γγ fit 1 ϕ P = 40.4(6) • |g| = 4.52(8) 0.34 (6)  given the discrepancies that we have found between different approaches trying to match the experimental data. However, our predictions are still acceptable and, in views of the dispersion of our results, we can do nothing but to encourage experimental groups to perform new measurements of this decay before claiming than a more refined treatment of this reaction is mandatory. This would help to avoid the controversy that lasted for years in the partner reaction η → π 0 γγ, where large discrepancies between different theoretical approaches and experimental data, and among the latter going from BR = 7.2(1.4) × 10 −4 in 1984 [14] to the recent values BR = 2.52(23) × 10 −4 [13] and BR = 2.21(24)(47) × 10 −4 [12], were reported. Also, a (first) measurement of the crossing-related process γγ → πη would be very welcome in this respect to extract further conclusions.
Finally, for η → ηγγ, the VMD contribution also dominates but the scalar meson effects seem to be sizable, in particular those related with the σ meson. In Fig. 6 [13] and Crystal Ball [12] measurements of the decay η → π 0 γγ as compared to our fit 1 (solid black curve) as presented in Table II. description of the scalar mesons effects in the decay η → ηγγ, one can consider a more sophisticated scalar amplitude for the scattering η η → π + π − , instead of the ones used in this work (cf. Eq. (14)), based on the analysis of Ref. [43] where the decay amplitude counterpart η → ηππ has been successfully proven against experimental data. However, we postpone it for a future analysis when a measurement becomes available.  Table II.

VI. CONCLUSIONS
In this paper, we have analyzed the invariant mass spectra and the branching ratio of the electromagnetic rare decays η → π 0 γγ and η → (π 0 , η)γγ in terms of scalar and vector meson exchange contributions using the frameworks of the linear sigma model and vector meson dominance, respectively. These decays are interesting for several reasons. First, due to the quantum numbers of the pseudoscalar mesons involved, if the decays proceed through resonances, which are mostly of scalar and vector nature, these give access to unravel its properties. Second, the presence of η and η in these reactions allows to study the mixing properties of both mesons. Third, and more generally, for the reasons exposed in section I, these decays allow to test Chiral Perturbation Theory (ChPT), and its natural extensions, processes in Table I, and provided a graphical account of the two-photon invariant mass distribution for η → π 0 γγ and η → π 0 γγ in Figs. 2 and 3, respectively, using different sets of representative values for the coupling |g| and for the η-η mixing angle ϕ P , and also for the vector ω-φ mixing angle ϕ V in some cases. It is seen that, while for the η ( ) → π 0 γγ decays the vector meson exchanges vastly dominate over the scalar contributions, for the η → ηγγ one, the VMD contribution also dominates but the scalar meson effects seem to be sizable. Our estimates for η → π 0 γγ, although are seen on the lower side with respect to the experimental measurements, agree well with the reported value. The reaction η → π 0 γγ deserves more attention since our predictions show some dispersion. This decay is vastly dominated by the production of the ω-meson, and the reason for the differing results resides mostly in the corresponding product of couplings g ωη γ g ωπ 0 γ which, in turn, depends strongly on |g| and in the ω-φ mixing angle ϕ V , and also in the η-η mixing angle ϕ P to less extent (cf. Eq. (7)).
These differences can be used to asses systematic uncertainties associated to our predictions but also to encourage the experimental collaborations to perform new measurements of this decay to see if more precise data demand a more refined theoretical treatment. Regarding the decay η → ηγγ, we have presented the two-photon decay spectrum, Fig. 6, and provided the very first estimate for the branching ratio which we consider is large enough to be measured in the near future by several experimental collaborations, such as BESIII, where measurements of various η and η decays are included in their physics programs [44]. This would help discriminate among different scalar models and see whether our predictions are corroborated.
On the second stage, in view of the sensitivity of the decay η → π 0 γγ to the numerical value of the coupling |g|, we have performed fits to the η → π 0 γγ and η → π 0 γγ experimental data in order to optimize its value to the measured spectra. The resulting fits are presented in Table II and a graphical account provided in Figs. 4 and 5 for η → π 0 γγ and η → π 0 γγ, respectively. We have found that while the value of |g| as determined from a fit to the η → π 0 γγ data from the A2 and Crystal Ball collaborations tends to lie on the upper side of the ballpark of determinations of this coupling, it is found to be on the lower side when it is extracted from the η → π 0 γγ BESIII experimental data. We have also seen that, for the η → π 0 γγ reaction, data seems to favor a scenario with a small contribution of the φ-meson.
As a final concluding remark, we would like to encourage experimental groups to measure these decays once again to see whether our predictions are corroborated or, on the contrary, more refined theoretical descriptions are demanded.