Search for the decay 
J/ψ→γ+
 invisible

M. Ablikim, M. N. Achasov, P. Adlarson, S. Ahmed, M. Albrecht, A. Amoroso, Q. An, Anita, Y. Bai, O. Bakina, R. Baldini Ferroli, I. Balossino, Y. Ban, K. Begzsuren, J. V. Bennett, N. Berger, M. Bertani, D. Bettoni, F. Bianchi, J. Biernat, J. Bloms, A. Bortone, I. Boyko, R. A. Briere, H. Cai, X. Cai, A. Calcaterra, G. F. Cao, N. Cao, S. A. Cetin, J. F. Chang, W. L. Chang, G. Chelkov, D. Y. Chen, G. Chen, H. S. Chen, M. L. Chen, S. J. Chen, X. R. Chen, Y. B. Chen, Z. J. Chen, W. S. Cheng, G. Cibinetto, F. Cossio, X. F. Cui, H. L. Dai, J. P. Dai, X. C. Dai, A. Dbeyssi, R. B. de Boer, 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, S. X. Du, J. Fang, S. S. Fang, Y. Fang, R. Farinelli, L. Fava, F. Feldbauer, G. Felici, C. Q. Feng, M. Fritsch, C. D. Fu, Y. Fu, X. L. Gao, Y. Gao, Y. Gao, Y. G. Gao, I. Garzia, E. M. Gersabeck, A. Gilman, K. Goetzen, L. Gong, W. X. Gong, W. Gradl, M. Greco, L. M. Gu, M. H. Gu, S. Gu, Y. T. Gu, C. Y. Guan, A. Q. Guo, L. B. Guo, R. P. Guo, Y. P. Guo, Y. P. Guo, A. Guskov, S. Han, T. T. Han, T. Z. Han, X. Q. Hao, F. A. Harris, K. L. He, F. H. Heinsius, T. Held, Y. K. Heng, M. Himmelreich, T. Holtmann, Y. R. Hou, Z. L. Hou, H. M. Hu, J. F. Hu, T. Hu, Y. Hu, G. S. Huang, L. Q. Huang, X. T. Huang, Y. P. Huang, Z. Huang, N. Huesken, T. Hussain, W. Ikegami Andersson, W. Imoehl, M. Irshad, S. Jaeger, S. Janchiv, Q. Ji, Q. P. Ji, X. B. Ji, X. L. Ji, H. B. Jiang, X. S. Jiang, X. Y. Jiang, J. B. Jiao, Z. Jiao, S. Jin, Y. Jin, T. Johansson, N. Kalantar-Nayestanaki, X. S. Kang, R. Kappert, M. Kavatsyuk, B. C. Ke, I. K. Keshk, A. Khoukaz, P. Kiese, R. Kiuchi, R. Kliemt, L. Koch, O. B. Kolcu, B. Kopf, M. Kuemmel, M. Kuessner, A. Kupsc, M. G. Kurth, W. Kühn, J. J. Lane, J. S. Lange, P. Larin, L. Lavezzi, H. Leithoff, M. Lellmann, T. Lenz, C. Li, C. H. Li, Cheng Li, D. M. Li, F. Li, G. Li, H. B. Li, H. J. Li, J. L. Li, J. Q. Li, Ke Li, L. K. Li, Lei Li, P. L. Li, P. R. Li, S. Y. Li, W. D. Li, W. G. Li, X. H. Li, X. L. Li, Z. B. Li, Z. Y. Li, H. Liang, H. Liang, Y. F. Liang, Y. T. Liang, L. Z. Liao, J. Libby, C. 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. M. Liu, Huanhuan Liu, Huihui Liu, J. B. Liu, J. Y. Liu, K. Liu, K. Y. Liu, Ke Liu, L. Liu, Q. Liu, S. B. Liu, Shuai Liu, T. Liu, X. Liu, Y. B. Liu, Z. A. Liu, Z. Q. Liu, Y. F. Long, X. C. Lou, F. X. Lu, H. J. Lu, J. D. Lu, J. G. Lu, X. L. Lu, Y. Lu, Y. P. Lu, C. L. Luo, M. X. Luo, P. W. Luo, T. Luo, X. L. Luo, S. Lusso, X. R. Lyu, F. C. Ma, H. L. Ma, L. L. Ma, M. M. Ma, Q. M. Ma, R. Q. Ma, R. T. Ma, X. N. Ma, X. X. Ma, X. Y. Ma, Y. M. Ma, F. E. Maas, M. Maggiora, S. Maldaner, S. Malde, Q. A. Malik, A. Mangoni, Y. J. Mao, Z. P. Mao, S. Marcello, Z. X. Meng, J. G. Messchendorp, G. Mezzadri, T. J. Min, R. E. Mitchell, X. H. Mo, Y. J. Mo, N. Yu. Muchnoi, H. Muramatsu, S. Nakhoul, Y. Nefedov, F. Nerling, I. B. Nikolaev, Z. Ning, S. Nisar, S. L. Olsen, Q. Ouyang, S. Pacetti, X. Pan, Y. Pan, A. Pathak, P. Patteri, M. Pelizaeus, H. P. Peng, K. Peters, J. Pettersson, J. L. Ping, R. G. Ping, A. Pitka, R. Poling, V. Prasad, H. Qi, H. R. Qi, M. Qi, T. Y. Qi, S. Qian, W.-B. Qian, Z. Qian, C. F. Qiao, L. Q. Qin, X. S. Qin, Z. H. Qin, J. F. Qiu, S. Q. Qu, K. H. Rashid, K. Ravindran, C. F. Redmer, A. Rivetti, V. Rodin, M. Rolo, G. Rong, Ch. Rosner, M. Rump, A. Sarantsev, Y. Schelhaas, C. Schnier, K. Schoenning, D. C. Shan, W. Shan, X. Y. Shan, M. Shao, C. P. Shen, P. X. Shen, X. Y. Shen, H. C. Shi, R. S. Shi, X. Shi, X. D. Shi , J. J. Song, Q. Q. Song, W.M. Song, Y. X. Song, S. Sosio, S. Spataro, F. F. Sui, G. X. Sun, J. F. Sun, L. Sun, S. S. Sun, T. Sun, W. Y. Sun, X. Sun, Y. J. Sun, Y. K. Sun, Y. Z. Sun, Z. T. Sun, Y. H. Tan, Y. X. Tan, C. J. Tang, G. Y. Tang, J. Tang, V. Thoren, B. Tsednee, I. Uman, B. Wang, B. L. Wang, C.W. Wang, D. Y. Wang, H. P. Wang, K. Wang, L. L. Wang, M. Wang, M. Z. Wang, Meng Wang, W. H. Wang, W. P. Wang, X. Wang, X. F. Wang, X. L. Wang, Y. Wang, Y. Wang, Y. D. Wang, Y. F. Wang, Y. Q. Wang, Z. Wang, Z. Y. Wang, Ziyi Wang, Zongyuan Wang, D. H. Wei, P. Weidenkaff, F. Weidner, S. P. Wen, D. J. White, U. Wiedner, G. Wilkinson, M. Wolke, L. Wollenberg, J. F. Wu, L. H. Wu, L. J. Wu, X. Wu, Z. Wu, L. Xia, H. Xiao, S. Y. Xiao, Y. J. Xiao, Z. J. Xiao, X. H. Xie, Y. G. Xie, Y. H. Xie, T. Y. Xing, X. A. Xiong, G. F. Xu, J. J. Xu, Q. J. Xu, W. Xu, X. P. Xu, L. Yan, L. Yan, W. B. Yan, W. C. Yan, Xu Yan, H. J. Yang, H. X. Yang, L. Yang, R. X. Yang, S. L. Yang, Y. H. Yang, Y. X. Yang, Yifan Yang, Zhi Yang, M. Ye, M. H. Ye, J. H. Yin, Z. Y. You, B. X. Yu, C. X. Yu, G. Yu, J. S. Yu, T. Yu, C. Z. Yuan, W. Yuan, X. Q. Yuan, Y. Yuan, Z. Y. Yuan, C. X. Yue, A. Yuncu, A. A. Zafar, Y. Zeng, B. X. Zhang, Guangyi Zhang, H. H. Zhang, H. Y. Zhang, J. L. Zhang, J. Q. Zhang, J. W. Zhang, J. Y. Zhang, J. Z. Zhang, Jianyu Zhang, Jiawei Zhang, L. Zhang, Lei Zhang, S. Zhang, S. F. Zhang, T. J. Zhang, X. Y. Zhang, Y. Zhang, PHYSICAL REVIEW D 101, 112005 (2020)


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][16]. 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 [17], 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].
In this paper, we search for the J=ψ → γ þ invisible decay using J=ψ produced through the process ψð3686Þ → π þ π − J=ψ in a data sample of ð448.1 AE 2.9Þ × 10 6 ψð3686Þ decays collected with the BESIII detector.

