Study of the Dalitz decay J / ψ → e + e − η

M. Ablikim, M. N. Achasov, S. Ahmed, M. Albrecht, M. Alekseev , A. Amoroso , F. F. An, Q. An, Y. Bai, O. Bakina, R. Baldini Ferroli, Y. Ban, K. Begzsuren, D. W. Bennett, J. V. Bennett, N. Berger, M. Bertani, D. Bettoni, F. Bianchi , E. Boger, I. Boyko, R. A. Briere, H. Cai, X. Cai, O. Cakir, A. Calcaterra, G. F. Cao, S. A. Cetin , J. Chai , J. F. Chang, W. L. Chang, G. Chelkov, G. Chen, H. S. Chen, J. C. Chen, M. L. Chen, P. L. Chen, S. J. Chen, X. R. Chen, Y. B. Chen, X. K. Chu, G. Cibinetto, F. Cossio , H. L. Dai, J. P. Dai, A. Dbeyssi, D. Dedovich, Z. Y. Deng, A. Denig, I. Denysenko, M. Destefanis , F. De Mori , Y. Ding, C. Dong, J. Dong, L. Y. Dong, M. Y. Dong, Z. L. Dou, S. X. Du, P. F. Duan, J. Fang, S. S. Fang, Y. Fang, R. Farinelli , L. Fava , S. Fegan, F. Feldbauer, G. Felici, C. Q. Feng, E. Fioravanti, M. Fritsch, C. D. Fu, Q. Gao, X. L. Gao, Y. Gao, Y. G. Gao, Z. Gao, B. Garillon, I. Garzia, A. Gilman, K. Goetzen, L. Gong, W. X. Gong, W. Gradl, M. Greco , L. M. Gu, M. H. Gu, Y. T. Gu, A. Q. Guo, L. B. Guo, R. P. Guo, Y. P. Guo, A. Guskov, Z. Haddadi, S. Han, X. Q. Hao, F. A. Harris, K. L. He, F. H. Heinsius, T. Held, Y. K. Heng, T. Holtmann, Z. L. Hou, H. M. Hu, J. F. Hu, T. Hu, Y. Hu, G. S. Huang, J. S. Huang, X. T. Huang, X. Z. Huang, Z. L. Huang, T. Hussain, W. Ikegami Andersson, M, Irshad, Q. Ji, Q. P. Ji, X. B. Ji, X. L. Ji, X. S. Jiang, X. Y. Jiang, J. B. Jiao, Z. Jiao, D. P. Jin, S. Jin, Y. Jin, T. Johansson, A. Julin, N. Kalantar-Nayestanaki, X. S. Kang, M. Kavatsyuk, B. C. Ke, T. Khan, A. Khoukaz, P. Kiese, R. Kliemt, L. Koch, O. B. Kolcu , B. Kopf, M. Kornicer, M. Kuemmel, M. Kuessner, A. Kupsc, M. Kurth, W. Kühn, J. S. Lange, M. Lara, P. Larin, L. Lavezzi , S. Leiber, H. Leithoff, C. Li, Cheng Li, D. M. Li, F. Li, F. Y. Li, G. Li, H. B. Li, H. J. Li, J. C. Li, J. W. Li, K. J. Li, Kang Li, Ke Li, Lei Li, P. L. Li, P. R. Li, Q. Y. Li, T. Li, W. D. Li, W. G. Li, X. L. Li, X. N. Li, X. Q. Li, Z. B. Li, H. Liang, Y. F. Liang, Y. T. Liang, G. R. Liao, L. Z. Liao, J. Libby, C. X. Lin, D. X. Lin, B. Liu, B. J. Liu, C. X. Liu, D. Liu, D. Y. Liu, F. H. Liu, Fang Liu, Feng Liu, H. B. Liu, H. L Liu, H. M. Liu, Huanhuan Liu, Huihui Liu, J. B. Liu, J. Y. Liu, K. Liu, K. Y. Liu, Ke Liu, L. D. Liu, Q. Liu, S. B. Liu, X. Liu, Y. B. Liu, Z. A. Liu, Zhiqing Liu, Y. F. Long, X. C. Lou, H. J. Lu, J. G. Lu, Y. Lu, Y. P. Lu, C. L. Luo, M. X. Luo, X. L. Luo, S. Lusso , X. R. Lyu, F. C. Ma, H. L. Ma, L. L. Ma, M. M. Ma, Q. M. Ma, X. N. Ma, X. Y. Ma, Y. M. Ma, F. E. Maas, M. Maggiora , Q. A. Malik, A. Mangoni , Y. J. Mao, Z. P. Mao, S. Marcello , Z. X. Meng, J. G. Messchendorp, G. Mezzadri, J. Min, T. J. Min, R. E. Mitchell, X. H. Mo, Y. J. Mo, C. Morales Morales, G. Morello, N. Yu. Muchnoi, H. Muramatsu, A. Mustafa, S. Nakhoul, Y. Nefedov, F. Nerling, I. B. Nikolaev, Z. Ning, S. Nisar, S. L. Niu, X. Y. Niu, S. L. Olsen , Q. Ouyang, S. Pacetti , Y. Pan, M. Papenbrock, P. Patteri, M. Pelizaeus, J. Pellegrino , H. P. Peng, Z. Y. Peng, K. Peters, J. Pettersson, J. L. Ping, R. G. Ping, A. Pitka, R. Poling, V. Prasad, H. R. Qi, M. Qi, T. Y. Qi, S. Qian, C. F. Qiao, N. Qin, X. S. Qin, Z. H. Qin, J. F. Qiu, K. H. Rashid, C. F. Redmer, M. Richter, M. Ripka, M. Rolo , G. Rong, Ch. Rosner, X. D. Ruan, A. Sarantsev, M. Savrié , C. Schnier, K. Schoenning, W. Shan, X. Y. Shan, M. Shao, C. P. Shen, P. X. Shen, X. Y. Shen, H. Y. Sheng, X. Shi, J. J. Song, W. M. Song, X. Y. Song, S. Sosio , C. Sowa, S. Spataro , G. X. Sun, J. F. Sun, L. Sun, S. S. Sun, X. H. Sun, Y. J. Sun, Y. K Sun, Y. Z. Sun, Z. J. Sun, Z. T. Sun, Y. T Tan, C. J. Tang, G. Y. Tang, X. Tang, I. Tapan , M. Tiemens, B. Tsednee, I. Uman, G. S. Varner, B. Wang, B. L. Wang, C. W. Wang, D. Wang, D. Y. Wang, Dan Wang, K. Wang, L. L. Wang, L. S. Wang, M. Wang, Meng Wang, P. Wang, P. L. Wang, W. P. Wang, X. F. Wang, Y. Wang, Y. F. Wang, Y. Q. Wang, Z. Wang, Z. G. Wang, Z. Y. Wang, Zongyuan Wang, T. Weber, D. H. Wei, P. Weidenkaff, S. P. Wen, U. Wiedner, M. Wolke, L. H. Wu, L. J. Wu, Z. Wu, L. Xia, X. Xia, Y. Xia, D. Xiao, Y. J. Xiao, Z. J. Xiao, Y. G. Xie, Y. H. Xie, X. A. Xiong, Q. L. Xiu, G. F. Xu, J. J. Xu, L. Xu, Q. J. Xu, Q. N. Xu, X. P. Xu, F. Yan, L. Yan , W. B. Yan, W. C. Yan, Y. H. Yan, H. J. Yang, H. X. Yang, L. Yang, S. L. Yang, Y. H. Yang, Y. X. Yang, Yifan Yang, M. Ye, M. H. Ye, J. H. Yin, Z. Y. You, B. X. Yu, C. X. Yu, J. S. Yu, C. Z. Yuan, Y. Yuan, A. Yuncu, A. A. Zafar, A. Zallo, Y. Zeng, Z. Zeng, B. X. Zhang, B. Y. Zhang, C. C. Zhang, D. H. Zhang, H. H. Zhang, H. Y. Zhang, J. Zhang, J. L. Zhang, J. Q. Zhang, J. W. Zhang, J. Y. Zhang, J. Z. Zhang, K. Zhang, L. Zhang, S. F. Zhang, T. J. Zhang, X. Y. Zhang, Y. Zhang, Y. H. Zhang, Y. T. Zhang, Yang Zhang, Yao Zhang, Yu Zhang, Z. H. Zhang, Z. P. Zhang, Z. Y. Zhang, G. Zhao, J. W. Zhao, J. Y. Zhao, J. Z. Zhao, Lei Zhao, Ling Zhao, M. G. Zhao, Q. Zhao, S. J. Zhao, T. C. Zhao, Y. B. Zhao, Z. G. Zhao, A. Zhemchugov, B. Zheng, J. P. Zheng, W. J. Zheng, Y. H. Zheng, B. Zhong, L. Zhou, Q. Zhou, X. Zhou, X. K. Zhou, X. R. Zhou, X. Y. Zhou, A. N. Zhu, J. Zhu, J. Zhu, K. Zhu, K. J. Zhu, S. Zhu, S. H. Zhu, X. L. Zhu, Y. C. Zhu, Y. S. Zhu, Z. A. Zhu, J. Zhuang, B. S. Zou, J. H. Zou

We study the electromagnetic Dalitz decay J/ψ → e + e − η and search for di-electron decays of a dark gauge boson (γ ) in J/ψ → γ η with the two η decay modes η → γγ and η → π + π − π 0 using (1310.6 ± 7.0) × 10 6 J/ψ events collected with the BESIII detector. The branching fraction of J/ψ → e + e − η is measured to be (1.42 ± 0.04(stat) ± 0.07(syst)) × 10 −5 , with a precision that is improved by a factor of 1.5 over the previous BESIII measurement. The corresponding di-electron invariant mass dependent modulus square of the transition form factor is explored for the first time, and the pole mass is determined to be Λ = 2.56 ± 0.04(stat) ± 0.13(syst) GeV/c 2 . We find no evidence of γ production and set 90% confidence level upper limits on the product branching fraction B(J/ψ → γ η) × B(γ → e + e − ) as well as the kinetic mixing strength between the Standard Model photon and γ in the mass range of 0.01 ≤ m γ ≤ 2.4 GeV/c 2 .

I. INTRODUCTION
The study of electromagnetic (EM) Dalitz decays of a vector meson (V = ρ, ω, φ, J/ψ) into a pseudoscalar meson (P = π 0 , η, η ) and a lepton-pair, V → + − P ( = e, µ), plays an important role in revealing the structure of hadrons and the interaction mechanism between photons and hadrons [1]. These decays proceed via V → γ * P in which the virtual photon γ * subsequently converts into a lepton pair. Assuming the mesons to be pointlike particles, the di-lepton invariant mass (m + − ) dependent decay rate of V → + − P can be described by quantum electrodynamics (QED) [2]. Any deviation from the QED prediction, caused by the dynamics of the EM structure arising at the V → P transition vertex, is formally described by a transition form factor (TFF) [1]. The dependence of the differential decay rate of V → P + − on the four-momentum transfer squared q 2 = m 2 + − is parameterized as [1] dΓ(V → P + − ) dqΓ(V → P γ) = 2α 3πq 1 − 4m 2 q 2 where QED(q 2 ) is the QED predicted q 2 dependent decay rate, α is the fine structure constant, F V P (q 2 ) is the q 2 dependent TFF, and m , m V and m P are the masses of leptons, V and P mesons, respectively. The TFFs of light mesons contribute to the hadronic light-by-light (HLBL) corrections [3] to the theoretical determination of the muon anomalous magnetic moment, a µ = (g µ − 2)/2, which provides a low-energy test of the completeness of the Standard Model (SM) [4,5]. Experimentally, it can directly be accessible by comparing the m + − spectrum of Dalitz decays V → + − P with that of the point-like QED prediction [1]. Within the vector meson dominance model (VMD) [6], the TFF is mainly governed by the coupling of the γ * to the V meson via an intermediate vector (V ) meson in the timelike region, and is commonly expressed as a multipole function in the charmonium mass region [7] where N is a normalization constant ensuring that F V P (0) = 1, V denotes the intermediate resonances ρ, ω, φ, and charmonium vector mesons, m V , Γ V and A V are the corresponding masses, widths and the coupling constants. The contribution of vector mesons with masses above (m V − m P ), non-resonant contribution, is often represented as where Λ is an effective pole mass. The inverse square value Λ −2 reflects the slope of the TFF at q 2 = 0. The EM Dalitz decays of light unflavored vector mesons ρ, ω and φ are well established by several collider and nuclear physics experiments [8][9][10][11][12]. The BE-SIII collaboration reported the first measurements of the branching fractions of J/ψ → e + e − P and the TFF of J/ψ → e + e − η using a data sample of 225 million J/ψ events [13]. The results agree well with the VMD predictions based on a simple pole approximation [14] within the statistical uncertainties. BESIII has recently accumulated 5 times more statistics of the J/ψ data-set [15], which can be used to improve the precision of these measurements and enable measurement of the TFFs of J/ψ → e + e − P .
The EM Dalitz decays can also be utilized to search for a hypothetical dark photon, γ , via the decay chain J/ψ → γ P , γ → + − [14,16]. The γ is a new type of force carrier in the simplest scenario of an Abelian U (1) interaction under which dark matter particles are considered to be charged [17][18][19]. A γ with mass below twice the proton mass can explain the features of the electron/positron excess observed by the cosmic ray experiments [20][21][22][23]. A dark photon with such a low mass can also explain the presently observed deviation of a µ up to the level of (3 − 4)σ between the measurement and SM prediction [19]. The γ couples with the SM photon through its kinetic mixing with the SM hypercharge field [24]. The coupling strength between the dark sector and the SM, , is parameterized as 2 = α /α, where α is the fine structure constant in the dark sector. A series of experiments have reported null results in γ searches, including the a µ favored region, and have constrained the values as a function of γ mass to be below 10 −3 [25][26][27]. More experimental information about the γ searches via new decay modes, such as J/ψ → γ P , might be helpful to understand some other possible scenarios of the γ coupling to the SM [28]. In this paper, we present a study of the EM Dalitz decays J/ψ → e + e − η and search for di-electron decays of a dark photon through J/ψ → γ η using (1310.6 ± 7.0) × 10 6 J/ψ events collected with the BESIII detector [15].

II. THE BESIII EXPERIMENT AND MONTE CARLO SIMULATION
The BESIII detector is a general purpose spectrometer containing four major detector sub-components with a geometric acceptance of 93% of the total solid angle as described in Ref. [29]. A helium-based (60% He, 40% C 3 H 8 ) multi-layer drift chamber (MDC), which contains 43 layers and operates in a 1.0 T (0.9 T) solenoidal magnetic field for the 2009 (2012) J/ψ data, is used to measure the momentum of the charged particles. Charged particle identification (PID) is based on the energy loss (dE/dx) in the tracking system and the time-of-flight (TOF) measured by a scintillation based TOF detector containing one barrel and two end-caps. A CsI(Tl) based electromagnetic calorimeter (EMC) is used to measure the energies of photons and electrons, while a muon counter containing nine (eight) layers of resistive plate chamber counters interleaved with steel in the barrel (end-cap) region is used for muon identification.
Monte Carlo (MC) simulated events are used to optimize the event selection criteria, to study the detection acceptance and to understand the potential backgrounds. The Geant4 [30] based simulation package contains the information about the detector geometry and material description, the detector response and signal digitization models, as well as the records of time dependent detector running conditions and performance. An MC sample of 1.225 billion inclusive J/ψ decays is generated for the background studies with the EvtGen generator [31] for the known J/ψ decay modes with the branching fractions set to their world average value taken from Ref. [32], and the lundcharm package [33] for the remaining unknown J/ψ decay modes. The kkmc event generator package [34] is used to simulate the production of the J/ψ resonance via e + e − annihilation, incorporating the effects of the beam energy spread and initial-state-radiation (ISR).
The angular distribution of the decay J/ψ → e + e − η is simulated according to a combined formula of Eqs. (4) and (6) of Ref. [14], where the dependence on the cosine of the η meson polar angle in the J/ψ rest frame (cos θ η ) is parameterized by (1 + α θ cos 2 θ η ) with α θ = 1.0 measured from the data as described in Section III B to take into account of the J/ψ polarization state in the e + e − annihilation system, and the TFF is assumed to follow Eq. (3) with Λ = 2.56 GeV/c 2 measured in this analysis also described in Sec. III B. The J/ψ → γ η decay is modeled by a helicity amplitude model and γ → e + e − decay by a model of a vector meson decaying to a leptonpair [31].

III. DATA ANALYSIS
In this analysis, the η meson candidates are reconstructed using the dominant decay modes η → γγ and η → π + π − π 0 , where the π 0 meson is reconstructed with a γγ pair. We select events of interest with two (four) charged tracks with zero net charge in the η → γγ (η → π + π − π 0 ) decay and at least two good photon candidates. The charged tracks are required to be measured in the active region of the MDC, | cos θ| < 0.93, where θ is the polar angle of the charged tracks. They must also have the points of closest approach to the beam line within ±10.0 cm from the interaction point in the beam direction and within 1.0 cm in the plane perpendicular to the beam. A PID algorithm, based on energy loss dE/dx in the MDC, TOF information, and energy deposited in the EMC, is performed to identify electrons. An electron-positron pair is required for the selected events. In the decay J/ψ → e + e − η, η → π + π − π 0 , the additional two charged tracks are assumed to be π candidates without any PID requirement.
The photon candidates are reconstructed from the clusters of energy deposits in the EMC that are separated from the extrapolated positions of any charged tracks by more than 10 degrees. The energy of each photon candidate is required to be larger than 25 MeV in the EMC barrel region (| cos θ γ | < 0.8) or 50 MeV in the EMC end-cap regions (0.86 < | cos θ γ | < 0.92), where θ γ is the polar angle of the photon. To improve the reconstruction efficiency and energy resolution, the energy deposited in nearby TOF counter is taken into account. The photons reconstructed poorly in the transition region between the barrel and the end-caps are discarded. The EMC timing is required to be within the range of [0, 700] ns to suppress electronic noise and energy deposits unrelated to the event.
The selected charged tracks are constrained to originate from a common vertex point by requiring a successful vertex fit. In order to improve the resolution and further suppress the background, a four-constraint (4C) kinematic fit that imposes overall momentum and energy conservation is implemented for the selected charged tracks and additional two photons under the hypothesis of J/ψ → e + e − (π + π − )γγ. The chi-square of the kinematic fit, χ 2 4C , is required to be less than 100, which rejects about 30% of the background events with a loss of the 10% of the signal events. If there are more than two good photons in an event, we try all the γγ combinations, and the one with the least χ 2 4C is chosen. The kinematic variables after the 4C kinematic fit are used in the further analysis. In the decay mode η → π + π − π 0 , the π 0 candidate is reconstructed with two selected photons by requiring m γγ within the range of [0.08, 0.16] GeV/c 2 . The η candidate is reconstructed with the selected γγ or π + π − π 0 , respectively, and the corresponding masses (m γγ and m π + π − π 0 ) are required to be within the range [0.45, 0.65] GeV/c 2 .
With the above selection criteria, the peaking background, which contains an η signal in the final state, is dominated by the events of the radiative decay J/ψ → γη followed by the conversion of the radiative photon into an e + e − pair in the detector material. In order to suppress this background, a photon-conversion finder algorithm [35] is exploited to reconstruct the photonconversion vertex point. The distance from the conversion vertex point to the origin in the x − y plane, δ xy = R 2 x + R 2 y , is used to separate the signal from the gamma conversion events, where R x and R y refer to the coordinates of the reconstructed vertex point along the x and y directions, respectively. The scatter plot of R y versus R x from the simulated γ conversion background MC sample J/ψ → γη and the signal MC sample J/ψ → e + e − η is shown in Fig. 1 (a), where the circles with radius of 3.5 cm and 6.5 cm correspond to the positions of the beam pipe and inner wall of the MDC, respectively. The corresponding distributions of δ xy from the signal and background MC samples, as well as data events are shown in Fig. 1 (b). The δ xy is then required to be less than 2 cm to remove around 98% of the γ conversion events from J/ψ → γη decay, while retaining about 80% of the signal events J/ψ → e + e − η.
In the decay mode η → γγ, the background is dominated by the non-peaking background from the QED processes e + e − → e + e − γ(γ) and e + e − → 3γ in which one of the photons converts into an e + e − pair. Since the η meson decays isotropically, the cosine of the helicity angle (cos θ heli ), defined as the angle between the direction of one of the photons and J/ψ direction in the η rest frame, is expected to be uniformly distributed for signal events and to peak near cos θ heli = ±1 for the background from QED processes. Thus a requirement | cos θ heli | < 0.9 is implemented in the decay mode η → γγ to suppress the non-peaking QED background.
After applying the above selection criteria, the distribution of the di-electron invariant mass m e + e − of surviving events (within the η signal region [0.51, 0.58] GeV/c 2 ) is shown in Fig. 2. Besides the EM Dalitz decay of interest, J/ψ → e + e − η, small signals of J/ψ → V η (V = ρ, ω, φ) with V subsequently decaying into e + e − pair are observed. Detailed MC studies indicate that the remaining peaking background is dominated by J/ψ → γη with γ converting into an e + e − pair, which accumulates in the low region of the m e + e − distribution. There are also small contributions of the peaking background of J/ψ → ρ/ωη with ρ/ω subsequently decaying into a π + π − pair and the direct three body decay J/ψ → π + π − η, in which the π + π − are mis-identified as an e + e − pair. The nonpeaking background, which is smoothly distributed in the high mass region of the m e + e − distribution, is almost negligible in the decay mode η → π + π − π 0 , but sizable in the decay mode η → γγ dominated by the radiative Bhabha e + e − → γe + e − process. The distributions of signal and individual background components are also depicted in Fig. 2. Here, the peaking backgrounds are estimated with the MC simulation normalized according to the branching fraction quoted from the PDG [32]; the three body decay J/ψ → π + π − η is simulated in accordance with the amplitude of J/ψ → π + π − η [36]; the non-peaking backgrounds are estimated with the events of data in the η sideband regions, which are defined as , and combined data of MC and side-band (cyan) for the decay modes (a) η → π + π − π 0 and (b) η → γγ.
A. Branching fraction measurement for the EM Dalitz decays J/ψ → e + e − η In order to suppress the peaking background from J/ψ → V η with meson V decaying into either the e + e − or the π + π − final state, the candidate events within regions of 0.65 < m e + e − < 0.90 GeV/c 2 or 0.96 < m e + e − < 1.08 GeV/c 2 are discarded. The number of remaining peaking background events, estimated by the MC simulation, for both η decay modes after this requirement is summarized in Table I. In the decay mode η → γγ, a sizable non-peaking background, which is dominated by the radiative Bhabha events e + e − → γe + e − and smoothly distributed in the high region of the m e + e − distribution, is suppressed by applying the further requirement p e ± < 1.45 GeV/c for m e + e − > 0.5 GeV/c 2 , where p e ± is the momentum of the e ± charged tracks. Other sources of peaking background are negligible in both η decay modes.
To determine the signal yields, we perform an unbinned Decay process η → γγ η → π + π − π 0 J/ψ → ρη, ρ → π + π − 2.3 0.6 extended maximum likelihood (ML) fit to the m γγ and m π + π − π 0 distributions, individually. In the fit, the probability density function (PDF) of the η signal is described with the corresponding signal MC simulated shape convolved with a Gaussian function with parameters that are left free during the fit to take into account the resolution difference between the data and MC simulation. The shape of the non-peaking background is described by a first order Chebyshev polynomial function with free parameters in the fit. The shape of the peaking background is described by that of MC simulation of the background J/ψ → γη, and the corresponding expected number of events is fixed during the fit. The ML fit yields N sig = 594.9 ± 25.3 and 1877.2 ± 76.1 events for the decay modes η → π + π − π 0 and η → γγ, respectively. The corresponding fit curves are shown in Fig. 3. The statistical uncertainty of the extracted signal yield in η → γγ slightly degrades compared to the previous BESIII measurement [13]; this is because the ML fit to the m γγ distribution is now performed in the full m e + e − range instead of m e + e − < 0.5 GeV/c 2 range as required by the previous measurement to avoid the large contamination from the radiative Bhabha background.

