Observation of $\psi(3686) \to \eta' e^+ e^- $

Using a data sample of $448.1 \times 10^6$ $\psi(3686)$ events collected with the BESIII detector at the BEPCII collider, we report the first observation of the electromagnetic Dalitz decay $\psi(3686) \to \eta' e^+ e^-$, with significances of 7.0$\sigma$ and 6.3$\sigma$ when reconstructing the $\eta'$ meson via its decay modes $\eta'\to\gamma \pi^+ \pi^-$ and $\eta'\to\pi^+\pi^-\eta$ ($\eta \to \gamma\gamma$), respectively. The weighted average branching fraction is determined to be $\mathcal{B}(\psi(3686) \to \eta' e^+ e^-)= (1.90 \pm 0.25 \pm 0.11) \times 10^{-6}$, where the first uncertainty is statistical and the second systematic.


Introduction
The electromagnetic (EM) Dalitz decays V → P ℓ + ℓ − , where V is a vector meson (V = ρ, ω, φ, ψ), P a pseudoscalar meson (P = π 0 , η, η ′ ) and ℓ a lepton (ℓ = e, µ), is of great interest for our understanding of both the intrinsic structure of hadrons and the fundamental mechanisms of the interactions between photons and hadrons [1]. These Dalitz decays proceed via a twobody radiative process of V decaying into P and an offshell photon, from which the lepton pair in the final state originates. The universal decay width of these Dalitz decays can be normalized to that of the corresponding radiative process V → P γ and can be parameterized as a product of the quantum electrodynamics prediction for a point-like particle and the transition form factor (TFF) F (q 2 ) at the V -P transition vertex [1], where q 2 = M 2 ℓ + ℓ − c 2 is the four-momentum transfer squared.
Knowledge of the q 2 -dependent TFF thus provides information about the EM structure arising at the V -P vertex.
EM Dalitz decays have been widely observed for light unflavored mesons, such as ω → π 0 e + e − [2,3], [5] and φ → ηe + e − [6,7]. The investigation of these decays motivated the authors of Ref. [8] to study the charmonium decays J/ψ → P ℓ + ℓ − and to calculate the branching fractions based on a monopole TFF F (q 2 ) = 1/(1 − q 2 /Λ 2 ) using a vector meson dominance model. Here Λ is an effective pole mass accounting for the overall effects from all possible resonance poles and scattering terms in the time-like kinematic region. The charmonium EM Dalitz decays J/ψ → P e + e − have been previously observed by the BESIII experiment using a data sample of 2.25 × 10 8 J/ψ events [9]. The results agree well with the theoretical predictions [8] for the P = η, η ′ cases. However, similar EM Dalitz decays have never been studied in ψ(3686) decays. The investigation of such processes will be important to understand the interaction of charmonium vector states with photons, and helpful for further studies on the ψ → V P process, including the related ρπ puzzle [10]. In this Letter, we report the first observation of the charmonium EM Dalitz decay ψ(3686) → η ′ e + e − using a data sample of 448.1×10 6 ψ(3686) events (107.0×10 6 [11] in 2009 and 341.1×10 6 [12] in 2012) collected with the BESIII detector [13]. Here, the intermediate η ′ meson is reconstructed via two decay modes, η ′ → γπ + π − and η ′ → π + π − η with η → γγ.

The BESIII experiment and Monte Carlo simulation
The BESIII detector [13] is a magnetic spectrometer operating at BEPCII, a double ring e + e − collider running at center-of-mass (c.m.) energies between 2.0 and 4.6 GeV with a peak luminosity of 1 × 10 33 cm −2 s −1 at a c.m. energy of 3.773 GeV. The cylindrical core of the BESIII detector comprises a helium-gas-based main drift chamber (MDC) to measure the momentum and the ionization energy loss (dE/dx) of charged particles, a plastic scintillator time-of-flight (TOF) system for particle identification (PID) information, a CsI(Tl) electromagnetic calorimeter (EMC) to measure photon and electron energies and a multilayer resistive plate chamber muon counter system (MUC) to identify muons. The MDC, TOF and EMC are enclosed in a superconducting solenoidal magnet providing a 1.0 T magnetic field. The geometrical acceptance is 93% of 4π for charged particles and photons. The momentum resolution is 0.5% for charged particles with transverse momentum of 1 GeV/c, and the energy resolution for photons is 2.5% (5%) at 1 GeV in the barrel (end cap) EMC.
Monte Carlo (MC) simulations are used to optimize the event selection criteria, to investigate potential backgrounds and to determine the detection efficiency. The GEANT4-based [14] simulation includes the description of geometry and material of the BE-SIII detector, the detector response, digitization models and tracking of the detector running conditions and its performance. An inclusive MC sample containing 506 × 10 6 generic ψ(3686) decays is used to study the potential backgrounds. The production of the ψ(3686) resonance is simulated by the MC generator KKMC [15], in which the beam energy spread and initial state radiation (ISR) effects are also included. The known decay modes of ψ(3686) are generated by EVTGEN [16] with branching fractions taken from the Particle Data Group (PDG) [17], while the remaining unknown decay modes are generated according to the LUNDCHARM [18] model. When generating the process ψ(3686) → η ′ e + e − , the TFF is parameterized as a monopole form factor with Λ = 3.773 GeV/c 2 . For the decay of η ′ → γπ + π − , the generator takes into account the ρ-ω interference and box anomaly [19]. The decays of η ′ → π + π − η and η → γγ are generated with a phase space model. The analysis is performed in the framework of the BESIII offline software system which takes care of the detector calibration and event reconstruction.

Data analysis
Charged tracks in BESIII are reconstructed from ionization signals of particles in the MDC. The point of closest approach of every charged track to the e + e − interaction point (IP) is required to be within ±10 cm in the beam direction and within 1 cm in the plane perpendicular to the beam direction. The polar angle θ between the direction of a charged track and that of the beam must satisfy | cos θ| < 0.93 for an effective measurement in the MDC. Four charged tracks are required with zero net charge for each candidate event. The combined information of the energy loss dE/dx and TOF is used to calculate PID confidence levels (C.L.) for the electron, pion and kaon hypotheses. Both the electron and positron require the highest PID C.L. for the electron hypothesis while the other two charged tracks are assumed to be pion candidates without any PID requirements.
Electromagnetic showers are reconstructed from clusters of energy depositions in the EMC. The shower energy of photon candidates in the EMC should be greater than 25 MeV in the barrel region (| cos θ| < 0.80) or 50 MeV in the endcap region (0.86 < | cos θ| < 0.92), whereas the showers located in the transition regions between the barrel and the endcaps are excluded due to bad reconstruction. The photon candidates are required to be separated from the extrapolated positions of any charged track by more than 10 • to exclude showers from charged particles. To suppress electronic noise and energy deposition unrelated to the event, the time at which the photon is recorded in the EMC with respect to the collision must be less than 700 ns. We require at least one photon in the decay mode η ′ → γπ + π − and at least two photons for the decay η ′ → π + π − η.
A vertex constraint is enforced on the four charged tracks π + π − e + e − to ensure they originate from the IP. To improve the resolution and suppress backgrounds, a kinematic fit with an energy-momentum constraint (4C) is performed. For events with more than the required number of photons, only the combination with the least χ 2 4C is retained. In all cases, events with χ 2 4C < 80 are kept for further analysis.
The radiative decay ψ(3686) → η ′ γ contributes as a peaking background to the distributions of the γπ + π − and γγπ + π − invariant masses (M (γπ + π − ) and M(γγπ + π − )), if the photon subsequently converts into an e + e − pair in the beam pipe or the inner wall of the MDC. The distance δ xy from the reconstructed vertex of the e + e − pair to the IP in the plane transverse to the beam axis (the x-y plane) is used to distinguish such γ conversion events from signal events [20], where δ xy = R 2 x + R 2 y and R x and R y refer to the coordinates of the reconstructed vertex position in the x and y directions. The scatter plot of R y versus R x from a simulated γ conversion MC sample ψ(3686) → η ′ γ, η ′ → γπ + π − is shown in Fig. 1 (a), where the inner and outer circles refer to the γ conversion occurs in the beam pipe and inner wall of the MDC, respectively. The distributions of δ xy for the data, γ conversion background, and signal from MC simulation are shown in Fig. 1 (b), where the two peaks around δ xy = 3 and δ xy = 6.5 cm match the positions of the beam pipe and inner wall of the MDC. From the MC study, requiring δ xy < 2 cm will remove more than 97% of the γ conver-sion background, and the number of remaining events is estimated to be 1.19 ± 0.06 and 0.43 ± 0.02 in the η ′ → γπ + π − and η ′ → π + π − η mode, respectively.
In an e + e − collider, a virtual photon can be emitted from each lepton. The interaction of these two virtual photons will produce even C-parity states such as pseudoscalar mesons, called two-photon process [21]. In the case of η ′ production, the two-photon process e + e − → e + e − η ′ leads to the same final state as signal if the outgoing e + and e − are both detected. It also contributes as a peaking background on the M (γπ + π − ) and M (γγπ + π − ) distributions. An independent ψ(3770) data sample taken at c.m. energy of 3.773 GeV, corresponding to an integrated luminosity of 2.93 fb −1 [22,23], is used to study this background. Scatter plots of the polar angle cos θ of e + and e − for the selected events from the signal MC sample and ψ(3770) data, dominated by two-photon events, are shown in Fig. 2 (a). For the signal events, in which the electron is mostly close to the positron in direction, they mainly accumulate in the diagonal band cos θ(e + ) = cos θ(e − ) in the scatter plot. For the two photon background evens, the outgoing direction of the e ± approaches its ingoing beam direction thus they mainly accumulate in the bands of cos θ(e + ) > 0.8 or cos θ(e − ) < −0.8, especially in the intersection part. The corresponding scatter plot of events from ψ(3686) data is shown in Fig. 2 (b). To suppress the background from two-photon process, cos θ(e + ) < 0.8 and cos θ(e − ) > −0.8 are further required. To estimate the number of reaming two-photon background events in the ψ(3686) data, we use ψ(3770) data as a normalization. After applying all above selection criteria, the number of survived two-photon events in ψ(3770) data is obtained by fitting the M (γπ + π − ) and M (γγπ + π − ) distributions. A scale factor f is defined as the ratio of the observed number of two-photon events N in ψ(3686) data to that in the ψ(3770) data where N , L, σ and ε refer to the observed number of two-photon events, integrated luminosity of data samples (L ψ(3686) = 668.55 pb −1 [12], L ψ(3770) = 2.93 fb −1 ), cross section and detection efficiency of two-photon process at the two c.m. energies. The details on the cross-section can be found in Ref. [21]. The detection efficiency ratios ε ψ(3686) ε ψ(3770) are determined to be 1.10 ± 0.01 and 1.19 ± 0.02 for the two modes by the simulation with generator EKHARA [24,25]. The scale factor is calculated to be 0.245 (0.265) and the normalized number of the remaining two-photon background events in the ψ(3686) data is 1.4 ± 1.7 (0.5 ± 0.4) for the decay mode η ′ → γπ + π − (η ′ → π + π − η).
After applying the above selection criteria, the studies with the inclusive MC sample indicate that the remaining background mainly arises from ψ(3686) → π + π − J/ψ, J/ψ → e + e − (γ) events, which contributes as a non-peaking background on the M (γπ + π − ) and M (γγπ + π − ) distributions. To determine the signal yield of ψ(3686) → η ′ e + e − , an unbinned maximum likelihood (ML) fit is performed to the M (γπ + π − ) and M (γγπ + π − ) distributions in the range of [0.85, 1.05] GeV/c 2 , as shown in Figs. 3 (a) and 3 (b). In the fit, the signal probability density function (PDF) is described by the signal MC shape convolved with a Gaussian function, which is used to compensate the resolution difference between data and MC simulation. The non-peaking background PDF is parameterized with a second order Chebychev polynomial function for the decay mode η ′ → γπ + π − and with an exponential function for the η ′ → π + π − η mode. The shape of the peaking background from ψ(3686) → η ′ γ due to γ conversion is derived from the MC simulation, and its magnitude is fixed to the value estimated by taking into account the corresponding branching fractions from PDG [17]. The peaking background from the twophoton process e + e − → e + e − η ′ is described using the shape obtained from ψ(3770) data and its magnitude is fixed at evaluated values. The corresponding distributions of e + e − invariant mass M (e + e − ) for the candidate events within η ′ mass region [0.93, 0.98] GeV/c 2 are shown in Figs. 3 (c) and 3 (d), where the number of signal MC events is normalized to the corresponding fitted yield. The signal MC sample generated with monopole TFF agrees well with ψ(3686) data.
The individual branching fractions for the two η ′ decay modes are calculated with where N sig is the signal yield obtained from fitting, Table 1: Signal and background yields, detection efficiency, significance and obtained branching fraction B of ψ(3686) → η ′ e + e − for η ′ → γπ + π − and η ′ → π + π − η modes. The first uncertainties of branching fractions are statistical while the second ones are systematic. 1.99 ± 0.33 ± 0.12 1.79 ± 0.38 ± 0.11 In (a) and (b), the black dots with error bars represent data, the blue solid line is the total fit result, the red dashed line shows the signal, the green dot-dashed line denotes the non-peaking background, the pink and green shaded areas indicate the peaking background from two-photon and γ conversion, respectively. In (c) and (d), the black dots with error bars represent data, the red solid and gray shaded histograms represent signal MC simulation and non-peaking background estimated from η ′ sideband, respectively, the insets show the M (e + e − ) distributions in a wider range.
N ψ(3686) = (448.1 ± 2.9) × 10 6 [12] is the total number of ψ(3686) events, B(η ′ → X) is the branching fraction of η ′ meson decaying to specific final state X and quoted from PDG [17], ǫ is the detection efficiency from signal MC simulation. The statistical significance, as determined by the ratio of maximum likelihood value and that with signal contribution set to zero, are 7.0σ and 6.3σ for the η ′ → γπ + π − and η ′ → π + π − η modes, respectively. The yields obtained from the fit, the detection efficiency, statistical significance, and the obtained branching fractions for each mode are listed in Table 1, individually.

Systematic uncertainties
Systematic uncertainties in the branching fraction measurement are summarized in Table 2. Most of them are determined by comparing the selection efficiency of control samples between data and MC simulations.
The tracking efficiency difference between data and MC simulation, both for electrons [26] and charged pions [27], is estimated to be 1% for each charged track, which results in a total systematic uncertainty 4% for both modes.
The uncertainty associated with the photon detection efficiency, derived from a control sample of J/ψ → π + π − π 0 , π 0 → γγ, is 1.5% for each photon in the endcap region and 0.5% for each photon in barrel region. The average value, weighted according to the ratio of numbers of photon in the endcap and barrel regions, is 0.6% for each photon. As a result, 0.6% is assigned as the photon uncertainty in η ′ → γπ + π − mode and 1.2% in η ′ → π + π − η mode.
The uncertainty on electron identification is studied with the control sample of radiative Bhabha scattering events e + e − → γe + e − . The average efficiency difference for electron identification between the data and MC simulation, weighted according to the polar angle and momentum distribution of signal MC samples, is determined to be 0.3% for electron and positron, individually. The average efficiency difference between data and MC simulation associated with the requirement E/p > 0.8 is estimated to be 0.2% with a similar method.
The systematic uncertainty related with the γ conversion veto criterion δ xy < 2 cm has been investigated with a control sample of J/ψ → π + π − π 0 , π 0 → γe + e − . The relative difference of efficiency associated with the γ conversion rejected criterion between data and MC simulation is 1% [9], which is taken as the systematic uncertainty.
In the 4C kinematic fit, the helix parameters of charged tracks are corrected for the signal MC samples to improve the consistency between data and MC simulation, as described in Ref. [28]. We compare the detection efficiencies obtained with and without helix parameters correction of signal MC samples. The relative change in results, 0.8% for η ′ → γπ + π − and 1.4% for η ′ → π + π − η modes, are taken as the systematic uncertainties associated with 4C kinematic fit.
The uncertainty for the η reconstruction using γγ pair is 1% based on a study of a control sample of J/ψ → ppη [29].
The uncertainty related to the RM (π + π − ) requirement is estimated by changing the selection criteria of it from 2.90 to 2.87 GeV/c 2 and from 3.20 to 3.17 GeV/c 2 for η ′ → γπ + π − and η ′ → π + π − η mode, respectively. The difference of branching fractions between the resulting and nominal requirement, 0.2% and 1.9%, are assigned as the systematic uncertainty for the two modes, respectively.
The nominal signal MC samples are generated based on the amplitude described in Ref. [8], where the parameter Λ for the monopole form factor F (q 2 ) is set to be 3.773 GeV/c 2 . Following the procedure used in Ref. [9], we adjust the Λ to a larger value 5.0 GeV/c 2 or a smaller value 3.2 GeV/c 2 and re-generate the alternative signal MC samples. The resultant largest efficiencies change, 0.9% and 0.2% for two individual η ′ decay modes, are regarded as systematic uncertainties associated with the uncertainty from the form factor.
In the nominal fit, an MC-based shape convolved with a Gaussian function is used to model the signal PDF. An alternative fit is performed in which the signal shape is described with the MC-simulated shape only. The changes of the signal yield result, 2.6% and 0.5% for the individual modes, are assigned as systematic uncertainties associated with the signal shape in the fit.
The systematic uncertainty due to non-peaking background is estimated by varying the fit range and changing its shape. In addition to the nominal fit range [0.85, 1.05] GeV/c 2 , two alternative ones are chosen by varying the edge of the fit range by ±20 MeV/c 2 . A third-order Chebychev polynomial function is selected as an alternative background shape for the η ′ → γπ + π − mode. For the η ′ → π + π − η mode, the MC shape of the major non-peaking background ψ(3686) → π + π − J/ψ, J/ψ → γe + e − is used to model the background shape. A series of alternative fits are performed for all possible combinations of fit ranges and modeling of non-peaking background. The resultant largest difference of signal yield with respective to the nominal values, 2.8% and 4.5% for each mode, are taken as the systematic uncertainties.
The uncertainty arising from peaking background due to the γ conversion process is negligible. For the two-photon process, the uncertainty associated with the scale factor is far less than the statistical uncertainties of the background events and can be ignored. We perform a series of alternative fits, varying the input normalized number of background events following a Gaussian function with a width of the statistical uncertainty. The standard deviation of the signal yields from these fit results, 1.3% and 0.7%, are taken as uncertainties for each mode.
The uncertainty from the total number of ψ(3686) events is 0.6% [12] and those of quoted branching fractions of B(η ′ → X) from PDG are 1.7% [17] for both modes.
Assuming all sources to be independent in a single mode and adding all individual contributions in quadrature, the total relative systematic uncertainties of the B(ψ(3686) → η ′ e + e − ), are determined to be 6.2% and 7.0% for the two η ′ modes, individually.

Results
The resulting B(ψ(3686) → η ′ e + e − ) from the two η ′ reconstructed modes η ′ → γπ + π − and η ′ → π + π − η with η → γγ are (1.99 ± 0.33 ± 0.12)×10 −6 and (1.79 ± 0.38 ± 0.11)×10 −6 , where the first uncertainties are statistical and second ones are systematic. The measured branching fractions from the two modes are consistent with each other within their uncertainties. Following the method described in Ref. [30], the measurements from the two modes are combined, taking into account the correlation between uncertainties among the two modes, as denoted with an asterisk in Table 2. The weighted averaged result for branching fraction of ψ(3686) → η ′ e + e − is calculated to be (1.90 ± 0.25 ± 0.11) × 10 −6 , where the first uncertainty is statistical and the second is systematic.

Summary
In summary, with a data sample of 448.1 × 10 6 ψ(3686) events collected with the BESIII detector, we observe the charmonium EM Dalitz decay ψ(3686) → η ′ e + e − for the first time by reconstructing η ′ meson via the two decay modes η ′ → γπ + π − and η ′ → π + π − η, with a statistical significance of 7.0σ and 6.3σ, respectively. The weighted average branching fraction of ψ(3686) → η ′ e + e − is measured to be (1.90 ± 0.25 ± 0.11) × 10 −6 , where the first uncertainty is statistical and second one is systematic. The observation of this process provides new information for the interaction of charmonium states with the EM field, although the statistics of current data does not allow for a precise TFF measurement.