# Search for $X(3872)\to\pi^{0}\chi_{c0}$ and $X(3872)\to\pi\pi\chi_{c0}$ at
BESIII

M. Ablikim1, M. N. Achasov10,b, P. Adlarson67, M. Albrecht4, R. Aliberti28, A.
Amoroso66A,66C, M. R. An32, Q. An63,50, X. H. Bai58, Y. Bai49, O. Bakina29, R.
Baldini Ferroli23A, I. Balossino24A, Y. Ban39,g, V. Batozskaya1,37, D.
Becker28, K. Begzsuren26, N. Berger28, M. Bertani23A, D. Bettoni24A, F.
Bianchi66A,66C, J. Bloms60, A. Bortone66A,66C, I. Boyko29, R. A. Briere5, A.
Brueggemann60, H. Cai68, X. Cai1,50, A. Calcaterra23A, G. F. Cao1,55, N.
Cao1,55, S. A. Cetin54A, J. F. Chang1,50, W. L. Chang1,55, G. Chelkov29,a, C.
Chen36, Chao Chen47, G. Chen1, H. S. Chen1,55, M. L. Chen1,50, S. J. Chen35,
S. M. Chen53, T. Chen1, X. R. Chen25,55, X. T. Chen1, Y. B. Chen1,50, Z. J.
Chen20,h, W. S. Cheng66C, X. Chu36, G. Cibinetto24A, F. Cossio66C, J. J.
Cui42, H. L. Dai1,50, J. P. Dai70, A. Dbeyssi14, R. E. de Boer4, D.
Dedovich29, Z. Y. Deng1, A. Denig28, I. Denysenko29, M. Destefanis66A,66C, F.
De Mori66A,66C, Y. Ding33, J. Dong1,50, L. Y. Dong1,55, M. Y. Dong1,50,55, X.
Dong68, S. X. Du72, P. Egorov29,a, Y. L. Fan68, J. Fang1,50, S. S. Fang1,55,
W. X. Fang1, Y. Fang1, R. Farinelli24A, L. Fava66B,66C, F. Feldbauer4, G.
Felici23A, C. Q. Feng63,50, J. H. Feng51, K Fischer61, M. Fritsch4, C.
Fritzsch60, C. D. Fu1, H. Gao55, Y. N. Gao39,g, Yang Gao63,50, S.
Garbolino66C, I. Garzia24A,24B, P. T. Ge68, Z. W. Ge35, C. Geng51, E. M.
Gersabeck59, A Gilman61, K. Goetzen11, L. Gong33, W. X. Gong1,50, W. Gradl28,
M. Greco66A,66C, L. M. Gu35, M. H. Gu1,50, C. Y Guan1,55, A. Q. Guo25,55, L.
B. Guo34, R. P. Guo41, Y. P. Guo9,f, A. Guskov29,a, T. T. Han42, W. Y. Han32,
X. Q. Hao15, F. A. Harris57, K. K. He47, K. L. He1,55, F. H. Heinsius4, C. H.
Heinz28, Y. K. Heng1,50,55, C. Herold52, M. Himmelreich11,d, G. Y. Hou1,55, Y.
R. Hou55, Z. L. Hou1, H. M. Hu1,55, J. F. Hu48,i, T. Hu1,50,55, Y. Hu1, G. S.
Huang63,50, K. X. Huang51, L. Q. Huang64, L. Q. Huang25,55, X. T. Huang42, Y.
P. Huang1, Z. Huang39,g, T. Hussain65, N Hüsken22,28, W. Imoehl22, M.
Irshad63,50, J. Jackson22, S. Jaeger4, S. Janchiv26, Q. Ji1, Q. P. Ji15, X. B.
Ji1,55, X. L. Ji1,50, Y. Y. Ji42, Z. K. Jia63,50, H. B. Jiang42, S. S.
Jiang32, X. S. Jiang1,50,55, Y. Jiang55, J. B. Jiao42, Z. Jiao18, S. Jin35, Y.
Jin58, M. Q. Jing1,55, T. Johansson67, N. Kalantar-Nayestanaki56, X. S.
Kang33, R. Kappert56, M. Kavatsyuk56, B. C. Ke72, I. K. Keshk4, A. Khoukaz60,
P. Kiese28, R. Kiuchi1, R. Kliemt11, L. Koch30, O. B. Kolcu54A, B. Kopf4, M.
Kuemmel4, M. Kuessner4, A. Kupsc37,67, W. Kühn30, J. J. Lane59, J. S. Lange30,
P. Larin14, A. Lavania21, L. Lavezzi66A,66C, Z. H. Lei63,50, H. Leithoff28, M.
Lellmann28, T. Lenz28, C. Li40, C. Li36, C. H. Li32, Cheng Li63,50, D. M.
Li72, F. Li1,50, G. Li1, H. Li63,50, H. Li44, H. B. Li1,55, H. J. Li15, H. N.
Li48,i, J. Q. Li4, J. S. Li51, J. W. Li42, Ke Li1, L. J Li1, L. K. Li1, Lei
Li3, M. H. Li36, P. R. Li31,j,k, S. X. Li9, S. Y. Li53, T. Li42, W. D. Li1,55,
W. G. Li1, X. H. Li63,50, X. L. Li42, Xiaoyu Li1,55, H. Liang63,50, H.
Liang27, H. Liang1,55, Y. F. Liang46, Y. T. Liang25,55, G. R. Liao12, L. Z.
Liao42, J. Libby21, A. Limphirat52, C. X. Lin51, D. X. Lin25,55, T. Lin1, B.
J. Liu1, C. X. Liu1, D. Liu14,63, F. H. Liu45, Fang Liu1, Feng Liu6, G. M.
Liu48,i, H. Liu31,j,k, H. M. Liu1,55, Huanhuan Liu1, Huihui Liu16, J. B.
Liu63,50, J. L. Liu64, J. Y. Liu1,55, K. Liu1, K. Y. Liu33, Ke Liu17, L.
Liu63,50, M. H. Liu9,f, P. L. Liu1, Q. Liu55, S. B. Liu63,50, T. Liu9,f, W. K.
Liu36, W. M. Liu63,50, X. Liu31,j,k, Y. Liu31,j,k, Y. B. Liu36, Z. A.
Liu1,50,55, Z. Q. Liu42, X. C. Lou1,50,55, F. X. Lu51, H. J. Lu18, J. G.
Lu1,50, X. L. Lu1, Y. Lu1, Y. P. Lu1,50, Z. H. Lu1, C. L. Luo34, M. X. Luo71,
T. Luo9,f, X. L. Luo1,50, X. R. Lyu55, Y. F. Lyu36, F. C. Ma33, H. L. Ma1, L.
L. Ma42, M. M. Ma1,55, Q. M. Ma1, R. Q. Ma1,55, R. T. Ma55, X. Y. Ma1,50, Y.
Ma39,g, F. E. Maas14, M. Maggiora66A,66C, S. Maldaner4, S. Malde61, Q. A.
Malik65, A. Mangoni23B, Y. J. Mao39,g, Z. P. Mao1, S. Marcello66A,66C, Z. X.
Meng58, J. G. Messchendorp56,11, G. Mezzadri24A, H. Miao1, T. J. Min35, R. E.
Mitchell22, X. H. Mo1,50,55, N. Yu. Muchnoi10,b, Y. Nefedov29, I. B.
Nikolaev10,b, Z. Ning1,50, S. Nisar8,l, Y. Niu 42, S. L. Olsen55, Q.
Ouyang1,50,55, S. Pacetti23B,23C, X. Pan9,f, Y. Pan49, A. Pathak1, A.
Pathak27, M. Pelizaeus4, H. P. Peng63,50, K. Peters11,d, J. Pettersson67, J.
L. Ping34, R. G. Ping1,55, S. Plura28, S. Pogodin29, V. Prasad63,50, F. Z.
Qi1, H. Qi63,50, H. R. Qi53, M. Qi35, T. Y. Qi9,f, S. Qian1,50, W. B. Qian55,
Z. Qian51, C. F. Qiao55, J. J. Qin64, L. Q. Qin12, X. P. Qin9,f, X. S. Qin42,
Z. H. Qin1,50, J. F. Qiu1, S. Q. Qu53, S. Q. Qu36, K. H. Rashid65, C. F.
Redmer28, K. J. Ren32, A. Rivetti66C, V. Rodin56, M. Rolo66C, G. Rong1,55, Ch.
Rosner14, S. N. Ruan36, H. S. Sang63, A. Sarantsev29,c, Y. Schelhaas28, C.
Schnier4, K. Schoenning67, M. Scodeggio24A,24B, K. Y. Shan9,f, W. Shan19, X.
Y. Shan63,50, J. F. Shangguan47, L. G. Shao1,55, M. Shao63,50, C. P. Shen9,f,
H. F. Shen1,55, X. Y. Shen1,55, B.-A. Shi55, H. C. Shi63,50, J. Y. Shi1, q. q.
Shi47, R. S. Shi1,55, X. Shi1,50, X. D Shi63,50, J. J. Song15, W. M. Song27,1,
Y. X. Song39,g, S. Sosio66A,66C, S. Spataro66A,66C, F. Stieler28, K. X. Su68,
P. P. Su47, Y.-J. Su55, G. X. Sun1, H. Sun55, H. K. Sun1, J. F. Sun15, L.
Sun68, S. S. Sun1,55, T. Sun1,55, W. Y. Sun27, X Sun20,h, Y. J. Sun63,50, Y.
Z. Sun1, Z. T. Sun42, Y. H. Tan68, Y. X. Tan63,50, C. J. Tang46, G. Y. Tang1,
J. Tang51, L. Y Tao64, Q. T. Tao20,h, M. Tat61, J. X. Teng63,50, V. Thoren67,
W. H. Tian44, Y. Tian25,55, I. Uman54B, B. Wang1, B. L. Wang55, C. W. Wang35,
D. Y. Wang39,g, F. Wang64, H. J. Wang31,j,k, H. P. Wang1,55, K. Wang1,50, L.
L. Wang1, M. Wang42, M. Z. Wang39,g, Meng Wang1,55, S. Wang9,f, T. Wang9,f, T.
J. Wang36, W. Wang51, W. H. Wang68, W. P. Wang63,50, X. Wang39,g, X. F.
Wang31,j,k, X. L. Wang9,f, Y. D. Wang38, Y. F. Wang1,50,55, Y. H. Wang40, Y.
Q. Wang1, Yi2020 Wang53, Z. Wang1,50, Z. Y. Wang1,55, Ziyi Wang55, D. H.
Wei12, F. Weidner60, S. P. Wen1, D. J. White59, U. Wiedner4, G. Wilkinson61,
M. Wolke67, L. Wollenberg4, J. F. Wu1,55, L. H. Wu1, L. J. Wu1,55, X. Wu9,f,
X. H. Wu27, Y. Wu63, Z. Wu1,50, L. Xia63,50, T. Xiang39,g, D. Xiao31,j,k, G.
Y. Xiao35, H. Xiao9,f, S. Y. Xiao1, Y. L. Xiao9,f, Z. J. Xiao34, C. Xie35, X.
H. Xie39,g, Y. Xie42, Y. G. Xie1,50, Y. H. Xie6, Z. P. Xie63,50, T. Y.
Xing1,55, C. F. Xu1, C. J. Xu51, G. F. Xu1, H. Y. Xu58, Q. J. Xu13, S. Y.
Xu62, X. P. Xu47, Y. C. Xu55, Z. P. Xu35, F. Yan9,f, L. Yan9,f, W. B.
Yan63,50, W. C. Yan72, H. J. Yang43,e, H. L. Yang27, H. X. Yang1, L. Yang44,
S. L. Yang55, Tao Yang1, Y. X. Yang1,55, Yifan Yang1,55, M. Ye1,50, M. H. Ye7,
J. H. Yin1, Z. Y. You51, B. X. Yu1,50,55, C. X. Yu36, G. Yu1,55, T. Yu64, C.
Z. Yuan1,55, L. Yuan2, S. C. Yuan1, X. Q. Yuan1, Y. Yuan1,55, Z. Y. Yuan51, C.
X. Yue32, A. A. Zafar65, F. R. Zeng42, X. Zeng Zeng6, Y. Zeng20,h, Y. H.
Zhan51, A. Q. Zhang1, B. L. Zhang1, B. X. Zhang1, D. H. Zhang36, G. Y.
Zhang15, H. Zhang63, H. H. Zhang27, H. H. Zhang51, H. Y. Zhang1,50, J. L.
Zhang69, J. Q. Zhang34, J. W. Zhang1,50,55, J. X. Zhang31,j,k, J. Y. Zhang1,
J. Z. Zhang1,55, Jianyu Zhang1,55, Jiawei Zhang1,55, L. M. Zhang53, L. Q.
Zhang51, Lei Zhang35, P. Zhang1, Q. Y. Zhang32,72, Shulei Zhang20,h, X. D.
Zhang38, X. M. Zhang1, X. Y. Zhang47, X. Y. Zhang42, Y. Zhang61, Y. T.
Zhang72, Y. H. Zhang1,50, Yan Zhang63,50, Yao Zhang1, Z. H. Zhang1, Z. Y.
Zhang36, Z. Y. Zhang68, G. Zhao1, J. Zhao32, J. Y. Zhao1,55, J. Z. Zhao1,50,
Lei Zhao63,50, Ling Zhao1, M. G. Zhao36, Q. Zhao1, S. J. Zhao72, Y. B.
Zhao1,50, Y. X. Zhao25,55, Z. G. Zhao63,50, A. Zhemchugov29,a, B. Zheng64, J.
P. Zheng1,50, Y. H. Zheng55, B. Zhong34, C. Zhong64, X. Zhong51, H. Zhou42, L.
P. Zhou1,55, X. Zhou68, X. K. Zhou55, X. R. Zhou63,50, X. Y. Zhou32, Y. Z.
Zhou9,f, J. Zhu36, K. Zhu1, K. J. Zhu1,50,55, L. X. Zhu55, S. H. Zhu62, S. Q.
Zhu35, T. J. Zhu69, W. J. Zhu9,f, Y. C. Zhu63,50, Z. A. Zhu1,55, B. S. Zou1,
J. H. Zou1 (BESIII Collaboration) 1 Institute of High Energy Physics, Beijing
100049, People’s Republic of China
2 Beihang University, Beijing 100191, People’s Republic of China
3 Beijing Institute of Petrochemical Technology, Beijing 102617, People’s
Republic of China
4 Bochum Ruhr-University, D-44780 Bochum, Germany
5 Carnegie Mellon University, Pittsburgh, Pennsylvania 15213, USA
6 Central China Normal University, Wuhan 430079, People’s Republic of China
7 China Center of Advanced Science and Technology, Beijing 100190, People’s
Republic of China
8 COMSATS University Islamabad, Lahore Campus, Defence Road, Off Raiwind Road,
54000 Lahore, Pakistan
9 Fudan University, Shanghai 200433, People’s Republic of China
10 G.I. Budker Institute of Nuclear Physics SB RAS (BINP), Novosibirsk 630090,
Russia
11 GSI Helmholtzcentre for Heavy Ion Research GmbH, D-64291 Darmstadt, Germany
12 Guangxi Normal University, Guilin 541004, People’s Republic of China
13 Hangzhou Normal University, Hangzhou 310036, People’s Republic of China
14 Helmholtz Institute Mainz, Staudinger Weg 18, D-55099 Mainz, Germany
15 Henan Normal University, Xinxiang 453007, People’s Republic of China
16 Henan University of Science and Technology, Luoyang 471003, People’s
Republic of China
17 Henan University of Technology, Zhengzhou 450001, People’s Republic of
China
18 Huangshan College, Huangshan 245000, People’s Republic of China
19 Hunan Normal University, Changsha 410081, People’s Republic of China
20 Hunan University, Changsha 410082, People’s Republic of China
21 Indian Institute of Technology Madras, Chennai 600036, India
22 Indiana University, Bloomington, Indiana 47405, USA
23 INFN Laboratori Nazionali di Frascati , (A)INFN Laboratori Nazionali di
Frascati, I-00044, Frascati, Italy; (B)INFN Sezione di Perugia, I-06100,
Perugia, Italy; (C)University of Perugia, I-06100, Perugia, Italy
24 INFN Sezione di Ferrara, (A)INFN Sezione di Ferrara, I-44122, Ferrara,
Italy; (B)University of Ferrara, I-44122, Ferrara, Italy
25 Institute of Modern Physics, Lanzhou 730000, People’s Republic of China
26 Institute of Physics and Technology, Peace Ave. 54B, Ulaanbaatar 13330,
Mongolia
27 Jilin University, Changchun 130012, People’s Republic of China
28 Johannes Gutenberg University of Mainz, Johann-Joachim-Becher-Weg 45,
D-55099 Mainz, Germany
29 Joint Institute for Nuclear Research, 141980 Dubna, Moscow region, Russia
30 Justus-Liebig-Universitaet Giessen, II. Physikalisches Institut, Heinrich-
Buff-Ring 16, D-35392 Giessen, Germany
31 Lanzhou University, Lanzhou 730000, People’s Republic of China
32 Liaoning Normal University, Dalian 116029, People’s Republic of China
33 Liaoning University, Shenyang 110036, People’s Republic of China
34 Nanjing Normal University, Nanjing 210023, People’s Republic of China
35 Nanjing University, Nanjing 210093, People’s Republic of China
36 Nankai University, Tianjin 300071, People’s Republic of China
37 National Centre for Nuclear Research, Warsaw 02-093, Poland
38 North China Electric Power University, Beijing 102206, People’s Republic of
China
39 Peking University, Beijing 100871, People’s Republic of China
40 Qufu Normal University, Qufu 273165, People’s Republic of China
41 Shandong Normal University, Jinan 250014, People’s Republic of China
42 Shandong University, Jinan 250100, People’s Republic of China
43 Shanghai Jiao Tong University, Shanghai 200240, People’s Republic of China
44 Shanxi Normal University, Linfen 041004, People’s Republic of China
45 Shanxi University, Taiyuan 030006, People’s Republic of China
46 Sichuan University, Chengdu 610064, People’s Republic of China
47 Soochow University, Suzhou 215006, People’s Republic of China
48 South China Normal University, Guangzhou 510006, People’s Republic of China
49 Southeast University, Nanjing 211100, People’s Republic of China
50 State Key Laboratory of Particle Detection and Electronics, Beijing 100049,
Hefei 230026, People’s Republic of China
51 Sun Yat-Sen University, Guangzhou 510275, People’s Republic of China
52 Suranaree University of Technology, University Avenue 111, Nakhon
Ratchasima 30000, Thailand
53 Tsinghua University, Beijing 100084, People’s Republic of China
54 Turkish Accelerator Center Particle Factory Group, (A)Istinye University,
34010, Istanbul, Turkey; (B)Near East University, Nicosia, North Cyprus,
Mersin 10, Turkey
55 University of Chinese Academy of Sciences, Beijing 100049, People’s
Republic of China
56 University of Groningen, NL-9747 AA Groningen, The Netherlands
57 University of Hawaii, Honolulu, Hawaii 96822, USA
58 University of Jinan, Jinan 250022, People’s Republic of China
59 University of Manchester, Oxford Road, Manchester, M13 9PL, United Kingdom
60 University of Muenster, Wilhelm-Klemm-Str. 9, 48149 Muenster, Germany
61 University of Oxford, Keble Rd, Oxford, UK OX13RH
62 University of Science and Technology Liaoning, Anshan 114051, People’s
Republic of China
63 University of Science and Technology of China, Hefei 230026, People’s
Republic of China
64 University of South China, Hengyang 421001, People’s Republic of China
65 University of the Punjab, Lahore-54590, Pakistan
66 University of Turin and INFN, (A)University of Turin, I-10125, Turin,
Italy; (B)University of Eastern Piedmont, I-15121, Alessandria, Italy;
(C)INFN, I-10125, Turin, Italy
67 Uppsala University, Box 516, SE-75120 Uppsala, Sweden
68 Wuhan University, Wuhan 430072, People’s Republic of China
69 Xinyang Normal University, Xinyang 464000, People’s Republic of China
70 Yunnan University, Kunming 650500, People’s Republic of China
71 Zhejiang University, Hangzhou 310027, People’s Republic of China
72 Zhengzhou University, Zhengzhou 450001, People’s Republic of China
a Also at the Moscow Institute of Physics and Technology, Moscow 141700,
Russia
b Also at the Novosibirsk State University, Novosibirsk, 630090, Russia
c Also at the NRC ”Kurchatov Institute”, PNPI, 188300, Gatchina, Russia
d Also at Goethe University Frankfurt, 60323 Frankfurt am Main, Germany
e Also at Key Laboratory for Particle Physics, Astrophysics and Cosmology,
Ministry of Education; Shanghai Key Laboratory for Particle Physics and
Cosmology; Institute of Nuclear and Particle Physics, Shanghai 200240,
People’s Republic of China
f Also at Key Laboratory of Nuclear Physics and Ion-beam Application (MOE) and
Institute of Modern Physics, Fudan University, Shanghai 200443, People’s
Republic of China
g Also at State Key Laboratory of Nuclear Physics and Technology, Peking
University, Beijing 100871, People’s Republic of China
h Also at School of Physics and Electronics, Hunan University, Changsha
410082, China
i Also at Guangdong Provincial Key Laboratory of Nuclear Science, Institute of
Quantum Matter, South China Normal University, Guangzhou 510006, China
j Also at Frontiers Science Center for Rare Isotopes, Lanzhou University,
Lanzhou 730000, People’s Republic of China
k Also at Lanzhou Center for Theoretical Physics, Lanzhou University, Lanzhou
730000, People’s Republic of China
l Also at the Department of Mathematical Sciences, IBA, Karachi , Pakistan

(February 1, 2022)

###### Abstract

Using 9.9 fb-1 of $e^{+}e^{-}$ collision data collected by the BESIII detector
at center-of-mass energies between 4.15 and 4.30 GeV, we search for the
processes $e^{+}e^{-}\to\gamma X(3872)$ with
$X(3872)\rightarrow\pi^{0}\chi_{c0}$ and $X(3872)\rightarrow\pi\pi\chi_{c0}$.
Depending on the fitting model, the statistical significance for
$X(3872)\to\pi^{0}\chi_{c0}$ ranges from 1.3$\sigma$ to 2.8$\sigma$. We set
upper limits (at 90% C.L.) of
$\frac{\mathcal{B}(X(3872)\rightarrow\pi^{0}\chi_{c0})}{\mathcal{B}(X(3872)\to\pi^{+}\pi^{-}J/\psi)}$$<3.6$,
$\frac{\mathcal{B}(X(3872)\rightarrow\pi^{+}\pi^{-}\chi_{c0})}{\mathcal{B}(X(3872)\to\pi^{+}\pi^{-}J/\psi)}$$<0.68$,
and
$\frac{\mathcal{B}(X(3872)\rightarrow\pi^{0}\pi^{0}\chi_{c0})}{\mathcal{B}(X(3872)\to\pi^{+}\pi^{-}J/\psi)}$$<1.7$.
Combined with the BESIII measurement of $X(3872)\to\pi^{0}\chi_{c1}$, we also
set an upper limit of
$\frac{\mathcal{B}(X(3872)\rightarrow\pi^{0}\chi_{c0})}{\mathcal{B}(X(3872)\to\pi^{0}\chi_{c1})}$$<4.4$.