B. Transition form factor
Due to large contamination from the radiative Bhabha process in the high m e + e − region in the decay mode η → γγ, only the events from the η → π + π − π 0 decay are used for the TFF study. The vicinities of ω and φ in the m e + e − distribution are also explored, and the resonant contribution of J/ψ → V η, V → e + e − is considered as a signal in the TFF measurement. Due to limited statistics in the high mass region as seen in Fig. 2, the TFF is extracted bin-by-bin from the efficiency and branching fractions corrected signal yields for the bin sizes of 0.10 GeV/c 2 between 2m e < m e + e − < 1.10 GeV/c 2 , 0.12 GeV/c 2 between 1.10 < m e + e − < 1.34 GeV/c 2 , 0.14 GeV/c 2 between 1.34 < m e + e − < 1.90 GeV/c 2 , 0.16 GeV/c 2 between 1.90 < m e + e − < 2.06 GeV/c 2 and 0.17 GeV/c 2 in the remaining m e + e − regions with a total 20 bins, where m e is the mass of the electron. The signal yield in each bin of m e + e − is extracted by performing ML fits to the m π + π − π 0 distribution as described in Sec. III A. The peaking background contribution from the J/ψ → γη exists only in the first and second bins of m e + e − . All the peaking background contributions including J/ψ → γη, J/ψ → ρ/ωη (ρ/ω → π + π − ) and J/ψ → π + π − η are estimated with the MC simulation and subtracted from the extracted signal yield from the ML fit in each bin.
The signal efficiency for the TFF measurement is calculated by the signal MC sample generated according to the method discussed in Sec. II, but with a constant TFF of F J/ψη (m 2 e + e − ) = 1.0. The angular distribution parameter α θ , used as an input parameter in this signal MC simulation, is evaluated after extracting the cos θ η dependent signal yield with a step size of 0.2 between −0.9 < cos θ η < 0.9 using a similar procedure of the ML fit mentioned above. Figure 4 shows the efficiency corrected signal yield versus cos θ η data and a fit with N (1 + α θ cos 2 θ η ), where N is a normalization constant and the efficiency for this study is evaluated after generating the simulated signal MC events with a flat distribution in cos θ η . The angular distribution parameter α θ is determined to be 1.0 +0.0 −0.2 with a condition of 0 ≤ |α θ | ≤ 1.0 to satisfy the theoretical constraints [37]. Table II summarizes the background subtracted fitted N i sig and branching fractions B(J/ψ → e + e − η) i for all 20 bins. The branching fraction of J/ψ → e + e − η is computed using where E is the signal selection efficiency, B(η → F ) is the branching fraction of subsequent η decays taken from the η θ cos FIG. 4. Fit to the efficiency corrected signal yield versus cos θη for data in the η → π + π − π 0 decay mode. The black dots with error bar are data, which include both statistical and systematic uncertainties, and the solid red curve shows the fit results.
PDG [32] and N J/ψ = (1310.6 ± 7.0) × 10 −6 is the number of J/ψ events from Ref. [15]. The distribution of B(J/ψ → e + e − η) i normalized to the m e + e − bin size superimposed with the QED predicted branching fractions, computed using the formula of Eq. (1), is shown in Fig. 5. The dark photon search is performed in the full m e + e − spectrum using the surviving event candidates within the TABLE II. Fitted N i sig , differential branching fraction B(J/ψ → e + e − η) i and the TFF |F (q 2 )| 2 , described in Section V, for all 20 bins. The first uncertainty is statistical and the second systematic discussed in Section IV.  We search for the γ signal in steps of 2 MeV/c 2 in the m e + e − distribution ranging from 10 MeV/c 2 to 2.4 GeV/c 2 excluding the vicinities of the ω and φ signals. The parameters of the signal PDF are kept fixed, while the parameters of the background PDF, the number of signal events (N sig ) and background events are determined by the fit. In order to address the fit problem associated with low-statistics, a lower bound of N sig is imposed with a requirement that the total signal and background PDF remains non-negative [38]. The statistical signal significance is computed as S = sign(N sig ) −2 ln(L 0 /L max ), where L max and L 0 are the likelihood values when N sig is left free and fixed at 0, respectively, and sign(N sig ) is the sign of N sig . The plots of N sig and signal significance as a function of m γ for both the η decay modes are shown in Fig. 6. The largest local significance is 2.92σ at m γ = 0.590 GeV/c 2 in the η → π + π − π 0 decay and 2.98σ at m γ = 2.144 GeV/c 2 in the η → γγ decay, which are less than 3σ. Therefore, we conclude that no evidence of γ production is found in both the η decay modes.  Table III summarizes the sources of additive and multiplicative systematic uncertainties considered in this analysis, where the additive systematic uncertainties arise from the fit procedure including the signal and background modeling, as well as the bias of the fit procedure.

