Search for the decay $J/\psi\to\gamma + \rm {invisible}$

We search for $J/\psi$ radiative decays into a weakly interacting neutral particle, namely an invisible particle, using the $J/\psi$ produced through the process $\psi(3686)\to\pi^+\pi^-J/\psi$ in a data sample of $(448.1\pm2.9)\times 10^6$ $\psi(3686)$ decays collected by the BESIII detector at BEPCII. No significant signal is observed. Using a modified frequentist method, upper limits on the branching fractions are set under different assumptions of invisible particle masses up to 1.2 $\mathrm{\ Ge\kern -0.1em V}/c^2$. The upper limit corresponding to an invisible particle with zero mass is 6.9$\times 10^{-7}$ at the 90\% confidence level.


I. INTRODUCTION
Understanding the nature of dark matter and finding direct evidence for its existence are among the primary goals of contemporary astronomy and particle physics [1,2]. Numerous experiments aim for the direct detection of dark matter, but no solid evidence has yet been found [3][4][5][6][7]. A series of Supersymmetric Standard Models [8], including the Next-to-Minimal Supersymmetric Model (NMSSM) [9,10], predict a light CP-odd pseudoscalar Higgs boson A 0 and a series of neutralinos. The light stable neutralino (χ 0 ), in particular, which is one possible explanation for the 511 keV γ ray feature observed by the INTEGRAL satellite [11], is one of the candidates for dark matter particles [12,13]. The χ 0 can couple with Standard Model particles via the A 0 boson, and the A 0 can be produced in the radiative decay of a quarkonium vector state, V [14,15]. The branching ratio of such a radiative decay is: where m A 0 , m V and m q are the masses of the A 0 , the quarkonium state, and the corresponding quark, respectively; α is the fine structure constant; G F is the Fermi coupling constant; C QCD is the combined QCD radiative and relativistic corrections [16], which depends on m A 0 ; and g q is the Yukawa coupling of the A 0 field to the quark-pair, and is g c = cosθ A / tan β for the charm quark and g b = cosθ A tan β for the bottom quark, where tan β is the usual ratio of vacuum expectation values and θ A is the Higgs mixing angle [13].
The CLEO-c [17], BaBar [18,19] and Belle [20] experiments have performed similar searches for J/ψ or Υ radiative decays into invisible particles, and no signal was observed. The upper limits at the 90% confidence level (C.L.) for the branching fraction of the decay J/ψ → γ + invisible, B(J/ψ → γ + invisible), are in the range (2.5 ∼ 6.3) × 10 −6 , depending on the mass of A 0 [17], where B(J/ψ → γ + invisible) is the product of B(J/ψ → γ + A 0 ) and B(A 0 → χ 0χ0 ). It is worth noting that the decay process J/ψ → γνν, which is allowed in the Standard Model, is an irreducible background in this analysis, but the predicted branching fraction is only 0.7 × 10 −10 , which is far below our experimental sensitivity [21]. Thus, this background is neglected.

