Radiative leptonic decay of heavy quarkonia

The radiative leptonic decays of heavy quarkonia, denoted as X, i.e., $ X\to \gamma l^+ l^- $ where $ l=e, \mu, \tau $, are crucial in tests of quantum electrodynamics (QED) due to the absence of hadronic final states. Moreover, in studies on quantum chromodynamics (QCD) and searches for physics beyond the Standard Model (SM) at colliders, where copious hadronic events are produced, an accurate simulation of QED radiative background events is essential. This necessitates the utilization of software tools, such as PHOTOS [1−3], which rely on appropriate theoretical inputs regarding not only the decay branching ratios of X but also the differential distributions of the final state particles. However, only the branching ratio for the decay $ J/\Psi\to \gamma e^+ e^- $ is available from the Particle Data Group (PDG) [4] when considering X states of quantum configuration $ I^G(J^{PC})=0^- (1^{--}) $, such as $ J/\Psi $, $ \Psi(2S) $, $ \Upsilon(1S) $, and $ \Upsilon(2S) $. This value was obtained based on the predicted ratio between $ {\Gamma_{J/\Psi \rightarrow \gamma e^{+} e^{-}}} $ and $ {\Gamma_{J/\Psi \rightarrow e^{+} e^{-}}} $ [5]. Consequently, it is imperative to encourage further experimental measurements and accurate theoretical predictions for the radiative leptonic decays of heavy quarkonia.

In Ref. [5], the ratio between $ {\Gamma_{X \rightarrow \gamma l^{+} l^{-}}} $ and $ {\Gamma_{X \rightarrow l^{+} l^{-}}} $ is calculated based on the QED formalism presented in Ref. [6]. The ratio is defined using quantities in the center-of-mass (c.m.) frame of the $ l^+ l^- $ system. By utilizing the measured value of $ {\Gamma_{X \rightarrow l^{+} l^{-}}} $ and aforementioned ratio, it becomes possible to determine the radiative leptonic decay width $ {\Gamma_{X \rightarrow \gamma l^{+} l^{-}}} $ and branching ratio without considering the wave functions of the bound state X at the origin explicitly. This ratio and results from Ref. [5] have also been employed in subsequent studies [7−10] to investigate the radiative leptonic decays of quarkonia. However, the c.m. frame of the $ l^+ l^- $ system varies with the energy and direction of the final state photon. As a result, the expression proposed in Ref. [5] does not directly provide information regarding the angular and energy/momentum distributions of the final state particles in the rest frame of X.

This study focuses on reexamination of the radiative leptonic decays of heavy quarkonia within the framework of non-relativistic QCD (NRQCD) [11]. For simplicity and to facilitate comparison with the results in literature and those of PDG, this study does not include the processes of $ X\to X^\prime \gamma \to \gamma l^+ l^- $, where $ X^\prime $ is another resonance lighter than X. The heavy quarkonia X is regarded as a bound state of a quark ($ Q=c,b $) and its antiquark ($ \bar{Q} $), and its wave function is analyzed in both momentum and position spaces. At the leading-order (LO) level, the decay width of the quarkonia X into specific final states is proportional to the square of the wave function at the origin, denoted as $ |\psi({\vec{0}})|^2 $. As the ratio between different decay widths of X decay modes is independent of $ |\psi({\vec{0}})|^2 $, we provide preditions for $ {\Gamma_{X \rightarrow \gamma l^{+} l^{-}}} $ in both the rest frame of X and c.m. frame of $ l^+ l^- $ by utilizing $ {\Gamma_{X \rightarrow l^{+} l^{-}}} $ and ratio $ {\Gamma_{X \rightarrow \gamma l^{+} l^{-}}}/{\Gamma_{X \rightarrow l^{+} l^{-}}} $. Numerical integration in the rest frame of X provides insights into the angular and energy/momentum distributions of the final state particles. On the other hand, the relation in the c.m. frame of $ l^+ l^- $ yields the ratio between $ {\Gamma_{X \rightarrow \gamma l^{+} l^{-}}} $ and $ {\Gamma_{X \rightarrow l^{+} l^{-}}} $, which can be compared to the corresponding ratio proposed in Ref. [5]. The results can be incorporated into some generators, similar to PHOTOS, as theoretical inputs for simulation of QED radiative processes during tests of the SM and search for new physics.