##  I  Introduction

The $X(3872)$, discovered by the Belle experiment in 2003 [1], was the first
state in the charmonium region that could not be easily explained by a simple
$c\bar{c}$ model. Despite having quantum numbers $J^{PC}=1^{++}$ [2] and a
mass near the predicted $\chi_{c1}(2P)$ mass of 3.95 GeV/$c^{2}$ [3], the
state has several properties that cannot be explained by a pure charmonium
state above open charm threshold. The state is exceptionally narrow, with a
measured width of $0.96^{+0.19}_{-0.18}\pm 0.21$ MeV [4] or $1.39\pm 0.24\pm
0.10$ MeV [5], depending on the assumed lineshape. In addition, the state has
large isospin violation effects in its decays, which is clearly seen in the
fact that the decays $X(3872)\to\rho J/\psi$ and $X(3872)\to\omega J/\psi$
occur at approximately the same rate. Many $X(3872)$ decay modes have been
observed, including $\rho J/\psi$ [6], $D^{0}\bar{D}^{*0}$ [7], $\gamma
J/\psi$ [8], $\pi^{0}\chi_{c1}$ [9] and $\omega J/\psi$ [10]. However the
nature of the $X(3872)$ remains unclear. Since the mass of this state is near
$D^{0}\bar{D}^{*0}$ threshold [11], one explanation of its exotic properties
is that the state has a $D^{0}\bar{D}^{*0}$ molecular component.

Table 1: Theoretical predictions on $\frac{\mathcal{B}(X(3872)\rightarrow\pi^{0}\chi_{c0})}{\mathcal{B}(X(3872)\to\pi^{+}\pi^{-}J/\psi)}$ and $\frac{\mathcal{B}(X(3872)\rightarrow\pi^{0}\chi_{c0})}{\mathcal{B}(X(3872)\to\pi^{0}\chi_{c1})}$ for different physical interpretations of the $X(3872)$ state.  Ref | Technique | Interpretation | $\frac{\mathcal{B}(X(3872)\rightarrow\pi^{0}\chi_{c0})}{\mathcal{B}(X(3872)\to\pi^{+}\pi^{-}J/\psi)}$ | $\frac{\mathcal{B}(X(3872)\rightarrow\pi^{0}\chi_{c0})}{\mathcal{B}(X(3872)\to\pi^{0}\chi_{c1})}$
---|---|---|---|---
[12] | Multipole expansion | Four-quark/molecule | NA | 2.97
[12] | Multipole expansion | $\chi_{c1}(2P)$ | 0.0 | 0.0
[13] | Effective field theory | $D^{0}D^{0*}$ | NA | 2.84 to 2.98
[14] | Effective field theory | $D^{0}\overline{D}^{0*}+D^{+}D^{-*}$ | 1.3 to 2.07 | 1.65 to 1.77
[15] | Effective field theory | $D^{0}\overline{D}^{0*}+D^{+}D^{-*}$ | NA | 3.72
[16] | Effective field theory |  $D^{0}\overline{D}^{0*}+D^{+}D^{-*}$+$\chi_{c1}(2P)$ | 0.094 | 1.15

By searching for new decay modes, we can learn about the quark configuration
of the $X(3872)$ state. The ratios
$\frac{\mathcal{B}(X(3872)\rightarrow\pi^{0}\chi_{c0})}{\mathcal{B}(X(3872)\to\pi^{+}\pi^{-}J/\psi)}$
and
$\frac{\mathcal{B}(X(3872)\rightarrow\pi^{0}\chi_{c0})}{\mathcal{B}(X(3872)\to\pi^{0}\chi_{c1})}$
are expected to be sensitive to different physical interpretations of the
$X(3872)$. Theoretical predictions under the hypothesis of a pure charmonium
state [12], a generic four quark state [12], a $D^{0}\bar{D}^{0*}$ molecule
[13], a $D\bar{D}^{*}$ molecule with charged and neutral components [14, 15],
or a combination of a molecular state and charmonium state [16] are summarized
in Table 1. Since BESIII recently measured
$\frac{\mathcal{B}(X(3872)\rightarrow\pi^{0}\chi_{c1})}{\mathcal{B}(X(3872)\to\pi^{+}\pi^{-}J/\psi)}$$=0.88^{+0.33}_{-0.27}\pm
0.10$ [9], we can measure both of these quantities at BESIII. The Belle
Collaboration has also performed a search for $X(3872)\to\pi^{0}\chi_{c1}$,
and they set an upper limit on
$\frac{\mathcal{B}(X(3872)\rightarrow\pi^{0}\chi_{c1})}{\mathcal{B}(X(3872)\to\pi^{+}\pi^{-}J/\psi)}$
of 0.97 at 90% confidence level [17], which is consistent with the measurement
of BESIII [9].

In this paper, we search for the process $e^{+}e^{-}\to\gamma X(3872)$ with
$X(3872)\to\pi^{0}\chi_{c0}$ and the $\chi_{c0}$ decaying hadronically to the
final states shown in Table 2. These final states are chosen because they have
large branching fractions and can be reconstructed with a high efficiency. In
addition, we search for the double pion transitions
$X(3872)\to\pi^{+}\pi^{-}\chi_{c0}$ and $X(3872)\to\pi^{0}\pi^{0}\chi_{c0}$
through the same hadronic decays of the $\chi_{c0}$. Based on a molecular
interpretation for the $X(3872)$, effective field theory (EFT) calculations
predict that
$\frac{\mathcal{B}(X(3872)\to\pi\pi\chi_{c0})}{\mathcal{B}(X(3872)\to\pi^{0}\chi_{c0})}\approx\mathcal{O}(10^{-3})$
[15] or $\mathcal{O}(10^{-5})$ [13], depending on the specific methods used.

Table 2:  Hadronic decay modes of the $\chi_{c0}$ reconstructed in this analysis, with branching fractions from Ref. [11]. Decay | Branching Fraction (%)
---|---
$\chi_{c0}\to\pi^{+}\pi^{-}$ | $0.567\pm 0.022$
$\chi_{c0}\to K^{+}K^{-}$ | $0.605\pm 0.031$
$\chi_{c0}\to\pi^{+}\pi^{-}\pi^{+}\pi^{-}$ | $2.34\pm 0.18$
$\chi_{c0}\to\pi^{+}\pi^{-}\pi^{0}\pi^{0}$ | $3.3\pm 0.4$
$\chi_{c0}\to\pi^{+}\pi^{-}K^{+}K^{-}$ | $1.81\pm 0.14$

We normalize our results to the process $e^{+}e^{-}\to\gamma X(3872)$ with
$X(3872)\to\pi^{+}\pi^{-}J/\psi$ with $J/\psi\to\ell^{+}\ell^{-}$, for
$\ell=e$ or $\ell=\mu$. This normalization channel is chosen because it has
large statistics and is easy to reconstruct. We measure the ratio
$\frac{\mathcal{B}(X(3872)\rightarrow\pi^{0}\chi_{c0})}{\mathcal{B}(X(3872)\to\pi^{0}\chi_{c1})}$
by using the efficiency and fit results for $X(3872)\to\pi^{0}\chi_{c1}$ from
Ref. [9]. In all cases, the $X(3872)$ is produced through $e^{+}e^{-}\to\gamma
X(3872)$. Since this is always the same for the search and normalization
channels, the production cross section cancels in the ratio, so we can combine
data sets from different energy points in our analysis. The normalizations
also cancel several systematic uncertainties.

##  II Experimental Background

The BESIII experiment, operating at the Beijing Electron Positron Collider
(BEPCII), has measured $\sigma(e^{+}e^{-}\to\gamma
X(3872))\mathcal{B}(X(3872)\to\pi^{+}\pi^{-}J/\psi)$ for center-of-mass
energies ($E_{\rm CM}$) between 4.008 and 4.60 GeV and found that the cross
section is the largest for 4.15$<E_{\rm CM}<$4.3 GeV [10]. BESIII has 9.9 fb-1
of data in this energy region, at the energies shown in Table 3. Center-of-
mass energies are measured in Refs. [18] and [19], while the luminosities are
measured in Refs. [20] and [21]. This data set makes it possible to search for
complicated decay modes of the $X(3872)$ state.

BESIII has also accumulated 5.3 fb-1 of data with $E_{\rm CM}>4.3$ GeV at the
energies shown in Table 4, which can be used to check that any potential
$X(3872)$ signal for $4.15<E_{\rm CM}<4.3$ GeV does not also occur outside of
the $e^{+}e^{-}\to\gamma X(3872)$ production region.

Table 3: Data with $4.15<E_{\rm CM}<4.3$ GeV, where the $X(3872)$ is produced via $e^{+}e^{-}\to\gamma X(3872)$. All of these energy points are used in our nominal fit. The total luminosity is 9.9 fb-1. Only statistical uncertainties are shown. Entries that do not include a reference are estimated values. The last line of the table shows the range of luminosity values for 29 different energy points that each have a much smaller luminosity than the other data points. Luminosity (pb-1) |  $E_{\rm CM}$ (MeV) | Year
---|---|---
401.5 | 4157.83$\pm 0.05$ [19] | 2019
3189.0 | 4178 | 2016
43.09 |  $\pm$0.03 [20] | 4188.59$\pm 0.15$ [18] | 2013
526.70 |  $\pm 2.16$ [19] | 4189.12$\pm 0.05$ [19] | 2017
526.60 |  $\pm 2.05$ [19] | 4199.15$\pm 0.06$ [19] | 2017
54.55 |  $\pm$0.03 [20] | 4207.73$\pm 0.14$ [18] | 2013
517.10 |  $\pm 1.81$ [19] | 4209.39$\pm 0.06$ [19] | 2017
54.13 |  $\pm$0.03 [20] | 4217.13$\pm 0.14$ [18] | 2013
514.60 |  $\pm 1.80$ [19] | 4218.83$\pm 0.06$ [19] | 2017
1047.34 |  $\pm$0.14 [20] | 4226.26$\pm 0.04$ [18] | 2013
44.40 |  $\pm$0.03 [20] | 4226.26$\pm 0.04$ [18] | 2013
530.30 |  $\pm 2.39$ [19] | 4235.77$\pm 0.04$ [19] | 2017
55.59 |  $\pm 0.04$ [20] | 4241.66$\pm 0.12$ [18] | 2013
538.10 |  $\pm 2.69$ [19] |  $4243.97\pm 0.04$ [19] | 2017
523.74 |  $\pm$0.10 [20] | 4257.97$\pm 0.04$ [18] | 2013
301.93 |  $\pm$0.08 [20] | 4257.97$\pm 0.04$ [18] | 2013
531.10 |  $\pm 3.13$ [19] | 4266.81$\pm 0.04$ [19] | 2017
175.70 |  $\pm 0.97$ [19] | 4277.78$\pm 0.11$ [19] | 2017
502.4 | 4288.43$\pm 0.06$ [19] | 2019
6.8 to 18.0 [21] | 29 energies | 2014
Table 4: Data with $E_{\rm CM}>4.3$ GeV, where there should be no $X(3872)$ signal. The total luminosity is 5.3 fb-1. Only statistical uncertainties are shown. Entries without a reference are estimated values. The last line of the table shows the range of luminosity values for 28 different energy points that each have a much smaller luminosity than the other data points. Luminosity (pb-1) |  $E_{\rm CM}$ (MeV) | Year
---|---|---
44.9 |  $\pm$0.03 [20] | 4307.89$\pm 0.17$ [18] | 2013
501.2 | 4312.68$\pm 0.06$ [19] | 2019
505.0 | 4337.93$\pm 0.06$ [19] | 2019
539.84 |  $\pm$0.10 [20] | 4358.26$\pm 0.05$ [18] | 2013
522.7 | 4377.88$\pm 0.06$ [19] | 2019
55.18 |  $\pm$0.04 [20] | 4387.40$\pm 0.17$ [18] | 2013
507.8 | 4396.83$\pm 0.06$ [19] | 2019
44.67 |  $\pm$0.03 [20] | 4415.58$\pm 0.04$ [18] | 2014
1028.89 |  $\pm$ 0.13 [20] | 4415.58$\pm 0.04$ [18] | 2014
569.9 | 4437.10$\pm 0.06$ [19] | 2019
109.94 |  $\pm$0.04 [20] | 4467.06$\pm 0.11$ [18] | 2014
109.98 |  $\pm$0.04 [20] | 4527.14$\pm 0.11$ [18] | 2014
47.67 |  $\pm$0.03 [20] | 4574.50$\pm 0.18$ [18] | 2014
566.93 |  $\pm$0.11 [20] | 4599.53$\pm 0.07$ [18] | 2014
6.4 to 9.3 [21] | 28 energies | 2014

The Beijing Spectrometer (BESIII) detector is described in detail in Ref.
[22]. A super-conducting solenoid provides a 1.0 T magnetic field. Inside the
magnet are the multi-layer drift chamber (MDC) for particle tracking, a CsI
(Tl) electromagnetic calorimeter (EMC) to measure the energy of
electromagnetic showers, and a time-of-flight system (TOF) using plastic
scintillators to help with particle identification. Charged particles with a
momentum of 1 GeV/c have a momentum resolution of 0.5%, and the $dE/dx$
resolution is 6% for electrons from Bhabha scattering. The EMC measures
photons with a resolution of 2.5% (5%) at 1 GeV in the barrel (end cap)
region. The time resolution of the TOF barrel region is 68 ps, while the end
cap has a resolution of 110 ps. The end cap TOF was upgraded in 2015 with
multi-gap resistive plate chambers, providing a time resolution of 60 ps [23].

Simulated data samples produced with a geant4-based [24] Monte Carlo (MC)
package, which includes the geometric description of the BESIII detector and
the detector response, are used to determine detection efficiencies and to
estimate backgrounds. The simulation models the beam energy spread and initial
state radiation (ISR) in the $e^{+}e^{-}$ annihilations with the generator
kkmc [25]. The inclusive MC simulation sample includes the production of open
charm processes, the ISR production of vector charmonium(-like) states, and
the continuum processes incorporated in KKMC [25]. The known decay modes are
modeled with evtgen [26] with the branching fractions taken from the Particle
Data Group (PDG) [11], and the remaining unknown charmonium decays are modeled
with lundcharm [27]. Final state radiation (FSR) from charged final state
particles is incorporated using the photos package [28].

Signal MC samples are generated for the search channels to estimate the
reconstruction efficiency. We assume the E1 transition dominates
$e^{+}e^{-}\to\gamma X(3872)$, so the angular distribution is given by
$1-\frac{1}{3}\cos^{2}\theta_{\gamma}$ [29], where $\theta_{\gamma}$ is the
helicity angle of the photon. The decays of
$X(3872)\to\pi^{+}\pi^{-}\chi_{c0}$ and $X(3872)\to\pi^{0}\pi^{0}\chi_{c0}$
are generated with a uniform distribution in phase space. The
$X(3872)\to\pi^{0}\chi_{c0}$ decay is generated as a P-wave transition.
Uniform distributions in phase space are also used to model both decays of
$\chi_{c0}\to\pi^{+}\pi^{-}$ and $\chi_{c0}\to K^{+}K^{-}$. For the four-body
$\chi_{c0}$ decays, we include a uniform phase space component as well as the
most common intermediate states through which the $\chi_{c0}$ can decay, with
sizes scaled according to the branching fractions measured in the PDG. This
means for $\chi_{c0}\to\pi^{+}\pi^{-}\pi^{+}\pi^{-}$ decays we include both
$\chi_{c0}\to\pi^{+}\pi^{-}\pi^{+}\pi^{-}$ and
$\chi_{c0}\to\rho^{0}\pi^{+}\pi^{-}$. The intermediate state included for
$\chi_{c0}\to\pi^{+}\pi^{-}K^{+}K^{-}$ is $\chi_{c0}\to
K_{1}^{\pm}(1270)K^{\mp}$, and for $\chi_{c0}\to\pi^{+}\pi^{-}\pi^{0}\pi^{0}$
it is $\chi_{c0}\to\rho^{\pm}\pi^{\mp}\pi^{0}$. To optimize the selection
criteria, we normalize the size of the $X(3872)\to\pi\pi\chi_{c0}$ MC sample
by setting
$\mathcal{B}(X(3872)\to\pi\pi\chi_{c0})=\mathcal{B}(X(3872)\to\pi^{+}\pi^{-}J/\psi)$,
while for $X(3872)\to\pi^{0}\chi_{c0}$, we scale the signal MC sample to match
the branching fraction ratio predicted by Ref. [12] for a four-quark state,
$\frac{\mathcal{B}(X(3872)\rightarrow\pi^{0}\chi_{c0})}{\mathcal{B}(X(3872)\to\pi^{0}\chi_{c1})}$$=2.97$.

We also generate background MC samples for the processes
$e^{+}e^{-}\to\omega\chi_{c0}$ [30], $e^{+}e^{-}\to\pi\pi\psi(2S)$ [31],
$e^{+}e^{-}\to\pi^{0}\psi(2S)$, $e^{+}e^{-}\to\eta\psi(2S)$, and
$e^{+}e^{-}\to\gamma_{ISR}\psi(2S)$, since all of these backgrounds could peak
at the $\chi_{c0}$ when $\psi(2S)\to\gamma\chi_{c0}$. We generate
$e^{+}e^{-}\to\pi^{0}\psi(2S)$ and $e^{+}e^{-}\to\eta\psi(2S)$ samples of the
same size as $e^{+}e^{-}\to\pi^{+}\pi^{-}\psi(2S)$ to get a conservative
estimation of the size of these backgrounds. The $\chi_{c0}$ decays are the
same as in the signal MC samples, and the other particles are allowed to decay
inclusively.

To ensure that any potential signal is not due to cross feed from other
$X(3872)$ decay channels, we also generate background MC samples for all known
$X(3872)$ decays. We simulate $e^{+}e^{-}\to\gamma X(3872)$, with the
$X(3872)$ state decaying to $\gamma J/\psi$, $\pi^{0}\chi_{c1}$,
$\rho^{0}J/\psi$, $\omega J/\psi$, and $D^{0}\overline{D}^{*0}+c.c.$, with all
particles decaying inclusively. We also simulate $X(3872)\to\gamma\psi(2S)$ at
the rate measured in Ref. [8], even though BESIII sets a more stringent upper
limit on this decay mode.

The inclusive MC sample is used to check for other possible backgrounds, and
no peaking backgrounds are found in the $X(3872)$ signal region.

##  III Event Selection

All the final state particles are required to be reconstructed in the
detector. Charged tracks detected in the MDC are required to be within a polar
angle ($\theta$) range of $|\rm{cos\theta}|<0.93$, where $\theta$ is defined
with respect to the $z$-axis. The distance of closest approach to the
interaction point (IP) must be less than 10 cm along the $z$-axis, $|V_{z}|$,
and less than 1 cm in the transverse plane, $|V_{xy}|$. No particle
identification is used.

Photon candidates are identified using showers in the EMC. The deposited
energy of each shower must be more than 25 MeV in the barrel region
($|\cos\theta|<0.80$) and more than 50 MeV in the end cap region
($0.86<|\cos\theta|<0.92$). To exclude showers that originate from charged
tracks, the angle between the position of each shower in the EMC and the
closest extrapolated charged track must be greater than 10 degrees. To
suppress electronic noise and showers unrelated to the event, the difference
between the EMC time and the event start time is required to be within (0,
700) ns.

We perform a (4+n)C kinematic fit, where 4 constraints conserve the total
four-momentum and the rest constrain the masses of the n $\pi^{0}$s in the
final state. The $\chi^{2}/DOF$, where $DOF$ is the number of degrees of
freedom in the kinematic fit, is optimized for each final state, and is the
selection criteria that removes the most background. We reconstruct
$\chi_{c0}$ candidates between 3.2 and 3.7 GeV so we can use $\chi_{c0}$
sidebands to estimate the shape and size of the non-peaking backgrounds. We
fit $X(3872)$ candidates between 3.75 and 4.0 GeV/$c^{2}$ after selecting
$\chi_{c0}$ candidates in a 50 MeV/$c^{2}$ window centered on the PDG mass of
the $\chi_{c0}$ resonance [11]. We only use data in the region $4.15<E_{\rm
CM}<4.3$ GeV for our fits, since that is where the production cross section is
the largest [10].

For the normalization channel, we use the same event selection as Ref. [9]. We
separate electrons and muons based on the energy deposited in the EMC divided
by the track momentum ($E/p$). Electrons are required to have $E/p>0.85$,
while muons must have $E/p<0.25$. Additionally, we require the kinematic fit
has $\chi^{2}/DOF<10$. The $J/\psi$ candidate is selected by using a 20
MeV$/c^{2}$ window in the $\ell^{+}\ell^{-}$ invariant mass distribution
centered on the PDG mass of the $J/\psi$ resonance. To suppress radiative
Bhabha events, the opening angle between the pions is required to satisfy
$\cos\theta_{\pi\pi}<0.98$. This background is further suppressed by requiring
the angle between any charged track and photon satisfies $\cos\theta_{\gamma
tk}<0.98$. There are additional backgrounds from $e^{+}e^{-}\to\eta J/\psi$
and $e^{+}e^{-}\to\eta^{\prime}J/\psi$ for the normalization channel, and they
are suppressed by requiring $M(\gamma\pi^{+}\pi^{-})>0.6$ GeV/$c^{2}$ and
$|M(\gamma\pi^{+}\pi^{-})-M(\eta^{\prime})|>0.02$ GeV/$c^{2}$, where
$M(\eta^{\prime})$ is the PDG mass of the $\eta^{\prime}$ resonance.

In order to refine the selection criteria for the search channels, we optimize
the figure of merit $\textrm{FOM}=\frac{S}{\sqrt{S+B_{1}+B_{2}}}$. The signal
$S$ is from a 50 MeV/$c^{2}$ box centered on the PDG masses of the $X(3872)$
and $\chi_{c0}$ states in the signal MC samples, scaled according to Ref.
[12]. The background $B_{1}$ is the estimated number of background events in
the signal region of the $X(3872)$ using $\chi_{c0}$ sidebands that are 50
MeV/$c^{2}$ wide in the $\chi_{c0}$ candidate masses on both the lower and
higher mass sides of the $\chi_{c0}$ signal region in data and extend from
3.75 to 4.0 GeV$/c^{2}$ in the $X(3872)$ candidate masses. The background
$B_{2}$ is the estimated number of peaking $\chi_{c0}$ background events,
which is determined using the peaking $\chi_{c0}$ background MC samples, which
are scaled according to previous measurements at BESIII [30, 31]. A plot
showing the signal region and sideband regions for
$X(3872)\to\pi^{0}\chi_{c0}$ decays with $\chi_{c0}\to\pi^{+}\pi^{-}$ is shown
in Figure 1. The data plots in the left column show that the $\chi_{c0}$
sidebands match the background in the signal region well, and the signal MC
plots in the right column show the selection window includes the majority of
our signal. Note that the upper right plot shows that part of the signal MC
falls in the sideband region. This $\chi_{c0}$ selection window was optimized
by maximizing the FOM, and includes approximately 87% to 90% of events for
most final states. Since there are two $\chi_{c0}$ sidebands that are each
five times as wide as the signal region in the $X(3872)$ invariant mass
distribution, this means we scale the sidebands down by a factor of ten for
the FOM calculation. These wide sidebands in $X(3872)$ candidate masses are
used to increase the data sample size for the FOM calculation. For final
states that have multiple possible $\chi_{c0}$ combinations, for instance,
$X(3872)\to\pi^{0}\chi_{c0}$ with $\chi_{c0}\to\pi^{+}\pi^{-}\pi^{0}\pi^{0}$,
an event is in the signal region if at least one combination of the
$\pi^{+}\pi^{-}\pi^{0}\pi^{0}$ invariant mass is in the signal region, and a
sideband event if at least one combination is in the sidebands and no
combinations fall in the signal region.

Figure 1: Plots showing the signal and sideband regions in data (left) and
signal MC (right) before selection criteria optimization. The top row shows
the $\chi_{c0}$ candidate mass projections, where the vertical red lines
denote the signal region and the two sideband regions are the areas between
the blue lines. The bottom row shows two-dimensional plots of the $\chi_{c0}$
candidate masses versus the $X(3872)$ candidate masses, where the z-axis is
the number of events in each bin. The central box is the signal region, and
the long rectangles above and below the $\chi_{c0}$ mass are the sidebands
used to estimate non-peaking $\chi_{c0}$ backgrounds.

In order to use a variable to suppress background, we require the FOM
increases by at least 5%. This is done to simplify the selection criteria as
well as to reduce the systematic uncertainties. We optimize the $\chi^{2}/DOF$
of the kinematic fit for each $\chi_{c0}$ decay mode, which is the requirement
that removes most of the background. In some cases, the FOM forms a broad
plateau for different $\chi^{2}/DOF$ values, so we choose the loosest
requirement that gives an FOM value within 1% of the maximum. By doing this,
we have essentially the same FOM value, but a much larger reconstruction
efficiency. The $\chi^{2}/DOF$ requirements for each final state range from
2.25 to 8.0, and the exact values are shown in Table 5. We also optimize $E/p$
to separate pions from electrons and positrons. This criteria is used for
three final states, $X(3872)\to\pi^{0}\chi_{c0}$ with
$\chi_{c0}\to\pi^{+}\pi^{-}$ requires $E/p<0.95$,
$X(3872)\to\pi^{+}\pi^{-}\chi_{c0}$ with $\chi_{c0}\to\pi^{+}\pi^{-}$ requires
$E/p<0.88$, and $X(3872)\to\pi^{+}\pi^{-}\chi_{c0}$ with $\chi_{c0}\to
K^{+}K^{-}$ requires $E/p<0.96$. This requirement does not improve the FOM
more than 5% for the other final states.

Table 5: Optimized $\chi^{2}/DOF$ values for all three $X(3872)$ decays. Decay | $X\to\pi^{0}\chi_{c0}$ | $X\to\pi^{+}\pi^{-}\chi_{c0}$ | $X\to\pi^{0}\pi^{0}\chi_{c0}$
---|---|---|---
$\chi_{c0}\to\pi^{+}\pi^{-}$ | 4.0 | 3.5 | 5.25
$\chi_{c0}\to K^{+}K^{-}$ | 3.75 | 4.25 | 2.5
$\chi_{c0}\to\pi^{+}\pi^{-}\pi^{+}\pi^{-}$ | 3.0 | 8.0 | 3.75
$\chi_{c0}\to\pi^{+}\pi^{-}K^{+}K^{-}$ | 3.75 | 6.25 | 3.75
$\chi_{c0}\to\pi^{+}\pi^{-}\pi^{0}\pi^{0}$ | 2.25 | 3.0 | 2.5

Several additional variables are investigated to reduce the background, but
none of them increase the FOM enough to be included. These include a
$\psi(2S)$ veto on the $\gamma\chi_{c0}$ invariant mass, an $\omega$ veto on
the $\gamma\pi^{0}$ system, and a veto on $D\to\pi K$. The background
$e^{+}e^{-}\to\gamma D^{*}\bar{D}$ with $D^{*}\to\gamma D$ with
$D\to\pi^{+}\pi^{0}K^{-}$ and $\bar{D}\to K^{+}\pi^{-}$ gives the same final
state as the search channel $X(3872)\to\pi^{0}\chi_{c0}$ with
$\chi_{c0}\to\pi^{+}\pi^{-}K^{+}K^{-}$, but inclusive MC samples show this
background does not cause any peaking backgrounds for $4.15<E_{\rm CM}<4.3$
GeV.

There is a small amount of peaking background at the $X(3872)$ mass for
several decay modes.Only two search final states have at least one predicted
$X(3872)$ background event. The $X(3872)\to\pi^{0}\pi^{0}\chi_{c0}$ channel
with $\chi_{c0}\to\pi^{+}\pi^{-}K^{+}K^{-}$ has a background from $X(3872)\to
D^{*0}\bar{D}^{0}$ with $D^{*0}\to\pi^{0}D^{0}$ and $D^{0}\to
K^{-}\pi^{+}\pi^{0}$ and $\bar{D}^{0}\to K^{+}\pi^{-}$, where the charge
conjugated mode is also implied. There will be $1.2\pm 0.1$ events due to
$X(3872)\to D^{0}\bar{D}^{*0}+c.c.$ for this search channel. The
$X(3872)\to\pi^{+}\pi^{-}\chi_{c0}$ mode with
$\chi_{c0}\to\pi^{+}\pi^{-}\pi^{0}\pi^{0}$ has a predicted rate of $1.0\pm
0.1$ background events from $X(3872)\to\omega J/\psi$, with
$\omega\to\pi^{+}\pi^{-}\pi^{0}$ and $J/\psi\to\pi^{+}\pi^{-}\pi^{0}$ decays.
Since the background levels are low, we do not veto them since they give a
more conservative upper limit and it keeps the selection criteria simple.

After optimizing the selection criteria, there could still be multiple
$X(3872)$ candidates due to the different photon combinations used to
reconstruct the $\pi^{0}$ candidates. To eliminate this double counting, we
rank the remaining combinations by their $\chi^{2}/DOF$ of the kinematic fit
and choose the best combination. This leaves at most one $X(3872)$ combination
per event, so there is no double counting when we fit the $X(3872)$ mass
spectrum. We measure the average number of combinations per event in the full
range of $X(3872)$ and $\chi_{c0}$ candidate masses before selecting the best
combination, and find that for all $\chi_{c0}$ decays except
$\chi_{c0}\to\pi^{+}\pi^{-}K^{+}K^{-}$ it ranges from 1.0 to 1.2 in both data
and signal MC. For the $\chi_{c0}\to\pi^{+}\pi^{-}K^{+}K^{-}$ mode, there can
be multiple combinations from swapping a pion and kaon during the
reconstruction. In this case, the average number of combinations before
selecting the best combination varies from 1.5 to 1.8 in data and signal MC
simulation.

##  IV Fitting

The production mechanism for the search and normalization channels are the
same, so we are able to combine data from different energies in our
measurements. Unbinned extended maximum likelihood fits are performed to the
$X(3872)$ candidate mass spectrum in the region 3.75 to 4.0 GeV/$c^{2}$ for
both the signal and normalization channels.

The $X(3872)$ signal shape for the normalization channel is well described by
a Voigtian function, the convolution of a Breit-Wigner and a Gaussian
function. The internal width of the Breit-Wigner function is set to be 0.96
MeV, the same value used for the signal MC simulation and measured in Ref [4].
We fit the signal MC sample to fix the mass and resolution parameters of the
Gaussian function. This fit is also used to determine the reconstruction
efficiency, which is calculated by integrating the fit function and dividing
by the number of events that were generated. To combine the efficiencies at
different energy points, we perform a weighted average, where the weights are
the luminosity times cross section for each $E_{CM}$. The default cross
section used is the $\sigma(e^{+}e^{-}\to\gamma X(3872))$ measured in Ref.
[10]. The total fit function for the normalization channel is the sum of the
Voigtian signal function and a first order polynomial function to describe the
background.

Each $X(3872)$ search channel includes five $\chi_{c0}$ decay modes, so we
perform a simultaneous fit to all five final states. The Voigtian signal
shapes and the reconstruction efficiencies are determined using the same
method as the normalization channel. We scale the relative sizes of the signal
yield for each $\chi_{c0}$ decay mode. To do so, we use scaling factors
$w_{i}=\epsilon_{i}\mathcal{B}_{i}(\chi_{c0})$, where $\epsilon_{i}$ is the
reconstruction efficiency and $\mathcal{B}_{i}(\chi_{c0})$ is the PDG
branching fraction for a specific $\chi_{c0}$ decay mode i. Using these
scales, we define

$N_{\textrm{tot},i}\equiv\frac{N_{i}}{w_{i}},$

where $N_{i}$ is the signal yield for a specific $\chi_{c0}$ decay mode i. We
then constrain the $N_{\textrm{tot},i}$ to be the same for all five
$\chi_{c0}$ decay modes, so we get a single yield value $N_{\textrm{tot}}$
from the fit.

The $X(3872)$ search channels all have a background process with the same
final state, including a $\chi_{c0}$. For $e^{+}e^{-}\to\gamma X(3872)$ with
$X(3872)\to\pi^{0}\chi_{c0}$, this is $e^{+}e^{-}\to\omega\chi_{c0}$ with
$\omega\to\gamma\pi^{0}$, and for $e^{+}e^{-}\to\gamma X(3872)$ with
$X(3872)\to\pi\pi\chi_{c0}$, it is $e^{+}e^{-}\to\pi\pi\psi(2S)$, with
$\psi(2S)\to\gamma\chi_{c0}$. All of these processes produce an asymmetric
background shape in the $X(3872)$ candidate mass spectrum. To account for this
background, we include a histogram of these background processes in the fit.
The size and shape of this histogram are fixed based on MC samples that were
generated using previously measured cross sections at BESIII. The
non-$\chi_{c0}$ peaking backgrounds are described using a first order
polynomial function for $X(3872)\to\pi^{0}\chi_{c0}$ and second order
polynomial functions for $X(3872)\to\pi\pi\chi_{c0}$. In order to more
accurately determine the shape of the non-$\chi_{c0}$ peaking backgrounds, we
constrain the shape of the background polynomial functions to be the same in
the $\chi_{c0}$ signal and sideband regions in data. To ensure no signal
contaminates the sideband fit, we exclude the $X(3872)$ signal region in the
fit to the sideband regions. We include a floating scale factor between the
polynomial function in the signal and sideband regions. The total fit
functions for the search channels are the sum of the polynomial function, the
peaking $\chi_{c0}$ background MC histogram, and the Voigtian signal function.

To measure the ratio of branching fractions, we use the formula

$\frac{\mathcal{B}(X(3872)\to\pi^{0}\chi_{c0})}{\mathcal{B}(X(3872)\to\pi^{+}\pi^{-}J/\psi)}=\frac{N_{\textrm{tot}}}{\mathcal{B}(\pi^{0})}\frac{\epsilon_{\pi\pi
J/\psi}\mathcal{B}(J/\psi)}{N_{\pi\pi J/\psi}},$

where $N_{\pi\pi J/\psi}$ and $\epsilon_{\pi\pi J/\psi}$ are the number of
events and reconstruction efficiency for $X(3872)\to\pi^{+}\pi^{-}J/\psi$,
respectively. The branching fraction for $\pi^{0}\to\gamma\gamma$ is denoted
by $\mathcal{B}(\pi^{0})$, and the branching fraction for
$J/\psi\to\ell^{+}\ell^{-}$ ($\ell=e$ or $\mu$) is $\mathcal{B}(J/\psi)$,
which are taken from the PDG. Note that the production cross section, ISR
correction factors, and the integrated luminosity are canceled in the ratio,
and the ratio of branching fractions is only sensitive to the ratio of
efficiencies.

To calculate the significance of the signal, we perform a fit with the signal
yield floating, as well as a fit with the signal yield fixed to zero. We then
use the likelihood ratio test to determine the statistical significance of the
fit. The systematic uncertainty due to the fitting model is determined using
216 alternative fit models, described in detail in Section V.

###  IV.1 Normalization Channel

The results of the fit to the normalization channel are shown in Figure 2.
There is a clear signal for the $X(3872)$ state. The reconstruction efficiency
is 32.5$\pm$0.2%, and the fit yield is $88.7^{+10.4}_{-9.7}$. These are both
consistent with Ref. [9]’s measured efficiency of $32.3\%$ and yield of
$84.1^{+10.1}_{-9.4}$, where the uncertainties on the yields are just the
statistical uncertainties from the fits. Note this analysis includes two more
data points not included in Ref. [9], which is why the yield is larger here.

Figure 2: Fit result to the $\pi^{+}\pi^{-}J/\psi$ mass spectrum with a first
order polynomial function to describe the background (dashed line) and a
Voigtian as a signal function (solid line). There is a clear signal for the
$X(3872)$ state.

###  IV.2 $X(3872)\to\pi^{0}\chi_{c0}$

The fit result to the $\pi^{0}\chi_{c0}$ mass spectrum is shown in Figure 3.
There is no obvious signal for any of the reconstructed $\chi_{c0}$ final
states. The efficiencies, scaling factors, and yields for each $\chi_{c0}$
decay mode are shown in Table 6. The signal has a total statistical
significance of $2.4\sigma$. In Section V, we perform 216 alternative fits
with different signal and background models. Figure 4 shows that the
significance for all these fit variations is always at least 1.3$\sigma$.
Figure 4 also shows the range of upper limit values measured in all the fit
variations.

Figure 3: Simultaneous fit to the $\pi^{0}\chi_{c0}$ mass spectrum for
$4.15<E_{\rm CM}<4.3$ GeV. The points are data in the signal region, the light
gray histogram is the background estimate from $\chi_{c0}$ sidebands in data,
and the dark gray histogram is the peaking $\chi_{c0}$ background estimated
from MC simulated with $e^{+}e^{-}\to\omega\chi_{c0}$. The total background
shape is the sum of the $e^{+}e^{-}\to\omega\chi_{c0}$ MC shape plus a first
order polynomial function. The solid lines show the fit with a signal
component, while the dashed lines are the background contributions. The
combination of $e^{+}e^{-}\to\omega\chi_{c0}$ contribution and $\chi_{c0}$
sidebands describes the size of the backgrounds well.

Figure 4: (Left) The distribution of upper limits on $\frac{\mathcal{B}(X(3872)\rightarrow\pi^{0}\chi_{c0})}{\mathcal{B}(X(3872)\to\pi^{+}\pi^{-}J/\psi)}$ from 216 alternative fits. (Right) Distribution of the significance (y-axis) vs the order of the background polynomial function used (x-axis) for the fit to the $\pi^{0}\chi_{c0}$ mass spectrum. Here we find the first order polynomial function gives a larger significance than both the second and third order polynomial functions. Table 6: Efficiencies, scaling factors, and yields for each decay mode of the $\chi_{c0}$ for $X(3872)\to\pi^{0}\chi_{c0}$. The uncertainties on the yields are statistical only. The scaling factors show the relative contribution of each $X(3872)$ decay in the simultaneous fit.  Decay | Efficiency | Scale ($w_{i}$) | Signal Yield
---|---|---|---
$\chi_{c0}\to\pi^{+}\pi^{-}$ | 23.4% | 0.00133 | $5.1\pm 2.4$
$\chi_{c0}\to K^{+}K^{-}$ | 21.6% | 0.00131 | $5.0\pm 2.3$
$\chi_{c0}\to\pi^{+}\pi^{-}\pi^{+}\pi^{-}$ | 13.5% | 0.00315 | $12.1\pm 5.6$
$\chi_{c0}\to\pi^{+}\pi^{-}K^{+}K^{-}$ | 12.9% | 0.00233 | $9.0\pm 4.1$
$\chi_{c0}\to\pi^{+}\pi^{-}\pi^{0}\pi^{0}$ | 5.83% | 0.00188 | $7.2\pm 3.3$

###  IV.3 $X(3872)\to\pi^{+}\pi^{-}\chi_{c0}$

The fit result for the $\pi^{+}\pi^{-}\chi_{c0}$ mass spectrum is shown in
Figure 5. There is no evidence of a signal for this decay mode. The
efficiencies, scaling factors, and yields are shown in Table 7.

Figure 5: Fit results to the $\pi^{+}\pi^{-}\chi_{c0}$ mass spectrum. The points are data in the signal region, the light gray histogram is $\chi_{c0}$ sidebands in data, and the dark gray histogram is the peaking $\chi_{c0}$ background estimated from MC simulated with $e^{+}e^{-}\to\pi^{+}\pi^{-}\psi(2S)$ with $\psi(2S)\to\gamma\chi_{c0}$. The solid lines show the fit with a signal component, while the dashed lines are the background functions. There is no evidence of a signal for the $X(3872)$ state. Table 7:  Efficiencies, scaling factors, and signal yields for each decay mode of the $\chi_{c0}$ for $X(3872)\to\pi^{+}\pi^{-}\chi_{c0}$. The scaling factors show the relative contribution of each component of the simultaneous fit. The uncertainties on the yields are statistical only.  Decay | Efficiency | Scale ($w_{i}$) | Signal Yield
---|---|---|---
$\chi_{c0}\to\pi^{+}\pi^{-}$ | 22.8% | 0.00130 | $0.18\pm 0.71$
$\chi_{c0}\to K^{+}K^{-}$ | 24.5% | 0.00148 | $0.20\pm 0.81$
$\chi_{c0}\to\pi^{+}\pi^{-}\pi^{+}\pi^{-}$ | 21.0% | 0.00492 | $0.7\pm 2.7$
$\chi_{c0}\to\pi^{+}\pi^{-}K^{+}K^{-}$ | 17.1% | 0.00309 | $0.4\pm 1.7$
$\chi_{c0}\to\pi^{+}\pi^{-}\pi^{0}\pi^{0}$ | 8.82% | 0.00284 | $0.4\pm 1.6$

###  IV.4 $X(3872)\to\pi^{0}\pi^{0}\chi_{c0}$

The fit to the $\pi^{0}\pi^{0}\chi_{c0}$ mass spectrum is shown in Figure 6.
There is no evidence for a signal. The efficiencies, scaling factors, and
signal yields for $X(3872)\to\pi^{0}\pi^{0}\chi_{c0}$ are summarized in Table
8. Note that $\chi_{c0}\to K^{+}K^{-}$ has a much smaller efficiency than
$\chi_{c0}\to\pi^{+}\pi^{-}$ in this case because of the more stringent
$\chi^{2}/DOF$ requirement shown in Table 5.

Table 8:  Efficiencies, scaling factors, and signal yields for each decay mode of the $\chi_{c0}$ for $X(3872)\to\pi^{0}\pi^{0}\chi_{c0}$. The uncertainties on the yields are statistical only. The scaling factors show the relative contribution of each component of the simultaneous fit. Decay | Efficiency | Scale ($w_{i}$) | Signal Yield
---|---|---|---
$\chi_{c0}\to\pi^{+}\pi^{-}$ | 11.2% | 0.000637 | $-0.8\pm 1.3$
$\chi_{c0}\to K^{+}K^{-}$ | 6.48% | 0.000392 | $-0.49\pm 0.80$
$\chi_{c0}\to\pi^{+}\pi^{-}\pi^{+}\pi^{-}$ | 5.66% | 0.00132 | $-1.6\pm 2.7$
$\chi_{c0}\to\pi^{+}\pi^{-}K^{+}K^{-}$ | 4.78% | 0.000865 | $-1.1\pm 1.8$
$\chi_{c0}\to\pi^{+}\pi^{-}\pi^{0}\pi^{0}$ | 2.31% | 0.000744 | $-0.9\pm 1.5$

Figure 6: Fit results to the $\pi^{0}\pi^{0}\chi_{c0}$ mass spectrum. The
points are data in the signal region, the light gray histogram is $\chi_{c0}$
sidebands in data, and the dark gray histogram is the peaking $\chi_{c0}$
background estimated from MC simulated with
$e^{+}e^{-}\to\pi^{0}\pi^{0}\psi(2S)$, with $\psi(2S)\to\gamma\chi_{c0}$. The
black histogram is from $X(3872)$ background MC simulation. The solid lines
show the fit with a signal component, while the dashed lines are the
background functions. There is no evidence of a signal for the $X(3872)$
state.

##  V Systematic Uncertainties

We include a 1% systematic uncertainty for each photon and charged track that
does not cancel in the ratio of branching fractions. The remaining systematic
uncertainties are discussed below.

###  V.1 Kinematic Fit

We use the $\chi^{2}/DOF$ of the kinematic fit to reduce the background for
all of the search channels. Since there are no high statistics final states
with similar kinematics, we use the procedure in Ref. [32] to determine this
systematic uncertainty. Ref [32] found that MC simulation has a significantly
narrower $\chi^{2}$ distribution than in data. Corrections to the track helix
parameters of charged particles are used to improve the agreement of the
$\chi^{2}$ distributions between data and MC simulation. We measure the ratio
of the search and normalization channel efficiencies before and after these
corrections are applied. By taking the ratio, we properly account for the fact
that some of the kinematic fit uncertainty cancels in the ratio. We take half
the percent difference as the systematic uncertainty due to the kinematic fit.
This results in a systematic uncertainty between 0.3% and 2.6% for all the
final states.

###  V.2 Signal MC Models

The nominal signal MC simulation includes the decays of the $\chi_{c0}$ into
final states both by uniform distributions in phase space and through the most
common intermediate states. To test how sensitive the kinematics are to the
intermediate states, we generate new signal MC samples in which the
$\chi_{c0}$ only decays through the most common intermediate states. The
angular distribution of the photon in $e^{+}e^{-}\to\gamma X(3872)$ is
modified from the nominal E1 transition to a uniform distribution in phase
space. We also modify the angular distribution of the decay
$X(3872)\to\pi^{0}\chi_{c0}$ from a P-wave transition to a uniform
distribution in phase space. There is no variation for the decays
$X(3872)\to\pi\pi\chi_{c0}$, since they are generated with a uniform
distribution in phase space in the nominal signal MC simulation, which results
in $X(3872)\to\pi^{0}\chi_{c0}$ having a larger systematic uncertainty.

These modifications are done for both the search and normalization channels,
so we compare the ratio of efficiencies for the nominal MC and for these
variations, and take the percent difference as a systematic uncertainty.
Normalizing the results to $X(3872)\to\pi^{+}\pi^{-}J/\psi$, the systematic
uncertainties range from 8.0% to 11.0% for $X(3872)\to\pi^{0}\chi_{c0}$ and
from 0.6% to 3.7% for $X(3872)\to\pi\pi\chi_{c0}$. When normalizing
$X(3872)\to\pi^{0}\chi_{c0}$ to $X(3872)\to\pi^{0}\chi_{c1}$, the uncertainty
in the simulated model of the denominator must also be taken into account. To
account for this, we include the uncertainty of 8.1% from Ref. [9], which is
added in quadrature with the total from $X(3872)\to\pi^{0}\chi_{c0}$.

###  V.3 $\chi_{c0}$ Mass Window

We select a 50 MeV/$c^{2}$ window centered on the $\chi_{c0}$ mass when we
perform our fits. To test the systematic uncertainty related to this
selection, we fit the $\chi_{c0}$ distribution in signal MC simulation using a
Voigtian function. For the initial fit, we fix the intrinsic width of the
Voigtian to be the intrinsic width of the $\chi_{c0}$ resonance from the PDG
[11], which is 10.8 MeV. As a variation, we widen the intrinsic width of the
$\chi_{c0}$ by its PDG uncertainty of 0.6 MeV, and we also widen the width of
the Gaussian function by 20%. A variation of 20% was chosen because a previous
analysis of the process $e^{+}e^{-}\to\gamma\eta_{c}$ at BESIII included a
study comparing hadronic final states in MC simulation and data, and it found
the resolutions can differ by up to 20% [33]. We take the percent difference
in the number of events from the two fits as the systematic uncertainty. These
vary from 1.9% to 6.0%.

###  V.4 $E_{\rm CM}$ Dependence on the Efficiency Ratio

The ratio of branching fractions depends on the ratio of efficiencies
$\frac{\epsilon(X(3872)\to\pi^{0}\chi_{c0})}{\epsilon(X(3872)\to\pi^{+}\pi^{-}J/\psi)}$.
To combine the efficiency ratio measurements at different energies, we perform
a weighted average, where the weights are the luminosity times cross section
for that energy. The default cross section we use is the $e^{+}e^{-}\to\gamma
X(3872)$ cross section measured in Ref. [10]. To probe the systematic
uncertainty due to the $E_{\rm CM}$ dependence of the efficiency ratio, we
also use the cross sections for the process
$e^{+}e^{-}\to\pi^{+}\pi^{-}J/\psi$ measured in Ref. [34], and one based on
the $\psi(4160)$, which is modeled as a Breit-Wigner function with parameters
taken from the PDG [11]. We take the largest deviation from the nominal ratio
of efficiencies as the systematic uncertainty. This systematic uncertainty
varies from 1.6% to 3.1%.

###  V.5 Input Branching Fractions

The branching fractions of the $\chi_{c0}$ decays are used to constrain the
relative sizes of the simultaneous fit components. The branching fractions for
the decays $\chi_{c0}\to\pi^{+}\pi^{-}$, $\chi_{c0}\to K^{+}K^{-}$,
$\chi_{c0}\to\pi^{+}\pi^{-}\pi^{+}\pi^{-}$, and
$\chi_{c0}\to\pi^{+}\pi^{-}K^{+}K^{-}$ are included in a constrained fit done
by the PDG that included 248 results from previous papers, so the correlated
uncertainties between these four branching fractions are known [11]. To
calculate this systematic uncertainty, we generate new sets of input branching
fractions using a multivariate Gaussian function, which uses the known values,
uncertainties, and correlations as input. For each set of branching fractions,
we refit the data with updated fit scales. This procedure is repeated 4000
times for each $X(3872)$ decay, and we measure a new value for the ratio of
branching fractions with respect to the decay $X(3872)\to\pi^{+}\pi^{-}J/\psi$
for each fit. The resulting distribution of the ratio of branching fractions
is then fit with a Gaussian function. The standard deviation of the Gaussian
distribution divided by the nominal ratio is the systematic uncertainty, which
is $4.7\%$ for $X(3872)\to\pi^{0}\chi_{c0}$, $14.9\%$ for
$X(3872)\to\pi^{+}\pi^{-}\chi_{c0}$, and $14.3\%$ for
$X(3872)\to\pi^{0}\pi^{0}\chi_{c0}$ decays.

###  V.6 Fit Model

The last systematic uncertainty is due to the fit model. We perform a total of
216 alternative fits, with the details shown in Table 9. For the upper limits,
we account for the fitting uncertainty by using the fit variation that results
in the largest upper limit. To determine the systematic uncertainty for the
central value regardless of what the central value is, we calculate

$\sigma^{2}=\sum_{i=1}^{216}\frac{(r_{\rm nom}-r_{i})^{2}}{216},$

where $r_{nom}$ is the nominal value for the ratio of branching fractions and
$r_{i}$ is the ratio value for each fit model variation. Using the first order
polynomial function for the background as the nominal for
$X(3872)\to\pi^{0}\chi_{c0}$ decays, we get an absolute fitting uncertainty of
$0.34$. The nominal fitting method for both
$X(3872)\to\pi^{+}\pi^{-}\chi_{c0}$ and $X(3872)\to\pi^{0}\pi^{0}\chi_{c0}$
decays is a second order polynomial function for the background. We measure
absolute fitting uncertainties of $0.12$ for
$X(3872)\to\pi^{+}\pi^{-}\chi_{c0}$ and $0.32$ for
$X(3872)\to\pi^{0}\pi^{0}\chi_{c0}$ modes. To get the total systematic
uncertainty for the central values, these uncertainties are added in
quadrature with the uncertainties from Table 10 and Table 11.

Table 9: List of all the fit variations used in this analysis. Variations to each of these components are done independently, so the total number of variations is $3\cdot 3\cdot 2\cdot 2\cdot 2\cdot 3=216$. Variation | Description | Variations
---|---|---
Polynomial function | 1st, 2nd, 3rd order | 3
Mass | Nominal, and $\pm 1$ MeV | 3
Internal width | 0.96 MeV, 1.39 MeV | 2
Resolution | Nominal, widen 20% | 2
Fit range | Nominal, $3.775-4.0$ GeV | 2
Beam energy | Nominal, $\pm 1$ MeV | 3

###  V.7 Total Systematic Uncertainties

Several of the systematic uncertainties are correlated between the different
$\chi_{c0}$ decay modes. For these uncertainties, we calculate a weighted
average to get the total systematic uncertainty from that source, using the
efficiency times branching fraction as the weights. This is done for the
tracking, photon, kinematic fit, MC simulation model, $\chi_{c0}$ selection,
and $\epsilon$ ratio systematic uncertainties. We summarize the total
systematic uncertainty values for $X(3872)\to\pi^{0}\chi_{c0}$ in Table 10,
and $X(3872)\to\pi^{+}\pi^{-}\chi_{c0}$ and
$X(3872)\to\pi^{0}\pi^{0}\chi_{c0}$ in Table 11. As previously mentioned, the
fitting uncertainty is included for the upper limits by using the variation
that results in the largest upper limit. The total systematic uncertainty for
the central values is obtained by adding the totals in Tables 10 and 11 with
the fitting uncertainties in the previous section.

Table 10: Total systematic uncertainties for the decay $X(3872)\to\pi^{0}\chi_{c0}$. Here we show the systematic uncertainty when the branching fraction is normalized to $X(3872)\to\pi^{+}\pi^{-}J/\psi$ as well as $X(3872)\to\pi^{0}\chi_{c1}$. This is done because some of the photon and charged track systematic uncertainties cancel in the ratio of branching fractions.  Source |  $R_{\pi^{+}\pi^{-}J/\psi}$ Total |  $R_{\pi^{0}\chi_{c1}}$ Total
---|---|---
Tracking | 5.1% | 3.7%
Photon Efficiency | 2.8% | 1.3%
$\chi^{2}/DOF$ Cut | 1.2% | 1.2%
Decay Models | 9.5% | 12.5%
Branching Fractions | 4.7% | 4.7%
$\chi_{c0}$ Selection | 3.3% | 3.3%
$\epsilon$ Ratio | 2.6% | 2.6%
Total | 12.8% | 14.5%
Table 11: Total systematic uncertainties for $X(3872)\to\pi^{+}\pi^{-}\chi_{c0}$ and $X(3872)\to\pi^{0}\pi^{0}\chi_{c0}$ modes. Source | $X(3872)\to\pi^{+}\pi^{-}\chi_{c0}$ | $X(3872)\to\pi^{0}\pi^{0}\chi_{c0}$
---|---|---
Tracking | 3.8% | 5.1%
Photons | 0.8% | 4.8%
$\chi^{2}/DOF$ Cut | 0.8% | 1.1%
Decay Models | 2.5% | 1.9%
$\mathcal{B}(\chi_{c0})$ | $14.9\%$ | $14.3\%$
$\chi_{c0}$ Selection | 2.9% | 3.3%
$\epsilon$ Ratio | 2.3% | 2.0%
Total | 16.0% | 16.5%

##  VI  Calculation of Statistical Uncertainties and Upper Limits

Our upper limits have to include the statistical uncertainty from both the
numerator and denominator of the ratio
$\frac{\mathcal{B}(X(3872)\rightarrow\pi^{0}\chi_{c0})}{\mathcal{B}(X(3872)\to\pi^{+}\pi^{-}J/\psi)}$.
To get the total statistical uncertainty correct, we perform a likelihood scan
for the search and normalization channels. This is done by performing several
hundred fits with the yield fixed to different values and with the background
parameters floating. Once we have both likelihood distributions, we randomly
sample them both to determine new yield values. These yield values are then
used to calculate a new value for the ratio of branching fractions. This
process is done a million times, and the resulting distribution gives us the
likelihood for the ratio of branching fractions. We plot the nominal fit
likelihood distributions for each ratio of branching fractions in Figure 7
(blue histogram). The nominal fit value is the solid blue line, and the dashed
blue lines show the statistical uncertainty, calculated by integrating to $\pm
33\%$ of the distribution. To determine the upper limit, we use the likelihood
function that results in the largest upper limit. We convolve this likelihood
function with a Gaussian function that has a width corresponding to the
systematic uncertainty. The upper limits are calculated by integrating the
resulting curve from 0 up to the point where 90% of the distribution is below
the upper limit. Figure 7 shows the likelihood curves for the largest upper
limit (red curve) after including the systematic uncertainties, as well as the
upper limit values (vertical red line).

The fit variation that results in the largest upper limit for
$X(3872)\to\pi^{0}\chi_{c0}$ is a fit with a first order polynomial function
where the signal function has a mass shifted $+1$ MeV$/c^{2}$ from the nominal
value, an internal width of 1.39 MeV, a resolution increased by 20%, and a
beam energy shifted $-1$ MeV. For $X(3872)\to\pi^{+}\pi^{-}\chi_{c0}$, the
variation that results in the largest upper limit is a fit with a second order
polynomial function where the signal function has a mass shifted $+1$
MeV$/c^{2}$, an internal width of 1.39 MeV, an increased resolution by 20%, a
fit range 50 MeV narrower than the nominal range, and a beam energy shifted
$+1$ MeV. For $X(3872)\to\pi^{0}\pi^{0}\chi_{c0}$, the variation that results
in the largest upper limit is a fit with a third order polynomial function
where the signal function has a mass shifted $-1$ MeV, and an increased
resolution by 20%.

Figure 7: Plots showing the likelihood scans for all four measured ratios of
branching fractions. The top left shows
$\frac{\mathcal{B}(X(3872)\rightarrow\pi^{0}\chi_{c0})}{\mathcal{B}(X(3872)\to\pi^{+}\pi^{-}J/\psi)}$,
top right shows
$\frac{\mathcal{B}(X(3872)\rightarrow\pi^{0}\chi_{c0})}{\mathcal{B}(X(3872)\to\pi^{0}\chi_{c1})}$,
bottom left shows
$\frac{\mathcal{B}(X(3872)\rightarrow\pi^{+}\pi^{-}\chi_{c0})}{\mathcal{B}(X(3872)\to\pi^{+}\pi^{-}J/\psi)}$
and bottom right shows
$\frac{\mathcal{B}(X(3872)\rightarrow\pi^{0}\pi^{0}\chi_{c0})}{\mathcal{B}(X(3872)\to\pi^{+}\pi^{-}J/\psi)}$.
The blue curve is the likelihood for the nominal fit model, while the red
curve is the likelihood curve that gives the largest upper limit after
convolving by the systematic uncertainty. The solid blue line is the nominal
fit value, and the dashed blue lines are the statistical uncertainties. The
vertical red line shows the 90% C.L. upper limit value.

##  VII Summary

In summary, we perform a search for the decays $X(3872)\to\pi^{0}\chi_{c0}$,
$X(3872)\to\pi^{+}\pi^{-}\chi_{c0}$, and $X(3872)\to\pi^{0}\pi^{0}\chi_{c0}$.
The significance for $X(3872)\to\pi^{0}\chi_{c0}$ is found to be between
1.3$\sigma$ and 2.8$\sigma$ . Since the fit results are all below 3$\sigma$,
we set an upper limit of
$\frac{\mathcal{B}(X(3872)\rightarrow\pi^{0}\chi_{c0})}{\mathcal{B}(X(3872)\to\pi^{+}\pi^{-}J/\psi)}$$<$
3.6 at 90% C.L, which is a significant improvement over the previous BESIII
upper limit of
$\frac{\mathcal{B}(X(3872)\rightarrow\pi^{0}\chi_{c0})}{\mathcal{B}(X(3872)\to\pi^{+}\pi^{-}J/\psi)}$$<19$
[9]. Combined with the previous BESIII result [9] of
$\frac{\mathcal{B}(X(3872)\rightarrow\pi^{0}\chi_{c1})}{\mathcal{B}(X(3872)\to\pi^{+}\pi^{-}J/\psi)}$$=0.88^{+0.33}_{-0.27}\pm
0.10$, we set an upper limit of
$\frac{\mathcal{B}(X(3872)\rightarrow\pi^{0}\chi_{c0})}{\mathcal{B}(X(3872)\to\pi^{0}\chi_{c1})}$$<$
4.4 at 90% C.L. This upper limit is too large to rule out any interpretation
of the $X(3872)$ state.

We find no significant signals for $X(3872)\to\pi^{+}\pi^{-}\chi_{c0}$ and
$X(3872)\to\pi^{0}\pi^{0}\chi_{c0}$, so we set upper limits of
$\frac{\mathcal{B}(X(3872)\rightarrow\pi^{+}\pi^{-}\chi_{c0})}{\mathcal{B}(X(3872)\to\pi^{+}\pi^{-}J/\psi)}$$<0.68$,
and
$\frac{\mathcal{B}(X(3872)\rightarrow\pi^{0}\pi^{0}\chi_{c0})}{\mathcal{B}(X(3872)\to\pi^{+}\pi^{-}J/\psi)}$$<1.7$
at 90% C.L. Both of them are consistent with theoretical predictions from Ref.
[12] that they should be suppressed regardless of whether the $X(3872)$ is a
four-quark or charmonium state. All of our results are summarized in Table 12.

Table 12: Central values and 90% C.L. upper limits for each of our measurements. The first uncertainty on the central values is statistical, and the second is systematic. Ratio | Central Value (Upper Limit)
---|---
$\frac{\mathcal{B}(X(3872)\rightarrow\pi^{0}\chi_{c0})}{\mathcal{B}(X(3872)\to\pi^{+}\pi^{-}J/\psi)}$ |  $1.70^{+0.55}_{-0.75}\pm 0.40$$\textrm{ }(3.6)$
$\frac{\mathcal{B}(X(3872)\rightarrow\pi^{0}\chi_{c0})}{\mathcal{B}(X(3872)\to\pi^{0}\chi_{c1})}$ |  $1.9^{+1.2}_{-0.7}\pm 0.5$$\textrm{ }(4.4)$
$\frac{\mathcal{B}(X(3872)\rightarrow\pi^{+}\pi^{-}\chi_{c0})}{\mathcal{B}(X(3872)\to\pi^{+}\pi^{-}J/\psi)}$ |  $0.06^{+0.24}_{-0.23}\pm 0.12$$\textrm{ }(0.68)$
$\frac{\mathcal{B}(X(3872)\rightarrow\pi^{0}\pi^{0}\chi_{c0})}{\mathcal{B}(X(3872)\to\pi^{+}\pi^{-}J/\psi)}$ |  $-0.55^{+0.53}_{-0.42}\pm 0.33$$\textrm{ }(1.7)$

##  VIII Acknowledgements

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 R&D Program of China under Contracts Nos. 2020YFA0406300, 2020YFA0406400;
National Natural Science Foundation of China (NSFC) under Contracts Nos.
11625523, 11635010, 11735014, 11822506, 11835012, 11935015, 11935016,
11935018, 11961141012, 12022510, 12025502, 12035009, 12035013, 12061131003;
the Chinese Academy of Sciences (CAS) Large-Scale Scientific Facility Program;
Joint Large-Scale Scientific Facility Funds of the NSFC and CAS under
Contracts Nos. U1732263, U1832207; CAS Key Research Program of Frontier
Sciences under Contract No. QYZDJ-SSW-SLH040; 100 Talents Program of CAS;
INPAC and Shanghai Key Laboratory for Particle Physics and Cosmology; ERC
under Contract No. 758462; European Union Horizon 2020 research and innovation
programme under Contract No. Marie Sklodowska-Curie grant agreement No 894790;
German Research Foundation DFG under Contracts Nos. 443159800, Collaborative
Research Center CRC 1044, FOR 2359, FOR 2359, GRK 214; Istituto Nazionale di
Fisica Nucleare, Italy; Ministry of Development of Turkey under Contract No.
DPT2006K-120470; National Science and Technology fund; Olle Engkvist
Foundation under Contract No. 200-0605; STFC (United Kingdom); The Knut and
Alice Wallenberg Foundation (Sweden) under Contract No. 2016.0157; The Royal
Society, UK under Contracts Nos. DH140054, DH160214; The Swedish Research
Council; U. S. Department of Energy under Contracts Nos. DE-FG02-05ER41374,
DE-SC-0012069

## References

  * [1] S. K. Choi et al. (Belle Collaboration), Phys. Rev. Lett. 91, 262001 (2003).

  * [2] R. Aaij et al. (LHCb Collaboration), Phys. Rev. Lett. 110, 222001 (2013).
  * [3] S. Godfrey and N. Isgur, Phys Rev. D 32 189 (1985).
  * [4] R. Aaij et al. (The LHCb Collaboration), JHEP 08 123 (2020).
  * [5] R. Aaij et al. (LHCb Collaboration), Phys. Rev. D 102, 092005 (2020).
  * [6] S.-K. Choi et al. (The Belle Collaboration) Phys. Rev. D 84, 052004 (2011).
  * [7] T. Aushev et al. (Belle Collaboration) Phys. Rev. D 81, 031103 (2010); B. Aubert et al. (BABAR Collaboration), Phys. Rev. D 77, 011102 (2008).
  * [8] R. Aaij et al. (LHCb Collaboration), Nucl. Phys. B886, 665 (2014); V. Bhardwaj et al. (Belle Collaboration), Phys. Rev. Lett. 107, 091803 (2011); B. Aubert et al. (BABAR Collaboration) Phys. Rev. Lett. 102, 132001 (2009).
  * [9] M. Ablikim et al. (BESIII Collaboration), Phys. Rev. Lett. 122, 202001 (2019).
  * [10] M. Ablikim et al. (BESIII Collaboration), Phys. Rev. Lett. 122 232002 (2019).
  * [11] P.A. Zyla et al. (Particle Data Group), Prog. Theor. Exp. Phys. 2020, 083C01 (2020).
  * [12] S. Dubynskiy, M. B. Voloshin, Phys. Rev. D. 77, 014013 (2008).
  * [13] S. Fleming, T. Mehen, Phys. Rev. D. 78, 094019 (2008).
  * [14] Q. Wu, D.Y. Chen, and T. Matsuki, Eur. Phys. J. C 81, 193 (2021).
  * [15] Y. Dong, A. Faessler, T. Gutsche, S. Koyalenko, V. Lyubovitskij, Phys Rev D 79 094013 (2009).
  * [16] Z. Zhou, M. Yu, Z. Xiao, Phys. Rev. D. 100, 094025 (2019).
  * [17] V. Bhardwaj et al. (Belle Collaboration), Phys. Rev. D, 99 111101(R) (2019).
  * [18] M. Ablikim et al. (BESIII Collaboration), Chinese Phys. C. 40 063001 (2016).
  * [19] M. Ablikim et al. (BESIII Collaboration), Chinese Phys. C. 45 103001 (2021).
  * [20] M. Ablikim et al. (BESIII Collaboration), Chinese Phys. C. 39 093001 (2015).
  * [21] M. Ablikim et al. (BESIII Collaboration), Chinese Phys. C. 41 063001 (2017).
  * [22] M. Ablikim et al. (BESIII Collaboration), Nucl. Instrum. Methods Phys. Res., Sect. A 614, 345, (2010).
  * [23] X. Li et al., Proceedings of IPAC2016, Busan, Korea, 2016.
  * [24] S. Agostinelli et al. (GEANT4 Collaboration), Nucl. Instrum. Methods Phys. Res., Sect. A 506, 250 (2003).
  * [25] S. Jadach, B. F. L. ward, and Z. Was, Phys. Rev. D 63, 113009 (2001).
  * [26] D. J. Lange, Nucl. Instrum. Methods Phys. Res., Sect. A 462, 152 (2001).
  * [27] J. C. Chen, G. S. Huang, X. R. Qi, D. H. Zhang and Y. S. Zhu, Phys. Rev. D 62, 034003 (2000); R. L. Yang, R. G. Ping and H. Chen, Chin. Phys. Lett. 31, 061301 (2014).
  * [28] P. Golonka and Z. Was, Eur. Phys. J. C 45, 97 (2006).
  * [29] E. Eichten, K. Gottfried, T. Kinoshita, J. B. Kogut, K. D. Lane and T. M. Yan, Phys. Rev. Lett. 34, 369-372 (1975) [erratum: Phys. Rev. Lett. 36, 1276 (1976)]
  * [30] M. Ablikim et al. (BESIII Collaboration), Phys. Rev. D. 99, 091103(R) (2019).
  * [31] M. Ablikim et al. (BESIII Collaboration), Phys. Rev. D. 96 032004 (2017).
  * [32] M. Ablikim et al. (BESIII Collaboration), Phys. Rev. D. 87, 012002 (2013).
  * [33] M. Ablikim et al. (BESIII Collaboration), Phys. Rev. D. 96, 051101(R) (2017).
  * [34] M. Ablikim et al. (BESIII Collaboration), Phys. Rev. Lett. 118, 092001 (2017).