Observation of the double Dalitz decay $\eta'\to e^+e^-e^+e^-$

Based on $(10087 \pm 44)\times10^6$ $J/\psi$ events collected with the BESIII detector at BEPCII, the double Dalitz decay $\eta'\to e^+e^-e^+e^-$ is observed for the first time via the $J/\psi\to\gamma\eta'$ decay process. The significance is found to be 5.7$\sigma$ with systematic uncertainties taken into consideration. Its branching fraction is determined to be $\mathcal{B}(\eta'\to e^+ e^- e^+ e^-) =(4.5\pm1.0(\mathrm{stat.})\pm0.5(\mathrm{sys.})) \times 10^{-6}$.


I. INTRODUCTION
The double Dalitz decays P → ℓ + ℓ − ℓ ′+ ℓ ′− , where P is a pseudoscalar meson (P = π 0 , η, or η ′ ) while ℓ and ℓ ′ are leptons (ℓ, ℓ ′ = e, µ), are expected to proceed through an intermediate state of two virtual photons. These processes are of great interest for understanding the pseudoscalar transition form factor (TFF) and the interactions between pseudoscalars and virtual photons [1,2]. These TFFs are necessary inputs to calculate the pseduoscalar-meson-pole contributions to the hadronic light-by-light scattering, which causes the second largest uncertainty in the Standard Model determi-nation of the muon anomalous magnetic moment [3][4][5][6]. Particularly, the double Dalitz decays of pseudoscalar mesons help to determine the TFFs in the small timelike momentum region, i.e. m 2 ll ≤ q 2 ≤ m 2 P , with m ll the invariant mass of the dilepton and m P the mass of the pseudoscalar meson, and thus are suitable to determine the slope of the TFFs at q 2 = 0 [3].
The pseudoscalar meson double Dalitz decays have been elucidated by many physics models, such as the hidden gauge model [2], the modified vector meson dominance model [2], the data driven approach [1], and the resonance chiral symmetric approach [7]. Experimental searches have been carried out in a variety of processes.

II. BESIII DETECTOR AND MONTE CARLO SIMULATION
The BESIII detector [13] records symmetric e + e − collisions provided by the BEPCII storage ring [14], which operates in the center-of-mass energy range between 2.0 and 4.95 GeV. BESIII has collected large data samples in this energy region [15]. The cylindrical core of the BESIII detector covers 93% of the full solid angle and consists of a helium-based multilayer drift chamber (MDC), a plastic scintillator time-of-flight system (TOF), and a CsI(Tl) electromagnetic calorimeter (EMC), which are all enclosed in a superconducting solenoidal magnet providing a 1.0 T (0.9 T in 2012) magnetic field. The solenoid is supported by an octagonal flux-return yoke with resistive plate counter muon identification modules interleaved with steel. The charged-particle momentum resolution at 1 GeV/c is 0.5%, and the dE/dx resolution is 6% for electrons from Bhabha scattering. The EMC measures photon energies with a resolution of 2.5% (5%) at 1 GeV in the barrel (end cap) region. The time resolution in the TOF barrel region is 68 ps, while that in the end cap region was initially 110 ps. The end cap TOF system was upgraded in 2015 using multigap resistive plate chamber technology, providing a time resolution of 60 ps [16].
Simulated data samples produced with geant4based [17] Monte Carlo (MC) software, which includes the geometric description of the BESIII detector and its response, are used to determine detection efficiencies and to estimate backgrounds. The simulation models the beam energy spread and initial state radiation in the e + e − annihilations with the generator KKMC [18]. The inclusive MC sample includes both the production of the J/ψ resonance and the continuum processes incorporated in KKMC. The known decay modes are modeled with evtgen [19] using branching fractions taken from the Particle Data Group (PDG) [8], and the remaining unknown charmonium decays are modeled with lundcharm [20]. Final state radiation from charged final state particles is incorporated using photos [21]. The decay J/ψ → γη ′ is generated with an angular distribution of 1 + cos 2 θγ , while the decay η ′ → e + e − e + e − is generated with a specific model based on the resonance chiral symmetric theory with SU(3) breaking [7].

III. DATA ANALYSIS
Four charged particles with zero net charge are required to be within the range of |cos θ| < 0.93, where θ is defined with respect to the z axis, which is the symmetry axis of the MDC. The distance of closest approach to the interaction point (IP) must be less than 10 cm along the z direction, and less than 1 cm in the plane transverse to the z axis. Particle identification (PID) for charged tracks combines measurements of the energy deposited in the MDC (dE/dx) and the flight time in the TOF to form likelihoods the electron [L(e)], pion [L(π)], and kaon [L(K)] hypotheses. Tracks are identified as electrons when the the electron hypotheses satisfy L(e) > 0.001 and has the highest PID likelihood among the three hypotheses, i.e. L(e) > L(π) and L(e) > L(K). To reduce background from hadrons and muons, the electron and positron candidates with higher momentum must satisfy E/cp > 0.8, where E is the deposited energy in the EMC and p is the momentum measured by the MDC [22][23][24].
Photon candidates are identified using showers in the EMC. The deposited energy of each shower must be more than 25 MeV in the barrel region (|cos θ| < 0.80) and more than 50 MeV in the end cap region (0.86 < |cos θ| < 0.92). To exclude showers that originate from charged tracks, the angle between the direction of each shower in the EMC and the direction of the closest extrapolated charged track must be greater than 10 • . To suppress electronic noise and showers unrelated to the event, the difference between the EMC time and the event start time is required to be within [0, 700] ns. At least one photon satisfying these selection criteria is required in the final state.
A vertex constraint on the four electrons is imposed to ensure they originate from the IP. To improve the resolution and suppress the background, a kinematic fit with an energy-momentum constraint (4C) is performed under the hypothesis of J/ψ → γe + e − e + e − . For a small portion of events with multiple photon candidates, all of them are looped in the 4C fit. According to the simulation study, the χ 2 γe + e − e + e − distribution of the 4C fit could identify the correct photon candidate effectively, so only the combination with the smallest χ 2 γe + e − e + e − is retained. The χ 2 γe + e − e + e − distributions from data and signal MC are in a good agreement, as shown in Fig. 1. Events with χ 2 γe + e − e + e − < 50 are kept for further analysis. To suppress background events with γπ + π − e + e − in the final state, a 4C kinematic fit is performed under the hypothesis of J/ψ → γπ + π − e + e − , and χ 2 γe + e − e + e − is required to be less than χ 2 γπ + π − e + e − . The processes J/ψ → γη ′ with η ′ → γγ, η ′ → γe + e − and η ′ → γρ 0 , ρ 0 → e + e − may contribute peaking backgrounds to the invariant mass M (e + e − e + e − ) distributions, if the photons subsequently convert into e + e − pairs in the beam pipe or the inner wall of the MDC. To distinguish the γ conversion events from signal events, the Photon Conversion Finder package [25,26] is used. Two variables are considered: (i) The distance δ xy = R 2 x + R 2 y from the reconstructed vertices of the e + e − pairs to the IP in the transverse plane. Here, R x and R y represent the coordinates of the reconstructed vertex position in the x and y directions. The distributions of δ xy for the data, γ conversion background, and signal from MC simulation are shown in Fig. 2. (ii) The angle θ eg between the momentum vector of the radiative photon and the direction from IP to conversion point. In each event, there are four possible e + e − pairs. The event with any e + e − pair satisfying δ xy ≥ 2 cm and cos θ eg ≥ −0.50 is rejected. This γ conversion veto requirement removes most of the γ conversion backgrounds with a relative detection efficiency loss less than 20%. The signal region is defined as 0.9 < M (e + e − e + e − ) < 1.0 GeV/c 2 . After applying the above selection criteria, the total detection efficiency is determined to be (12.67 ± 0.02)% and the number of survived events is 86 as shown in the Fig. 3.
The possible background contaminations are studied with an inclusive MC sample of 1.0011 × 10 10 simulated events and data samples taken at √ s = 3.08 and √ s = 3.773 GeV. Since the cross section of e + e − → γη ′ is less than 1 nb [7], the possible background from this process is negligible. For the dominant background channels, dedicated exclusive MC samples are generated to estimate their contributions. After applying all selection criteria, most backgrounds are negligible. The only remaining peaking background in the e + e − e + e − mass spectrum The distribution of δxy. The black dots with error bars represent data, the red dashed histogram shows the signal MC simulation. The blue dash-dotted and magenta dashdot-doted histograms show the simulated γ conversion samples from J/ψ → γη ′ , η ′ → γe + e − and J/ψ → γη ′ , η ′ → γγ, respectively.
After the above selection, the resulting M (e + e − e + e − ) distribution is shown in Fig. 3, where a clear η ′ signal is visible. An unbinned extended maximum likelihood fit is performed to determine the η ′ signal yield. The signal probability density function (PDF) is represented by the signal MC shape convolved with a Gaussian function to account for the resolution difference between data and MC simulation. The shape for the nonpeaking background is described by a linear function. The background yield and its PDF parameters are allowed to float in the fit. The peaking background from the γ conversion process of J/ψ → γη ′ , η ′ → γe + e − is described by the MCsimulated shape [27] with the yield fixed at 2.48.
The BF of η ′ → e + e − e + e − is calculated as where N signal = 30.1 ± 7.0 is the signal yield obtained from fitting, N J/ψ is the number of J/ψ events, ǫ is the MC-determined detection efficiency, and B(J/ψ → γη ′ ) is the BF of J/ψ → γη ′ from the PDG [8]. The significance is determined by evaluating −2ln(L 0 /L max ), where L 0 (L max ) is the smeared likelihood without (with maximum) signal events. When the systematic uncertainties are taken into account by convolving the likelihood function of the signal yield with a Gaussian function, the significance is estimated to 5.7σ  Table I lists sources of systematic uncertainties considered in the BF measurement. Most systematic uncertainties are determined from comparisons of low-background, high-statistics data samples and MC simulation samples.

IV. SYSTEMATIC UNCERTAINTIES
The PID and tracking efficiencies of electrons are determined from a control sample of radiative Bhabha scattering e + e − → γe + e − (including J/ψ → γe + e − ) corresponding to the center-of-mass energy of the J/ψ resonance. For the electron PID study, the same PID requirements as used in this analysis are applied to the control samples. Similarly, for the electron-tracking study, the same conditions on the polar angle and the distances of closest approach are used. Differences in PID (tracking) efficiencies between the data and MC simulations are obtained for each bin of a two-dimensional distribution of the momentum (transverse momentum) versus the polar angle of the electron tracks. These results are subsequently used to determine the overall weighted differences per track for PID and tracking, and the obtained PID and tracking uncertainties for electrons are 3.8% and 2.2%, respectively. The photon detection efficiency is studied with the control sample based on J/ψ → π + π − π 0 , π 0 → γγ events [28]. The difference between data and MC simulation is 0.5% (1.5%) for a photon in the EMC barrel (end cap) region. Since most photon candidates are located in the barrel region, the difference 0.5% is taken as systematic uncertainty.
In the 4C kinematic fit, the track helix parameters are corrected to reduce the differences between data and MC simulation [27,29]. The 2.7% difference between the efficiencies with and without helix parameter corrections is taken as the systematic uncertainty. The sys-tematic uncertainty related with the γ conversion veto criterion has been investigated with a control sample of e + e − → γe + e − at √ s = 3.08 GeV, where the invariant mass of e + e − is limited to a range of extremely small mass (< 10 MeV/c 2 ). The relative difference of efficiency associated with the γ conversion rejection criterion between data and MC simulation is 1.4%. Since the γ conversion veto criterion is used for two e + e − pairs in the final state, the total systematic uncertainty is 5.6%.
In the nominal signal MC model, the parameters describing η − η ′ mixing and the π 0 → γγ decay width are fixed. To estimate the systematic uncertainty due to the form-factor parametrizations, signal MC events are also generated using the global fit results from Ref. [7]. The relative difference 7.1% in the detection efficiency is taken as the uncertainty associated with the signal model.
The uncertainty due to the nonpeaking background shape is estimated by varying the PDF shape and fitting range. The fit range is changed to (0.89, 1.01) or (0.91, 0.99) GeV/c 2 . A second-order Chebyshev polynomial is selected as an alternative background shape. A series of alternative fits are performed for all combinations of fit ranges and nonpeaking background shapes to account for possible correlations. The largest relative difference of the signal yield with respect to the nominal values, 2.7%, is taken as the systematic uncertainty. In the nominal fit, the number of peaking background events from J/ψ → γη ′ , η ′ → γe + e − is fixed. To determine the systematic uncertainty, the number of peaking background events is varied within its uncertainty, and the resulting difference is negligible.
The uncertainty in the quoted BF of J/ψ → γη ′ is 1.3% [8]. The uncertainty in the number of J/ψ events is determined to be 0.4% [12].
Assuming that all systematic uncertainties in Table I are independent, the total systematic uncertainty, obtained from their quadratic sum, is 10.8%. Using a sample of about 1 × 10 10 J/ψ events collected with the BESIII detector, we observe the double Dalitz decay process η ′ → e + e − e + e − via J/ψ → γη ′ , with a significance of 5.7σ after taking the systematic uncertainty into consideration. The BF of η ′ → e + e − e + e − is measured to be (4.5 ± 1.0(stat.) ± 0.5(sys.)) × 10 −6 . The measured BF is consistent with the theoretical predictions within the uncertainties and provides new information for the studies about η ′ TFF and the interactions between η ′ and virtual photons [1,2].