This paper is organized as follows. In Section II, a brief introduction to the production and decay of heavy quarkonia within the framework of NRQCD is provided. Section III presents an analytical description of the kinematics of the final state particles in the decay $ X \to \gamma l^{+} l^{-} $, focusing on the X rest frame. In Section IV, to illustrate the analysis, the decay $ J/\Psi \to \gamma \mu^+ \mu^- $ is used as an example, and the angular and energy/momentum distributions of the final state particles are derived at the LO level within the NRQCD framework. In Section V, the Lorentz transformation from the X rest frame to the $ l^+ l^- $ c.m. frame is performed. This section presents the derivation of the ratio between $ {\Gamma_{X \rightarrow \gamma l^{+} l^{-}}} $ and $ {\Gamma_{X \rightarrow l^{+} l^{-}}} $, and the proposed expression is compared with that in Ref. [5]. Finally, the conclusions are presented in Section VI.

In this study, we use the symbol X to represent heavy quarkonia of the $ Q\bar{Q} $ ($ Q=c,b $) state. Within the framework of NRQCD [11], the physics of X production can be divided into a perturbative part and non-perturbative part based on the energy scale. The perturbative part involves the production of a heavy-quark pair $ Q\bar{Q}[n] $ with quantum configuration n at a hard scale $ \Lambda\sim m_Q $, where $ n=\; ^{2S+1}L_J^{[1,8]} $ (S, L, and J represent the total spin, orbital angular momentum, and total angular momentum of the $ Q\bar{Q} $ pair, respectively; the $ Q\bar{Q} $ pair is in either a color-singlet or color-octet state indicated by the superscript [1, 9]). The non-perturbative part refers to the hadronization of the $ Q\bar{Q}[n] $ pair into the bound state X at a significantly lower scale $ \Lambda^\prime<m_Q $. This factorization can be expressed as [12]

where ${\rm d}\hat{\sigma}_{ab}(Q\bar{Q}[n]+\{k\})$ describes the production of the heavy-quark pair $ Q\bar{Q} $ of quantum configuration n with additional final state partons $ \{k\} $. $ \langle \mathcal{O}_X^n\rangle $represents the non-perturbative long-distance matrix element (LDME), which describes the hadronization of the $ Q\bar{Q}[n] $ state into the bound state X.

Similarly, within the NRQCD framework, the decay of X into light hadrons can be factorized into two parts as [11]

where $ \delta\mathcal{L}_\text{4-fermion} $ represents the four-fermion correction terms in the NRQCD Lagrangian. This includes the summation of different local four-fermion operators $ \mathcal{O}_n $ of quantum configuration n multiplied by corresponding coefficients $ f_n $. $ d_n $ represents the naive scaling dimension of the operator $ \mathcal{O}_n $. The coefficient $ \text{Im}f_n $, which can be computed perturbatively, is proportional to the rates for on-shell Q and $ \bar{Q} $ to annihilate from the initial quantum configuration n into the final states. The non-perturbative matrix elements $ \langle X |\mathcal{O}_n|X \rangle $ indicate the probability of finding Q and $ \bar{Q} $ in the quantum configuration n within the bound state X, which is suitable for annihilation.

In this study, we consider only the (radiative) leptonic decay of X with the quantum configuration $ I^G(J^{PC})= 0^- (1^{--}) $, denoted as $ X\longrightarrow (\gamma) l^{+} l^{-} $. Within the NRQCD framework, the decay width can be expressed as [11]