II. BESIII DETECTOR AND MONTE CARLO SIMULATION
The BESIII detector is a magnetic spectrometer [22] located at the Beijing Electron Positron Collider (BEPCII) [23]. The cylindrical core of the BE-SIII 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.
The performance of the BESIII detector is evaluated using a geant4-based [24] Monte Carlo (MC) program that includes the description of the detector geometry and response. To check for potential backgrounds, an inclusive MC sample of ψ(3686) decays is used. The sample includes approximately the same number of ψ(3686) decays as in data. The production of the ψ(3686) resonance is simulated by the MC event generator kkmc [25], taking into account the beam energy spread; the known decay modes are generated using evtgen [26] with the branching fractions as given by the particle data group (PDG) [3]; the unknown decay modes are modeled with the lundcharm model [27]. Signal MC samples, corresponding to ψ(3686) → π + π − J/ψ with the subsequent decay J/ψ → γ + invisible, are used to evaluate the detection efficiencies and model the line shapes of variables of interest. The samples are generated under different assumptions for m A 0 . In these signal MC samples, the decay J/ψ → γ + invisible is modeled with an angular distribution of 1 + cos 2 θ γ (θ γ is the angle of the radiative photon relative to the positron beam direction in the J/ψ rest frame). Throughout the text, the decay ψ(3686) → π + π − J/ψ is modeled according to the formulas and measurement in Ref. [28]. In this analysis, detailed MC studies indicate that the dominant backgrounds are from ψ(3686) → π + π − J/ψ with subsequent decays J/ψ → γπ 0 , γη and γK L K L . These backgrounds are each generated exclusively with more than 100 times the statistics in data, where the decays of J/ψ → γπ 0 and γη are generated with the angular distribution of 1 + cos 2 θ γ , and J/ψ → γK L K L is modeled with the partial wave analysis (PWA) results of J/ψ → γK S K S [29] by assuming isospin symmetry. Many potential backgrounds of the form ψ(3686) → π + π − J/ψ with J/ψ decaying into purely neutral particles in the final states, or with large branching fractions, are generated exclu-sively with different generators, i.e. J/ψ → γη ′ , γη(1405) and γη c with the angular distribution of 1 + cos 2 θ γ ; J/ψ → γπ 0 π 0 and γπ + π − according to PWA results of J/ψ → γπ 0 π 0 [30] with isospin symmetry assumption; J/ψ → γK + K − and γK S K S (with K S → π 0 π 0 ) according to PWA results of J/ψ → γK S K S [29], as well as J/ψ → γπ 0 η, γγγ, K S K L , π 0 nn and ηnn with phase space distribution. The above MC samples with much larger statistics than in data are helpful to check potential backgrounds.

A. ANALYSIS METHOD
In this analysis, the J/ψ sample originates from the decay ψ(3686) → π + π − J/ψ. The analysis strategy is to first tag J/ψ events by selecting two oppositely charged pions, and then to search for the decay J/ψ → γ + invisible within the tagged J/ψ sample. The branching fraction of the decay J/ψ → γ+ invisible is calculated using: where N sig and N J/ψ are the yields of the signal candidates of J/ψ → γ + invisible and ψ(3686) → π + π − J/ψ, respectively, and ǫ sig and ǫ J/ψ are the corresponding detection efficiencies, evaluated with the corresponding MC samples. A semi-blind analysis is performed to avoid possible bias, where only one quarter of the full data sample is used to optimize the event selection criteria and to decide upon the upper limit calculation approach. The final results are obtained with the full data sample by repeating the analysis only after all the analysis methods are frozen. In this paper, only the results based on the full data sample are presented.
B. J/ψ TAG PROCEDURE J/ψ events are tagged using the two oppositely charged pions produced in the process ψ(3686) → π + π − J/ψ. For each charged pion candidate, the point of closest approach to the e + e − interaction point must be within ±10 cm in the beam direction and 1 cm in the plane perpendicular to the beam, and the polar angle θ with respect to the axis of the drift chamber must satisfy the condition |cosθ| < 0.93. The charged pions are identified by combining the information of the flight time measured from TOF and the dE/dx measured in MDC. The corresponding likelihood for the pion hypothesis is required to be larger than that of the kaon hypothesis and 0.001. To suppress pions not from the decay ψ(3686) → π + π − J/ψ, the momentum of a pion is required to be less than 0.45 GeV/c. Additionally, to further suppress the background from γ conversion occurring in the inner detector, the angle between the two selected pions (θ 1 ) is required to satisfy cosθ 1 <0.95. To veto γγ fusion events, the polar angle (θ 2 ) of the total momentum vector of the pion pair should fullfill |cosθ 2 | <0.95.
To identify ψ(3686) → π + π − J/ψ candidate events, the recoiling mass of the π + π − system, M rec where E CMS is the center-of-mass energy of the initial e + e − system, and E π + π − and p π + π − are the sum of the energies and momenta of the pions in the rest frame of the initial e + e − system, respectively. The distribution of M rec π + π − in the range [3.06, 3.14] GeV/c 2 is shown in Fig. 1, where multiple entries per event are allowed. A clear J/ψ peak with low level of background events is observed. To extract the signal yield, a binned maximum likelihood fit to the M rec π + π − distribution is performed. To better model the J/ψ signal shape, a control sample of ψ(3686) → π + π − J/ψ with the subsequent decay J/ψ → e + e − , which has almost no background, is selected. In the fit, the signal shape is modeled using the M rec π + π − distribution of the control sample convoluted with a Gaussian function, which represents the resolution difference between J/ψ → e + e − and the J/ψ inclusive decay. The background is described by a 2 nd order Chebychev polynomial function. The fit results are shown in Fig. 1, and the resolution difference of the M rec π + π − distribution between J/ψ → e + e − and the inclusive decay is found to be small, i.e., the width of the Gaussian function is close to zero. Candidate events in the J/ψ signal region [3.082, 3.112] GeV/c 2 , which is roughly three times the M rec π + π − resolution, are used for further analysis. The number of tagged J/ψ events in the signal region is (8848 ± 1) × 10 4 , obtained by integrating the fitted signal curve in the J/ψ signal region.

