Four-lepton decays of neutral vector mesons

We investigate the rare electromagnetic decays of neutral vector mesons (exemplified by $J/\psi$, $\Upsilon$, $\rho^0$, $\omega,\phi$, \ldots) into four charged leptons ($2(e^+e^-)$, $2(\mu^+\mu^-)$, $2(\tau^+\tau^-)$, $e^+e^-\mu^+\mu^-$, $e^+e^-\tau^+\tau^-$, \ldots) at lowest order in QED. In contrast to the case of vector meson decay into a single pair of leptons, the lepton mass must be retained in these four-lepton decay channels to cutoff the mass singularity. Owing to more pronounced collinear enhancement, the branching fractions of vector mesons decays into $ 2(e^+e^-)$ are considerably greater than those for decays into $2(\mu^+\mu^-)$. The decay channels $J/\psi(\Upsilon)\to 2(e^+e^-), e^+e^-\mu^+\mu^-$ are predicted to have branching fractions of order $10^{-5}$, which appear to have bright observation prospect at {\tt BESIII}, {\tt Belle 2} and {\tt LHC} experiments.

The leptonic decays of neutral vector mesons are arguably one of the cleanest and most exhaustively studied hadronic processes. To dematerialize into a lepton pair, the valence quark and antiquark inside the vector meson have to first annihilate into a virtual photon, which then subsequently transition into a lepton pair. Since the hadronic decay part and the electromagnetic production part of the amplitude simply factorizes, one can write down the well-known formula for vector meson leptonic decay: where α denotes the fine structure constant, e q signifies the electric charge of the valence quark of V . M V denotes the vector meson mass, and f V denotes the decay constant associated with V 1 : where ε µ signifies the polarization vector of the V meson.
In (1), the lepton mass has been neglected for simplicity. It is obviously a perfect approximation for the vector meson decay into electrons, and still a decent one for heavy vector quarkonia such as J/ψ and Υ decay into a muon pair.  As one can tell from Table I and II, the branching fractions of various leptonic decays of vector mesons have been measured quite precisely. From these experimental inputs, one sees that discarding the lepton mass is indeed a satisfactory approximation in (1), and one can also deduce the values of the nonperturbative factor f V for different vector mesons.
At present, BESIII [4], Belle 2 [5], KLOE experiment [6], and LHC experiments [7] have accumulated a gigantic number of J/ψ and Υ, as well as light vector mesons ρ, ω and φ, which make it feasible to search some truly rare decay channels. The aim of this paper is to explore a special kind of rare electromagnetic decay processes, that is, neutral vector meson decays into four charged leptons. Concretely speaking, we are interested in considering two types of electromagnetic decays, V → 2(l + l − ) and l + l − l ′+ l ′− with l, l ′ = e, µ, τ . Such multilepton rare decay of quarkonia represents a novel testing bed for QED predictions, and we hope they can be seen at BESIII, Belle 2, KLOE and LHC experiments in near future.
The partial width of neutral vector meson decay to four leptons reads where dΠ 4 is the four-body phase space with possible symmetry factors due to identical leptons. We will dedicate Appendix A for a detailed account of the technicalities related to 4-body phase space integration. In the second equality, we utilize the fact that, entirely analogous to the process V → 2l, the decay amplitude simply factorizes into the product of the V → γ * hadronic part and the γ * → 4l leptonic part.
We employ the packages FeynArts [8] to generate the Feynman diagrams for V → 4l. The corresponding four Feynman diagrams for V → l + l − l ′+ l ′− have been displayed in Fig. 1, while the corresponding eight diagrams for V → 2(l + l − ) V → 2(e + e − ) are illustrated in Fig. 2. Subsequently we use FeynCalc [9] to handle Dirac trace/Lorentz tensor algebra.
For phenomenological purpose, it is more convenient to introduce the following dimensionless ratio R, rather than directly compute the partial width: As is evident in (1) and (3), the advantage of introducing the R ratio is that the nonpertubative factor f V , as well as the electric charge of valence quark e q , exactly cancel between the numerator and the denominator. Therefore, the prediction of the R ratio is dictated entirely by QED, without contamination by any nonperturbative effects. Since R is dimensionless, it can only be a function of the mass ratio m l /M V . We envisage that the leading order QED prediction should suffice to provide a decent account for the R values for various V → 4l channels.
In (1) we have safely discarded the lepton mass in V → l + l − since m e,µ ≪ M V . One may wonder whether the same simplification can be taken or not for V → 4l. As a matter of fact, for the four-body decays, it is possible that three particles can move collinear and some are simultaneously soft, in which the collinear and soft singularity may be developed. In fact, to assess the degreee of mass singularity, we invoke the reverse unitarity method [10] to analytically compute the R ratio for vector meson V decay to four massless leptons (as depicted in Fig. 2) in spacetime dimension d = 4 − 2ǫ: where µ is the 't Hooft unit mass. Note that after supplemented with proper color factor, this piece may be identified with the Abelian subset of the double-real part of the next-tonext-leading order (NNLO) QCD corrections for e + e − → γ * → qq, which has already been known analytically long ago [11,12]. Note the occurrence of the triple infrared pole in (5) unambiguously indicates that the R ratio in (4) does not have a safe m l → 0 limit, in sharp contrast to V → e + e − . In our numerical study we work in d = 4 spacetime dimension and keep m l finite to regularize the emerging mass singularity. It is intuitively envisaged the leading 1/ǫ 3 pole in (5) would be literally translated into the triple-logarithm ln 3 (M 2 V /m 2 l ) in our case. Expanding the µ-dependent prefactor in (5) to order ǫ 3 , multiplying with the 1/ǫ 3 pole, replacing µ by m l , we readily identify the leading triple-logarithmic term: It is illuminating to trace the origin of the triple-logarithm in (6) for the case of realistic lepton, or equivalently, the 1/ǫ 3 pole in (5) for the case of massless lepton. Near the tricollinear limit, i.e., when one of the l + and two l − become nearly collinear, the squared amplitude for γ * → 4l can be factorized into the squared amplitude for γ * → l + l − times the universal NNLO splitting function for l − → l − + (l − l + ) [13]. Double collinear poles would immediately arise in the tri-collinear limit. Furthermore, when the invariant mass of the l − l + -pair also becomes simultaneously soft, the adjacent photon propagator also becomes soft, which would develop an additional soft pole. We have explicitly verified that, upon integrating the NNLO splitting function over three-lepton phase space, the leading triple IR pole in (5) can indeed be recovered [14].
For numerical analysis, we adopt the following values for various input parameters [1]: FIG. 3. The ratio R(V → 4l) as a function of ξ = m l M V . Note that this is a log-log plot, and the logarithmical divergence is reflected by the negative slope of the curve in the ξ → 0 limit.  III. Predictions for R ratios and the affiliated branching fractions for light vector mesons → 4l. Note the latter is obtained by multiplying R with the corresponding measured branching fraction of V → l + l − , with our theoretical uncertainty entirely inherited from experimental error on di-leptonic decay of V .
We use the numerical package CUBA [15] to conduct the 4-body phase space integration, which implements the adaptive Monte Carlo integration algorithm.
The dependence of the R ratio for V → 4l on ξ ≡ m l /M V is shown in Fig. 3. It can be readily visualized that the curve becomes singular as ξ → 0. A numerical fit reveals that the small-m l enhancement is compatible with the triply-logarithmical scaling behavior.
Our detailed predictions for ρ 0 , ω, φ → 4l and J/ψ, Υ → 4l are tabulated in Table III  and Table IV. By multiplying the measured branching fractions for V → l + l − as given in Table I and II, we then predict the branching fractions for V → 4l. As one can see, the branching fractions for ρ 0 , ω → 2(e + e − ), e + e − µ + µ − are of order 10 −8 , while those for φ → 2(e + e − ), e + e − µ + µ − may reach 10 −7 . The rather tiny branching fractions make the observation of these rare decays of light vector mesons a great challenge. On the other hand, the branching fractions of J/ψ → 2(e + e − ), e + e − µ + µ − can reach the order of 10 −5 2 . Concerning the billions of J/ψ events collected at BESIII, it is hopeful to discover these two channels in near future. In contrast, because muon is much heavier than electron, the triple logarithmic enhancement for J/ψ → 2(µ + µ − ) is much less pronounced. From Table IV, we observe B(J/ψ → 2(µ + µ − )) is only about 2% of B(J/ψ → 2(e + e − )),  Predictions for R ratios and the affiliated branching fractions for J/ψ, Υ → 4l. Note the latter is obtained by multiplying R with the measured branching fraction for V → l + l − , with our theoretical uncertainty entirely propagating from experimental error.
which renders the observation prospect in near future obscure. The pattern of predictions for Υ → 4l is similar. B(Υ → 2(e + e − )), e + e − µ + µ − ), e + e − τ + τ − ) are of order 10 −5 , while B(Υ → 2(µ + µ − )), µ + µ − τ + τ − ) are 1 ∼ 2 orders of magnitude smaller. Compared with these channels, B(Υ → 2(τ + τ − )) is too tiny to have any chance for observation. It is interesting to note that, the branching fractions associated with three decay channels Υ → 2(e + e − ), e + e − µ + µ − and Υ → e + e − τ + τ − do not exhibit considerable hierarchy, all of which is of order 10 −5 . We hope the dedicated analysis in Belle 2 and perhaps LHC will observe these channels in the foreseeable future. The four-body leptonic final state allows one to probe a variety of differential spectra, which can sharpen the test of our QED predictions in a greater detail. From Fig. 4 to Fig. 6, we present the invariant mass spectra of various combination of lepton pairs in the final state. It is desirable if the future experiments can confront our predictions for invariant mass distributions.
In summary, in this work we have investigated a specific type of rare electromagnetic decays of vector mesons, where the final state is composed of four charged leptons. To the best of our knowledge, the current work represents the first theoretical study of these processes. For simplicity, we have only considered the lowest-order QED contribution, which should already constitute decent predictions according to our experiences. In sharp contrast to the predicted width for vector meson into a lepton pair, where the lepton mass is usually discard, we must explicitly retain the finite lepton mass to regularize the potential collinear/soft divergences. The branching fractions of V → 2(e + e − ) is orders of magnitude greater than that of V → 2(µ + µ − ), because of more pronounced triple-logarithmic enhancement received by the former channel. We predict that the branching fractions associated with the J/ψ → 2(e + e − ), e + e − µ + µ − and Υ → 2(e + e − ), e + e − µ + µ − , e + e − τ + τ − channels all reach the order 10 −5 , which may indicate bright observation opportunity at BESIII, Belle 2 and LHC experiments in near future.