where v represents the velocity of Q ($ \bar{Q} $) in the bound state X. The average value of $ v^2 $ is approximately 0.3 for charmonium and around 0.1 for bottomonium [13]. The second term in Eq. (3) is proportinal to $ v^2 $. As the processes of $ X\longrightarrow l^{+} l^{-} $ and $ X\longrightarrow (\gamma) l^{+} l^{-} $ share almost the same Feynman diagrams according to the optical theorem [11],

where α and $ \alpha_S $ are coupling constants of electromagnetic interaction and strong interaction, respectively. Thus,

where $ \text{Im}f^{\text{LO}}_{l^+ l^-} $ and $ \text{Im}f^{\text{LO}}_{\gamma l^+ l^-} $ are calculated from the diagrams in Fig. 1 and Fig. 2, respectively. Eq. (6) demonstrates that, by considering only the diagrams in Fig. 2 and Fig. 1 to calculate the ratio between $ \Gamma(X\longrightarrow \gamma l^{+} l^{-}) $ and $ \Gamma(X\longrightarrow l^{+} l^{-}) $, the next-to-leading order (NLO) corrections are of $ O(\alpha) $ and $ O(v^4) $.

Figure 1.  Diagrams for decay $ X\longrightarrow l^{+} l^{-} $.

Figure 2.  Diagrams for decay $ X\longrightarrow \gamma l^{+} l^{-} $.

For the radiative leptonic decay of a quarkonium, denoted as $ X\longrightarrow \gamma l^{+} l^{-} $, the four-momenta of the final state particles in the rest frame of X can be parameterized as follows:

where $ p_X^\mu=(m_X,\; 0,\; 0,\; 0) $. Utilizing the on-shell condition of $ l^+ $, we obtain the relation

The variables $ p_{l^-} $, $ E_\gamma $, $ \theta_{l^-} $, and $ \theta_\gamma $ are complete to describe the kinematics of the final state particles. The sensible regions of the parameter space should ensure that $ E_{l^+} $ obtained from Eq. (9) is larger than $ m_l $, and $ \cos\phi_\gamma $ obtained by Eq. (10) ranges from $ -1 $ to $ 1 $. Consequently, the accessible ranges of the parameters are determined as follows:

where $ E_{\gamma\text{cut}} $ represents the photon energy cut (considering only photons with energies higher than $ E_{\gamma\text{cut}} $), and

At the LO level ($ v=0 $) within the NRQCD, we can consider X with the quantum configuration $ I^G(J^{PC})= 0^- (1^{--}) $ as a bound state of $ Q \bar{Q}[n] $ with $ n=\; ^{3}S_1^{[1]} $. The momentum-space quantum state of X can be expressed as [14]

where $ \widetilde{\psi}(\vec{k}) $ represents the wave function in momentum space, and $ \vec{k} $ ($ -\vec{k} $) is the three-momentum of Q ($ \bar{Q} $) in the rest frame of X, satisfying the condition $ |\vec{k}|\ll m_Q $. The amplitude for X to decay into specific final states is expressed by [14]

To obtain Eq. (19), we approximately treat Q and $ \bar{Q} $ as static since $ v=0 $ ($ |\vec{k}|=0 $) at the LO level within the NRQCD. We also use the approximation $ 2m_Q\approx m_X $. In Eq. (19), $ \psi(\vec{0}) $ represents the wave function of X at the origin ($ \vec{0} $), obtained from the Fourier transformation of $ \widetilde{\psi}(\vec{k}) $ :

Consequently, the decay width of X into specific final states can be approximated as [14]

where $ \Pi_f $ represents the phase space of the final state particles. The correlation between Eq. (21) and Eq. (2) is established through

Considering the described parameters,

where $ \cos\phi_\gamma $ is expressed by Eq. (10), with $ \theta_l $ and $ \theta_\gamma $ replaced by $ \dfrac{\theta_{+}+\theta_{-}}{2} $ and $ \dfrac{\theta_{+}-\theta_{-}}{2} $, respectively.