C. SIGNAL SEARCH PROCEDURE
We search for the decay J/ψ → γ + invisible in the remaining J/ψ candidate events by requiring no additional charged track is present and there is exactly one photon candidate. Photon candidates are reconstructed from EMC and must satisfy the following requirements. The minimum energy is 25 MeV for barrel showers (| cos θ| <0.80) or 50 MeV for end-cap showers (0.86< | cos θ| <0.92). To eliminate showers associated with charged particles, the photon candidates must be separated by at least 20 degrees from any charged tracks in EMC. To suppress electronic noise or the showers unrelated to the events, the time of the cluster measured from EMC is required to be within 0 and 700 ns after the event start time. To further suppress background with multiple photons in the final state, the total energy of the remaining showers in the EMC, not satisfying the requirements on photon candidates, is required to be less than 0.1 GeV. In order to improve the resolution, to further suppress background, and to make sure the invisible particle is within the detector volume, the directions of the signal photon and the missing particle (calculated as the recoiling momentum against the system of π + π − pair and signal photon) are required to be within the EMC barrel region.
After the above selection criteria, detailed MC studies indicate that the dominant backgrounds are from ψ(3686) → π + π − J/ψ with J/ψ decays into final states including neutral hadrons, e.g., nn, γK L K L , π 0 nn. To further suppress these backgrounds, a series of requirements on the shower shape variables, i.e., the second moment should be larger than 5 cm 2 and less than 25 cm 2 , the lateral moment should be larger than 0.1 and less than 0.4, the ratio of energy in 3 × 3 and 5 × 5 crystals should be larger than 0.95 due to the narrow shower shape for γ, as well as the number of crystals(N crystals ) and energy(E shower ) of the shower should satisfy 4< N crystals − 10 × E shower (GeV)<20 due to the strong relation between these two variables for γ, are implemented, where these selection criteria are optimized with the control samples of γ,n/n and K L selected from the decay processes J/ψ → π + π − π 0 (π 0 → γγ), J/ψ → pπ −n + c.c. and J/ψ → KπK L , J/ψ → π + π − φ(φ → K S K L ), respectively.
The variable E * γ , which is defined as the energy of the selected photon in the J/ψ rest frame, is used to identify the signal. For the signal process J/ψ → γ + invisible with a given mass and zero width for the invisible particle, the E * γ is expected to be convoluted with the corresponding detector resolution function. The distribution of E * γ above 1.25 GeV for the selected events is shown in Fig. 2. The dominant backgrounds are from ψ(3686) → π + π − J/ψ with subsequent decays J/ψ → γK L K L , γη and γπ 0 , where the latter two produce the peak in the E * γ distribution. The above three backgrounds, depicted in Fig. 2, are estimated with the corresponding exclu-sive MC samples and normalized according to the PDG branching fractions [3]. The contribution from the non-J/ψ process is found to be small and is estimated by the normalized data sample in the J/ψ sideband region (on the M rec π + π − distribution), also shown in Fig. 2. To better model the peaking backgrounds from J/ψ → γη and J/ψ → γπ 0 in the follow up procedure, a binned maximum likelihood fit is performed on the two corresponding exclusive MC samples, individually. In the fit, the peaking component, where the detected photon is from the J/ψ radiative decay, is described by a Crystal Ball function [31], while the others, which distribute relatively uniformly and correspond to the case that the detected photon is not from the J/ψ radiative decay, is described by a second order Chebychev polynomial function. The Crystal Ball functions obtained are used to represent the peaking background from J/ψ → γη/π 0 in the following analysis. The number of events are normalized according to the PDG [3] and the yield of tagged J/ψ in data.
Unbinned likelihood fits are performed on the E * γ range from 1.25 to 1.65 GeV/c 2 , where the signal is described by the shape from signal MC and the background is described by an exponential function and two fixed crystal ball functions, with different m invisible (mass of invisible particle) assumptions up to 1.2 GeV/c 2 . (For heavier invisible particle assumption, the signal shape is broken.) As no strong peaks are observed in all fits, the upper limits are calculated by using the modified frequentist method known as CL s [32,33] combined with the asymptotic approximation [34]. In this approach, the test statistic is the profile likelihood ratio, where the likelihood is given with the Possion function: where s i represents the signal probability in the i-th bin, B is the branching fraction B(J/ψ → γ + invisible), b exp ij is the expected background number in the i-th bin for the j-th source. Here the background is modeled with the exponential function and the two fixed crystal functions from zero-signal assumption fit result. Additionally, systematic uncertainties are included assumimg Gaussian distributions for nuisance parameters. The upper limit is determined by integrating the test statistic in the range of positive assumed branching fractions.

