A theoretical analysis of the semileptonic decays $\eta ^{(\prime )}\rightarrow \pi ^0l^+l^-$ and $\eta ^\prime \rightarrow \eta l^+l^-$

Abstract

A complete theoretical analysis of the C- conserving semileptonic decays $\eta ^{(\prime )}\rightarrow \pi ^0l^+l^-$ and $\eta ^\prime \rightarrow \eta l^+l^-$ ($l=e$ or $\mu$) is carried out within the framework of the Vector Meson Dominance (VMD) model. An existing phenomenological model is used to parametrise the VMD coupling constants and the associated numerical values are obtained from an optimisation fit to $V\rightarrow P\gamma$ and $P\rightarrow V\gamma$ radiative decays ($V=\rho ^0$, $\omega$, $\phi$ and $P=\pi ^0$, $\eta$, $\eta ^{\prime }$). The decay widths and dilepton energy spectra for the two $\eta \rightarrow \pi ^0l^+l^-$ processes obtained using this approach are compared and found to be in good agreement with other results available in the published literature. Theoretical predictions for the four $\eta ^{\prime }\rightarrow \pi ^0l^+l^-$ and $\eta ^\prime \rightarrow \eta l^+l^-$ decay widths and dilepton energy spectra are calculated and presented for the first time in this work.

1 Introduction

The electromagnetic and strong interactions conserve parity (P) and charge conjugation (C) within the well-established and well-tested Standard Model of particle physics (SM). In this framework, the $\eta$ and $\eta ^{\prime }$ pseudoscalar mesons are specially suited for the study of rare decay processes, for instance, in search of C, P and CP violations, as these mesons are C and P eigenstates of the electromagnetic and strong interactions [1].

Specifically, the semileptonic decays $\eta ^{(\prime )}\rightarrow \pi ^0l^+l^-$ and $\eta ^\prime \rightarrow \eta l^+l^-$ ($l=e$ or $\mu$) are of special interest given that they can be used as fine probes to assess if new physics beyond the Standard Model (BSM) is at play. This is because any contribution from BSM physics ought to be relatively small¹ and the above decay processes only get a  contribution from the SM through the C-conserving exchange of two photons that is highly suppressed, as there is no contribution at tree-level but only corrections at one-loop and higher orders. This small SM contribution would presumably be of the same order of magnitude as that of physics BSM, which, in turn, means that the $\eta$-$\eta ^{\prime }$ phenomenology might play an interesting role and be an excellent arena for stress testing SM predictions [1, 2]. As an example, the $\eta ^{(\prime )}\rightarrow \pi ^0l^+l^-$ and $\eta ^\prime \rightarrow \eta l^+l^-$ decays could be mediated by a single intermediate virtual photon, but this would entail that the electromagnetic interactions violate C-invariance (e.g. [3, 4]) and, therefore, would represent a departure from the SM.

Early theoretical studies of semileptonic decays of pseudoscalar mesons date back to the late 1960s. A very significant contribution was made by Cheng in Ref. [5] where he analysed the $\eta \rightarrow \pi ^0e^+e^-$ decay mediated by a C-conserving, two-photon intermediate state within the Vector Meson Dominance (VMD) framework. By setting the electron mass to $m_e=0$ and neglecting in the numerator of the amplitude terms that were second or higher order in the electron or positron 4-momenta, he found theoretical estimations for the decay width eV, the relative branching ratio , as well as the associated decay energy spectrum. This was an enormous endeavour given the very limited access to computer algebra systems at the time. For this reason, a number of strong assumptions had to be made, as pointed out above, which may have had an undesired effect on the accuracy of Cheng’s estimates. A different approach was followed by Smith [6] also in the late 1960s, whereby an S-wave $\eta \pi ^0\gamma \gamma$ coupling and unitary bounds² were used for the calculation of the C-conserving modes associated to both $\eta \rightarrow \pi ^0l^+l^-$ decay processes. By neglecting p-wave contributing terms to simplify the calculations and noting that the unknown $\eta \pi ^0\gamma \gamma$ coupling constant cancels out when calculating relative branching ratios, Smith was able to find and , after estimating the real part of the matrix element from a single dispersion relation and employing a cut-off $\varLambda = 2m_{\eta }$. The calculation of the latter ratio was possible due to the fact that Smith did not approximate the lepton mass to zero.