In Fig. 2, we present the diagrams for the decay $ X\longrightarrow \gamma l^{+} l^{-} $. As $ m_X\ll m_Z, m_H $ (where $ m_Z $ and $ m_H $ represent the masses of the Z boson and Higgs boson, respectively), we consider only the Feynman diagrams where the photon serves as the intermediate line for the decay $ X\longrightarrow \gamma l^{+} l^{-} $. Additionally, due to charge-conjugation invariance [5], the photon in the final state can be emitted only by one of the charged leptons in the final state. Therefore, we consider the diagrams shown in Fig. 2. A direct calculation yields the following result:

where $ \mathcal{Q}_{Q/\bar{Q}} $ and $ \lambda_{Q/\bar{Q}} $ represent the electric charge and helicity of $ Q/\bar{Q} $ from the quarkonia X, respectively. Experimental findings show that, for the vector particle $ J/\Psi $, the spin projection along the z-axis is restricted to values of $ \pm 1 $, which corresponds to the state being represented as $ |S,S_z \rangle=|1,\pm 1 \rangle $. This requirement ensures that $ \lambda_Q=-\lambda_{\bar{Q}} $ at the LO within the NRQCD framework. Consequently, the squared amplitude where $ \lambda_Q=-\lambda_{\bar{Q}} $ exclusively contributes to the decay process of $ J/\Psi $. Using the relation

we obtain

As an example, we consider the decay $ J/\Psi \to \gamma \mu^+ \mu^- $. By considering the squared amplitude with $ \lambda_c=-\lambda_{\bar{c}} $ and phase space in the rest frame of $ J/\Psi $, numerical integration yields the angular and energy/momentum distributions of the final state particles, as shown in Fig. 3 with $ E_{\gamma\text{cut}}=0.025,\; 0.1 $ GeV and in Fig. 4 with $ E_{\gamma\text{cut}}=0.1 $ GeV. In the left middle panel of Fig. 3, noticeable kinks are observed in the lepton momentum ($ p_{\mu^-} $) distributions. These kinks arise due to a shift in the permissible range of $ E_\gamma $ values when $ p_{\mu^-} $ surpasses $ p_{l\text{mid}} $, as detailed in Eq. (12). The separation between the kink location ($ p_{l\text{mid}} $) and maximum value of $ p_{\mu^-} $ ($ p_{l\text{max}} $) is $ E_{\gamma\text{cut}} $, in accordance with Eq. (15). Regarding the photon energy distribution depicted in the right middle panel of Fig. 3, variations in the photon energy cutoff ($ E_{\gamma\text{cut}} $) merely alter the position of the left endpoint of the curve. The distribution tends toward divergence as the photon energy approaches zero. In the lowest panel of Fig. 3, radiative photons predominantly align along the collinear ($ \theta_{\gamma\mu^-}= 0 $) and antiparallel ($ \theta_{\gamma\mu^-}= \pi $) directions relative to the $ \mu^- $ movement, a phenomenon attributed to collinear singularity. Notably, the distribution exhibits asymmetry, with a higher propensity for photons to align in the antiparallel ($ \theta_{\gamma\mu^-}= \pi $) direction to the $ \mu^- $ trajectory. This asymmetry is due to the uneven phase space available at $ \theta_{\gamma\mu^-}= 0 $ and $ \theta_{\gamma\mu^-}= \pi $. With $ E_{\gamma\text{cut}} $ set at 0.1 GeV, the allowed ranges for $ p_{\mu^-} $ are [0, 1.44] GeV and [0, 1.54] GeV at $ \theta_{\gamma\mu^-}= 0 $ and $ \theta_{\gamma\mu^-}= \pi $, respectively.

Figure 3.  (color online) Angular and energy/momentum distributions of γ and $ \mu^- $ from the decay $ J/\Psi \to \gamma \mu^+ \mu^- $ in the $ J/\Psi $ rest frame. $ \theta_{\gamma} $, $ \theta_{\mu^-} $, and $ \theta_{\gamma\mu^-} $ are the polar angle of γ, polar angle of $ \mu^- $, and angle between the moving directions of γ and $ \mu^- $, respectively. $ E_\gamma $ and $ p_{\mu^-} $ are the photon energy and $ \mu^- $ momentum, respectively. Two kinds of photon energy cut are chosen.