ACKNOWLEDGMENTS
We thank Haibo Li for interesting discussion that stimulate us to initiate this work. We are also grateful to Xiaohui Liu for useful discussions on triple logarithms in the massless lepton limit. W. C. wishes to acknowledge the warm hospitality extended to him by Theory In this Appendix, we present some technical details about the four-body phase space integration. To be specifical, let us consider J/ψ(P ) → e − (k 1 )e + (k 2 )µ − (k 3 )µ + (k 4 ). The four-body phase-space integration can be carried out recursively, owing to its factorization property [1]. To our purpose, it is most convenient to lump the e + e − into a compound particle, and lump µ + µ − into another compound particle. The four-body phase space can then be factorized into the convolution of three two-body phase space integration measures: [dk i ](2π) 4 δ (4) (P −  6. e + e − (µ + µ − , τ + τ − ) invariant mass spectra in Υ → e + e − τ + τ − (µ + µ − τ + τ − ).
where dΠ q 1 q 2 ≡ [dq 1 ][dq 2 ](2π) 4 δ (4) (Q − q 1 − q 2 ) denotes the two-body phase-space integration measure dΠ 2 , with [dq] ≡ d 3 q (2π) 3 2Eq . In (A1), we define the momentum of the compound system as K ab ≡ k a + k b , with its invariant mass M ab ≡ K 2 ab . Since the two-body phase space dΠ 2 is Lorentz invariant, we can choose to calculate in any reference frame. We prefer to work in the center-of-mass frame for each individual two-body phase space integral: