Observation of $J/\psi$ decays to $e^{+}e^{-}e^{+}e^{-}$ and $e^{+}e^{-}\mu^{+}\mu^{-}$

Using a data sample of $4.481\times 10^8 \psi^\prime$ events collected with the BESIII detector, we report the first observation of the four-lepton-decays $J/\psi\to e^+e^-e^+e^-$ and $J/\psi\to e^+e^-\mu^+\mu^-$ utilizing the process $\psi^\prime\to \pi^+\pi^- J/\psi$. The branching fractions are determined to be $[5.48\pm0.31~(\rm stat)\pm0.45~(\rm syst)]\times 10^{-5}$ and $[3.53~\pm0.22~(\rm stat)\pm0.13~(\rm syst)]\times 10^{-5}$, respectively. The results are consistent with theoretical predictions. No significant signal is observed for $J/\psi\to \mu^+\mu^-\mu^+\mu^-$, and an upper limit on the branching fraction is set at $1.6\times 10^{-6}$ at the 90$\%$ confidence level. A CP asymmetry observable is constructed for the first two channels, which is measured to be $(-0.012\pm0.054\pm0.010)$ and $(0.062\pm0.059\pm0.006)$, respectively. No evidence for CP violation is observed in this process.

The standard model (SM) [1][2][3] of particle physics has been remarkably successful in explaining almost all experimental results in this field and has precisely predicted a wide variety of phenomena.However, it is clear that the model is incomplete.Many questions still remain unanswered, such as the origin of neutrino masses, the mechanism behind the observed matter-antimatter asymmetry, the nature of dark matter, etc.These facts have prompted particle physicists to put forward new theoretical frameworks, and to search experimentally for new physics beyond the SM.On the one hand, lepton flavor universality (LFU) is expected to be obeyed in SM.In recent years, however, indications for violation of LFU have been reported in semileptonic decays of the kind b → s ℓ + ℓ − .In 2014, LHCb measured the ratio of branching fractions (BFs) R K = B(B + → K + µ + µ − )/B(B + → K + e + e − ), and found a deviation from the SM prediction by 2.6σ [4,5].The measurements have continuously been updated by LHCb [6] and Belle [7].Very recently, LHCb reported their latest result with full Run I and Run II data, which deviates from SM prediction by more than 3σ [8].On the other hand, recently, the FNAL Muon g − 2 experiment [9] released their first measurement of the anomalous magnetic moment of the muon.The result confirmed the long-standing discrepancy with the SM prediction, which was previously observed by the E821 experiment at Brookhaven National Laboratory (BNL) [10].After taking into account the contribution from BNL-E821, the deviation between measurement and theoretical prediction amounts to around 4.2σ.These findings indicate there possibly exist new particles or new couplings to leptons, especially for the muon [11][12][13][14][15][16][17][18][19][20][21][22].It is therefore urgent to investigate the validity of LFU in other experiments.
In addition, to date, the strength of charge-parity (CP) violation in the quark-sector weak interaction has been found to be insufficient to account for the observed matter-antimatter asymmetry in the universe, motivating further searches for CP violation.Any new CP violation mechanism is usually constrained by the neutron electric dipole moment (nEDM) [25][26][27].Recently, Sanchez-Puertas [28] proposed a new test for CP violation in the electromagnetic decay η → µ + µ − e + e − , in which the sources of CP violation are derived from dimension-six terms in the Standard Model Effective Field Theory [29].This purely leptonic decay avoids the strong constraints from the nEDM and could be studied at the proposed η facility experiment REDTOP [30].Similarly, a test can be performed in leptonic decays of J/ψ → ℓ + 1 ℓ − 1 ℓ + 2 ℓ − 2 at BESIII, especially for the decay of J/ψ → e + e − µ + µ − .
Details about the design and performance of the BESIII detector are given in Refs.[31,38,39].Monte Carlo (MC) simulated data samples produced with a Geant4-based [40] software package, which includes the geometric description of the BESIII detector and the detector response, are used to determine the detection efficiencies and to estimate backgrounds.The simulation includes the beam energy spread and initial state radiation in the e + e − annihilations modeled with the generator kkmc [41,42].The data sample consists of (448.1 ± 2.9) × 10 6 ψ(3686) events [43] collected with the BESIII detector.Comparable amounts of inclusive MC simulated events are used to study the backgrounds from ψ(3686) decays.The production of the ψ(3686) resonance is simulated by the MC event generator kkmc [41,42].The known decay modes are generated by evtgen [44,45] with branching fractions taken from the Particle Data Group (PDG) [46], and the remaining unknown decays are generated with the lundcharm model [47,48].For the signal MC, the ψ(3686) → π + π − J/ψ channel is generated by jpipi [44,45], and the simulation of the J/ψ → ℓ + 1 ℓ − 1 ℓ + 2 ℓ − 2 decay includes the polarization of the J/ψ and the angular distributions of the final states.Since the π + π − system is dominated by the S−wave [32,49], the produced J/ψ fully inherits the state of the ψ(3686) polarization.
Charged particle tracks in the polar angle range |cos θ| < 0.93 are reconstructed from hits in the main drift chamber (MDC).Tracks with their point of closest approach to the beam line within ±10 cm of the interaction point (IP) in the beam direction, and within 1 cm in the plane perpendicular to the beam, are selected.At least six charged tracks fulfilling these criteria are required.The time-of-flight and specific energy loss (dE/dx) information are used to calculate particle identification (PID) probabilities (prob) for the electron, pion, muon and kaon hypotheses.A track is considered to be an electron if it satisfies prob (e) > prob (π) and prob (e) > prob (K).For the two channels with muons in the final states, muon candidates must satisfy prob (µ) > prob (e) and prob (µ) > prob (K), and the deposited energy in the calorimeter (E EMC ) for the muon candidate must be in the range of [0.1, 0.3] GeV.Any track not identified as an electron or muon is assigned as a pion.
To identify π + π − J/ψ candidates, two oppositely charged tracks with momentum less than 0.45 GeV/c are required and are used to calculate the mass recoiling against them, M rec.
π + π − .Using double Gaussian functions, fits are performed on the M rec.
π + π − distributions of data, as shown in Fig. 1 (a) and (c).The combinations with M rec.
π + π − within a 5σ mass window around the J/ψ peak are retained, where σ is the resolution of the fitted double Gaussian function (σ = 4.0 ± 0.2 MeV/c 2 for the J/ψ → e + e − e + e − channel; σ = 2.8 ± 0.5 MeV/c 2 for the J/ψ → e + e − µ + µ − channel).An energy-momentum constraint (4C) kinematic fit is performed on the selected For the channel J/ψ → e + e − e + e − , if more than one combination satisfies the 4C fit, the combination with the smallest χ 2 4C is retained.For the last two channels, the combination with M rec.
π + π − closest to the nominal J/ψ mass is selected.The χ 2 4C values of the candidate events are required to be less than 200.
Possible background contributions are studied with data taken at √ s = 3.773 GeV [51] and the ψ(3686) inclusive MC sample.The former indicates that the QED background is negligible.Examination of the latter with an event topology analysis tool, TopoAna [52], shows that the dominant backgrounds are from the channels J/ψ → γe + e − and J/ψ → γ FSR e + e − for the J/ψ → e + e − e + e − decay, J/ψ → γµ + µ − and J/ψ → γ FSR µ + µ − for J/ψ → e + e − µ + µ − and J/ψ → µ + µ − µ + µ − decay modes, where γ FSR is a photon from final state radiation which converts to an e + e − pair in the detector material.For the J/ψ → µ + µ − µ + µ − channel, two background events from J/ψ → π + π − π + π − in the inclusive MC sample survive the above selection criteria and are considered as peaking background.
A photon-conversion finder [53] is used to reconstruct the photon conversion point.The distance from the IP to the reconstructed conversion point of the lower momentum e + e − pair, R xy , is used to separate the signal events from the photon conversion background events.For the J/ψ → e + e − e + e − and J/ψ → e + e − µ + µ − signal events, the R xy distribution accumulates around 0 cm, because the e + e − pair comes directly from the J/ψ decay.For the photon conversion background, which is usually associated with the lower momentum e + e − pair, most of the photon conversions occur in the beam pipe and inner wall of the MDC, and the R xy distribution accumulates beyond 2 cm.However, there can be some contamination under the peak at R xy = 0 from events where another e + e − pair occurs from photon conversion.Usually the conversion pairs in these events have similar momentum with the lower momentum pair.In order to remove this background, the differences of the momenta between the two electrons and the two positrons in the final state are required to be greater than 1.0 GeV/c, i.e. ∆p = (p e ± h − p e ± l ) > 1.0 GeV/c (e ± l : e ± with lower momentum, e ± h : e ± with higher momentum).In order to determine the signal yields, unbinned maximum likelihood fits are performed on the R xy distributions for the first two decay modes and to the π + π − recoil mass (M rec. π + π − ) spectrum for the J/ψ → µ + µ − µ + µ − channel, as shown in Fig. 1 (b), (d) and (e), respectively.In the fits to the R xy distributions for the first two channels, the signal is described by signal MC shape, while the background is described by a combination of the inclusive MC, the exclusive J/ψ → γe + e − (J/ψ → γµ + µ − ) MC, and a 2 nd order polynomial function.For the J/ψ → µ + µ − µ + µ − channel, the signal is described by a MC-simulated shape convolved with a Gaussian function, while the background is described by a combination of the J/ψ → π + π − π + π − MC shape and a 1 st order polynomial function.The number of events of the peaking background from J/ψ → π + π − π + π − is fixed to the value estimated using the world average branching fraction [46].The fitted signal yields, detection efficiencies and measured branching fractions are listed in Table 1.The measured branching fractions are determined by where N sig is the number of observed signal events, N ψ(3686) is the number of ψ(3686) events [43], and ǫ is the detection efficiency determined by MC simulation.The branching fraction B(ψ(3686)) ≡ B(ψ(3686) → π + π − J/ψ) is taken from the PDG [46].To obtain a more accurate of detecting the decay of J/ψ → e + e − e + e − , the events of the signal MC sample are weighted according to the e + l momentum versus e − l momentum distribution in data.The weight factors are the ǫ i data /ǫ i MC ratios which are obtained in different momentum bins, where ǫ i data and ǫ i MC are the efficiencies of e + e − e + e − candidates in the ith bin from data and MC simulation, respectively.The resultant detection efficiency is (8.22 ± 0.02)%.π + π − distribution for the J/ψ → µ + µ − µ + µ − channel.(f) Likelihood distribution versus the branching fraction of the J/ψ → µ + µ − µ + µ − decay mode.The red arrow indicates the UL at the 90% confidence level.
The following sources of systematic uncertainty for the branching fraction measurement are considered: tracking efficiency, PID efficiency, total number of ψ(3686) events, branching fraction of ψ(3686) → π + π − J/ψ, kinematic fit, π + π − recoil mass window, fit range, signal PDF shape, background PDF shape, and the ∆p cut.
To study the systematic uncertainties related to the tracking and PID efficiencies of electrons and positrons, control samples of radiative Bhabha events are selected, i.e.
e + e − → γe + e − , at the ψ(3686) energy for both the experimental data and MC simulation.The differences between the data and MC simulation in different momentum and polar angle regions are obtained from the control samples.The systematic uncertainties Table 1.Observed signal yields, detection efficiencies, numbers of observed events with CT > 0 and CT < 0, measured branching fraction and AT values for different decay modes.The first uncertainty is statistical and the second is systematic.
In the 4C kinematic fits, the helix parameters of simulated charged tracks are corrected to reduce the discrepancy between data and MC simulation as described in Ref. [55].
The correction factors are obtained by studying a control sample of ψ(3686) → π + π − J/ψ, J/ψ → ℓ + ℓ − .The MC samples with the corrected track helix parameters are taken as the nominal ones.The differences between detection efficiencies obtained from MC samples with and without the correction are taken as the uncertainties, which are 1.0%, 0.7% and 0.9% for the three channels, respectively.
The systematic uncertainties due to the maximum likelihood fits are studied as follows: For the first two channels, the signal shape is replaced with the signal MC convolved with a gaussian function, and the uncertainties from the signal shape are evaluated to be 1.7% and 1.2%, respectively.The 2 nd order polynomial function is changed to a 3 rd order one to estimate the uncertainties from the background PDF shape, which are assigned to be 0.8% and 0.1%, respectively.To evaluate the systematic effects from the requirement on the recoil mass of π + π − , an alternative shape from the signal MC sample is used.The relative differences in the number of signal events, 0.4% and 0.1%, respectively, are taken as the systematic uncertainties.To estimate the systematic uncertainties from the fit range, the range of R xy is changed from [0.0, 8.0] to [0.0, 9.0], and the changes of the measured branching fractions are taken as the systematic uncertainties, which are evaluated to be 1.0% and 0.1%, respectively.The systematic uncertainty from the selection on ∆p is estimated by using a different requirement ((p e ± h − p e ± l ) > 1.10 GeV/c) on this variable, set to be 3.0%.To obtain reliable signal efficiency, the signal MC sample J/ψ → e + e − e + e − is weighted to match the data.To estimate the associated systematic uncertainty, the weight factors are randomly changed within one standard deviation in each bin one thousand times to re-obtain the signal efficiency.The distribution of the resulting signal efficiency is fit with a Gaussian function, and the standard deviation of 6.5% is taken as the systematic uncertainty.For the J/ψ → µ + µ − µ + µ − mode, the π + π − recoil mass window is changed from [3.04, 3.16] GeV/c 2 to [3.06, 3.16] GeV/c 2 , and 0.9% is taken as the systematic uncertainty from this source.The 1 st order polynomial function is changed to a 2 nd order one to estimate the uncertainty from the background PDF shape, which is assigned to be 0.3%.In estimating the uncertainty from the background events of J/ψ → π + π − π + π − , the fits are performed twenty thousand times with the number of the background events sampled from a Gaussian distribution based on its statistical error.The distribution of the relative difference on the signal yield is fitted by a Gaussian function, and its width of 0.2% is assigned as the systematic uncertainty from this source.The total systematic uncertainties of the BF measurement are 8.2%, 3.6%, and 4.2% for the three channels, respectively, obtained by adding the above effects quadratically.
As indicated in Fig. 1 (e), no significant signal is observed for the J/ψ → µ + µ − µ + µ − channel, so the branching fraction upper limit (UL) is determined according to Eq. 1 based on a Bayesian method [56].A series of fits of the M rec.
π + π − distribution are carried out fixing the BF at different values, and the resultant curve of likelihood values as a function of the BF is convolved with a Gaussian function to take into account the overall systematic uncertainty, as shown in Fig. 1 (f).The UL on the BF is obtained when the integral of the likelihood curve in the positive domain reaches 90% of its total value.
Following Refs.[33,34], we measure correlations between the final state leptons.In the rest frame of J/ψ, we define , where ℓ a and ℓ b are final-state leptons: ℓ a = ℓ b = e for J/ψ → e + e − e + e − , while ℓ a = e and ℓ b = µ for J/ψ → e + e − µ + µ − .In the former case, ℓ + a and ℓ − a are taken as the leptons with lower momenta, and ℓ + b those with higher momenta [33,57].A CP asymmetry observable based on CPT invariance [33,34], A T , can be constructed with the number of events of positive and negative C T values: The numbers of events with C T > 0 and C T < 0 are evaluated by a simultaneous fit to the R xy distributions of data using the same fit method as that used in the branching fraction measurements.The parameters of the signal shape are shared in the fits, while the other parameters are free.Taking into account the contributions of the same systematic uncertainties as in the branching fraction measurement, the respective numbers of events and measured asymmetries of the J/ψ → e + e − e + e − and J/ψ → e + e − µ + µ − channels are listed in Table 1.
In summary, based on a data sample of 4.481 × 10 8 ψ(3686) events collected with the BESIII detector, the decay processes of J/ψ → ℓ + 1 ℓ − 1 ℓ + 2 ℓ − 2 are investigated.The decay channels of J/ψ → e + e − e + e − and J/ψ → e + e − µ + µ − are observed for the first time, and the corresponding branching fractions are measured.The theoretical predictions based on NRQCD [24] are consistent with our measurements of the branching fractions for these two channels.No signal is observed for the J/ψ → µ + µ − µ + µ − channel, and the UL of the branching fraction is determined.No CP violation is observed.All measured results are summarized in Table 1.The ratio B eeee /B eeµµ is calculated to be (1.55 ± 0.13 ± 0.14), which agrees with theory within 1σ [24].The UL on the ratio B µµµµ /B eeee is calculated with Here N 4µ UL is the 90% UL on the number of observed events for J/ψ → µ + µ − µ + µ − decay; ǫ 4µ is the MCdetermined efficiency for the channel; N 4e is the number of events for the J/ψ → e + e − e + e − decay; ǫ 4e is the MC-determined efficiency for this channel; σ 13 = (σ stat 13 ) 2 + (σ sys 13 ) 2 = 10.9%,where σ stat 13 is the relative statistical error of N 4e (5.7%) and σ sys 13 is the total relative systematic error for J/ψ → e + e − e + e − and J/ψ → µ + µ − µ + µ − decay channels.The UL on the ratio B µµµµ /B eeee is determined to be 0.033.Similarly, the UL on the ratio B µµµµ /B eeµµ is determined to be 0.050.
The BESIII collaboration thanks the staff of BEPCII and the IHEP computing center for their strong support.The authors would like to extend thanks to Prof. Xian-Wei Kang for useful discussions and helpful advice.
This work is supported in part by National Key Research and Development Program