Figure 4.  (color online) Left panel: polar angle-energy distribution of γ from the decay $ J/\Psi \to \gamma \mu^+ \mu^- $ in the $ J/\Psi $ rest frame at $ E_{\gamma \text{cut}}=0.1 $ GeV. Right panel: polar angle-momentum distribution of $ \mu^- $ from the decay $ J/\Psi \to \gamma \mu^+ \mu^- $ in the $ J/\Psi $ rest frame at $ E_{\gamma \text{cut}}=0.1 $ GeV.

Although the final states consist of three particles ($ l^+ $, $ l^- $, and γ), it is advantageous to analyze the decay of X in the c.m. frame of the $ l^+ l^- $ system. In this frame, we can employ a set of variables to describe the angular distributions of the photon and leptons, as well as their relative motion. We introduce $ \Omega_\gamma^{\prime} $ ($ \theta_\gamma^\prime $ and $ \phi_\gamma^\prime $) and $ \Omega_l^{\prime} $ ($ \theta_l^\prime $ and $ \phi_l^\prime $) to represent the angles of the photon and leptons, respectively, in the $ l^+ l^- $ c.m. frame. Furthermore, we use $ \theta^\prime_{\gamma l} $ to denote the angle between the directions of motions of the photon and lepton in this frame. To quantify the kinematics, we introduce $ \beta^\prime $, $ E_\gamma^\prime $, and $ s^\prime $ to represent the velocity of the leptons, energy of the photon, and square of the invariant mass of the $ l^+ l^- $ system, respectively. Using the given quantities, we can parameterize the four-momentum of the final states in the c.m. frame of $ l^+ l^- $ as follows:

where we have chosen the direction of motion of $ l^- $ as the z-axis direction, resulting in $ \theta_\gamma^\prime=\theta^\prime_{\gamma l} $. Using the relation $ p^{\prime 2}_X=(p_{l^-}^{\prime}+p_{l^+}^{\prime}+p_{\gamma}^{\prime})^2=m_X^2 $, we obtain

To ensure $ s^\prime > 4m_l^2 $, we obtain the upper limit for $ E_\gamma^\prime $ by

The photon energy cut $ E_{\gamma\text{cut}} $ in the X rest frame corresponds to the lower limit $ E^\prime_{\gamma\text{cut}} $ for $ E_\gamma^\prime $ in the $ l^+ l^- $ c.m. frame, expressed by

The parameter of the Lorentz transformation from the X rest frame to the $ l^+ l^- $ c.m. frame is expressed by

where $ \vec{p}^\prime_\gamma $ is the three-momentum of the photon in the $ l^+ l^- $ c.m. frame. Thus, the $ l^+ l^- $ c.m. frame changes with the photon energy and angles. We need to determine the Jacobi factor J appearing in the following equation:

where $ \vec{p}_{l^-} $, $ \vec{p}_{l^+} $, and $ \vec{p}_{\gamma} $ ($ \vec{p}^\prime_{l^-} $, $ \vec{p}^\prime_{l^+} $, and $ \vec{p}^\prime_{\gamma} $) are the three-momenta of $ l^- $, $ l^+ $, and γ in the X rest frame ($ l^+ l^- $ c.m. frame), respectively. Direct derivations demonstrate that the Jacobi factor J is related to the partial derivatives of the four-momentum of $ l^- $ and γ in the X rest frame with respect to those in the $ l^+ l^- $ c.m. frame, which yields

Therefore,

which indicates that the production of high-energy photons in the $ l^+ l^- $ c.m. frame is suppressed.

According to Eq. (21), the decay width of quarkonia X decaying into specific final states is directly proportional to the squared amplitude of the wave function at the origin, represented as $ |\psi(\vec{0})|^2 $. Consequently, the ratio between the decay widths of different decay modes of X is independent of the specific value of $ |\psi(\vec{0})|^2 $. Considering the phase space, Jacobi factor resulting from the Lorentz transformation, and squared amplitude in the center-of-mass frame of the $ l^+ l^- $ pair, by integration over the final state angles, we derive the following expression:

where

The analytical solution for the integration outlined in Eq. (40) is presented in Appendix A.

The predicted ratios of $ {\Gamma_{X \rightarrow \gamma l^{+} l^{-}}}/{\Gamma_{X \rightarrow l^{+} l^{-}}} $, calculated using $ g^\prime (b,c) $, are depicted in Fig. 5. These ratios exhibit a decrease as either the lepton mass $ m_l $ or photon energy cut $ E_{\gamma\text{cut}} $ increases. Utilizing the experimentally measured values of $ R_{J/\Psi \to e^+ e^-} $ = (5.971±0.032)% and $ R_{J/\Psi \to \mu^+ \mu^-} $ = (5.961 ± 0.033)% obtained from the PDG [4], and employing the ratio $ {\Gamma_{X \rightarrow \gamma l^{+} l^{-}}}/{\Gamma_{X \rightarrow l^{+} l^{-}}} $ derived from Eq. (40), we can calculate the branching ratios of $ J/\Psi \to \gamma e^+ e^- $ and $ J/\Psi \to \gamma \mu^+ \mu^- $ at various $ E_{\gamma \text{cut}} $ values, as depicted in Fig. 5. Notably, the ratio between $ R_{J/\Psi \to \gamma e^+ e^-} $ and $ R_{J/\Psi \to \gamma \mu^+ \mu^-} $ consistently remains around 3.0 within the range of $ E_{\gamma\text{cut}} $ spanning from $ \mathcal{O}(0.01) $ to $ \mathcal{O}(0.1) $ GeV. The future experimental results of this ratio will serve as a direct test of lepton flavor universality (LFU). A more comprehensive analysis of this ratio, encompassing the robustness of theoretical predictions in the face of radiative and power corrections, susceptibility to possible breaches of LFU, and practicality of experimentally detecting such BSM phenomena via this decay rate ratio, is beyond the scope of this study. These facets are designated for subsequent research endeavors.

Figure 5.  (color online) Upper panel: predicted values of $ {\Gamma_{X \rightarrow \gamma l^{+} l^{-}}}/{\Gamma_{X \rightarrow l^{+} l^{-}}} $ using $ g^\prime (b,c) $. Lower left panel: branching ratios of $ J/\Psi \to \gamma e^+ e^- $ and $ J/\Psi \to \gamma \mu^+ \mu^- $ at different $ E_{\gamma \text{cut}} $. Lower right panel: ratio between $ R_{J/\Psi \to \gamma e^+ e^-} $ and $ R_{J/\Psi \to \gamma \mu^+ \mu^-} $ at different $ E_{\gamma \text{cut}} $.

In Ref. [5], the differential decay width for the process $ X \rightarrow \gamma l^{+} l^{-} $ is expressed by

where

By integrating out the final state angles ($ \Omega^\prime_\gamma $ and $ \Omega^\prime_l $), we obtain the ratio of the decay widths:

The analytical result for the integration in Eq. (45) is provided in Appendix A.

The comparison between $ g(b,c) $ derived from the relation proposed in Ref. [5] and $ g^\prime(b,c) $ obtained in this study reveals differences arising from the prefactors and integrated functions of $ \beta^\prime $.

The dependencies of $ \dfrac{g(b,c)}{2\pi g^\prime(b,c)} $ on b and c are illustrated in Fig. 6. The ratio $ \dfrac{g(b,c)}{2\pi g^\prime(b,c)} $ is larger than one for small values of b (e.g., $ b<0.03 $ when $ c=0.02 $) and decreases as b increases.