D. SYSTEMATIC UNCERTAINTIES
Three categories of systematic uncertainties, which are associated with the number of tagged J/ψ events(N J/ψ ), the signal efficiency and the estimated numbers of peaking backgrounds, are considered individually.
The systematic uncertainty related to N J/ψ comes from the binned fit procedure and includes the fit range, bin size, and the shapes of the signal and background. The uncertainties from the fit range and bin size are estimated to be 0.6% by varying the fit range by ± 5 MeV and 0.3% by changing the bin size from 0.4 to 0.2 MeV, respectively. The uncertainties from the signal and background shapes are determined as 0.1%, individually, estimated by the alternative fits without convoluting the Gaussian function on the signal shape or using a 3-rd order Chebychev function for background. The total uncertainty related to N J/ψ is 0.7%, obtained by adding the above components in quadrature.
To estimate the uncertainty related to the signal efficiency, two control samples, e + e − → γe + e − and J/ψ → π + π − π 0 (π 0 → γγ), are selected. The former is used to estimate the uncertainty associated with the event topology requirement, i.e., no extra photons or charged tracks, as well as the remaining energy requirement. And the latter is used to estimate the uncertainty associated with the shower shape requirements. The resulting differences on the efficiency between the data and MC simulation are assigned to be the systematic uncertainty, individually. The numerical results are 0.6% and 0.9% for the "no extra photons or charged tracks" requirement and the shower shape requirements, respectively. The uncertainty due to the energy cut on the remaining showers in the EMC is less than 0.1% and negligible. For the photon reconstruction efficiency, the uncertainty is 1% [35]. By adding all the above uncertainties in quadrature, the systematic uncertainty from the signal efficiency is 1.5%.
The uncertainties due to the estimated numbers of two peaking backgrounds come from the J/ψ yield, the decay branching fractions, and the selection efficiency (or fake rate) for the process J/ψ → γη/π 0 . The uncertainty of J/ψ yield is discussed above, 0.7%. The uncertainties of decay branching fractions are quoted from the PDG [3], 3.0% for J/ψ → γη and 4.8% for J/ψ → γπ 0 . The uncertainties associated with the selection efficiency include those of γ selection (including photon reconstruction and shower shape requirements) and the event topology requirement (including charged tracks number, photon number and extra showers' energy requirements). The uncertainty associated with the γ selection is discussed above. The uncertainty associated with the event topology requirement is investigated by studying a control sample of ψ(3686) → π + π − J/ψ, J/ψ → φη. For the decay of J/ψ → γη, the control sample is selected by tagging a π + π − pair and a K + K − pair as well as the J/ψ and φ mass window requirements on the π + π − recoiling system and K + K − system, respectively. The corresponding efficiency is computed for both data and MC samples by fitting to the η signal on the recoiling mass of π + π − K + K − system before and after implementing the event topology requirements. The resulting difference in the efficiencies is taken as the systematic uncertainty. For the decay J/ψ → γπ 0 , no extra charged tracks is required, since the π 0 decays into the γγ final state dominantly. Then the same procedure is applied. Since the efficiency of the event topology requirement is extremely low, ∼ 0.2%/0.3% for the peaking backgrounds of J/ψ → γη/π 0 , the resulting uncertainties, 16% for both J/ψ → γη/π 0 , are dominated by the statistical uncertainty of the data control sample, and are conservatively taken as the systematic uncertainties in this analysis. By adding all uncertainties in quadrature, the systematic uncertainties for the number of peaking backgrounds are 17% for both ψ(3686) → π + π − J/ψ and J/ψ → γη/π 0 .

E. UPPER LIMIT RESULT
Taking into account all systematic uncertainties and the signal detection efficiencies obtained from MC simulation with different m invisible assumptions, the expected upper limits on the branching fraction of J/ψ → γ + invisible at the 90% C.L. are calculated with the CL s approach and are shown in Fig. 3. The expected upper limits as well as their uncertainties are also obtained using toy MC sample, which is generated using the background model from no signal assumption fit with the same luminosity as data set. The result from data is consistent with the zero-signal assumption in the 2σ region. And for the zero mass assumption of the invisible particle the upper limit is 6.9×10 −7 . The signal significances with different mass assumptions are also shown in Fig. 3, where the signal significance is calculated by 2ln( Lsig L0 ) incorporating the maximum likelihood with floating signal yield L sig and with zero-signal yield L 0 .

IV. SUMMARY AND DISCUSSION
In summary, we search for the J/ψ radiative decay into a weakly interacting neutral particle in the process ψ(3686) → π + π − J/ψ by using a ψ(3686) sample of (448.1 ± 2.9) × 10 6 events collected with the BESIII detector. No significant signal is observed, and the upper limits at the 90% C.L. on the decay branching fraction of J/ψ → γ + invisible are obtained for different m invisible assumptions up to 1.2 GeV/c 2 . The observed upper limit for a zero mass of the invisible particle is improved by a factor 6.2 compared to the previous results [17].
To further investigate the physical parameters in NMSSM, and to better compare the physical results from the different quarkonium decays, according to Ref. [15] and Eq. (1), the upper limits of g c × tan 2 β × B(A 0 → χ 0χ0 ) based on the measured upper limits of the J/ψ → γ + invisible decay branching fractions are extracted for tan β = 0.7, 0.8 and 0.9, individually, as presented in Fig. 4 (a). The extracted results are directly compared to g b × B(A 0 → χ 0χ0 )(= g c × tan 2 β × B(A 0 → χ 0χ0 )), which is obtained based on the Belle results [20] and also presented in Fig. 4 (a).
We obtain better sensitivity in the range tan β ≤ 0.6 compared to the Belle results. Combining the results from Belle [20], we also extract upper limits on cos θ A (= √ g b g c ) × B(A 0 → χ 0χ0 ), as presented in Fig. 4 (b). V. ACKNOWLEDGMENTS