Study of electromagnetic Dalitz decays $\chi_{cJ} \rightarrow \mu^{+}\mu^{-}J/\psi$

Using $4.48 \times 10^{8}$ $\psi(3686)$ events collected with the BESIII detector, we search for the decays $\chi_{cJ} \rightarrow \mu^{+}\mu^{-}J/\psi$ through the radiative decays $\psi(3686) \rightarrow \gamma\chi_{cJ}$, where $J=0,1,2$. The decays $\chi_{c1,2} \rightarrow \mu^{+}\mu^{-}J/\psi$ are observed, and the corresponding branching fractions are measured to be $\mathcal{B}(\chi_{c1} \rightarrow \mu^{+}\mu^{-}J/\psi) = (2.51 \pm 0.18 \pm 0.20)\times10^{-4}$ and $\mathcal{B}(\chi_{c2} \rightarrow \mu^{+}\mu^{-}J/\psi) = (2.33 \pm 0.18 \pm 0.29)\times10^{-4}$, where the first uncertainty is statistical and the second one systematic. No significant $\chi_{c0} \rightarrow \mu^{+}\mu^{-}J/\psi$ decay is observed, and the upper limit on the branching fraction is determined to be $2.0\times10^{-5}$ at 90% confidence level. Also, we present a study of di-muon invariant mass dependent transition form factor for the decays $\chi_{c1,2} \rightarrow \mu^{+}\mu^{-}J/\psi$.


I. INTRODUCTION
The electromagnetic (EM) Dalitz decays M 1 → M 2 ℓ + ℓ − (M for meson, ℓ = e or µ) provide information on the internal structure of the mesons and the interactions of the mesons with the electromagnetic field [1][2][3][4]. Such decays are well studied in light-quark meson sector [5], but very rare in charm sector, let alone in bottom sector. The q-dependent transition form factor (TFF), where q is the invariant mass of the lepton pair, serves as a sensitive probe to the inner structure of the mesons involved, thus provides crucial tests to the theoretical models developed to describe the nature of the mesons, especially the charmonium-like states which manifested exotic properties compared with conventional charmonium states. One example is the X(3872); while it is a candidate for the radial excitation of the P -wave charmonium state χ c1 , it is also a good candidate of DD * molecule. Precision measurement of its EM Dalitz decays and comparison with those of χ c1 decays and the relevant theoretical models may eventually reveal its nature.
The branching fractions for χ cJ → µ + µ − J/ψ are predicted in Ref. [6] (Throughout this paper, χ cJ refers to χ c0 , χ c1 , and χ c2 ), and it is demonstrated that the µ + µ − decay channels are more suitable for the investigation of χ cJ → γ * J/ψ decay vertices than e + e − decay channels, which have been observed at BESIII [7]. Very recently, LHCb reported the observation of χ c1,2 → µ + µ − J/ψ [8] and measured the χ c1,2 resonance parameters. At BESIII, since the branching fractions of ψ(3686) → γχ cJ can be calculated very precisely, we can measure the absolute branching fractions of χ cJ → µ + µ − J/ψ. The branching fractions in theoretical calculations are related to the TFF, so the measurements can provide more constraints on theoretical calculations about TFF correction [6].

II. BESIII DETECTOR AND MONTE CARLO SIMULATION
The BESIII detector is a magnetic spectrometer [10] located at the BEPCII [11]. 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 counter muon identifier modules interleaved with steel. The acceptance of charged particles and photons is 93% over 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 part is 68 ps, while that of the end cap part is 110 ps.
Simulated samples produced with the geant4based [13] 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 backgrounds. The simulation includes the beam energy spread and initial state radiation (ISR) in the e + e − annihilations modelled with the generator kkmc [14]. The signal MC samples are generated using evtgen [15] with a q-dependent decay amplitude based on the assumption of a point-like meson, as described in Refs. [6,16]. The inclusive MC sample consists of the production of the ψ(3686) resonance, the ISR production of the J/ψ, and the continuum processes incorporated in kkmc [14]. The known decay modes are modelled with evtgen [15] using branching fractions taken from the Particle Data Group [5], and the remaining unknown decays from the charmonium states with lundcharm [17]. The final state radiations (FSR) from charged final state particles are incorporated with the photos package [18].