Ng et al. [8] also found in the early 1990s lower limits for the decay widths of the two $\eta \rightarrow \pi ^0l^+l^-$ processes by making use of unitary bounds and the decay chain $\eta \rightarrow \pi ^0\gamma \gamma \rightarrow \pi ^0l^+l^-$. The transition form factors associated to the $\eta \rightarrow \pi ^0\gamma \gamma$ decay, which are required to perform the above calculation, were obtained using the VMD model supplemented by the exchange of an $a_0$ scalar meson. The lower bounds that they found are eV and eV, making use of VMD only. By adding the $a_0$ exchange³ to the latter process, they obtained eV and eV for a constructive and destructive interference, respectively. The real parts of the amplitudes were estimated by means of a cut-off dispersive relation and the authors argued that the expected dispersive contribution should be no larger than 30% of the absorptive one. Only a few months later, Ng and Peters provided in Ref. [9] new estimations for the unitary bounds of the $\eta \rightarrow \pi ^0l^+l^-$ decay widths. This new contribution was two-fold; on one hand, they calculated the $\eta \rightarrow \pi ^0\gamma \gamma$ decay width within a constituent quark model framework; on the other hand, they recalculated the VMD transition form factors from Ref. [8] by performing a Taylor expansion and keeping terms linear in $M_{\eta }^2/M_V^2$, $x_1$ and $x_2$ ($x_i\equiv P_{\eta }\cdot q_{\gamma _i}/M_{\eta }^2$), which had been neglected in their previous work. Their new findings were: (i) eV and eV for a constituent quark mass $m=330$ MeV$/c^2$; and (ii) eV and eV. It is important to highlight that their estimations using the quark-box mechanism were strongly dependent on the specific constituent quark mass selected, especially for the electron mode.

Fig. 1

Feynman diagrams contributing to the C-conserving semileptonic decays $\eta ^{(\prime )}\rightarrow \pi ^0l^+l^-$ and $\eta ^{\prime }\rightarrow \eta l^+l^-$ ($l=e$ or $\mu$). Note that $q=p_{+}+p_{-}$ and

On the experimental front, new upper limits have recently been established by the WASA-at-COSY collaboration for the $\eta \rightarrow \pi ^0e^+e^-$ decay width [10]. This is a very useful contribution, as the previous available empirical measurements date back to the 1970s which provided an upper limit for the relative branching ratio of the above process that was many orders of magnitude larger than the corresponding theoretical estimations at the time. In particular, Adlarson et al. [10] found from the analysis of a total of $3\times 10^7$ events of the reaction $\text {pd}\rightarrow \!^3\!\text {He}\eta$, with a recorded excess energy of $Q=59.8$ MeV, that the results are consistent with no C-violating single-photon intermediate state event being recorded. Based on their analysis, the new upper limits and (CL $=90\%$) have been established for the C-violating $\eta \rightarrow \pi ^0\gamma ^*\rightarrow \pi ^0e^+e^-$ decay. In addition, the WASA-at-COSY Collaboration is currently analysing additional data from the $\text {pp}\rightarrow \text {pp}\eta$ reaction collected over three periods in 2008, 2010 and 2012 which should put more stringent upper limits on the $\eta \rightarrow \pi ^0e^+e^-$ branching ratio. The experimental state of play is expected to be further improved in the near future with the advent of new experiments such as the REDTOP Collaboration, which will focus on rare decays of the $\eta$ and $\eta ^{\prime }$ mesons, providing increased sensitivity in the search for violations of SM symmetries by several orders of magnitude beyond the current experimental state of the art [11].

The present work is structured as follows: In Sect. 2, we present the detailed calculations for the decay widths associated to the six $\eta ^{(\prime )}\rightarrow \pi ^0l^+l^-$ and $\eta ^{\prime }\rightarrow \eta l^+l^-$ processes. In Sect. 3, numerical results from theory for the decay widths and the corresponding dilepton energy spectra are presented and discussed for the six reactions. Some final remarks and conclusions are given in Sect. 4.

2 Calculations of $\pmb {\eta ^{(\prime )}\rightarrow \pi ^0 l^+ l^-}$ and $\pmb {\eta ^{\prime }\rightarrow \eta l^+ l^-}$

The calculations in this work assume that the $\eta ^{(\prime )}\rightarrow \pi ^0l^+l^-$ and $\eta ^{\prime }\rightarrow \eta l^+l^-$ decays processes are dominated by the exchange of vector resonances⁴ [5]; that is, they proceed through the C-conserving virtual transition $\eta ^{(\prime )}\rightarrow V\gamma ^*$ (with $V = \rho ^0$, $\omega$ or $\phi$), followed by $V\rightarrow \pi ^0\gamma ^*$ (or $V\rightarrow \eta \gamma ^*$) and $2\gamma ^{*}\rightarrow l^+l^-$ (see Fig. 1 for details).⁵

In order to perform the calculations, one first needs to select an effective vertex that contains the appropriate interacting terms. The $VP\gamma$ interaction amplitude consistent with Lorentz, P, C and electromagnetic gauge invariance can be written as [13]

$$\begin{aligned} {\mathscr {M}}(V\rightarrow P\gamma ) = g_{VP\gamma }\,\epsilon _{\mu \nu \alpha \beta } \epsilon ^\mu _{(V)}p_V^\nu \epsilon ^{*\alpha }_{(\gamma )}q^\beta {\hat{F}}_{VP\gamma }(q^2)\ , \end{aligned}$$

(1)