Fig. 1 .
Fig. 1.(a) Fit to the M rec.π + π − distribution of the ψ(3686) data using a double Gaussian function for J/ψ → e + e − e + e − channel.(b) Fit to the Rxy distribution of for J/ψ → e + e − e + e − channel.The black error bars represent data; the red solid curve indicates the overall fit; the black dashed line shows the signal; the blue dotted line is the inclusive MC background; the pink dashed line is the J/ψ → γe + e − background and the green dash-dotted line is the continuum background.(c) Fit to the M rec.π + π − distribution using a double Gaussian function for the J/ψ → e + e − µ + µ − channel.(d) Fit to the Rxy distribution for the J/ψ → e + e − µ + µ − channel.(e) Fit to the M rec.π + π − distribution for the J/ψ → µ + µ − µ + µ − channel.(f) Likelihood distribution versus the branching fraction of the J/ψ → µ + µ − µ + µ − decay mode.The red arrow indicates the UL at the 90% confidence level.
of Modern Physics, Fudan University, Shanghai 200443, People's Republic of China i Also at Harvard University, Department of Physics, Cambridge, MA, 02138, USA j Currently at: Institute of Physics and Technology, Peace Ave.54B, Ulaanbaatar 13330, Mongolia k Also at State Key Laboratory of Nuclear Physics and Technology, Peking University, Beijing 100871, People's Republic of China l School of Physics and Electronics, Hunan University, Changsha 410082, China m Also at Guangdong Provincial Key Laboratory of Nuclear Science, Institute of Quantum Matter, South China Normal University, Guangzhou 510006, China