III. EVENT SELECTION
Candidate events are required to have four charged tracks, with zero net charge, and at least one photon. For each charged track, the distance of the closest approach to the interaction point (IP) is required to be smaller than 1 cm on the radial direction and smaller than 10 cm along the beam axis. The polar angle (θ) of the tracks 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 20 • away from the nearest charged track. The photon energy is required to be at least 25 MeV in the barrel region (| cos θ| < 0.8) or 50 MeV in the endcap region (0.86 < | cos θ| < 0.92). In order to suppress electronic noise and energy depositions which are unrelated to the event, the time after the collision at which the photon is recorded in the EMC is required to satisfy 0 ≤ t ≤ 700 ns.
According to the study of signal MC, the tracks with momentum larger than 1 GeV/c are assumed to be leptons from J/ψ decay, otherwise they are considered as muons from χ cJ decay. The EMC deposited energy is used to separate electrons and muons from J/ψ, leptons from the J/ψ decay with energy deposited in EMC larger than 1.0 GeV are identified as electrons, less than 0.3 GeV as muons. The J/ψ signal is selected by requiring the invariant mass of the lepton pair to be in the mass region [3.085, 3.110] GeV/c 2 . A vertex fit is performed on the four charged tracks to restrict the tracks originated from the IP. In order to reduce backgrounds and improve the mass resolution, a four-constraint (4C) kinematic fit is performed by constraining the total four momentum to that of the initial beams. All the photons are looped with the four tracks in the kinematic fit and only those with a χ 2 < 40 are retained. If there is more than one photon candidate in an event, only the one with the least χ 2 is retained for further analysis.
A study of the ψ(3686) inclusive MC sample shows that, after applying the above selection criteria, the main backgrounds come from the four processes: Cat. I: To suppress the backgrounds from Cats. I and II, where one photon is converted into two electrons, a photon-conversion finder [19] is used to reconstruct the photon-conversion vertex. There are no additional requirements in photon-conversion finder. The distance from the reconstructed conversion vertex to the z axis, R xy , is used to distinguish the photon conversion background from the signal. Figure 1 shows the R xy distribution of the decay χ c1 → µ + µ − J/ψ as an example. By studying the MC samples of Cats. I and II, the peaks around R xy = 3 and 6 cm match the positions of the beam pipe and the inner wall of MDC [10], respectively. For the background from Cat. III, it enhances around R xy = 0 cm. In order to remove all these backgrounds, a requirement R xy > 8.5 cm is applied. For the signal events, the reconstructed R xy is almost proportional to the opening angle of the two µ tracks, so if the angle of the two tracks is large, the variable R xy is also large.
where Rxy is the distance from the reconstructed conversion vertex to the z axis calculated from photo-conversion finder [19]. The points with error bars are data, the red histograms are for the signal MC simulations and the blue dotted lines are for the background MC simulations.
To remove the background from Cat. IV, which has the same final state as the signal event, a requirement M (γµ + µ − ) < 0.535 or > 0.560 GeV/c 2 is applied, where M (γµ + µ − ) is the invariant mass of γµ + µ − . The requirement removes almost all the Cat. IV background events (it accounts for about 60% of the remaining background, is about 140 events), with an efficiency loss of about 15% for signal events.
After applying all the above criteria for ψ(3686) inclusive MC sample, which does not include the signal processes, only a few events are left, and the overall contribution from ψ(3686) decays in the M (µ + µ − J/ψ) distribution is found to be smooth. Here M (µ + µ − J/ψ) = M (µ + µ − ℓ + ℓ − ) − M (ℓ + ℓ − ) + m(J/ψ) is used to reduce the resolution effect of the lepton pairs, and m(J/ψ) is the nominal mass of J/ψ [5]. The continuum background is studied by using the data collected at √ s = 3.65 GeV, and the contribution is found to be negligible. Figure 2 shows the M (µ + µ − J/ψ) distribution for selected events from data. Clear enhancements at the masses of χ c1,2 are seen, corresponding to the decays χ c1,2 → µ + µ − J/ψ, while no significant signals for the χ c0 → µ + µ − J/ψ decay are found. An unbinned maximum likelihood fit is performed to the M (µ + µ − J/ψ) distribution to extract the signal yields. We use the MCdetermined shapes to describe the χ cJ signals, where the magnitudes are free parameters. The background is described by a linear function with the number of events as free parameter. The fit result is shown in Fig. 2 and the corresponding signal yields are summarized in Table I. The significances for χ c1,2 are larger than 10σ by comparing the likelihood values for the fits with or without χ c1,2 signals and taking the change of the number of degreesof-freedom into account. Since no significant signal is observed for χ c0 → µ + µ − J/ψ decay, we give the upper limit at 90% confidence level (C.L.) using Bayesian method. With the fit function described before, we scan the number of χ c0 signal yield to obtain the likelihood distribution, and smear it with the systematic uncertainty. The upper limit of the number of χ c0 signal yield N up χc0 at 90% C.L. is obtained via  Distribution of M (µ + µ − J/ψ) in data (dots with error bars). The solid curve is the overall fit result, the dashed curve is for background contribution.