Figure 6.  (color online) Ratio between the predicted values of $ {\Gamma_{X \rightarrow \gamma l^{+} l^{-}}}/{\Gamma_{X \rightarrow l^{+} l^{-}}} $ obtained using $ \dfrac{g(c,b)}{2\pi} $ and $ g^\prime (b,c) $.

In the specific case where $ E_{\gamma \text{cut}}=0.1 $ GeV and X includes $ J/\Psi $, $ \Psi(2S) $, $ \Upsilon(1S) $, and $ \Upsilon(2S) $, the ratio between $ {\Gamma_{X \rightarrow \gamma l^{+} l^{-}}} $ and $ {\Gamma_{X \rightarrow l^{+} l^{-}}} $, as determined by $ g(b,c) $ and $ g^\prime (b,c) $, is presented in Table 1. For modes with b approximately equal to or smaller than $ \mathcal{O}(10^{-2}) $, $ \dfrac{g(c,b)}{2\pi} $, rather than $ g(c,b) $, is comparable to $ g^\prime (c,b) $, with a relative difference within 10%, which is comparable to the NLO corrections of $ O(\alpha) $ and $ O(v^4) $ according to Eq. (6). The largest discrepancy between $ \dfrac{g(c,b)}{2\pi} $ and $ g^\prime (c,b) $ is observed for the mode $ {{\Psi(2S) \rightarrow (\gamma) \tau^{+} \tau^{-}}} $, where $ b=0.232 $. $ g^\prime (c,b) $ is four times larger than $ \dfrac{g(c,b)}{2\pi} $. Compared to $ \dfrac{g(c,b)}{2\pi} $, the factor $ \dfrac{1}{(1+2b)\sqrt{1-4b}} $ in $ g^\prime (c,b) $ yields a 2.5-fold enhancement for the decay mode $ {{\Psi(2S) \rightarrow (\gamma) \tau^{+} \tau^{-}}} $ when $ b=0.232 $. By disregarding the differences in the prefactors between ${g(c,b)}/{2\pi}$ and $ g^\prime (c,b) $, which are negligible when $ b\ll 1 $, the remaining enhancement in $ g^\prime (c,b) $ for the $ {{\Psi(2S) \rightarrow (\gamma) \tau^{+} \tau^{-}}} $ decay mode can be attributed to the variations in the integrated functions.

Mode	b	c ($ E_{\gamma\text{cut}}=0.1 $ GeV)	$ g^\prime(c,b) $	$ \frac{g(c,b)}{2\pi g^\prime(c,b)} $
$ {{J/\Psi \rightarrow (\gamma) e^{+} e^{-}}} $	$ 2.72\times 10^{-8} $	$ 3.34\times 10^{-2} $	0.152	1.11
$ {{J/\Psi \rightarrow (\gamma) \mu^{+} \mu^{-}}} $	$ 1.16\times 10^{-3} $	$ 3.34\times 10^{-2} $	$ 5.04\times 10^{-2} $	1.10
$ {{\Psi(2S) \rightarrow (\gamma) e^{+} e^{-}}} $	$ 1.92\times 10^{-8} $	$ 2.79\times 10^{-2} $	0.168	1.10
$ {{\Psi(2S) \rightarrow (\gamma) \mu^{+} \mu^{-}}} $	$ 8.22\times 10^{-4} $	$ 2.79\times 10^{-2} $	$ 5.84\times 10^{-2} $	1.10
$ {{\Psi(2S) \rightarrow (\gamma) \tau^{+} \tau^{-}}} $	0.232	$ 2.79\times 10^{-2} $	$ 6.09\times 10^{-6} $	0.25
$ {{\Upsilon(1S) \rightarrow (\gamma) e^{+} e^{-}}} $	$ 2.92\times 10^{-9} $	$ 1.07\times 10^{-2} $	0.266	1.08
$ {{\Upsilon(1S) \rightarrow (\gamma) \mu^{+} \mu^{-}}} $	$ 1.25\times 10^{-4} $	$ 1.07\times 10^{-2} $	0.111	1.07
$ {{\Upsilon(1S) \rightarrow (\gamma) \tau^{+} \tau^{-}}} $	$ 3.53\times 10^{-2} $	$ 1.07\times 10^{-2} $	$ 2.80\times 10^{-2} $	0.98
$ {{\Upsilon(2S) \rightarrow (\gamma) e^{+} e^{-}}} $	$ 2.60\times 10^{-9} $	$ 1.01\times 10^{-2} $	0.273	1.07
$ {{\Upsilon(2S) \rightarrow (\gamma) \mu^{+} \mu^{-}}} $	$ 1.11\times 10^{-4} $	$ 1.01\times 10^{-2} $	0.115	1.07
$ {{\Upsilon(2S) \rightarrow (\gamma) \tau^{+} \tau^{-}}} $	$ 3.14\times 10^{-2} $	$ 1.01\times 10^{-2} $	$ 3.04\times 10^{-2} $	0.99