IV. SYSTEMATIC UNCERTAINTY
The multiplicative systematic uncertainty arises from the systematic uncertainty on the number of J/ψ events, the branching fractions in the cascade decay and the event reconstruction and selection efficiencies.
In the measurements of the branching fraction of J/ψ → e + e − η, the signal yields are determined by fitting the corresponding m γγ and m π + π − π 0 distributions, while in the TFF studies, the signal yields are extracted with the same fit procedure in the m π + π − π 0 distribution only in different m e + e − bins. The uncertainty associated with the signal model in the fit is studied by replacing the corresponding PDF to be the sum of two CB functions convolved with a Gaussian function, where the parameters of the CB functions are extracted from fits to the signal MC samples. The uncertainty associated with the peaking background is studied by varying its expected number of events within ±1σ of the uncertainties in the fit, and observed to be negligible. The uncertainty associated with the non-peaking background is studied by replacing the corresponding PDF to be a second order Chebyshev polynomial function in the fit. In the fit to search for the γ boson, the signal is modeled with the sum of two CB function whose parameters are extracted and extrapolated from the simulated MC samples at 27 m γ points. The corresponding uncertainty is studied by changing the parameters of the CB functions within ±1σ of their uncertainties, taking into account the correlation between the different parameters. The uncertainty due to the background model is studied by changing the order of polynomial functions in the fit. The changes in the signal yields due to the PDF parameters are considered as the uncertainties. To validate the reliability of fits, we produce a large number of pseudo-experiments, which are of the same statistics of data, and perform the same fit procedure in each pseudo-experiment. The resultant average difference between the input and output signal yields is found to be very small and considered as one of the systematic uncertainties.
The tracking efficiency for charged pions is studied with the control sample of J/ψ → π + π − π 0 [39]. The difference between data and MC simulation is found to be 1%, and is considered as the systematic uncertainty of the charged pion. The efficiencies of tracking and PID for e ± is explored with the control sample of radiative Bhabha events e + e − → γe + e − in 2-dimensional bins of momentum versus polar angle. The resultant average differences on efficiency between data and MC simulation, 1.2% for tracking and 0.6% for the PID, weighted according to the momentum and polar angle distribution of the MC samples, are considered as the systematic uncertainties.
The photon reconstruction efficiency is studied with a control sample of e + e − → γµ + µ − , in which the momentum of the ISR photon is inferred from the four-momenta of the µ + µ − pair [40]. The difference in the efficiency between data and MC simulation is smaller than 1%, which is taken as the systematic uncertainty. In the decay mode η → π + π − π 0 , the uncertainty related with the π 0 mass TABLE III. Summary of systematic uncertainties. The systematic uncertainties correlated between the decay modes η → π + π − π 0 and η → γγ are denoted by asterisks. Here 'Negl.' means negligible, and '-' means the corresponding source of systematic uncertainty is not applicable in a particular decay process.
Additive systematic uncertainties (events) window requirement is studied with a high statistics control sample of J/ψ → ppπ 0 and is assigned to be 1%. In the γ search, the uncertainty associated with the η mass window requirement is studied with a control sample of J/ψ → ppη, and is assigned to be 1%, too.
The uncertainty associated with the 4C kinematic fit is explored by utilizing a control sample of J/ψ → π + π − π 0 in which the π 0 dominant decay modes of π 0 → γγ and π 0 → γe + e − are utilized to mimic the J/ψ → e + e − η signal with subsequent decay modes η → γγ and η → π + π − π 0 , respectively. The relative difference in efficiencies between data and MC simulation in the corresponding control samples is observed to be up to the level of 0.9%, and considered as the systematic uncertainty.
The control sample of J/ψ → π + π − π 0 , π 0 → γe + e − is also utilized to evaluate the systematic uncertainty for the δ xy < 2 cm requirement used to suppress the γconversion background. The simulated MC events for π 0 → γe + e − are generated with a simple monopole ap-proximation TFF, F (m 2 e + e − ) = 1+a π m 2 e + e − /m 2 π 0 , where m π 0 is the nominal π 0 mass and a π = 0.032 ± 0.004 is the slope parameter [32]. We extract the π 0 → γe + e − signal from the data by performing a ML fit to the m e + e − distribution before and after the selection of δ xy < 2 cm requirement. The corresponding differences in efficiencies, 1.0% in the measurement of branching fraction of J/ψ → e + e − η and (0.0 − 1.5)% depending on m e + e − in the TFF measurement and γ search, are taken as the systematic uncertainties.
We similarly utilize the control sample J/ψ → π + π − π 0 , π 0 → γγ to evaluate the systematic uncertainty due to the photon helicity angle requirement | cos θ heli | < 0.9 in the η → γγ decay. The background in this control sample, π 0 → γe + e − , has a flat shape in m γγ , and is eliminated by performing a ML fit to the m γγ distribution. The uncertainty is evaluated to be up to the level of 1.9% by comparing the efficiencies between the data and MC simulation, where the efficiency is the ratio of signal yields with and without this requirement applied. We extract the signal yield in η → γγ decay by varying the requirement of e ± momentum within one standard deviation of its statistical uncertainty, and one of the largest values of the relative difference between the signal yields is considered as the systematic uncertainty and found to be negligible.
In the branching fraction measurement of the J/ψ → e + e − η Dalitz decay, the signal MC samples used to evaluate the detection efficiency are generated by following the TFF of Eq. (3) with measured Λ value of 2.56 GeV/c 2 as described in Sec. III. Two alternative MC samples with values of the pole mass Λ differing by ±1σ are generated, and the resulting largest relative difference in efficiencies, 1.5% for the decay mode η → π + π − π 0 and 1.0% for the decay mode η → γγ, are considered as the systematic uncertainty.
The systematic uncertainties associated with the decay branching fractions of η → π + π − π 0 and η → γγ are taken from the PDG [32]. In the measurements of TFF and coupling strength between the dark sector and the SM , the related uncertainty associated with the branching fraction of J/ψ → γη is taken from the PDG [32], and in the measurement, the uncertainty of the theoretical branching fraction B(γ → e + e − ), dominated by the uncertainty of the R value [32], varies in the range of (0 − 14)% depending upon the m γ [41]. The uncertainty of the number of J/ψ events is determined to be 0.5% using the inclusive hadronic events of the J/ψ decays [15].
The TFF in each m e + e − bin is determined by dividing B(J/ψ → e + e − η) i by the integrated QED prediction in each m e + e − interval (Table II). Figure 7 shows a plot of the resultant TFF versus m e + e − . A chi-square fit to the TFF versus m e + e − data is performed using a modified multipole function of Eq. (2), in which the contributions of the ρ resonance and non-resonance are included, and the interference is neglected: where the mass and width of the ρ resonance are fixed to the values in the PDG [32]. The statistical uncertainties and uncorrelated systematic uncertainty (between the different m e + e − bins) are considered when building the chi-square function. The fit curve is depicted in Fig. 7, too. The statistical significance of the ρ signal is 4.3σ estimated with the change of chi-square values with and without ρ signal included in the fit. We fit the TFF of data once again by including the interference between ρ and non-resonant components. The resultant change on Λ, 0.02 GeV/c 2 , is taken to be one of the systematic uncertainties. We also fit the TFF of data without including the systematic uncertainty, the resultant change on Λ, 0.13 GeV/c 2 , is taken as another systematic uncertainty. Finally, the pole mass is determined to be Λ = 2.56 ± 0.04(stat) ± 0.13(syst) GeV/c 2 . We compute the upper limits on the product branching fraction B(J/ψ → γ η) × B(γ → e + e − ) at the 90% confidence level (C.L.) as a function of m γ using a Bayesian method after incorporating the systematic uncertainty by smearing the likelihood curve with a Gaussian function with a width of the systematic uncertainty. The combined result is obtained by adding the logarithm likeli-hoods of two η decays by taking into account their correlated and uncorrelated systematic uncertainties. As shown in Fig. 8(a), the combined limits on product branching fraction B(J/ψ → γ η) × B(γ → e + e − ) vary in the range of (1.9 − 91.1) × 10 −8 for 0.01 ≤ m γ ≤ 2.4 GeV/c 2 depending on m γ points. The upper limit on B(J/ψ → γ η) at the 90% C.L. at each m γ point is computed by dividing the combined upper limit on the product branching fraction B(J/ψ → γ η) × B(γ → e + e − ) by the expected dark photon decay branching fraction of γ → e + e − obtained from Ref. [41]. We then compute the upper limits of the coupling strength between the dark sector and the SM at the 90% C.L. as a function of m γ using the Eq. (4.6) of Ref. [16], where the TFF is given by Eq. (3) with Λ = 2.56 GeV/c 2 . As shown in Fig. 8(b), the upper limits on at the 90% C.L. vary in the range of 10 −2 − 10 −3 for 0.01 ≤ m γ ≤ 2.4 GeV/c 2 depending on m γ .

VI. SUMMARY
In summary, with a data sample of (1310.6 ± 7.0) million J/ψ events collected with the BESIII detector, we study the EM Dalitz decay of J/ψ → e + e − η and search for a dark photon in J/ψ → γ η decay using two different η decay modes η → γγ and η → π + π − π 0 . The branching fraction of J/ψ → e + e − η is measured to be (1.42 ± 0.04(stat) ± 0.07(syst)) × 10 −5 , which supersedes the previous BESIII measurement [13]. We present the first measurement of TFF as a function of m e + e − for the decay J/ψ → e + e − η. The corresponding pole mass of the TFF is determined to be Λ = 2.56 ± 0.04(stat) ± 0.13(syst) GeV/c 2 by fitting the TFF versus m e + e − data with a modified TFF function. No evidence of dark photon γ production is observed, and we set upper limits on the product branching fraction B(J/ψ → γ η) × B(γ → e + e − ) at the 90% C.L. to be in the range of (1.9−91.1)×10 −8 for 0.01 ≤ m γ ≤ 2.4 GeV/c 2 depending on m γ . The upper limits on the coupling strength between the dark sector and the SM at the 90% C.L. are also set at the level of 10 −2 − 10 −3 , which are above the existing stringent experimental results [25][26][27].