IV. BRANCHING FRACTION MEASUREMENT
The branching fractions B(χ cJ → µ + µ − J/ψ) are cal-culated according to where N is the signal yields obtained from the fit, N ψ(3686) is the number of ψ(3686) events [9], ǫ is the average selection efficiency of the decays J/ψ → e + e − and J/ψ → µ + µ − determined from the signal MC samples, B rad is the branching fraction of the radiative transitions ψ(3686) → γχ cJ , and B J/ψ→ℓ + ℓ − is the sum of branching fractions of J/ψ → e + e − and J/ψ → µ + µ − . All the branching fractions used are taken from Ref. [5]. The results of χ cJ → µ + µ − J/ψ are listed in Table I. Figure 3 shows comparisons of the observed q distributions without efficiency correction in data and MC simulation for the decays χ c1,2 → µ + µ − J/ψ, where the χ c1 and χ c2 signals are extracted by requiring a mass within [3.500, 3.520] and [3.545, 3.565] GeV/c 2 , respectively; with these criteria the backgrounds are expected to be about 5%. The data are in reasonable agreement with the MC simulation generated by using the point-like model described in Refs. [6,16]. To measure the TFF, the q distributions in the decays χ c1,2 → µ + µ − J/ψ are divided into 4 and 5 regions, respectively. The bin-by-bin signal yields and corresponding branching fractions are listed in Table II. The quantum electrodynamics (QED) predicted branching fraction results of χ c1,2 → µ + µ − J/ψ are obtained from Eq.(2) in Ref. [6], and the uncertainty is from the branching fractions of χ c1,2 → γJ/ψ. The TFFs are the ratios of measured branching fractions and QED predicted branching fractions in each bin, which are also listed in Table II. Figure 4 shows the TFF distributions for the decays χ c1,2 → µ + µ − J/ψ. If we use the parametrization F (q) = 1 1−q 2 /Λ 2 [3] to fit TFF distributions, the fit  and TFF |F (q)| 2 for the decays χc1,2 → µ + µ − J/ψ in each bin. Here the first uncertainty is statistical and the second systematic.

VI. SYSTEMATIC UNCERTAINTY
The systematic uncertainties for the branching fraction measurement arise from the following sources: track reconstruction, photon detection, kinematic fit, J/ψ mass criteria, M (γµ + µ − ) requirement, R xy requirement, fit procedure, angular distribution, number of ψ(3686) events, and the branching fractions of the cascade decays. All uncertainties are discussed in detail below.
The difference between data and MC simulation on the tracking efficiency of high momentum tracks is estimated to be 1% [20] using control sample ψ(3686) → π + π − J/ψ, J/ψ → ℓ + ℓ − . To study the difference on the low momentum muon tracking efficiency between data and MC simulation, we select a sample of ψ(3686) → π + π − J/ψ, J/ψ → µ + µ − γ. The weighted difference between data and MC simulation is about 4% for the low momentum µ + µ − pair. We also checked cosθ dependence of low momentum tracking efficiency using control sample J/ψ → ppπ + π − . The π tracking efficiency is cosθ dependent, and we use these results to correct the efficiency for µ + µ − pair, while the weighted difference between data and MC simulation is also about 4%. Totally, a 6% systematic uncertainty on tracking efficiency is attributed to all channels. The uncertainty on the photon detection efficiency is derived from a control sample of J/ψ → ρ 0 π 0 and is 1.0% per photon [21].
In the 4C kinematic fit, the helix parameters of charged tracks are corrected to reduce the discrepancy between data and MC simulation as described in Ref. [22]. The correction factors are obtained by studying a control sample of ψ(3686) → π + π − J/ψ, J/ψ → ℓ + ℓ − . To determine the systematic uncertainty from this source, we determine the efficiencies from the MC samples without the helix correction; the resulting differences with respect to the nominal values are taken as systematic uncertainties.
The uncertainty associated with the J/ψ mass requirement is 1.0%, which is determined by studying a control sample of ψ(3686) → ηJ/ψ, η → γγ (where one γ undergoes conversion to an e + e − pair) or η → γe + e − decays. The systematic uncertainty related to the M (γµ + µ − ) requirement is studied by removing the requirement and then repeat the analysis to get the result. The difference from the nominal result is taken as systematic uncertainty from this source. Likewise to estimate the systematic uncertainty from R xy requirement, we also remove the requirement to get the result and the difference is taken as systematic uncertainty. Due to the absence of χ c0 signal, the uncertainties for χ c0 channel on M (γµ + µ − ) requirement and R xy requirement are taken from the larger one in the χ c1 and χ c2 channels.
The sources of uncertainty in the fit procedure include the fit range, the signal shape, and the background shape. The uncertainty related to the fit range is obtained by varying the limits of the fit range by ±5 MeV/c 2 . The largest difference in the signal yields with respect to the nominal values is taken as systematic uncertainty. In the nominal fit, the signal shapes are described with the MC simulated signal shapes. An alternative fit is performed with the signal MC simulated shapes convolved with a Gaussian function. The resulting change in the signal yields is taken as systematic uncertainty. The uncertainty associated with the background shape is estimated by an alternative fit replacing the first order polynomial function with a second order polynomial function. The change in the signal yields is taken as systematic uncertainty. About the uncertainty from fit procedure for χ c0 channel, we try to use different combinations of fit range, signal shape and background shape to get the upper limits, and choose the largest one as nominal upper limit.
The number of ψ(3686) events is measured with an uncertainty of 0.7% by using the inclusive hadronic events [9]. The uncertainties of the branching fractions in the cascade decays are taken from Ref. [5]. Table III summarizes all individual systematic uncertainties, and the overall uncertainties are the quadrature sums of the individual ones, assuming they are independent. The overall uncertainties are also taken as the systematic uncertainties of the branching fraction measurements for the decays χ c1,2 → µ + µ − J/ψ in each bin in Table II.

ACKNOWLEDGMENTS
The BESIII collaboration thanks the staff of BEPCII and the IHEP computing center for their strong support.
This work is supported in part by National Key