Table 1.  $ {\Gamma_{X \rightarrow \gamma l^{+} l^{-}}}/{\Gamma_{X \rightarrow l^{+} l^{-}}} $ determined by $ g^\prime (b,c) $ and ratio between the predicted values of $ {\Gamma_{X \rightarrow \gamma l^{+} l^{-}}}/{\Gamma_{X \rightarrow l^{+} l^{-}}} $ using $ \frac{g(c,b)}{2\pi} $ and $ g^\prime (b,c) $. X represents $ J/\Psi $, $ \Psi(2S) $, $ \Upsilon(1S) $, and $ \Upsilon(2S) $.

In this research, we delved into the leptonic decays of heavy quarkonia X, conceptualized as a zero-momentum bound state of a quark (Q) and an antiquark ($ \bar{Q} $) within the NRQCD framework at the LO level. Through a detailed analysis of the squared amplitude for the decay process $ X \to \gamma l^+ l^- $and an examination of the kinematics of the final state particles in the rest frame of X, we calculated the decay widths and branching ratios for the three-body leptonic decays of X. Our study was also extended to the determination of the angular and energy/momentum distributions of the final state particles in the rest frame of X via numerical integration.

To explore the Lorentz transformation from the rest frame of X to the c.m. frame of $ l^+ l^- $, we derived the Jacobian of this transformation and analyzed the phase space in the $ l^+ l^- $ c.m. frame. Through these analytical approaches, we propose a prediction for the ratio between $ {\Gamma_{X \rightarrow \gamma l^{+} l^{-}}} $ and $ {\Gamma_{X \rightarrow l^{+} l^{-}}} $, highlighting discrepancies from previously reported expressions, particularly from that in Ref. [5]. Our analysis identifies a missing factor of $ 1/2\pi $ in the earlier expression and suggests that a detailed comparison of the equations can elucidate the differences. The discrepancies become pronounced when the ratio $ m_l^2/m_X^2 $ crosses a threshold of approximately $ \mathcal{O}(0.1) $.

Our investigation spans various heavy quarkonia, including $ J/\Psi $, $ \Psi(2S) $, $ \Upsilon(1S) $, and $ \Upsilon(2S) $. The ratios derived from our proposed relation align closely with those from literature (after adjustment by $ 2\pi $), which showcases relative differences generally within or below 10%. This level of agreement is on par with the expected NLO corrections of $ O(\alpha) $ and $ O(v^4) $. An outlier is observed in the ratio between $ {\Gamma_{\Psi(2S) \to \gamma \tau^+ \tau^-}} $ and $ {\Gamma_{\Psi(2S) \to \tau^+ \tau^-}} $, where our proposed expression predicts a value four times larger than expected. The findings of this research offer a solid foundation for comparison with future experimental results from facilities such as BESIII, B-factory experiments, and CEPC, thereby enhancing our understanding of heavy quarkonia decays and contributing to tests of the SM and search for new physics.

The analytical solutions for the integrations described in Eqs. (41) and (46) are

and

respectively, where

and $ \Re(z) $ represents the real part of z.