where $g_{VP\gamma }$ is the coupling constant for the $VP\gamma$ transition involving on-shell photons, $\epsilon _{\mu \nu \alpha \beta }$ is the totally antisymmetric Levi-Civita tensor, $\epsilon _{(V)}$ and $p_V$ are the polarisation and 4-momentum vectors of the initial V, $\epsilon _{(\gamma )}^*$ and q are the corresponding ones for the final $\gamma$, and ${\hat{F}}_{VP\gamma }(q^2)\equiv F_{VP\gamma }(q^2)/F_{VP\gamma }(0)$ is a normalised form factor to account for off-shell photons mediating the transition.⁶ In addition to this, the usual QED vertex is used to describe the subsequent $2\gamma ^{*}\rightarrow l^+l^-$ transition. Accordingly, there are six diagrams contributing to each one of the six semileptonic decay processes and the corresponding Feynman diagrams are shown in Fig. 1.

The invariant decay amplitude in momentum space can, therefore, be written as follows

(2)

where $q=p_{+}+p_{-}$ is the sum of lepton-antilepton pair 4-momenta, e is the electron charge, and $g_{V\eta ^{(\prime )}\gamma }$ and $g_{V\pi ^0(\eta )\gamma }$ are the corresponding coupling constants in Eq. (1). Noting that the Levi-Civita tensors are antisymmetric under the substitutions and , whilst the products of loop momenta $k^{\mu }k^{\alpha }$ and $k^{\rho }k^{\delta }$ are symmetric under these substitutions, one finds that the terms in Eq. (2) containing these combinations vanish and that the superficial degree of divergence for the loop integrals of the two diagrams in Fig. 1 is $-1$. Accordingly, both diagrams are convergent individually.

The numerator of ${\mathcal {M}}$ can be simplified using the usual Dirac algebra manipulations and the equations of motion. For these calculations, the mass of the leptons are not approximated to zero, as we are interested in both the electron and muon modes for the six decay processes; as a result, the task of manipulating and simplifying the algebraic expressions would be daunting should computer algebra packages not be available. In the present work, use of the Mathematica package FeynCalc 9.2.0 [14, 15] is made for this purpose.

Let us now proceed to calculate the loop integral. As usual, one first introduces the Feynman parametrisation and completes the square in the new denominators $\varDelta _{iV}$ ($i=1,2$ and ) by shifting to a new loop momentum variable $\ell$ [16]. Hence, the denominators become

$$\begin{aligned} \varDelta _{1V}&= 2yz(P\cdot q)+2xy(p_{+}\cdot q)+(y-1)yq^2\nonumber \\&\quad +2x z(P\cdot p_{+})+x^2m_{l}^2+z\big [(z-1)m_{\eta ^{(\prime )}}^2 \nonumber \\&\quad +m_V(m_V-i\varGamma _V)\big ]\ , \nonumber \\ \varDelta _{2V}&\equiv \varDelta _{1V} \ \ \mathrm {with} \ \ p_{+}\leftrightarrow p_{-}\ . \end{aligned}$$

(3)

Rewriting the numerators of the Feynman diagrams 1 and 2 (i.e. t-channel and u-channel diagrams, respectively, in Fig. 1) in terms of the new momentum variable $\ell$, one finds

(4)

where the explicit expressions for the parameters $A_i$, $B_i$, $C_i$ and $D_i$ ($i=1,2$) are provided in Appendix A. Finally, we perform a Wick rotation and change to four-dimensional spherical coordinates [16, 17] to carry out the momentum integral. The following expressions for the amplitudes of the Feynman diagrams are found

(5)

where $m_l$ is the corresponding lepton mass, and the parameters , in Eq. (5) are defined as

(6)

with x, y and z being the Feynman integration parameters. Therefore, the full amplitude can now be expressed as

(7)

where $\varOmega$ and $\Sigma$ are defined as follows

$$\begin{aligned} \begin{aligned} \varOmega= & {} \sum _{V=\rho ^0,\omega ,\phi }\alpha _V + \sigma _V\ , \\ \Sigma= & {} \sum _{V=\rho ^0,\omega ,\phi }\beta _V + \tau _V\ , \end{aligned} \end{aligned}$$

(8)

and the unpolarised squared amplitude is

$$\begin{aligned} \overline{\vert {\mathcal {M}}\vert ^2}&= 4\Big \{2(P\cdot p_{+})(P\cdot p_{-})-m_{\eta ^{(\prime )}}^2\big [(p_{+}\cdot p_{-})+m_l^2\big ]\Big \}\nonumber \\&\quad \times \mathrm {Abs}(\varOmega )^2+8m_l^2\big [(P\cdot p_{+})-(P\cdot p_{-})\big ]\mathrm {Re}(\varOmega \Sigma ^{*})\nonumber \\&\quad + 4m_l^2\big [(p_{+}\cdot p_{-})-m_l^2\big ]\mathrm {Abs}(\Sigma )^2\ . \end{aligned}$$

(9)

Finally, the differential decay rate for a three-body decay can be written as [18]

$$\begin{aligned} d\varGamma = \frac{1}{(2\pi )^3}\frac{1}{32\ \!m^3_{\eta ^{(\prime )}}} \overline{\vert {\mathcal {M}}\vert ^2}\ dm_{l^+l^-}^2 dm_{l^-\pi ^0(\eta )}^2\ , \end{aligned}$$

(10)

where $m_{ij}^2=(p_i+p_j)^2$.

3 Theoretical results

Making use of the theoretical expressions that have been presented in Sect. 2, one can find numerical predictions for the decay widths of the $\eta ^{(\prime )}\rightarrow \pi ^0l^+l^-$ and $\eta ^{\prime }\rightarrow \eta l^+l^-$ decay processes, as well as their associated dilepton energy spectra. Both, the integral over the Feynman parameters as well as the integral over phase space, must be carried out numerically, as algebraic expressions cannot be obtained. In addition, the numerical integrals over the Feynman parameters are to be performed using adaptive Monte Carlo methods [19]; this is driven by the complexity of the expressions to be integrated and their multidimensional nature.

In the conventional VMD model, pseudoscalar mesons do not couple directly to photons but through the exchange of intermediate vectors. Thus, in this framework, a particular $VP\gamma$ coupling constant times its normalised form factor, cf. Eq. (1), is given by⁷

$$\begin{aligned} g_{VP\gamma }\,{\hat{F}}_{VP\gamma }(q^2)= \sum _{V^\prime }\frac{g_{VV^\prime P}\,g_{V^\prime \gamma }}{M_{V^\prime }^2-q^2}\ , \end{aligned}$$

(11)

where $g_{VV^\prime P}$ are the vector-vector-pseudoscalar couplings, $g_{V^\prime \gamma }$ the vector-photon conversion couplings, and $M_{V^\prime }$ the intermediate vector masses. In the SU(3)-flavour symmetry and OZI-rule respecting limits, one could express all the $g_{VP\gamma }$ in terms of a single coupling constant and SU(3)-group factors [20]. However, to account for the unavoidable SU(3)-flavour symmetry-breaking and OZI-rule violating effects, we make use of the simple, yet powerful, phenomenological quark-based model first presented in Ref. [21], which was developed to describe $V\rightarrow P\gamma$ and $P\rightarrow V\gamma$ radiative decays. According to this model, the decay couplings can be expressed as

$$\begin{aligned} \begin{aligned} g_{\rho ^0\pi ^0\gamma }&= \frac{1}{3}g\ ,\\ g_{\rho ^0\eta \gamma }&= gz_{\text {NS}}\cos {\phi _P}\ ,\\ g_{\rho ^0\eta ^{\prime }\gamma }&= gz_{\text {NS}}\sin {\phi _P}\ ,\\ g_{\omega \pi ^0 \gamma }&= g\cos {\phi _V}\ ,\\ g_{\omega \eta \gamma }&= \frac{1}{3}g\Big (z_{\text {NS}}\cos {\phi _P}\cos {\phi _V} - 2 \frac{{\overline{m}}}{m_s}z_{\text {S}}\sin {\phi _P}\sin {\phi _V}\Big )\ ,\\ g_{\omega \eta ^{\prime }\gamma }&= \frac{1}{3}g\Big (z_{\text {NS}}\sin {\phi _P}\cos {\phi _V} + 2 \frac{{\overline{m}}}{m_s}z_{\text {S}}\cos {\phi _P}\sin {\phi _V}\Big )\ ,\\ g_{\phi \pi ^0\gamma }&= g\sin {\phi _V}\ ,\\ g_{\phi \eta \gamma }&= \frac{1}{3}g\Big (z_{\text {NS}}\cos {\phi _P}\sin {\phi _V} + 2 \frac{{\overline{m}}}{m_s}z_{\text {S}}\sin {\phi _P}\cos {\phi _V}\Big )\ ,\\ g_{\phi \eta ^{\prime }\gamma }&= \frac{1}{3}g\Big (z_{\text {NS}}\sin {\phi _P}\sin {\phi _V} - 2 \frac{{\overline{m}}}{m_s}z_{\text {S}}\cos {\phi _P}\cos {\phi _V}\Big )\ , \end{aligned}\nonumber \\ \end{aligned}$$

(12)

where g is a generic electromagnetic constant, $\phi _P$ is the pseudoscalar $\eta$-$\eta ^{\prime }$ mixing angle in the quark-flavour basis, $\phi _V$ is the vector $\omega$-$\phi$ mixing angle in the same basis, ${\overline{m}}/m_s$ is the quotient of constituent quark masses, and $z_{\text {NS}}$ and $z_{\text {S}}$ are the non-strange and strange multiplicative factors accounting for the relative meson wavefunction overlaps [21, 22]. By performing an optimisation fit to the most up-to-date $VP\gamma$ experimental data [18], one can find values for the above parameters

$$\begin{aligned}&g = 0.70 \pm 0.01 \ \text {GeV}^{-1}, \quad z_{\text {S}}{\overline{m}}/m_s = 0.65 \pm 0.01 \ ,\nonumber \\&\quad \phi _{P} = (41.4 \pm 0.5)^\circ , \quad \phi _{V} = (3.3 \pm 0.1)^\circ \ ,\nonumber \\&\quad z_{\text {NS}} = 0.83 \pm 0.02 \ . \end{aligned}$$

(13)

Given the very wide decay width associated to the $\rho ^0$ resonance, which, in turn, is linked to its very short lifetime, the use of the usual Breit-Wigner approximation for the $\rho ^0$ propagator is not justified. Instead, an energy-dependent approximation for the vector propagator must be used. This can be written for the t-channel process (cf. t-channel diagram in Fig. 1) as follows

$$\begin{aligned} \varGamma _{\rho ^0}(t) = \varGamma _{\rho ^0}\times \Bigg (\frac{t-4m_{\pi ^{\pm }}^2}{m_{\rho ^0}^2 -4m_{\pi ^{\pm }}^2}\Bigg )^{\frac{3}{2}}\times \theta (t-4m_{\pi ^{\pm }}^2)\ , \end{aligned}$$

(14)

where $\theta (x)$ is the Heaviside step function. Likewise, for the u-channel process (cf. u-channel diagram in Fig. 1), one only needs to substitute $t\leftrightarrow u$ in Eq. (14). The energy dependent propagator is not needed, though, for the $\omega$ and $\phi$ resonances, as their associated decay widths are narrow and, therefore, use of the usual Breit-Wigner approximation suffices.

Using the most recent empirical data for the meson masses and total decay widths from Ref. [18], together with all the above considerations, one arrives at the decay width results shown in Table 1 for the six $\eta ^{(\prime )}\rightarrow \pi ^0l^+l^-$ and $\eta ^{\prime }\rightarrow \eta l^+l^-$ processes. The total decay widths associated to the electron modes turn out to be larger than the ones corresponding to the muon modes despite the second and third terms in the unpolarised squared amplitude (cf. Eq. (9)) being helicity suppressed for the electron modes. This suppression, though, does not overcome the phase space suppression for the muon modes, yielding and .

Table 1 Decay widths and branching ratios for the six C-conserving decays $\eta ^{(\prime )}\rightarrow \pi ^0l^+l^-$ and $\eta ^{\prime }\rightarrow \eta l^+l^-$ ($l=e$ or $\mu$). First error is experimental, second is down to numerical integration and third is due to model dependency

Let us now look at the contributions from the different vector meson exchanges to the total decay widths. For the first decay, i.e. $\eta \rightarrow \pi ^0e^+e^-$, we find that the contribution from the $\rho ^0$ exchange is $\sim 25\%$, the contribution from the $\omega$ is $\sim 22\%$, whilst the one from the $\phi$ is negligible, i.e. $\sim 0\%$. The interference between the $\rho ^0$ and the $\omega$ is constructive, accounting for the $\sim 47\%$; similarly, the interference between the $\rho ^0$ and the $\omega$ with the $\phi$ is constructive and about $\sim 6\%$. The contributions to the second decay, i.e. $\eta \rightarrow \pi ^0\mu ^+\mu ^-$, are $\sim 26\%$, $\sim 22\%$ and $\sim 0\%$ from the $\rho ^0$, $\omega$, and $\phi$ exchanges, respectively. As before, the interference between the $\rho ^0$ and the $\omega$ is constructive, weighing $\sim 49\%$, and the interference between the $\rho ^0$ and the $\omega$ with the $\phi$ is constructive and accounts for approximately the $\sim 3\%$. For the third decay, i.e. $\eta ^{\prime }\rightarrow \pi ^0e^+e^-$, the contributions from the $\rho ^0$, $\omega$ and $\phi$ turn out to be $\sim 17\%$, $\sim 36\%$ and $\sim 0\%$, respectively; the interference between the $\rho ^0$ and the $\omega$ exchanges is constructive and accounts for the $\sim 51\%$, whilst the interference between the $\rho ^0$ and $\omega$ with the $\phi$ is destructive and weighs approximately $\sim 4\%$. The contributions to the fourth decay, i.e. $\eta ^{\prime }\rightarrow \pi ^0\mu ^+\mu ^-$, from the $\rho ^0$, $\omega$ and $\phi$ exchanges are $\sim 21\%$, $\sim 32\%$ and $\sim 0\%$, respectively. The interference between the $\rho ^0$ and the $\omega$ is constructive, representing a $\sim 52\%$ contribution, whilst the interference between the $\rho ^0$ and the $\omega$ with the $\phi$ is destructive and accounts for the $\sim 5\%$. The fifth decay, i.e. $\eta ^{\prime }\rightarrow \eta e^+e^-$, gets contributions from the exchange of $\rho ^0$, $\omega$ and $\phi$ resonances of approximately $\sim 79\%$, $\sim 1\%$ and $\sim 2\%$, respectively; the interference between the $\rho ^0$ and the $\omega$ is constructive weighing $\sim 26\%$, and the interference between the $\rho ^0$ and the $\omega$ with the $\phi$ is destructive and contributes with roughly the $\sim 8\%$. Finally, for the sixth decay, i.e. $\eta ^{\prime }\rightarrow \eta \mu ^+\mu ^-$, we find that the contribution from the $\rho ^0$ exchange is $\sim 93\%$, the contribution from the $\omega$ is $\sim 2\%$ and the one from the $\phi$ is $\sim 3\%$; the interference between the $\rho ^0$ and the $\omega$ is constructive and accounts for the $\sim 26\%$, whilst the interference between the $\rho ^0$ and the $\omega$ with the $\phi$ is destructive weighing close to $\sim 24\%$. The tiny contribution from the $\phi$ exchange to the decay widths of the six processes is explained by the relatively small product of VMD $VP\gamma$ coupling constants. Likewise, the comparatively minute contribution from the $\omega$ exchange to the decay widths of the last two reactions is down to the significantly smaller product of coupling constants, if compared to that of the $\rho ^0$ exchange.

In order to assess a systematic error associated to the model dependency of our predictions, we repeat all the above calculations in the context of Resonance Chiral Theory (RChT). In this framework, the $VP \gamma$ effective vertex is made of two contributions, a local $VP\gamma$ vertex weighted by a coupling constant, $h_V$, and a non-local one built from the exchange of an intermediate vector which, again, is weighted by a second coupling constant, $\sigma _V$, times the vector-photon conversion factor $f_V$. For a given $VP\gamma$ transition, this effective vertex can be written in the SU(3)-flavour symmetry limit as [13]

$$\begin{aligned} g_{VP\gamma }\,{\hat{F}}_{VP\gamma }(q^2)= C_{VP\gamma }\frac{4\sqrt{2}\,h_V}{f_\pi } \left( 1+\frac{\sigma _V f_V}{\sqrt{2}\,h_V} \frac{q^2}{M_{V^\prime }^2-q^2}\right) \ ,\nonumber \\ \end{aligned}$$

(15)

where $C_{VP\gamma }$ are SU(3)-group factors and, depending on the process, the exchanged vector is or is not the same as the initial vector (see Refs. [13, 23] for each particular case). To fix the $VP\gamma$ couplings in this second approach, we make use of the extended Nambu–Jona-Lasinio (ENJL) model, where $h_V$ is found to be $h_V=0.033$ [13]. The VVP coupling $\sigma _V$ obtained using the ENJL model turns out to be $\sigma _V=0.28$. However, $\sigma _V$ can also be obtained from the analysis of the dilepton mass spectrum in $\omega \rightarrow \pi ^0\mu ^+\mu ^-$ decays, where one finds $\sigma _V \approx 0.58$ [24]. Due to the fact that $\sigma _V$ is poorly known and the dispersion of the above estimations is large, we do not consider the $q^2$ dependence of the form factors in the subsequent calculations. An alternative model to fix $g_{VP\gamma }$, the normalisation of the form factors, is the Hidden Gauge Symmetry (HGS) model [25], where the vector mesons are considered as gauge bosons of a hidden symmetry. Within this model, a $VP\gamma$ transition proceeds uniquely through the exchange of intermediate vector mesons. In this sense, it is equivalent to the conventional VMD model with the relevant exception of including direct $\gamma P^3$ terms (P being a pseudoscalar meson), which are forbidden in VMD [26]. Due to this similarity, we will not make use of the HGS model to assess the systematic model error and refer the interested reader to Ref. [20] for a detailed calculation of the $g_{VP\gamma }$ couplings in this model.

Next, our results for the semileptonic decays $\eta ^{(\prime )}\rightarrow \pi ^0l^+l^-$ and $\eta ^\prime \rightarrow \eta l^+l^-$ in the conventional VMD framework using the $VP\gamma$ couplings from the phenomenological quark-based model in Eq. (12) are discussed and, if available, compared with previous literature. These predictions include a first experimental error ascribed to the propagation of the parametric errors in Eq. (13), a second error down to the numerical integration, and a third systematic error associated to the model dependence of our approach. The latter is calculated as the absolute difference between the predicted central values obtained from the VMD and RChT frameworks (cf. Table 1).

Fig. 2

Dilepton energy spectra corresponding to the six C-conserving semileptonic decay processes $\eta ^{(\prime )}\rightarrow \pi ^0 l^+l^-$ and $\eta ^{\prime }\rightarrow \eta l^+l^-$ ($l=e$ or $\mu$) as a function of the dilepton invariant mass $q^2$

