Study of $e^+e^- \rightarrow \pi^{0}X(3872)\gamma$ and search for $Z_c(4020)^{0}\rightarrow X(3872)\gamma$

Using data samples collected with the BESIII detector operating at the BEPCII storage ring at center-of-mass energies from 4.178 to 4.600 GeV, we study the process $e^+e^-\rightarrow\pi^{0}X(3872)\gamma$ and search for $Z_c(4020)^{0}\rightarrow X(3872)\gamma$. We find no significant signal and set upper limits on $\sigma(e^+e^-\rightarrow\pi^{0}X(3872)\gamma)\cdot\mathcal{B}(X(3872)\rightarrow\pi^{+}\pi^{-}J/\psi)$ and $\sigma(e^+e^-\rightarrow\pi^{0}Z_c(4020)^{0})\cdot\mathcal{B}(Z_c(4020)^{0}\rightarrow X(3872)\gamma)\cdot\mathcal{B}(X(3872)\rightarrow\pi^{+}\pi^{-}J/\psi)$ for each energy point at $90\%$ confidence level, which is of the order of several tenths pb.


I. INTRODUCTION
The recent discovery of several charmonium-like states has attracted great experimental and theoretical interests [1]. The charmonium-like states are also called XY Z states where X is the isospin-singlet state with J P C = 1 −− , Y is the isospin-singlet state with J P C = 1 −− , and Z is the isospin-triplet state [2]. The masses of these states are above the open-charm thresholds, and due to the unexpected resonance parameters and decay channels, these states can not be described by conventional quark models. Therefore, they are good candidates for exotic states, such as hybrids, tetraquarks, molecules, etc. [3][4][5].
The X(3872) is a rather narrow state with a mass that is consistent with D 0D * 0 threshold. It decays through open-charm, radiative and isospin-violating pion emission decays, and is found to be an isospin singlet with J P C = 1 ++ [1]. Among these features, the extremely small mass difference between the X(3872) and D 0D * 0 threshold which we will denote as δ, is of particular interest. Taking the values for the D 0 , D * 0 and X(3872) masses from the PDG [1], δ is calculated to be (−10±180) keV/c 2 . Very recently, the LHCb reported a new measurement yielding δ = (70 ± 120) keV/c 2 [7,8]. However the improved precision is still insufficient to tell whether the X(3872) mass is above or below the D 0D * 0 threshold. Better knowledge of δ will be an important step towards a deeper understanding of the nature of the X(3872) [9,10], and eventually of other related XY Z states. A completely new method to measure the δ value by measuring the X(3872)γ line shape, which is sensitive to the δ value due to a triangle singularity, is proposed by Ref. [11][12][13]. Here, the X(3872)γ needs to be produced associated with another positive C-parity neutral meson, e.g. e + e − → π 0 X(3872)γ. In principle, this method could be applied at the BESIII experiment, based on the sizable data samples taken for XY Z studies. According to the theoretical prediction in Ref. [13], the cross section of the process e + e − → π 0 X(3872)γ is expected to be small. However, there could be other scenarios where the expected cross section is large.
Recently, the BESIII Collaboration reported an enhancement around 4.2 GeV for the e + e − → γX(3872) production cross sections [14], which suggests a connection between Y and X states. BESIII also reported another connection, now between Y and Z states, with the observation of a Y (4220) resonance in the process e + e − → π 0 Z c (3900) 0 [15]. Those observations may indicate some common nature among the XY Z states. Therefore, it is important to search for possible connections between Z and X states. Establishing connections among XY Z states may be a clue that can facilitate a better theoretical interpretation of these. One such connection [16] could be a transition Z c (4020) 0 → X(3872)γ in the scenario where the X(3872) is dominantly an S-wave D 0D * 0 molecule and the Z c (4020) 0 is an isotopic triplet of near-threshold S-wave D * D * resonances. Therefore, the search for the transition Z c (4020) 0 → X(3872)γ will help to quantitatively study the molecular picture of the X(3872). The Z c (4020) is observed in the e + e − → πZ c (4020) process, so the study of e + e − → π 0 X(3872)γ allows one to search for the transition Z c (4020) 0 → X(3872)γ.

II. BESIII DETECTOR AND MONTE CARLO SIMULATION
The BESIII detector is a magnetic spectrometer [17] located at the Beijing Electron Positron Collider (BEPCII) [18]. The cylindrical core of the BESIII detector 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 magnetic field. The solenoid is supported by an octagonal flux-return yoke with resistive plate chamber muon identifier modules interleaved with steel. The acceptance of charged particles and photons is 93% over the 4π solid angle. The chargedparticle momentum resolution at 1 GeV/c is 0.5%, and the dE/dx resolution is 6% for the 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 of the TOF barrel section is 68 ps, while that of the end cap section is 110 ps. The end cap TOF system was upgraded in 2015 with multigap resistive plate chamber technology, providing a time resolution of 60 ps [19]. About 70% of the data sample used here was taken after this upgrade.
Simulated data samples produced with the geant4based [20] Monte Carlo (MC) package, which includes the geometric description of the BESIII detector and the detector response, are used to determine the detection efficiency and to estimate the background contributions. The simulation includes the beam energy spread and initial-state radiation (ISR) in the e + e − annihilations modeled with the generator kkmc [21]. The ISR production of vector charmonium(-like) states and the continuum processes are incorporated also in kkmc [21]. The known decay modes are modeled with evtgen [22], using branching fractions summarized and averaged by the PDG [1], and the remaining unknown decays from the charmonium states are generated with lundcharm [23]. Final state radiation from charged final state particles is incorporated with the photos package [24].
Signal MC samples for e + e − → π 0 X(3872)γ and e + e − → π 0 Z c (4020) 0 → π 0 X(3872)γ are generated according to phase space at each center-of-mass energy point, assuming that the cross section follows the function fit for the e + e − → π + π − h c line shape in Ref. [25]. The event selection criteria and the detection efficiencies are determined and studied based on signal MC samples of 1 × 10 5 signal events generated for each value of √ s. Detection efficiencies are determined by the ratio of the reconstructed event yields (after the selection criteria) and the number of generated events. Inclusive MC samples consisting of open charm production processes are employed to investigate potential backgrounds.

III. EVENT SELECTION
For each charged track, the distance of closest approach to the 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 angles (θ) of the tracks with respect to the beam axis (ignoring the small crossing angle), must be within the fiducial volume of the MDC (| cos θ| < 0.93). Photons are reconstructed from isolated showers in the EMC, which are at least 10 • away from the nearest charged track. The photon energy is required to be at least 25 MeV in the barrel region (| cos θ| < 0.80) or 50 MeV in the end cap region (0.86 < | cos θ| < 0.92). To suppress electronic noise and energy depositions unrelated to the event, the EMC cluster timing from the reconstructed event start time is further required to satisfy 0 ≤ t ≤ 700 ns.
Since the reaction e + e − → π 0 X(3872)γ results in the final states γγγπ + π − ℓ + ℓ − , candidate events are required to have four tracks with zero net charge and at least three photons. Tracks with momenta larger than 1.0 GeV/c are assigned as leptons from the J/ψ decay; otherwise, they are regarded as pions. Leptons from the J/ψ decay with energy deposited in the EMC larger than 1.0 GeV are identified as electrons, while those with less than 0.4 GeV as muons. The π 0 candidates are reconstructed from photon pairs with invariant mass in the range 110 < M γγ < 150 MeV/c 2 .
To reduce the background contributions and to improve the mass resolution, a five-constraint (5C) kinematic fit is performed. Four constraints come from the total initial four momentum of the colliding beams; the last one is from constraining the M γγ invariant mass to the nominal π 0 value [1]. If there is more than one combination in an event, the one with the smallest χ 2 5C is chosen. Furthermore, the χ 2 5C is required to be less than 60. The J/ψ is reconstructed via ℓ + ℓ − decays, and the invariant mass of lepton pairs is required to be in the J/ψ mass window [3.080, 3.120] GeV/c 2 .

IV. BORN CROSS SECTION MEASUREMENT
IV.I e + e − → π 0 X(3872)γ After applying the above requirements, the remaining background is mainly coming from e + e − → γ ISR ηJ/ψ, η → π + π − π 0 and e + e − → γωJ/ψ, ω → π + π − π 0 events. In order to veto these events, the π + π − π 0 invariant mass is required to be outside the η and ω mass regions [0.535, 0.560] GeV/c 2 and [0.750, 0.810] GeV/c 2 , respectively. Besides the η and ω backgrounds, the π + π − ψ(3686), ψ(3686) → π 0 π 0 J/ψ background is removed by requiring the π + π − recoil mass to be outside the ψ(3686) mass region of [3.670, 3.700] GeV/c 2 . Figure 1 shows distributions of the π + π − J/ψ invariant mass M (π + π − J/ψ) for data and the MC samples of e + e − → π 0 X(3872)γ. The X(3872) signal region is taken as [3.860, 3.885] GeV/c 2 , while the sideband regions are set to be [3.825, 3.850] GeV/c 2 and [3.895, 3.920] GeV/c 2 . No significant X(3872) signals are seen at any energies. The signal yield is determined from the event yields in the X(3872) signal and sideband regions. The sideband yields are scaled by the ratio of the relevant mass-window widths in order to predict the background expected in the signal region. Upper limits on the number of signal events at the 90% C.L. are calculated by using a frequentist method [26] with unbounded profile likelihood treatment of systematic uncertainties, which is implemented by the package trolke [27] in the root framework [28], where the signal and background obey Poisson statistics, and the efficiencies are Gaussian-distributed. The numerical results are summarized in Table I. The Born cross section multiplied by the branching fraction σ(e + e − → π 0 X(3872)γ) · B(X(3872) → π + π − J/ψ) is calculated as: where N X(3872) is the number of X(3872) signal events, ǫ is the detection efficiency (excluding intermediate branching fractions), L int is the integrated luminosity [29], 1+δ(s) is the ISR correction factor obtained from a quantum electrodynamics calculations [21,30], 1 |1−Π| 2 is vacuum polarization factor [31]. The corresponding upper limits for this cross section at the 90% C.L. for each energy point are listed in Table I and shown in Fig. 2 (a).
Assuming the Born cross section of e + e − → π 0 X(3872)γ is constant at 4.178 ≤ √ s ≤ 4.600 GeV, the average Born cross section multiplied by the branching fractionσ(e + e − → π 0 X(3872)γ) · B(X(3872) → π + π − J/ψ) for data is calculated as: is the total number of X(3872) signal events, and i denotes each energy point. The corresponding upper limit for the average Born cross section multiplied by the branching fraction is determined to be 21.9 fb at the 90% C.L.
The possible connection between X and Z charmonium-like states can be studied via the decay Z c (4020) 0 → X(3872)γ. In order to search for the process, the X(3872) signal region is set to be [3.860, 3.885] GeV/c 2 , which is the same as for the e + e − → π 0 X(3872)γ study. As we do not observe any X(3872) signal, there cannot be any Z c (4020) 0 being produced in the given channel.
Assuming the Born cross section of e + e − → π 0 Z c (4020) 0 is constant at 4.178 ≤ √ s ≤ 4.600 GeV, the average Born cross section multiplied by branching fractionsσ(e + e − → π 0 Z c (4020) 0 ) · B(Z c (4020) 0 → X(3872)γ) · B(X(3872) → π + π − J/ψ) for data is calculated with the following formula: where N total Zc(4020) 0 is the total number of Z c (4020) 0 signal events. The corresponding upper limit for the average Born cross section multiplied by the branching fraction is determined to be 1.6 fb at the 90% C.L.
The integrated luminosity at each point has been measured with a precision of 1.0% using the Bhabha process [29].
The uncertainty from the tracking efficiency is 1.0% per track [32] and the uncertainty in photon detection efficiency is 1.0% per photon [33].
The uncertainty due to the kinematic fit requirements is estimated by correcting the helix parameters of charged tracks according to the method described in Ref. [34]. The difference between detection efficiencies obtained from MC samples with and without this correction is taken as the uncertainty.
The uncertainty for the J/ψ mass window is estimated using the control sample of e + e − → γ ISR ψ(3686), ψ(3686) → π + π − J/ψ. The difference of the efficiency between data and MC simulation is found to be 1.6% [35], which is taken as the uncertainty.
The uncertainty from the X(3872) mass window is estimated by changing the window range by ±10%, and the largest efficiency change is taken as the uncertainty.
The uncertainties arising from the Z c (4020) 0 mass and width are estimated by changing them by one standard deviation values [1] while generating the signal MC. The largest efficiency difference relative to the nominal one is taken as the uncertainty.
The line shape affects the ISR correction factor and the efficiency. No obvious signal was found for our Z c (4020) 0 search, so we use the line shape from e + e − → π + π − h c in Ref. [25] as the input line shape to get the nominal results. To get the uncertainty introduced by the line shape, we change it to a Breit-Wigner function describing the ψ(4230) or ψ(4415), with the masses and widths fixed to the values from PDG [1]. The largest difference of the final result is taken as a systematic uncertainty.
For the systematic uncertainty from the MC simulation describing the process e + e − → π 0 X(3872)γ, we use the three-body phase space MC simulation to get the nominal efficiency, then change to the e + e − → π 0 Z c (4020) 0 → π 0 X(3872)γ. The difference on the detection efficiency with and without the intermediate resonant state is taken as the uncertainty due to the MC generator model. The systematic uncertainty from the MC simulation describing the process e + e − → π 0 Z c (4020) 0 → π 0 X(3872)γ is estimated by varying the distribution of the π 0 polar angle θ. The nominal efficiency is determined assuming a flat distribution in cos θ. A conservative estimate of the systematic uncertainty is obtained using alternative MC samples with angular distributions of 1 ± cos 2 θ. The largest change of efficiency is taken as the uncertainty due to the MC generator model.
The ISR correction factor is obtained from quantum electrodynamics calculations [21,30]. We also analyze MC samples with and without ISR effects considered to get the ISR correction factor, the difference of the two results is taken as the systematic uncertainty on the ISR correction factor.