II. BESIII DETECTOR AND MONTE CARLO SIMULATION
The BESIII detector is a magnetic spectrometer [23] located at the Beijing Electron Positron Collider (BEPCII) [24]. 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 charged-particle 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 [25] 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 [26], taking into account the beam energy spread; the known decay modes are generated using EVTGEN [27] with the branching fractions as given by the particle data group (PDG) [3]; the unknown decay modes are modeled with the LUNDCHARM model [28]. 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. [29]. 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 [30] 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 exclusively with different generators, i.e., J=ψ → γη 0 , γηð1405Þ and γη c with the angular distribution of 1 þ cos 2 θ γ ; J=ψ → γπ 0 π 0 and γπ þ π − according to PWA results of J=ψ → γπ 0 π 0 [31] 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 [30], 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 semiblind 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 AE10 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 j cos θj < 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 fulfill j cos θ 2 j < 0.95.
To identify ψð3686Þ → π þ π − J=ψ candidate events, the recoiling mass of the π þ π − system, M rec , is used, 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 2nd 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 AE 1Þ × 10 4 , obtained by integrating the fitted signal curve in the J=ψ signal region. By performing same procedure on the inclusive MC sample, the efficiency for tagging J=ψ is determined as ð56.80 AE 0.01Þ%.

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 (j cos θj < 0.80) or 50 MeV for endcap showers (0.86 < j cos θj < 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 exclusive 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 [32], 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 , corresponding to a mass from 0 up to 1.2 GeV=c 2 for the invisible particle. In the fit, the signal shape is taken from the signal MC simulation convoluted with a Gaussian function representing the resolution difference between data and MC, where the parameters of the Gaussian function are obtained by studying a clean control sample of ψð3686Þ → π þ π − J=ψ; J=ψ → γηðη → γγÞ. The background shape is described by the sum of an exponential function and two crystal ball functions with fixed amplitudes and shapes presenting for the background of ψð3686Þ → π þ π − J=ψ with subsequent decay J=ψ → γη and γπ 0 , respectively, where amplitudes and shapes are estimated by the MC simulation, and the same correction on shape as the signal description is implemented. (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 [33,34] combined with the asymptotic approximation [35]. In this approach, the test statistic is the profile likelihood ratio, where the likelihood is given with the Poisson function: where s i represents the signal probability in the ith bin, B is the branching fraction BðJ=ψ → γ þ invisibleÞ, b exp ij is the expected background number in the ith bin for the jth source. Here background is modeled with the exponential function and the two fixed crystal functions from zerosignal assumption fit result. Additionally, systematic uncertainties are included assuming 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 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 AE5 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 3rd 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% [36]. 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 .
The uncertainties due to the continuum background, representing by the exponential function, are also included. Both the shape and magnitude are considered, and the corresponding uncertainties are evaluated by performing a fit on E Ã γ distribution with zero-signal assumption. The all discussed systematic uncertainties are listed in the Table I. Fit procedure Number of ψð3686Þ → π þ π − J=ψ; J=ψ → γη 17% Number of ψð3686Þ → π þ π − J=ψ; J=ψ → γπ 0 17% Number of continuum background 4.4%

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 with most mass assumptions. And for the zero mass assumption of the invisible particle the upper limit is 7.0 × 10 −7 . The local signal significances with different mass assumptions are also shown in Fig. 3, where the local signal significance is calculated by 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 AE 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 [18].
To further investigate the physical parameters in NMSSM, and to better compare the physical results from the different quarkonium decays, according to Ref. [16] and Eq. (1), the upper limits of g c × tan 2 β × ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi BðA 0 → invisibleÞ p 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 × ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi BðA 0 → invisibleÞ p ð¼g c × tan 2 β × ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi BðA 0 → invisibleÞ p Þ, which is obtained based on the Belle results [21] 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 [21], we also extract upper limits on cos θ A ð¼ ffiffiffiffiffiffiffiffi ffi g b g c p Þ × ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi BðA 0 → invisibleÞ p , as presented in Fig. 4(b).

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