Our prediction for the decay width eV is about an order of magnitude smaller than the one provided by Cheng in Ref. [5] (cf. Sect. 1), i.e.  eV; however, by plugging into our expressions the couplings that Cheng used in his work, we find a decay width eV, which is less than a factor of two larger than Cheng’s result, and the difference might be down to the simplifications that he had to carry out in his calculations, as well as the more sophisticated propagators that have been employed in the present work.⁸ In addition, from our calculations one can also get a prediction for the ratio of branching ratios⁹, which is very approximate to Cheng’s model independent estimation of , but more than two orders of magnitude larger than Smith’s result¹⁰ [6]; as well as this, for the muon mode we find the relative branching ratio¹¹, which is roughly an order of magnitude smaller than Smith’s estimation [6]. Moreover, our decay widths for both the $\eta \rightarrow \pi ^0e^+e^-$ and $\eta \rightarrow \pi ^0\mu ^+\mu ^-$ processes are consistent with the lower bounds provided by Ng et al. in Ref. [8], i.e. eV and eV making use of VMD, and eV and eV using VMD supplemented by the exchange of an $a_0$ scalar meson. Using the quark-box diagram and a constituent quark mass $m=330$ MeV$/c^2$, Ng et al. provided in Ref. [9] an estimation for the electron mode, eV, which is in accordance with our result, and an estimate for the muon mode, eV, which in this case is incompatible with our calculation.¹² Additionally, Ng et al. also presented in Ref. [9] a recalculation of their previous VMD results from Ref. [8], yielding eV and eV, which are consistent with our results if one considers the associated errors. Our decay width calculations for the other four processes, i.e. $\eta ^{\prime }\rightarrow \pi ^0l^+l^-$ and $\eta ^{\prime }\rightarrow \eta l^+l^-$, cannot be compared with any previously published theoretical results, as these decays have been calculated, to the best of our knowledge, for the first time in the present work. Likewise, comparison with the most up-to-date empirical data provides limited value given that the corresponding current experimental upper bounds, though consistent with our theoretical predictions, are many orders of magnitude larger (cf. Table 1).

Finally, theoretical results for the dilepton energy spectra of the six C-conserving semileptonic decays are presented in Fig. 2. The energy spectra for the three electron modes, which are displayed in Fig. 2a, c, e, take off as the dilepton invariant mass $q^2$ approaches zero. This is in line with Cheng’s and Ng et al.’s energy spectra for the $\eta \rightarrow \pi ^0 e^+e^-$ (Refs. [5] and [8], respectively), which exhibit the same behaviour at low $q^2$. It appears as though the electron modes prefer to proceed through the emission of (relativistic) collinear electron-positron pairs (i.e. $\theta _{e^+e^-}\simeq 0$, where $\theta _{e^+e^-}$ is the electron-positron angle). The reason for this can easily be understood from dynamics¹³ if one assumes the electron and positron to be massless, $m_e\approx 0$; then, by inspection of Eq. (9), one can determine that the unpolarised squared amplitude is maximised when $q^2\rightarrow 0$, which occurs when $\theta _{e^+e^-}\simeq 0$¹⁴. Physically, it may be explained to some extent by the fact that the diphoton invariant spectra for the three $\eta ^{(\prime )}\rightarrow \pi ^0\gamma \gamma$ and $\eta ^{\prime }\rightarrow \eta \gamma \gamma$ peak at low $m_{\gamma \gamma }^2$ (cf., e.g., Ref. [12] and references therein). On the other hand, the dilepton energy spectra for the muon modes, shown in Fig. 2b, d, f, are bell-shaped, which is driven by the kinematics of the processes. This, once more, seems to be consistent with Ng et al.’s [8] energy spectrum for the $\eta \rightarrow \pi ^0 \mu ^+\mu ^-$. It is interesting to observe that the energy spectra for the $\eta \rightarrow \pi ^0\mu ^+\mu ^-$ and $\eta ^{\prime }\rightarrow \pi ^0 \mu ^+\mu ^-$ are skewed to the left (i.e. small values for $\theta _{\mu ^+\mu ^-}$ are favoured, where $\theta _{\mu ^+\mu ^-}$ is the muon-antimuon angle), whilst the energy spectrum for the $\eta ^{\prime }\rightarrow \eta \mu ^+\mu ^-$ is skewed to the right (i.e. somewhat larger values for $\theta _{\mu ^+\mu ^-}$ are preferred). This is more difficult to explain¹⁵ given that this effect, which is connected to the fact that $m_{\eta }>m_{\pi ^0}$, is a consequence of the complex dynamical interplay between the different terms in Eq. (9). Surprisingly, the kinematics of the reactions do not seem to play a significant role in this difference in skewness.

4 Conclusions

In this work, the C-conserving decay modes $\eta ^{(\prime )}\rightarrow \pi ^0l^+l^-$ and $\eta ^{\prime }\rightarrow \eta l^+l^-$ ($l=e$ or $\mu$) have been analysed within the theoretical framework of the VMD model. The associated decay widths and dilepton energy spectra have been calculated and presented for the six decay processes. To the best of our knowledge, the theoretical predictions for the four $\eta ^{\prime }\rightarrow \pi ^0l^+l^-$ and $\eta ^{\prime }\rightarrow \eta l^+l^-$ reactions that we have provided in this work are the first predictions from theory that have been published.

The decay width results that we have obtained from our calculations, which are summarised in Table 1, have been compared with those available in the published literature. In general, the agreement is reasonably good considering that the previous analysis either contain important approximations or consist of unitary lower bounds. Predictions for the dilepton energy spectra have also been presented for all the above processes, cf. Fig. 2.

Experimental measurements to date have provided upper limits to the decay processes studied in this work. These upper limits, though, are still many orders of magnitude larger than the theoretical results that we have presented. For this reason, we would like to encourage experimental groups, such as the WASA-at-COSY and REDTOP Collaborations, to study these semileptonic decays processes, as we believe that they can represent a fruitful arena in the search for new physics beyond the Standard Model.

Appendix A: Definition of parameters \(\pmb {A_i}\), \(\pmb {B_i}\), \(\pmb {C_i}\) and \(\pmb {D_i}\)

The parameters $A_i$, $B_i$, $C_i$ and $D_i$ ($i=1,2$) from Eq. (4) are defined as follows

$$\begin{aligned} A_1&= -(x+y+2z)(P\cdot p_{+})-(y+2z-1)(P\cdot p_{-}) \nonumber \\&\quad +(-x+y+1)(p_{-}\cdot p_{+})+m_{l}^2(2x+y+1) \nonumber \\&\quad +\frac{3 z m_{\eta ^{(\prime )}}^2}{2}\ , \end{aligned}$$

(A.1)

$$\begin{aligned} A_2&\equiv -A_1\ \ \mathrm {with}\ \ p_{+}\leftrightarrow p_{-}\ , \end{aligned}$$

(A.2)

$$\begin{aligned} B_1&= -2z(x+y)(x+y+z-1)(P\cdot p_{+})^2-2y z(y+z-1) \nonumber \\&\quad \times (P\cdot p_{-})^2-2x y(x+y-1)(p_{-}\cdot p_{+})^2+(P\cdot p_{-}) \nonumber \\&\quad \times \Big \{-2z\big [x(2y+z-1)+2y(y+z-1)\big ](P\cdot p_{+})\nonumber \\&\quad +2y(p_{-}\cdot p_{+})(-2x z+x+y+z-1) \nonumber \\&\quad +m_l^2\big [x^2(1-2z)+2x y+2y(y+z-1)\big ]\Big \}+ (P\cdot p_{+})\nonumber \\&\quad \times \Big \{m_l^2\big [x^2(2z-1)+2x y(2z-1)-2y(y+z-1)\big ] \nonumber \\&\quad -2y(p_{-}\cdot p_{+})(x+y+z-1)\Big \}-x^3 m_l^2(p_{-}\cdot p_{+}) \nonumber \\&\quad +z m_{\eta ^{(\prime )}}^2 \Big \{(p_{-}\cdot p_{+})\big [x(4y+z-1)+4y(y+z-1)\big ] \nonumber \\&\quad +m_l^2\big [2x^2+x(4y+3z-3)+4y(y+z-1)\big ]\Big \} \nonumber \\&\quad +x m_l^4\big [x^2+2x y+2(y-1)y\big ], \end{aligned}$$

(A.3)

$$\begin{aligned} B_2&\equiv -B_1\ \ \mathrm {with}\ \ p_{+}\leftrightarrow p_{-}\ , \end{aligned}$$

(A.4)

$$\begin{aligned} C_1&= (x+1)(P\cdot p_{-})+(4x+1)(P\cdot p_{+}) \nonumber \\&\quad -\Big (\frac{5}{2} x+1\Big ) m_{\eta ^{(\prime )}}^2\ , \end{aligned}$$

(A.5)

$$\begin{aligned} C_2&\equiv C_1\ \ \mathrm {with}\ \ p_{+}\leftrightarrow p_{-}\ , \end{aligned}$$

(A.6)

$$\begin{aligned} D_1&= 2\big [x^3+x^2(2y+2z-1)+x y(y+2z)+y(y+z-1)\big ] \nonumber \\&\quad \times (P\cdot p_{+})^2+2y(x y+y+z-1)(P\cdot p_{-})^2+x(P\cdot p_{+}) \nonumber \\&\quad \times \Big \{-2(x^2-y^2+y)(p_{-}\cdot p_{+})+m_l^2\big [x^2+2(y-1)y\big ] \nonumber \\&\quad +(z-1)z m_{\eta ^{(\prime )}}^2\Big \}+(P\cdot p_{-})\Big \{x\big [-2(y-1)y \nonumber \\&\quad \times (p_{-}\cdot p_{+}+m_l^2)+x^2m_l^2-(z-1)z m_{\eta ^{(\prime )}}^2\big ]+4(x+1) \nonumber \\&\quad \times y(x+y+z-1)(P\cdot p_{+})\Big \}-m_{\eta ^{(\prime )}}^2 \Big \{(p_{-}\cdot p_{+}+m_l^2) \nonumber \\&\quad \times \big [x^2(4y+2z-1)+4 x y(y+z)+4y(y+z-1)\big ] \nonumber \\&\quad +2x^3m_l^2\Big \}\ , \end{aligned}$$

(A.7)

$$\begin{aligned} D_2&\equiv D_1\ \ \mathrm {with}\ \ p_{+}\leftrightarrow p_{-}\ . \end{aligned}$$

(A.8)