# Study of the $h_{c}$($1^{1}P_{1}$) meson via
$\psi(2S)\rightarrow\pi^{0}h_{c}$ decays at BESIII

M. Ablikim1, M. N. Achasov10,b, P. Adlarson68, S. Ahmed14, M. Albrecht4, R.
Aliberti28, A. Amoroso67A,67C, M. R. An32, Q. An64,50, X. H. Bai58, Y. Bai49,
O. Bakina29, R. Baldini Ferroli23A, I. Balossino24A, Y. Ban39,h, V.
Batozskaya1,37, D. Becker28, K. Begzsuren26, N. Berger28, M. Bertani23A, D.
Bettoni24A, F. Bianchi67A,67C, J. Bloms61, A. Bortone67A,67C, I. Boyko29, R.
A. Briere5, A. Brueggemann61, H. Cai69, 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, G. Chen1, H. S. Chen1,55, M. L. Chen1,50, S. J.
Chen35, T. Chen1, X. R. Chen25, X. T. Chen1, Y. B. Chen1,50, Z. J. Chen20,i,
W. S. Cheng67C, G. Cibinetto24A, F. Cossio67C, J. J. Cui42, H. L. Dai1,50, J.
P. Dai71, A. Dbeyssi14, R. E. de Boer4, D. Dedovich29, Z. Y. Deng1, A.
Denig28, I. Denysenko29, M. Destefanis67A,67C, F. De Mori67A,67C, Y. Ding33,
J. Dong1,50, L. Y. Dong1,55, M. Y. Dong1,50,55, X. Dong69, S. X. Du73, P.
Egorov29,a, Y. L. Fan69, J. Fang1,50, S. S. Fang1,55, Y. Fang1, R.
Farinelli24A, L. Fava67B,67C, F. Feldbauer4, G. Felici23A, C. Q. Feng64,50, J.
H. Feng51, K Fischer62, M. Fritsch4, C. D. Fu1, Y. N. Gao39,h, Yang Gao64,50,
I. Garzia24A,24B, P. T. Ge69, C. Geng51, E. M. Gersabeck59, A Gilman62, K.
Goetzen11, L. Gong33, W. X. Gong1,50, W. Gradl28, M. Greco67A,67C, M. H.
Gu1,50, C. Y Guan1,55, A. Q. Guo25, A. Q. Guo22, L. B. Guo34, R. P. Guo41, Y.
P. Guo9,g, 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,e, T. Holtmann4, G. Y. Hou1,55, Y.
R. Hou55, Z. L. Hou1, H. M. Hu1,55, J. F. Hu48,j, T. Hu1,50,55, Y. Hu1, G. S.
Huang64,50, K. X. Huang51, L. Q. Huang65, X. T. Huang42, Y. P. Huang1, Z.
Huang39,h, T. Hussain66, N Hüsken22,28, W. Imoehl22, M. Irshad64,50, S.
Jaeger4, S. Janchiv26, Q. Ji1, Q. P. Ji15, X. B. Ji1,55, X. L. Ji1,50, Y. Y.
Ji42, H. B. Jiang42, S. S. Jiang32, X. S. Jiang1,50,55, J. B. Jiao42, Z.
Jiao18, S. Jin35, Y. Jin58, M. Q. Jing1,55, T. Johansson68, N. Kalantar-
Nayestanaki56, X. S. Kang33, R. Kappert56, M. Kavatsyuk56, B. C. Ke73, I. K.
Keshk4, A. Khoukaz61, P. Kiese28, R. Kiuchi1, R. Kliemt11, L. Koch30, O. B.
Kolcu54A, B. Kopf4, M. Kuemmel4, M. Kuessner4, A. Kupsc37,68, W. Kühn30, J. J.
Lane59, J. S. Lange30, P. Larin14, A. Lavania21, L. Lavezzi67A,67C, Z. H.
Lei64,50, H. Leithoff28, M. Lellmann28, T. Lenz28, C. Li40, C. Li36, C. H.
Li32, Cheng Li64,50, D. M. Li73, F. Li1,50, G. Li1, H. Li64,50, H. Li44, H. B.
Li1,55, H. J. Li15, H. N. Li48,j, J. L. Li42, J. Q. Li4, J. S. Li51, Ke Li1,
L. J Li1, L. K. Li1, Lei Li3, M. H. Li36, P. R. Li31,k,l, S. X. Li9, S. Y.
Li53, T. Li42, W. D. Li1,55, W. G. Li1, X. H. Li64,50, X. L. Li42, Xiaoyu
Li1,55, Z. Y. Li51, H. Liang27, H. Liang64,50, H. Liang1,55, Y. F. Liang46, Y.
T. Liang25, G. R. Liao12, J. Libby21, A. Limphirat52, C. X. Lin51, D. X.
Lin25, T. Lin1, B. J. Liu1, C. X. Liu1, D. Liu14,64, F. H. Liu45, Fang Liu1,
Feng Liu6, G. M. Liu48,j, H. M. Liu1,55, Huanhuan Liu1, Huihui Liu16, J. B.
Liu64,50, J. L. Liu65, J. Y. Liu1,55, K. Liu1, K. Y. Liu33, Ke Liu17, L.
Liu64,50, M. H. Liu9,g, P. L. Liu1, Q. Liu55, S. B. Liu64,50, T. Liu9,g, W. M.
Liu64,50, X. Liu31,k,l, Y. Liu31,k,l, 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. Luo72, T. Luo9,g, 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,h, F. E.
Maas14, M. Maggiora67A,67C, S. Maldaner4, S. Malde62, Q. A. Malik66, A.
Mangoni23B, Y. J. Mao39,h, Z. P. Mao1, S. Marcello67A,67C, Z. X. Meng58, J. G.
Messchendorp56,d, G. Mezzadri24A, H. Miao1, T. J. Min35, R. E. Mitchell22, X.
H. Mo1,50,55, N. Yu. Muchnoi10,b, H. Muramatsu60, S. Nakhoul11,e, Y.
Nefedov29, F. Nerling11,e, I. B. Nikolaev10,b, Z. Ning1,50, S. Nisar8,m, S. L.
Olsen55, Q. Ouyang1,50,55, S. Pacetti23B,23C, X. Pan9,g, Y. Pan59, A. Pathak1,
A. Pathak27, M. Pelizaeus4, H. P. Peng64,50, K. Peters11,e, J. Pettersson68,
J. L. Ping34, R. G. Ping1,55, S. Plura28, S. Pogodin29, R. Poling60, V.
Prasad64,50, H. Qi64,50, H. R. Qi53, M. Qi35, T. Y. Qi9,g, S. Qian1,50, W. B.
Qian55, Z. Qian51, C. F. Qiao55, J. J. Qin65, L. Q. Qin12, X. P. Qin9,g, X. S.
Qin42, Z. H. Qin1,50, J. F. Qiu1, S. Q. Qu53, S. Q. Qu36, K. H. Rashid66, K.
Ravindran21, C. F. Redmer28, K. J. Ren32, A. Rivetti67C, V. Rodin56, M.
Rolo67C, G. Rong1,55, Ch. Rosner14, M. Rump61, H. S. Sang64, A. Sarantsev29,c,
Y. Schelhaas28, C. Schnier4, K. Schoenning68, M. Scodeggio24A,24B, K. Y.
Shan9,g, W. Shan19, X. Y. Shan64,50, J. F. Shangguan47, L. G. Shao1,55, M.
Shao64,50, C. P. Shen9,g, H. F. Shen1,55, X. Y. Shen1,55, B.-A. Shi55, H. C.
Shi64,50, R. S. Shi1,55, X. Shi1,50, X. D Shi64,50, J. J. Song15, W. M.
Song27,1, Y. X. Song39,h, S. Sosio67A,67C, S. Spataro67A,67C, F. Stieler28, K.
X. Su69, P. P. Su47, Y.-J. Su55, G. X. Sun1, H. K. Sun1, J. F. Sun15, L.
Sun69, S. S. Sun1,55, T. Sun1,55, W. Y. Sun27, X Sun20,i, Y. J. Sun64,50, Y.
Z. Sun1, Z. T. Sun42, Y. H. Tan69, Y. X. Tan64,50, C. J. Tang46, G. Y. Tang1,
J. Tang51, L. Y Tao65, Q. T. Tao20,i, J. X. Teng64,50, V. Thoren68, W. H.
Tian44, Y. T. Tian25, I. Uman54B, B. Wang1, D. Y. Wang39,h, F. Wang65, H. J.
Wang31,k,l, H. P. Wang1,55, K. Wang1,50, L. L. Wang1, M. Wang42, M. Z.
Wang39,h, Meng Wang1,55, S. Wang9,g, T. J. Wang36, W. Wang51, W. H. Wang69, W.
P. Wang64,50, X. Wang39,h, X. F. Wang31,k,l, X. L. Wang9,g, Y. D. Wang38, Y.
F. Wang1,50,55, Y. Q. Wang1, Y. Y. Wang31,k,l, Ying Wang51, Z. Wang1,50, Z. Y.
Wang1, Ziyi Wang55, D. H. Wei12, F. Weidner61, S. P. Wen1, D. J. White59, U.
Wiedner4, G. Wilkinson62, M. Wolke68, L. Wollenberg4, J. F. Wu1,55, L. H. Wu1,
L. J. Wu1,55, X. Wu9,g, X. H. Wu27, Y. Wu64, Z. Wu1,50, L. Xia64,50, T.
Xiang39,h, H. Xiao9,g, S. Y. Xiao1, Y. L. Xiao9,g, Z. J. Xiao34, X. H.
Xie39,h, Y. G. Xie1,50, Y. H. Xie6, Z. P. Xie64,50, T. Y. Xing1,55, C. F. Xu1,
C. J. Xu51, G. F. Xu1, Q. J. Xu13, S. Y. Xu63, X. P. Xu47, Y. C. Xu55, F.
Yan9,g, L. Yan9,g, W. B. Yan64,50, W. C. Yan73, H. J. Yang43,f, H. X. Yang1,
L. Yang44, S. L. Yang55, Y. X. Yang1,55, Yifan Yang1,55, Zhi Yang25, M.
Ye1,50, M. H. Ye7, J. H. Yin1, Z. Y. You51, B. X. Yu1,50,55, C. X. Yu36, G.
Yu1,55, J. S. Yu20,i, T. Yu65, C. Z. Yuan1,55, L. Yuan2, S. C. Yuan1, X. Q.
Yuan1, Y. Yuan1, Z. Y. Yuan51, C. X. Yue32, A. A. Zafar66, X. Zeng6, Y.
Zeng20,i, Y. H. Zhan51, A. Q. Zhang1, B. L. Zhang1, B. X. Zhang1, G. Y.
Zhang15, H. Zhang64, H. H. Zhang27, H. H. Zhang51, H. Y. Zhang1,50, J. L.
Zhang70, J. Q. Zhang34, J. W. Zhang1,50,55, J. Y. Zhang1, J. Z. Zhang1,55,
Jianyu Zhang1,55, Jiawei Zhang1,55, L. M. Zhang53, L. Q. Zhang51, Lei Zhang35,
P. Zhang1, Shulei Zhang20,i, X. D. Zhang38, X. M. Zhang1, X. Y. Zhang47, X. Y.
Zhang42, Y. Zhang62, Y. T. Zhang73, Y. H. Zhang1,50, Yan Zhang64,50, Yao
Zhang1, Z. H. Zhang1, Z. Y. Zhang69, Z. Y. Zhang36, G. Zhao1, J. Zhao32, J. Y.
Zhao1,55, J. Z. Zhao1,50, Lei Zhao64,50, Ling Zhao1, M. G. Zhao36, Q. Zhao1,
S. J. Zhao73, Y. B. Zhao1,50, Y. X. Zhao25, Z. G. Zhao64,50, A.
Zhemchugov29,a, B. Zheng65, J. P. Zheng1,50, Y. H. Zheng55, B. Zhong34, C.
Zhong65, X. Zhong51, L. P. Zhou1,55, X. Zhou69, X. K. Zhou55, X. R. Zhou64,50,
X. Y. Zhou32, Y. Z. Zhou9,g, J. Zhu36, K. Zhu1, K. J. Zhu1,50,55, S. H. Zhu63,
T. J. Zhu70, W. J. Zhu36, W. J. Zhu9,g, Y. C. Zhu64,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 Minnesota, Minneapolis, Minnesota 55455, USA
61 University of Muenster, Wilhelm-Klemm-Str. 9, 48149 Muenster, Germany
62 University of Oxford, Keble Rd, Oxford, UK OX13RH
63 University of Science and Technology Liaoning, Anshan 114051, People’s
Republic of China
64 University of Science and Technology of China, Hefei 230026, People’s
Republic of China
65 University of South China, Hengyang 421001, People’s Republic of China
66 University of the Punjab, Lahore-54590, Pakistan
67 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
68 Uppsala University, Box 516, SE-75120 Uppsala, Sweden
69 Wuhan University, Wuhan 430072, People’s Republic of China
70 Xinyang Normal University, Xinyang 464000, People’s Republic of China
71 Yunnan University, Kunming 650500, People’s Republic of China
72 Zhejiang University, Hangzhou 310027, People’s Republic of China
73 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 Currently at Istanbul Arel University, 34295 Istanbul, Turkey
e Also at Goethe University Frankfurt, 60323 Frankfurt am Main, Germany
f 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
g 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
h Also at State Key Laboratory of Nuclear Physics and Technology, Peking
University, Beijing 100871, People’s Republic of China
i Also at School of Physics and Electronics, Hunan University, Changsha
410082, China
j Also at Guangdong Provincial Key Laboratory of Nuclear Science, Institute of
Quantum Matter, South China Normal University, Guangzhou 510006, China
k Also at Frontiers Science Center for Rare Isotopes, Lanzhou University,
Lanzhou 730000, People’s Republic of China
l Also at Lanzhou Center for Theoretical Physics, Lanzhou University, Lanzhou
730000, People’s Republic of China
m Also at the Department of Mathematical Sciences, IBA, Karachi , Pakistan

(July 29, 2022)

###### Abstract

Using 448 million $\psi(2S)$ events, the spin-singlet $P$-wave charmonium
state $h_{c}$(1${}^{1}P_{1}$) is studied via the
$\psi(2S)\rightarrow\pi^{0}h_{c}$ decay followed by the
$h_{c}\to\gamma\eta_{c}$ transition. The branching fractions are measured to
be $\mathcal{B}_{\rm
Inc}(\psi(2S)\rightarrow\pi^{0}h_{c})\times\mathcal{B}_{\rm
Tag}(h_{c}\rightarrow\gamma\eta_{c})=(4.22^{+0.27}_{-0.26}\pm 0.19)\times
10^{-4}$ , $\mathcal{B}_{\rm Inc}(\psi(2S)\rightarrow\pi^{0}h_{c})=(7.32\pm
0.34\pm 0.41)\times 10^{-4}$, and $\mathcal{B}_{\rm
Tag}(h_{c}\rightarrow\gamma\eta_{c})=(57.66^{+3.62}_{-3.50}\pm 0.58)\%$, where
the uncertainties are statistical and systematic, respectively. The
$h_{c}$(1${}^{1}P_{1}$) mass and width are determined to be $M=(3525.32\pm
0.06\pm 0.15)\leavevmode\nobreak\ \rm{MeV}/c^{\rm{2}}$ and
${\mathit{\Gamma}}=(0.78^{+0.27}_{-0.24}\pm 0.12)\leavevmode\nobreak\
\rm{MeV}$. Using the center of gravity mass of the three
$\chi_{cJ}$(1${}^{3}P_{J}$) mesons ($M(\rm c.o.g.)$), the $1P$ hyperfine mass
splitting is estimated to be $\Delta_{\rm hyp}=M(h_{c})-M(\rm c.o.g.)=(0.03\pm
0.06\pm 0.15)\leavevmode\nobreak\ \rm{MeV}/c^{\rm{2}}$, which is consistent
with the expectation that the $1P$ hyperfine splitting is zero at the lowest-
order.

BESIII, charmonium spectroscopy, QCD

##  I INTRODUCTION

Despite extensive studies of the charmonium system since its discovery in 1974
Aubert  _et al._ (1974); Augustin  _et al._ (1974); Bacci  _et al._ (1974);
Abrams  _et al._ (1974), knowledge of the singlet state
$h_{c}$(1${}^{1}P_{1}$) is sparse. Only nine decay modes have been observed,
with the $h_{c}\rightarrow\gamma\eta_{c}$ predominant, with a branching
fraction of (50$\leavevmode\nobreak\ \pm\leavevmode\nobreak\ $9)% Zyla  _et
al._ (2020). The $\psi(2S)$ to $h_{c}$(1${}^{1}P_{1}$) hadronic transition is
also known with a relatively large uncertainty, its branching ratio being
(8.6$\leavevmode\nobreak\ \pm\leavevmode\nobreak\ $1.3)$\times 10^{-4}$ Zyla
_et al._ (2020).
Although several measurements of its mass have been performed by the BESIII
Ablikim  _et al._ (2010a, 2012), the CLEO Dobbs  _et al._ (2008), and the E835
Andreotti  _et al._ (2005) collaborations, a better precision is desirable.
This would allow further tests of the hypothesis of zero spin-spin hyperfine
mass splitting relative to the center-of-gravity mass of the three
$\chi_{cJ}$(1${}^{3}P_{J}$) states, defined as Kwong  _et al._ (1987)

$M({\rm c.o.g.})=\frac{M(\chi_{c0})+3M(\chi_{c1})+5M(\chi_{c2})}{9},$

where the $M(\chi_{cJ})$ are the masses of the $\chi_{cJ}$(1${}^{3}P_{J}$)
states. Finally, only one measurement of the $h_{c}$(1${}^{1}P_{1}$) width
exists, provided by the BESIII experiment based on 106 million $\psi(2S)$
events Ablikim  _et al._ (2012). The dataset used in the current analysis
includes that sample.

A more precise knowledge of the $h_{c}$(1${}^{1}P_{1}$) resonance parameters
is also important given the recent discoveries in the $XYZ$ (charmonium-like)
sector with the $h_{c}$ resonance as an intermediate state. Indeed, BESIII
observed the $Z^{\pm}_{c}$(4020) decaying to $\pi^{\pm}h_{c}$ Ablikim  _et
al._ (2013) as well as two resonant structures in the cross section for
$e^{+}e^{-}\rightarrow\pi^{+}\pi^{-}h_{c}$ Ablikim  _et al._ (2017). The
signal decays of these exotic states all employ the tagged channel
$h_{c}\rightarrow\gamma\eta_{c},\eta_{c}\rightarrow{\rm hadrons}$. Thus,
decreasing the uncertainty on the $h_{c}$ branching fractions and resonance
parameters is of the utmost importance.

This article reports an improved determination of the $h_{c}$(1${}^{1}P_{1}$)
resonance parameters taking advantage of $(448.1\leavevmode\nobreak\
\pm\leavevmode\nobreak\ 2.9)\times\leavevmode\nobreak\ 10^{6}$ $\psi(2S)$
collected by the BESIII detector in 2009 and 2012 Ablikim  _et al._ (2018).
The mass and width are extracted by studying the
$\psi(2S)\rightarrow\pi^{0}h_{c}\rightarrow(\gamma\gamma)(\gamma\eta_{c})$
process following the same approach as Ref. Ablikim  _et al._ (2010a). In
particular, this study uses the $\pi^{0}$ recoil mass distribution to
reconstruct the $h_{c}$(1${}^{1}P_{1}$) mass, both inclusively
($\psi(2S)\rightarrow\pi^{0}h_{c}$) and by tagging the decay via the electric
dipole (E1) transition photon from $h_{c}\rightarrow\gamma\eta_{c}$. The
$\pi^{0}$ recoil mass ($RM(\pi^{0})$) is defined as follows:

$RM(\pi^{0})=\sqrt{(E_{\psi(2S)}-E_{\pi^{0}})^{2}-(\bar{p}_{\psi(2S)}-\bar{p}_{\pi^{0}})^{2}},$

where $E_{\pi^{0}}$ ($E_{\psi(2S)}$) and $\bar{p}_{\pi^{0}}$
($\bar{p}_{\psi(2S)}$) are the $\pi^{0}$ ($\psi(2S)$) energy and the momentum
in the reference frame of the laboratory.

For the rest of this article, the two data samples will be referred as
inclusive and tagged (Inc or Tag as a subscript), respectively. The choice of
using two channels is motivated by the necessity of measuring the branching
fractions $\mathcal{B}_{\rm Inc}(\psi(2S)\rightarrow\pi^{0}h_{c})$ and
$\mathcal{B}_{\rm Tag}(h_{c}\rightarrow\gamma\eta_{c})$, the uncertainties of
which are still large Zyla  _et al._ (2020).

##  II BESIII DETECTOR AND DATASETS

The BESIII detector Ablikim  _et al._ (2010b) records symmetric $e^{+}e^{-}$
collisions provided by the BEPCII storage ring Yu  _et al._ (2016), which
operates in the center-of-mass energy range from 2.0 to 4.9 GeV. BESIII has
collected large data samples in this energy region Ablikim  _et al._ (2020).
The cylindrical core of the BESIII detector covers 93% of the full solid angle
and 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 identification
modules interleaved with steel. The charged-particle momentum resolution at
1$\,{\rm GeV}/c$ is $0.5\%$, and the $dE/dx$ resolution is $6\%$ for 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
in the TOF barrel region is 68 ps, while in the end cap is 110 ps.

Simulated data samples produced with a geant4-based Agostinelli  _et al._
(2003) 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 Jadach  _et al._ (2001, 2000). The
inclusive MC sample consists of the production of the $\psi$(2S) resonance,
the ISR production of the $J/\psi$, and the continuum processes
($e^{+}e^{-}\rightarrow e^{+}e^{-}$, $e^{+}e^{-}\rightarrow$ hadrons, and
$e^{+}e^{-}\rightarrow\gamma\gamma$) incorporated in kkmc. The known decay
modes are modelled with evtgen Lange (2001); Ping (2008) using branching
fractions taken from the Particle Data Group (PDG) Zyla  _et al._ (2020), and
the remaining unknown charmonium decays are modelled with lundcharm Chen  _et
al._ (2000); Yang  _et al._ (2014). Final state radiation (FSR) from charged
final state particles is incorporated using the photos Richter-Was (1993)
package. In the signal MC samples, the $\psi(2S)\rightarrow\pi^{0}h_{c}$ decay
and the E1 transition $h_{c}\rightarrow\gamma\eta_{c}$ are both generated by
evtgen.

##  III EVENT SELECTION AND ANALYSIS PROCEDURE

Although charged tracks are not directly used in the reconstruction, candidate
events must have at least two good tracks to suppress background. Good tracks
reconstructed in the MDC must pass the following fiducial and production
vertex cuts. Only tracks with momenta less than 2.0 GeV/$c$ are considered and
they are required to satisfy $|\rm{cos\theta}|$ < 0.93, where $\theta$ is the
angle between the momentum and beam axis. The distance of closest approach to
the interaction point must be less than 10 cm along the beam axis, and less
than 1 cm in the transverse plane. 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 ($|\rm{cos\theta}|$ < 0.80) and more than 50 MeV in
the end cap region (0.86 < $|\rm{cos\theta}|$ < 0.92). 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. 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∘. Finally, to suppress the background contribution from the
continuum processes, such as Bhabha scattering, the total energy deposit in
the EMC ($E^{\rm Tot}_{\rm EMC}$) is required to satisfy 0.6 < $E^{\rm
Tot}_{\rm EMC}$ < 3.2 GeV. Candidate events are required to have at least two
(three) photon candidates passing the fiducial and energy cuts for the
inclusive (tagged) channel. The $\pi^{0}$ candidate is reconstructed via its
decay to $\gamma\gamma$, with both photon candidates laying in the barrel
($|\rm{cos\theta}|$ < 0.80) and having an energy $\geq$ 40 MeV. The di-photon
pair is accepted as a $\pi^{0}$ candidate if its invariant mass satisfies 120
< $M_{\gamma\gamma}$ < 145 MeV/$c^{2}$. Multiple $\pi^{0}$ candidates are
allowed in one event. A one-constraint (1C) kinematic fit fixing the $\pi^{0}$
mass to its known value Zyla  _et al._ (2020) is used to improve the energy
resolution. Candidates with a fit $\chi^{2}>200$ are rejected. The dominant
background for the inclusive sample comes from combining photons from two
different $\pi^{0}$s; the shape of this background is discussed below. The
photon coming from the E1 decay ($\gamma_{\rm Tag}$) is expected to peak
around 500 MeV. Taking into consideration that $\Gamma_{\eta_{c}}\approx
30\,$MeV Zyla  _et al._ (2020), the E1 photon is required to satisfy 465 <
$E\gamma_{\rm{Tag}}$ < 535 MeV and must not form a $\pi^{0}$ with any other
photons in the event. In the tagged channel, if more than one $\pi^{0}$ is
found in the $\pi^{0}$ recoil mass signal region (3.500–3.550 GeV/$c^{2}$) the
$\pi^{0}$ with the minimum 1C fit $\chi^{2}$ is kept.

The three main background sources are $\psi(2S)\rightarrow
J/\psi\,\pi^{+}\pi^{-}$, $\psi(2S)\rightarrow J/\psi\,\pi^{0}\pi^{0}$, and
$\psi(2S)\rightarrow\gamma\chi_{c0}$. The first two are suppressed by
requiring all combinations of the $\pi^{+}\pi^{-}$ ($\pi^{0}\pi^{0}$) recoil
mass to be outside the range $M_{J/\psi}\pm 4$ MeV/$c^{2}$
($M_{J/\psi}\,^{+38}{}_{-8}$ MeV/$c^{2}$), where $M_{J/\psi}$ is the nominal
${J/\psi}$ mass Zyla  _et al._ (2020). This veto window is optimized based on
the figure of merit $\frac{\rm{S}}{\sqrt{\rm{S+B}}}$, where S and B are the
numbers of signal and background events estimated from the MC. This selection
allows to remove $\sim 76\%$ ($\sim 96\%$) of the $\psi(2S)\rightarrow
J/\psi\,\pi^{+}\pi^{-}$ ($\psi(2S)\rightarrow J/\psi\,\pi^{0}\pi^{0}$)
background events, while cutting less than $5\%$ ($1\%$) of the signal events.
No suppression of the $\gamma\chi_{c0}$ decay is needed since the photon
energy is not in the tagged photon energy range, and thus this decay
contributes to the background with a typical combinatorial shape.

The signal efficiencies, $\epsilon_{\rm{Tag}}$ = 12.37% and
$\epsilon_{\rm{Inc}}$ = 14.25%, are estimated for the tagged and inclusive
channels, respectively, based on two signal MC samples of 300000 events each.
The signal shape is modelled with a resolution function (shown in Fig. 1)
based on MC. It is a sum of Gaussian and Crystal Ball functions; both
contribute to symmetric smearing while the latter also includes the low-energy
tail due to $\pi^{0}$ reconstruction. The parameters of this resolution
function are obtained from a signal MC dataset where the
$h_{c}$(1${}^{1}P_{1}$) width is set to 0. The convolution of the resolution
function and a Breit-Wigner distribution is used to describe the signal in the
final dataset to extract $N\rm{{}_{\rm Inc}}$ and $N\rm{{}_{\rm Tag}}$ (i.e.,
the number of the inclusive and tagged events, respectively). Defining
$N(\psi$(2S)) as the total number of $\psi$(2S) events, and
$\mathcal{B}(\pi^{0}\rightarrow\gamma\gamma)$ as the branching fraction of the
$\pi^{0}\rightarrow\gamma\gamma$ decay (taken from PDG Zyla  _et al._ (2020)),
the various branching fractions, $\mathcal{B}$, are obtained as

  * •

$\mathcal{B}_{\rm Inc}(\psi(2S)\rightarrow\pi^{0}h_{c})\times\mathcal{B}_{\rm
Tag}(h_{c}\rightarrow\gamma\eta_{c})=\frac{N_{{\rm Tag}}}{\epsilon_{{\rm
Tag}}\times N(\psi(2S))\times\mathcal{B}(\pi^{0}\rightarrow\gamma\gamma)}$;

  * •

$\mathcal{B}_{\rm Inc}(\psi(2S)\rightarrow\pi^{0}h_{c})=\frac{{N_{\rm
Inc}}}{\epsilon_{{\rm Inc}}\times
N(\psi(2S))\times\mathcal{B}(\pi^{0}\rightarrow\gamma\gamma)}$;

  * •

$\mathcal{B}_{\rm Tag}(h_{c}\rightarrow\gamma\eta_{c})=\frac{\mathcal{B}_{\rm
Inc}(\psi(2S)\rightarrow\pi^{0}h_{c})\times\mathcal{B}_{\rm
Tag}(h_{c}\rightarrow\gamma\eta_{c})}{\mathcal{B}_{\rm
Inc}(\psi(2S)\rightarrow\pi^{0}h_{c})}=\frac{N_{\rm Tag}\times\epsilon_{\rm
Inc}}{N_{\rm Inc}\times\epsilon_{\rm Tag}}$.

FIG. 1: Signal MC $\pi_{0}$ recoil mass (referred as RM($\pi_{0}$))
distributions for the $h_{c}$(1${}^{1}P_{1}$) with its width set to 0. The
solid blue line refers to the global fit result, while the black dots with
error bars are the MC simulation. The blue dashed line represents the Crystal
Ball component, while the red dotted one is the Gaussian component. Fig. (a)
shows the resolution function for the inclusive channel, while the one for the
tagged chain is presented in Fig. (b).

To assess the shape of the irreducible background, a MC dataset (described in
Sec. II) of 400 million inclusive $\psi(2S)$ decays with the signal
contributions removed is studied. This shows that a 5th-order Chebychev
polynomial satisfactorily describes the irreducible background in both
channels.

Fits to the distributions of the $\pi^{0}$ recoil mass are performed
minimizing the negative log-likelihood on the inclusive and tagged channels
separately. The total probability function is constructed from a 5th-order
Chebychev polynomial, to describe the background, added to a convolution of
the resolution function and a Breit-Wigner distribution for the signal. In the
tagged channel, the parameters of the Chebychev polynomial are allowed to
float, as well as the mass and the width of the $h_{c}$(1${}^{1}P_{1}$)
resonance. But for the inclusive channel, the signal shape parameters are
fixed to the values found in the tagged channel, while the parameters of the
Chebychev polynomial are left floating. Input-output tests on MC samples are
used to validate both the model and the fitting procedure; no bias is found.

The fit results are summarized in Table 1. A comparison between this analysis
and the PDGZyla  _et al._ (2020) shows the compatibility between the central
values with improved precision. Fig. 2 presents the fits for the $\pi^{0}$
recoil mass spectra, RM($\pi_{0}$), for both the inclusive and tagged
channels. Separate fits to the 2009 and 2012 datasets yield parameters
compatible with each other within 1.0$\sigma$. The systematic uncertainties
reported in Table 1 are described in detail in the next section.

TABLE 1: Results of the fits to the $\pi^{0}$ recoiling mass spectra with statistical (first) and systematic (second) uncertainties. The last column provides current PDG values. Variable | Value | PDG Value Zyla  _et al._ (2020)
---|---|---
$M(h_{c})$ (MeV/$c^{\rm{2}}$) | $3525.32\pm 0.06\pm 0.15$ | 3525.38 $\pm$ 0.11
${\mathit{\Gamma}}(h_{c})$ (MeV) | $0.78^{+0.27}_{-0.24}\pm 0.12$ | 0.70 $\pm$ 0.28 $\pm$ 0.22
(BESIII Ablikim  _et al._ (2012))
$N_{\rm{Tag}}(h_{c})$ | 23118 ${}^{+1500}_{-1398}$ | —
$\mathcal{B}_{\rm Inc}\times\mathcal{B}_{\rm Tag}$ (10-4) | $4.22^{+0.27}_{-0.26}\pm 0.19$ | $4.58\pm 0.64$
(BESIII Ablikim  _et al._ (2010a))
$4.16\pm 0.48$
(CLEO Dobbs  _et al._ (2008))
$N_{\rm{Inc}}(h_{c})$ | 46187 $\pm$ 2123 | —
$\mathcal{B}_{\rm Inc}$ (10-4) | $7.32\pm 0.34\pm 0.41$ | $8.40\pm 1.30\pm 1.00$
(BESIII Ablikim  _et al._ (2010a))
$9.00\pm 1.5\pm 1.3$
(CLEO Ge  _et al._ (2011))
$\mathcal{B}_{\rm Tag}$ (%) | $57.66^{+3.62}_{-3.50}\pm 0.58$ | $53\pm 7\pm 8$
(BESIII Ablikim  _et al._ (2010a))
$48\pm 6\pm 7$
(CLEO Dobbs  _et al._ (2008))

FIG. 2: Fits to the $\pi^{0}$ recoil mass spectra for the (a) inclusive and
(b) tagged samples. Red solid lines denote the fit results, while the black
dots with error bars are the data. The blue dashed lines represent the signal
component and the black dashed lines are the background. Insets show the
background-subtracted data with the signal shape overlaid.

##  IV SYSTEMATIC UNCERTAINTIES

Sources of systematic uncertainties on the measurement of the
$h_{c}$(1${}^{1}P_{1}$) resonance parameters and branching fractions include
the background line-shape, the mass ranges for the veto of $J/\psi$ from
$\psi(2S)$, the photon energy calibration and reconstruction efficiency, the
$\pi^{0}$ reconstruction efficiency, and the luminosity. The contributions
from each source are shown in Table 2. For each measurement, the total
systematic uncertainty corresponds to a quadrature sum of all individual
sources, which are discussed in detail next.

### Background line-shape.

To determine the systematic uncertainties associated to the fit function, a
4th-order Chebychev polynomial to describe the background behaviour in both
the inclusive and tagged channels is tested. The discrepancy with respect to
the nominal fit result is taken as a systematic uncertainty. Assuming that at
the $\pi^{0}$ reconstruction level (due to the fact that the $\eta_{c}$ is
left to decay inclusively) some background events might be miscounted as
signal, a first-order polynomial is added to the resolution function,
maintaining a 4th order Chebychev polynomial to describe the main background
contribution. The discrepancy between this method and the nominal one is taken
as a systematic uncertainty and summed in quadrature with the result from the
first variation.

### Mass ranges for the veto of $J/\psi$ from $\psi(2S)$.

Mass veto ranges for the $\pi^{0}\pi^{0}$ and $\pi^{+}\pi^{-}$ recoiling
masses are changed according to the observed variations of the figure of
merit. Two different scenarios are tested,

  * •

for the $\pi^{+}$ $\pi^{-}J/\psi$ background, $M_{J/\psi}\pm$ 3
MeV/$c^{\rm{2}}$ and $M_{J/\psi}\pm$ 5 MeV/$c^{\rm{2}}$;

  * •

for the $\pi^{0}$ $\pi^{0}J/\psi$ background, $M_{J/\psi}$${}^{+37}_{-7}$
MeV/$c^{\rm{2}}$ and $M_{J/\psi}$${}^{+39}_{-9}$ MeV/$c^{\rm{2}}$ .

In both cases, smaller veto ranges give the largest changes with respect to
the nominal vetoes; these changes are taken as systematic uncertainties.

### Photon reconstruction efficiency.

The systematic uncertainties arising from potential inconsistencies of the
photon-energy measurements between data and MC simulation are obtained from
Ref. Ablikim  _et al._ (2012).

### Photon energy calibration.

The uncertainty due to the photon energy distribution is obtained from Ref.
Ablikim  _et al._ (2012).

### Signal shape.

This uncertainty includes contributions from the signal line-shape and the 1-C
kinematic fit, and is estimated in Ref. Ablikim  _et al._ (2012).

###  $\pi^{0}$ reconstruction efficiency.

The uncertainty on the branching fractions due to the $\pi^{0}$ reconstruction
efficiency is estimated in Ref. Ablikim  _et al._ (2010c).

### Number of $\psi(2S)$ events.

This uncertainty is estimated in Ref. Ablikim  _et al._ (2018), and included
in the branching fraction uncertainties.

Other sources of the possible systematic uncertainties are studied, but found
to be negligible. These include the bin size, the $\pi^{0}$ recoil mass signal
region, the trigger efficiency, the numbers of $\pi^{0}$ and charged tracks,
the masses and widths of the $\psi(2S)$ and $\eta_{c}$, and the sample sizes
of the MC simulations.

TABLE 2: Summary of the relative systematic uncertainties on the $h_{c}$(1${}^{1}P_{1}$) resonance parameters and branching fractions. The — symbol is used when the systematic sources are found to be negligible or not applicable. For each measurement, the total systematic uncertainty corresponds to a quadrature sum of all individual contributions. Source | $M(h_{c})$ | ${\mathit{\Gamma}}(h_{c})$ | $\mathcal{B}_{\rm Inc}\times\mathcal{B}_{\rm Tag}$ | $\mathcal{B}_{\rm Inc}$ | $\mathcal{B}_{\rm Tag}$
---|---|---|---|---|---
($10^{-3}$) | ($10^{-2}$) | ($10^{-2}$) | ($10^{-2}$) | ($10^{-2}$)
Background | 0.02 | — | 0.44 | 2.86 | —
shape
$J/\psi$ veto | 0.01 | 9.33 | 3.17 | 4.15 | —
$\gamma$ reconstruction | — | — | 3.00 | 2.00 | 1.00
efficiency
$\gamma$ energy | 0.04 | 8.97 | — | — | —
calibration
Signal shape | — | 7.69 | — | — | —
$\pi^{0}$ reconstruction | — | — | 1.00 | 1.00 | —
efficiency
Number of | — | — | 0.65 | 0.65 | —
$\psi(2S)$ events
Total systematic | 0.04 | 15.06 | 4.55 | 5.55 | 1.00
uncertainty

##  V RESULTS AND SUMMARY

In this article, two decay chains involving the $h_{c}(1^{1}P_{1})$ charmonium
state are studied,

  * •

inclusive: $\psi(2S)\rightarrow\pi^{0}h_{c}$ with $h_{c}\rightarrow{\rm
anything}$,

  * •

tagged: $\psi(2S)\rightarrow\pi^{0}h_{c}$ with
$h_{c}\rightarrow\gamma\eta_{c}$.

The measurements of the $h_{c}(1^{1}P_{1})$ resonance parameters and branching
fractions are performed with the world’s largest $\psi(2S)$ data sample,
collected by the BESIII experiment Ablikim  _et al._ (2018). This work
provides the second estimate of the $h_{c}(1^{1}P_{1})$ width, given with a
similar uncertainty as a previous BESIII measurement which used 16 exclusive
decays to reconstruct the $\eta_{c}$ state Ablikim  _et al._ (2012). The
novelty provided by this analysis resides in the increased branching fractions
precision, which is improved between two and three times with respect to the
PDG measurements Zyla  _et al._ (2020). Despite the decreased uncertainties,
this work’s branching ratios remain compatible with both CLEO Dobbs  _et al._
(2008); Ge  _et al._ (2011) and BESIII Ablikim  _et al._ (2010a) estimates.
All the discussed features, summarised in Table 1 altogether with the PDG ones
Zyla  _et al._ (2020), are compatible with the world values Zyla  _et al._
(2020) within one standard deviation.

Furthermore with the $h_{c}(1^{1}P_{1})$ mass estimate of this work, using
current PDG Zyla  _et al._ (2020) values for the center-of-gravity mass of the
three $\chi_{cJ}$(1 ${}^{3}P_{J}$) states, $M(\rm c.o.g.)$, an updated value
for the 1$P$ hyperfine mass splitting ($\Delta_{\rm hyp}=M(h_{c})-M({\rm
c.o.g.})$) may be obtained. This value is predicted by QCD to be identical to
0 at leading-order Kwong  _et al._ (1987). No mass splitting is observed with
this measurement,

$\Delta_{\rm hyp}=0.03\pm 0.06\leavevmode\nobreak\
(\rm{stat.})\pm\leavevmode\nobreak\ 0.15\leavevmode\nobreak\
(\rm{syst.})\leavevmode\nobreak\ \rm{MeV}/c^{\rm{2}}.$

These results on $\Delta_{\rm hyp}$ show that, with the foreseen increase in
$\psi(2S)$ statistics by the BESIII experiment Ablikim  _et al._ (2020),
significant efforts to reduce systematic uncertainties will be necessary.

## 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

  * Aubert  _et al._ (1974) J. J. Aubert _et al._ (E598 Collaboration), Phys. Rev. Lett. 33, 1404 (1974).
  * Augustin  _et al._ (1974) J. E. Augustin _et al._ (SLAC-SP-017), Phys. Rev. Lett. 33, 1406 (1974).
  * Bacci  _et al._ (1974) C. Bacci _et al._ , Phys. Rev. Lett. 33, 1408 (1974), [Erratum: Phys.Rev.Lett. 33, 1649 (1974)].
  * Abrams  _et al._ (1974) G. S. Abrams _et al._ , Phys. Rev. Lett. 33, 1453 (1974).
  * Zyla  _et al._ (2020) P. Zyla _et al._ (Particle Data Group), PTEP 2020, 083C01 (2020).
  * Ablikim  _et al._ (2010a) M. Ablikim _et al._ (BESIII Collaboration), Phys. Rev. Lett. 104, 132002 (2010a).
  * Ablikim  _et al._ (2012) M. Ablikim _et al._ (BESIII Collaboration), Phys. Rev. D 86, 092009 (2012).
  * Dobbs  _et al._ (2008) S. Dobbs _et al._ (CLEO Collaboration), Phys. Rev. Lett. 101, 182003 (2008).
  * Andreotti  _et al._ (2005) M. Andreotti _et al._ (Fermilab E835 Collaboration), Phys. Rev. D 72, 032001 (2005).
  * Kwong  _et al._ (1987) W. Kwong, J. L. Rosner, and C. Quigg, Ann. Rev. Nucl. Part. Sci. 37, 325 (1987).
  * Ablikim  _et al._ (2013) M. Ablikim _et al._ (BESIII Collaboration), Phys. Rev. Lett. 111, 242001 (2013).
  * Ablikim  _et al._ (2017) M. Ablikim _et al._ (BESIII Collaboration), Phys. Rev. Lett. 118, 092002 (2017).
  * Ablikim  _et al._ (2018) M. Ablikim _et al._ (BESIII Collaboration), Chin. Phys. C 42, 023001 (2018).
  * Ablikim  _et al._ (2010b) M. Ablikim _et al._ (BESIII Collaboration), Nucl. Instrum. Methods Phys. Res., Sect. A 614, 345 (2010b).
  * Yu  _et al._ (2016) C. H. Yu _et al._ , Proceedings of IPAC2016, Busan, Korea (2016), 10.18429/JACoW-IPAC2016-TUYA01.
  * Ablikim  _et al._ (2020) M. Ablikim _et al._ (BESIII Collaboration), Chin. Phys. C 44, 040001 (2020).
  * Agostinelli  _et al._ (2003) S. Agostinelli _et al._ (GEANT4), Nucl. Instrum. Methods Phys. Res., Sect. A 506, 250 (2003).
  * Jadach  _et al._ (2001) S. Jadach, B. F. L. Ward, and Z. Wa¸s, Phys. Rev. D 63, 113009 (2001).
  * Jadach  _et al._ (2000) S. Jadach, B. F. L. Ward, and Z. Was, Comput. Phys. Commun. 130, 260 (2000).
  * Lange (2001) D. J. Lange, Nucl. Instrum. Methods Phys. Res., Sect. A 462, 152 (2001).
  * Ping (2008) R.-G. Ping, Chin. Phys. C 32, 599 (2008).
  * Chen  _et al._ (2000) J. C. Chen, G. S. Huang, X. R. Qi, D. H. Zhang, and Y. S. Zhu, Phys. Rev. D 62, 034003 (2000).
  * Yang  _et al._ (2014) R.-L. Yang, R.-G. Ping, and H. Chen, Chin. Phys. Lett. 31, 061301 (2014).
  * Richter-Was (1993) E. Richter-Was, Phys. Lett. B 303, 163 (1993).
  * Ge  _et al._ (2011) J. Y. Ge _et al._ (CLEO Collaboration), Phys. Rev. D 84, 032008 (2011).
  * Ablikim  _et al._ (2010c) M. Ablikim _et al._ (BESIII Collaboration), Phys. Rev. D 81, 052005 (2010c).


# Singular Spectrum Analysis of Exoplanetary Transits

Hossein Fatheddin Leiden Observatory, Leiden University, PO Box 9513, NL-2300
RA Leiden, the Netherlands<EMAIL_ADDRESS>Sedighe Sajadian
Department of Physics, Isfahan University of Technology, Isfahan 84156-83111,
Iran

###### Abstract

Transit photometry is currently the most efficient and sensitive method for
detecting extrasolar planets (exoplanets) and a large majority of confirmed
exoplanets have been detected with this method. The substantial success of
space-based missions such as NASA’s Kepler/K2 and Transiting Exoplanet Survey
Satellite (TESS) has generated a large and diverse sample of confirmed and
candidate exoplanets. Singular Spectrum Analysis (SSA) provides a useful tool
for studying photometric time series and exoplanetary transits. SSA is a
technique for decomposing a time series into a sum of its main components,
where each component is a separate time series that incorporates specific
information from the behavior of the initial time series. SSA can be
implemented for extracting important information (such as main trends and
signals) from the photometry data or reducing the noise factors. The
detectability and accurate characterization of an exoplanetary transit signal
is principally determined by its signal-to-noise ratio (SNR). Stellar
variability of the host star, small planet to star radius ratio, background
noises from other sources in the field of observations and instrumental noise
can cause lower SNRs and consequently, more complexities or inaccuracies in
the modeling of the transit signals, which in turn leads to the inaccurate
inference of the astrophysical parameters of the planetary object. Therefore,
implementing SSA leads to a more accurate characterization of exoplanetary
transits and is also capable of detecting transits with low SNRs ($SNR<10$).
In this paper, after discussing the principles and properties of SSA, we
investigate its applications for studying photometric transit data and
detecting low SNR exoplanet candidates.

Astrostatistics techniques (1886) — Time series analysis (1916) — Exoplanet
detection methods (489) — Transits (1711)

##  1 Introduction

Since the first discovery of a Jovian-mass planet around a solar-type star in
1995 (Mayor & Queloz, 1995), the theoretical and observational efforts for
detecting and characterizing exoplanets have evolved tremendously. Since
exoplanets are extremely faint objects, several indirect methods have been
developed for detecting them. So far, the major exoplanet detection methods
include (Wright & Gaudi, 2013; Fatheddin & Sajadian, 2023): Radial Velocity
measurements (which is based on observing the Doppler shifts in the spectrum
of the planet’s host star), Gravitational Microlensing (which is based on the
temporal magnification in the light of a source star due to a gravitational
binary lens system consisting of a star and a planet) and Transits (which is
based on the periodic dimming in the photometric light curve of a star, caused
by a passing planet).

Currently, the transit photometry has been the most efficient method; out of
the 5535 confirmed exoplanet detections to
date111https://exoplanetarchive.ipac.caltech.edu/docs/counts_detail.html,
about 4100 planets were discovered via the transit method. Data from space-
based missions, such as Kepler/K2 (Borucki et al., 2010; Howell et al., 2014)
and the Transiting Exoplanet Survey Satellite (TESS) (Ricker et al., 2014),
has provided a large and diverse sample of transiting exoplanets. These
missions implement high precision and high cadence time-domain photometric
observations of a large sample of stars in our Galactic neighborhood from the
stable environment of space (Stapelfeldt et al., 2021).

For detecting the periodic dimming caused by a transiting exoplanet in the
time series photometry of a star, the Box-fitting Least Squares (BLS)
algorithm (Kovács et al., 2002) is implemented in most cases. In the BLS
formalism, the transit signals are modelled with periodic box-shaped functions
(upside down boxcar function) with the following parameters: Period, Depth,
Duration and Mid-Transit Time (as Reference Time). For determining the best
fit of these parameters to the transit signal, a range of periods is
considered and the time series is folded to each trial period in that range.
The best least squares fit (maximized over depth) for a period determines its
corresponding signal power in a resulting periodogram. Finally, the parameters
corresponding to the highest power in the periodogram, are considered as the
best transit model.

The most important planet detection metric in the transit method is the
signal-to-noise ratio (SNR) of the transit signal (Kipping, 2023). The SNR is
a precise metric based on the properties and characteristics of the transit
signal and noise in the photometric observations. In this paper, we estimate
the SNR of a transit as the ratio of the measured transit depth ($\delta$) to
the averaged noise estimate from the Combined Differential Photometric
Precision (CDPP) metric ($\sigma_{cdpp}$) (Gaudi, 2000; Kunimoto et al.,
2018):

$\displaystyle S/N=\sqrt{N_{tran}}\frac{\delta}{\sigma_{cdpp}}$ (1)

where $N_{tran}$ is the number of transits and $\sqrt{N_{tran}}$ acts as the
constant of proportionality. The depth of a transit ($\delta$ in equation 1)
is proportional to the ratio of squared planet to star radius (Heller, 2019)
and can be easily measured from the best fit model from the BLS algorithm.

The main idea behind the Combined Differential Photometric Precision (CDPP)
metric is that after removing all long term trends in the light curve, an
estimate of the noise can be measured as the averaged fluctuations in data
(Gilliland et al., 2011). In order to measure the averaged noise estimate
($\sigma_{cdpp}$), we start by detrending the light curve and removing the
outliers. Then, we remove low frequency signals using a Savitzky-Golay filter
(Savitzky & Golay, 1964) with window length equal to the duration of the
transit. Finally, $\sigma_{cdpp}$ is calculated as the standard deviation of a
running mean over the data. The window length of the running mean is also set
to the duration of a single transit.

Singular Spectrum Analysis (SSA) is a powerful tool for time series analysis
with many applications (Broomhead & King, 1986; Vautard & Ghil, 1989; Vautard
et al., 1992). The SSA can be used to decompose a time series into a sum of
its components with separate periodicities occurring on different time scales
(Hassani, 2007). These components can be grouped together as trend,
periodicity and noise signals. The original time series can be reconstructed
by summing all of its components.

The applications of SSA in various fields such as Climate Physics, Economics
and Oceanology have been studied extensively (for example, see Groth & Ghil,
2015; Coussin, 2022; Majumder & Kanjilal, 2019). But, contrary to its great
potential, there are currently very few instances of its use in astronomical
data analysis (Greco et al., 2016; Weinberg & Petersen, 2021; Thekkeppattu et
al., 2023).

In this paper, after reviewing the principles and properties of the SSA
algorithm in Section 2, we will study its applications to the analysis of the
photometric time series data and exoplanetary transits in section 3. In
Section 4, we investigate the detectability and characterization of low SNR
transits with implementing the SSA algorithm and report our results. Finally,
in the last section, we discuss and summarize the results.

##  2 Singular Spectrum Analysis (SSA)

The Singular Spectrum Analysis (SSA) is a model-free nonparametric technique
for time series analysis with the main assumption that the time series can be
decomposed into a sum of different components such as trend, modulated
periodicities and noise (Agarwal et al., 2020). The name of ”singular spectrum
analysis” refers to the spectrum of eigenvalues in a singular value
decomposition (SVD) of the covariance matrix.

By decomposing a time series via SSA, one can:

  * •

Extract the main trends (which can be defined as the long-term change in the
level of the time series)

  * •

Extract the periodicities that form the overall composition of the time series

  * •

Reconstruct a smoothed version of the time series by denoising the data

  * •

Fill in the missing observations or forecast the future trends

  * •

Perform a linear frequency filter

The first step in performing SSA is embedding the time series into a
”Trajectory matrix” (which is explained in section 2.1). Then we can decompose
that trajectory matrix into a collection of eigenvalues and eigenvectors with
the SVD (Golyandina & Zhigljavsky, 2013).

The SSA algorithms for decomposing and reconstructing the time series based on
its trajectory matrix are explained in sections 2.2 and 2.3, respectively.

###  2.1 The Trajectory Matrix

For performing SSA, we need a procedure that takes a univariate time series
and turn it into a multivariate set of observations (Elsner & Tsonis, 1996).

We start by considering a time series of length $N$ as
$P_{N}=(p_{1},p_{2},...,p_{N})$. We also consider a fixed integer $L$ as the
”window length”, where: $2<L<N/2$. Then we form column vectors (each with the
length L) from subsets of the time series as:

$\displaystyle\textbf{X}_{i}=(p_{i},...,p_{i+L-1})^{T}$ (2)

where $i=1,2,...,K$ and $K=N-L+1$.

Now, we can use these vectors to construct the $L\times K$ dimensional
trajectory matrix:

$\displaystyle\textbf{X}=[\textbf{X}_{1}:...:\textbf{X}_{K}]=\begin{bmatrix}p_{1}&p_{2}&p_{3}&\dots&p_{K}\\\
p_{2}&p_{3}&p_{4}&\dots&p_{K+1}\\\ p_{3}&p_{4}&p_{5}&\dots&p_{K+2}\\\
\vdots&\vdots&\vdots&\ddots&\vdots\\\
p_{L}&p_{L+1}&p_{L+2}&\dots&p_{N}\end{bmatrix}_{L\times K}$ (3)

We notice that this matrix has equal anti-diagonal arrays which is a
characteristic of the Hankel matrices (Aoki, 1983).

In the SSA formalism, it is necessary that the window length ($L$) is
sufficiently large, so that each $\textbf{X}_{i}$ column vector (equation 2)
would incorporate an essential part of the initial time series and we can
separate underlying periodicities from the main trends. If $L$ is large
enough, we are able to consider each $\textbf{X}_{i}$ vector as a separate
time series which embodies the dynamics of certain characteristics from the
original time series ($P_{N}$).

It is interesting to note that this notion of embedding and constructing the
trajectory matrix is also very useful in the context of data analysis for the
non-linear dynamical systems and chaos theory (Packard et al., 1980; Takens,
1981). In this context, the variables that describe the evolution of the
dynamical system (which are a set of $N$ first-order differential equations)
are reduced to a single differential equation of $N^{th}$ order by successive
differentiation. So, a difficult multi-variable description of the system is
replaced with a single-variable description (Elsner & Tsonis, 1996).

###  2.2 Decomposition

At this point, we can decompose the trajectory matrix X (equation 3) with the
SVD algorithm. In the SVD formalism (Klema & Laub, 1980), factorization of the
matrix X is written as:

$\displaystyle\textbf{X}=\textbf{U}\mathbf{\Sigma}\textbf{V}^{T}$ (4)

where U and V (referred to as ”left” and ”right” singular vectors) are
orthogonal unitary matrices with dimensions $L\times L$ and $K\times K$,
respectively. $\textbf{V}^{T}$ is the transpose matrix of V and
$\mathbf{\Sigma}$ is a $L\times K$ rectangular diagonal matrix.

Furthermore, we can use equation 4 and some basic linear algebra to show the
following relations:

$\displaystyle\textbf{X}\textbf{X}^{T}=\textbf{U}\mathbf{\Sigma}\textbf{V}^{T}\textbf{V}\mathbf{\Sigma}^{T}\textbf{U}^{T}=\textbf{U}(\mathbf{\Sigma}\mathbf{\Sigma}^{T})\textbf{U}^{T}$
(5)
$\displaystyle\textbf{X}^{T}\textbf{X}=\textbf{V}\mathbf{\Sigma}^{T}\textbf{U}^{T}\textbf{U}\mathbf{\Sigma}\textbf{V}^{T}=\textbf{V}(\mathbf{\Sigma}^{T}\mathbf{\Sigma})\textbf{V}^{T}$
(6)

From equation 5, one can simply show that the columns of the U are
eigenvectors of $\textbf{X}\textbf{X}^{T}$ and that they form an orthonormal
basis set spanning the subsets of the time series in the columns of the
trajectory matrix.

Similarly, from equation 6, we can show that the columns of V are eigenvectors
of $\textbf{X}^{T}\textbf{X}$ and that they form an orthonormal basis set
spanning the subsets of the time series in the rows of the trajectory matrix.

Additionally, from both equations 5 and 6, we see that the elements of
$\mathbf{\Sigma}$ (referred to as singular values) are the square roots of the
non-zero eigenvalues of $\textbf{X}\textbf{X}^{T}$ or
$\textbf{X}^{T}\textbf{X}$.

The trajectory space of the time series (defined by U and V) is at most
L-dimensional. This occurs when all of the columns in the trajectory matrix
are linearly independent. But, if there are any columns that are linearly
dependent, we would have arrays with zero values in the singular values matrix
($\mathbf{\Sigma}$) which consequently means that the dimension of the
trajectory space would be less than $L$. To avoid this, we only consider the
trajectory space as d-dimensional, where d is the rank of the trajectory
matrix. The rank of a matrix is defined as the maximal number of linearly
independent columns of that matrix (Banerjee & Roy, 2014).

Finally, the decomposition of the trajectory matrix can be written as:

$\displaystyle\textbf{X}=\textbf{X}_{1}+...+\textbf{X}_{d}=\sum^{d}_{i=1}\sqrt{\lambda_{i}}U_{i}V_{i}^{T}$
(7)

where $(\lambda_{1},\lambda_{2},...,\lambda_{d})$ are the eigenvalues of
$\textbf{X}\textbf{X}^{T}$ or $\textbf{X}^{T}\textbf{X}$. In the SSA
formalism, $X_{i}=\sqrt{\lambda_{i}}U_{i}V_{i}^{T}$ are called ”elementary
matrices” and $(\sqrt{\lambda_{i}},U_{i},V_{i}^{T})$ are called ”Eigen-
triples”.

It is also worth noting that $\sqrt{\lambda_{i}}$ can be interpreted as a
scaling factor that indicates the relative importance of each eigentriple in
the expansion of X (for more information, see: Golyandina et al. (2001)).

###  2.3 Reconstruction

Now that we have decomposed the initial time series ($P_{N}$) into its
elementary matrices (equation 7), we need a procedure that enables us to
reproduce $P_{N}$ as the sum of $m$ reconstructed time series:

$\displaystyle P_{N}=\sum^{m}_{k=1}\tilde{p}^{(k)}_{n}$ (8)

where $\tilde{p}_{n}^{(k)}$ are separate component time series from $P_{N}$.

We start by partitioning the indices in equation 7 $(1,...,d)$ into $m$
disjoint subsets $I_{1},...,I_{m}$. If $I=(i_{1},...,i_{p})$, then the
corresponding $\textbf{X}_{I}$ matrix is defined as:

$\displaystyle\textbf{X}_{I}=\textbf{X}_{i1}+...+\textbf{X}_{ip}$ (9)

If we compute the resulting matrices $\textbf{X}_{I}$ for each
$I\in(I_{1},...,I_{m})$, the grouped expansion of X can be written as:

$\displaystyle\textbf{X}=\textbf{X}_{I1}+...+\textbf{X}_{Im}$ (10)

Since our intention is to reconstruct the component time series with length
$N$ ($\tilde{p}_{n}^{(k)}$) from Each matrix $\textbf{X}_{Ij}$, we implement a
method called ”Diagonal averaging” (Hassani, 2007). In diagonal averaging, the
values of the reconstructed time series $\tilde{p}_{n}^{(k)}$ are defined as
averages of the corresponding anti-diagonals of the matrices
$\textbf{X}_{Ij}$:

$\displaystyle
A_{k}=\begin{cases}\frac{1}{k}\sum^{k}_{m=1}X_{m,k-m+1}&\text{for $1\leq
k<L$}\\\ \frac{1}{L}\sum^{L}_{m=1}X_{m,k-m+1}&\text{for $L\leq k\leq K$}\\\
\frac{1}{N-k+1}\sum^{N-K+1}_{m=k-K+1}X_{m,k-m+1}&\text{for $K<k\leq
N$}\end{cases}$ (11)

where $A_{k}$ are the elements of time series $\tilde{p}_{n}^{(k)}$. The
resulting component time series from diagonal averaging are the main and final
result of the SSA algorithm and by adding them together, we get the original
time series (equation 8).

In order to understand the relationship between each component time series and
whether if they should be grouped together or separately (as Trend,
Periodicities or Noise), we use the ”W-Correlation Matrix” (Bógalo et al.,
2017). Let us consider two component time series $\tilde{p}_{i}$ and
$\tilde{p}_{j}$ derived from matrices $\tilde{X}_{Ii}$ and $\tilde{X}_{Ij}$
respectively. The weighted inner product for these two time series is:

$\displaystyle(\tilde{p}_{i},\tilde{p}_{j})_{w}=\sum_{k=1}^{N}w_{k}\tilde{p}_{i,k}\tilde{p}_{j,k}$
(12)

where $w_{k}$ are the weights, defined as:

$\displaystyle w_{k}=\begin{cases}k&\text{for $1\leq k<L$}\\\ L&\text{for
$L\leq k\leq K$}\\\ N-k+1&\text{for $K<k\leq N$}\end{cases}$ (13)

Now, we can define the $L\times L$ W-correlation matrix ($W_{corr}$) for every
component time series. The arrays of this matrix are:

$\displaystyle
W_{i,j}=\frac{(\tilde{p}_{i},\tilde{p}_{j})_{w}}{||\tilde{p}_{i}||_{w}||\tilde{p}_{j}||_{w}}$
(14)

From this equation, it can be easily shown that if $W_{ij}=0$, the two
component time series are orthogonal and completely separable (i.e. should not
be grouped together). Consequently, if $W_{ij}$ is close to 1, the two
components are correlated and should be grouped together.

##  3 SSA of Photometric Time Series

Until now, we have discussed and reviewed the SSA algorithm for an arbitrary
time series. A very important category of time series in astronomy, is the
time-domain photometry; which is the study of how the light (in numerous
passbands or filters) from astronomical objects change over time.

The variations in a photometric time series (or light curve) of an star can
reveal many useful information, including but not limited to: periodic and
semi-regular stellar pulsations (Cox, 1974), stars with outbursts (Forgan &
Rice, 2010), asteroseismology studies of variable stars (Di Mauro, 2016),
eclipses and occultation (in binary stars or planetary transits) (Willems et
al., 2006) and gravitational microlensing (Mao, 2012; Sajadian & Fatheddin,
2023).

After the immense success of NASA’s Kepler/K2 space-based missions, TESS is
providing a large and diverse amount of photometric time series. Unlike
Kepler/K2 which covered a fixed region of sky (Borucki et al., 2010), TESS is
continuously surveying $85\%$ of the sky divided into sectors. Each sector is
observed for two orbits of the telescope around the Earth (about 27 days on
average) with 2 seconds cadence (Ricker et al., 2014).

In this section, we investigate the implementation and applications of the SSA
algorithm (section 2) for simulated photometric time series and an actual
light curve from TESS (TIC $94727308$) as our case studies.

###  3.1 Case Study 1: Simulated Photometric Time Series

Figure 1: Left Panels: Simulated light curves of low SNR planetary transits
with different depths ($\delta$) and noise levels ($\sigma_{cdpp}$). In each
simulation, the period is set to one day. The red curve in each figure
corresponds to the best transit model derived from the BLS algorithm. The blue
curves correspond to the actual transits without the noise factors. Right
Panels: The resulting light curves after the SSA de-noising. Before SSA,
differences can be noted between the BLS models and the actual transits. These
differences are resolved after SSA. On top of each figure, Depth,
$\sigma_{cdpp}$ and SNR of the transits are reported.

In left panels of figure 1, we have generated three artificial light curves
for stars with different exoplanetary transit signals. Each light curve
contains a noise-like background feature which corresponds to the variability
of the host star and the photometric measurement error of the telescope. The
periods are all set equal to one day, but we have considered a different
planetary radius in each simulation in order to simulate different transit
depths.

As we can see in each light curve in the left panel, there is a discrepancy
between the actual transit signals (blue curves) and the BLS transit model
(red curves) derived from the ”noisy” data. Specifically, simulations 1 and 1
show considerable differences between actual transit depths and depths derived
from the model. As it was mentioned in the introduction, the depth of a
transit is directly used for inferring the planetary radius, hence even a
slight discrepancy between the model and actual transit depth may yield
inaccurate results for the astrophysical characteristics of the planet. For
simulation 1, where the signal has a very low SNR (4.63), there is also a
difference in mid-transit reference times of the model and signal.

Figure 2: Left Panels: Visualizations of W-correlation matrices for SSA
components of the simulations. It can be noted that some component time series
show weak correlations with most of the other components which represents
their noise-like behaviour. Right Panels: W-correlation matrices of the de-
noised light curves. The number of noise components have decreased after SSA.
The colorbar at the right of each figure shows the mapping from scalar
weighted correlations in the W-correlation matrix to colors.

For each simulation, we perform SSA with window length $L=100$. We note that
since the length of each time series is $N=8000$, our choice of L is well
below N/2. The resulting $100\times 100$ W-correlation matrix visualizations
for $100$ component time series of each simulation are shown in the left
panels of figure 2. In these visualizations, we see that some components are
slightly correlated with most of the other component time series. This can be
explained by considering the fact that the noise factors have no discernible
”separable” structures (regardless of how large the window length is) and as a
result, they have a small non-zero weighted inner products with most of the
other components.

In order to separate the noise components, we loop over the W-correlation
matrix while counting the number of weekly correlated ($W_{ij}<0.2$)
components for each component. If this number is larger than a critical value,
we omit the corresponding component. The idea behind this procedure is that we
are keeping the components that were either completely separated from the
other components or have strong correlations, while deleting the components
with mostly negligible correlations (indication of their noise-like
behaviour). Finally, we can group the remaining components to reconstruct a
de-noised version of the time series (section 2.3).

For each simulated light curve, the de-noised version is illustrated in the
right panels of figure 1. As we can see, after performing the de-noising
procedure, we have less fluctuations (lower $\sigma_{cdpp}$) in the light
curves and the new BLS transit models, which are derived from the de-noised
data, match the real transit signals. In each case, the SNR (equation 1) has
also increased substantially. In particular, the transit signal from
simulation 1 was undetectable before the SSA, as the SNR of its signal was
well below the detection criteria ($SNR>10$) for a transiting exoplanet. This
implies that SSA can be implemented to detect low SNR exoplanetary signals in
photometric time series (which will be discussed in section 4 ).

For comparison, we can now perform SSA on the de-noised versions of the
simulated data and calculate the corresponding W-correlation matrices (right
panels of figure 2). As we can see in the new visualizations, the number of
components with noise-like behaviours has decreased, indicating a de-noised
version of the original time series.

###  3.2 Case Study 2: TIC 94727308

Figure 3: Figure (a) and (c) show the initial light curve and the river plot
for TOI-$4622$, respectively. The light curve is folded over seven transits
and the river plot shows the evolution of the noisy periodic transit signals
over the course of the entire light curve. Figure (b) and (d) show the
denoised versions of these plots after the implementation of SSA. The SNR of
the transit has increased from 13.37 to 15.59.

TIC $94727308$ is the photometric time series for the TESS Object of Interest
(TOI) $4622$ (Gaia DR3 $170033017904053376$). TOI-$4622$ is a very faint star
with I-band apparent magnitude of $16.1779\pm 0.009$ (from TESS) and G-band
mean magnitude of $17.621546$ (from Gaia). Additionally, data from
Astrophysical parameters table of the Gaia third data release (Gaia
Collaboration et al., 2023) indicate that this star is a faint young M dwarf.

The TESS photometry data for TOI-$4622$ includes two $600$ seconds exposure
time observations from sectors $43$ and $44$ conducted in $2021$. After
normalizing and stitching (via the LightKurve package) these time series, the
BLS periodogram reveals a relatively low SNR transit signal with period
$P=1.297$ days. We fold the detrended light curve over a time interval of
$7\times P$ in order to inspect and compare seven transits in figure 3. As we
can see, the transit signals seem very weak and they almost resemble noise
rather than signal.

For an additional and more insightful visualization, we also make the ”river
plot” for this time series (figure 3). River plots depict the relative
fluctuations in flux values by splitting the light curve into parts of equal
length (i.e. the transit period), and displaying them side by side (Lightkurve
Collaboration et al., 2018). This allows us to compare these variations over
the course of an entire light curve. The river plot in figure 3 shows a noisy
periodic dimming in the light of the TOI-$4622$.

As described in the previous case study, We perform SSA decomposition on the
time series in order to reconstruct a de-noised version of the light curve, by
excluding the noise-like components.

After reconstructing the de-noised version of the time series, we once again
fold the data and make the river plot (figures 3 and 3 respectively). These
figures indicate that by omitting the noise components in SSA, we have
successfully retrieved the exoplanetary transit signal which was not
accurately detectable due to the faintness of the host star. We discuss the
further implementation of this de-noising algorithm for detecting and modeling
low SNR exoplanets in the following section.

##  4 Low SNR Candidates from TESS

Figure 4: Folded light curves of six encouraging TESS candidates. The left
side light curve in each panel corresponds to the initial time series, whereas
the right side ones belong to the resulting reconstructed data from SSA. The
red curve in each figure corresponds to the best transit model derived from
the BLS algorithm. On top of each figure, Period, Depth, $\sigma_{cdpp}$ and
SNR of the transits are reported. A decrease in fluctuations, improvements in
the SNRs and changes in the best fit BLS models can be noted in the folded
light curves after the SSA.

As mentioned in the introduction, the depth of an exoplanetary transit
($\delta$) is proportional to the ratio of squared planet to star radius
(Heller, 2019). Consequently, if an exoplanet is relatively too small compared
to its host star (as is the case with the terrestrial exoplanets), its transit
signal would have a low SNR and may even remain unnoticeable in the
fluctuations of the raw data (as demonstrated in the third simulation of
section 3.1).

Additionally, the variability or the faintness of a host star or background
and instrumental noises can cause a higher range of fluctuations in the light
curve which would result in a higher $\sigma_{cdpp}$ and a lower SNR,
respectively. As the SNR is also dependant on the number of transits (N in
equation 1), in exoplanetary transits with relatively long periods and the
same amount of observations, the SNR would be lower as well.

One of the main applications of SSA is time series de-noising and smoothing
and it has been implemented as an effective noise removal technique in
different contexts (Zabalza et al., 2014; Yuan-feng & Mei, 2005). As we have
discussed in our case study in the previous section, SSA can also be used for
de-noising a photometric time series in order to retrieve an exoplanetary
transit signal. There are currently more than $250$ TESS project candidates in
the Exoplanet Follow-up Observing
Program222https://exofop.ipac.caltech.edu/tess/view_toi.php (ExoFOP) (Akeson &
Christiansen, 2019) that show very low SNR transit-like signals.

This indicates a great potential for SSA to investigate the true natures of
these signals and provides a preliminary insight into whether there should be
further photometric or spectroscopic observations. we implemented our SSA de-
noising algorithm and compared the SNRs of these signals before and after this
procedure. If the transit-like signal disappears (which manifests in a lower
SNR) it can be deduced that the signal was due to minor background
periodicities or a noise-like feature. Whereas, if the signal becomes more
prominent (a higher SNR), the target star would make an appropriate candidate
for follow-up photometric or spectroscopic observations to validate the
transit.

Furthermore, as low SNR transits have very weak signals that are infused with
noise, the BLS model derived from the raw data may yield inaccurate planetary
radii and other astrophysical parameters. In these situations, SSA can be used
to resolve these issues.

Table 1: Photometric accuracy and Transit data for the six TESS candidates displayed in figure 4. By comparing columns corresponding to the initial and SSA time series, we can note the changes in transit depths, decreases in fluctuations ($\sigma_{cdpp}$) and improvements in SNRs.  TIC ID | TOI | P [days] |  $\delta$ [ppm] (initial) |  $\sigma_{cdpp}$ [ppm] (initial) | SNR (initial) |  $\delta$ [ppm] (SSA) |  $\sigma_{cdpp}$ (SSA) | SNR (SSA)
---|---|---|---|---|---|---|---|---
65416038 | 4190 | 15.616 | 664.914 | 166.295 | 7.996 | 429.464 | 88.179 | 9.741
122617317 | 4311 | 15.065 | 814.866 | 180.733 | 9.017 | 719.016 | 131.796 | 10.911
88114177 | 5531 | 4.535 | 209.169 | 63.321 | 10.956 | 184.599 | 51.509 | 11.886
354173595 | 5974 | 1.726 | 618.874 | 315.647 | 7.593 | 450.852 | 200.061 | 8.728
39018923 | 2216 | 5.681 | 1822.132 | 394.611 | 13.060 | 1693.483 | 356.496 | 14.436
195187769 | 5550 | 13.956 | 266.438 | 96.917 | 6.734 | 223.957 | 69.411 | 8.176

Figure 4 displays light curves of the most promising TESS candidates that
showed an improved SNR after the SSA. The left side light curve in each panel
corresponds to the initial time series, whereas the right side ones belong to
the resulting reconstructed data. More information for each candidate
(including the period, depth ($\delta$), averaged noise estimate
($\sigma_{cdpp}$) and SNR of the transit signal) is given in table 1.

Figure 5: Visualization of w-correlation matrix for TIC-$65416038$ with
$L=100$.

As an example of how this process works, let us consider the folded time
series of TIC-$65416038$ (figure 4). There is a transit signal with period
$P=15.616$ and $SNR=7.996$. The w-correlation matrix for a SSA decomposition
of this time series is shown in figure 5. After removing the components with a
noise-like behaviour, the resulting reconstructed light curve is shown in the
right side of figure 4. The SNR of the transit signal has increased to
$9.741$, which means that it is a genuine signal and not part of the
underlying noise. Additionally, we can see differences in the durations and
depths of the transit between the initial and the de-noised BLS models of the
time series.

In figure 4, the SNRs of transit signals were improved by the SSA de-noising,
but these improvements were less substantial compared to the simulations of
case study 1 (section 3.1). This can be explained by considering the fact that
as a result of longer transit periods and shorter observation, these TESS
light curves have very small number of transits (N) in their available
photometric time series. AS SSA increases the SNR by decreasing the noise-like
fluctuations ($\sigma_{cdpp}$) in data, future observations of the same
targets leads to a higher N and an increased SNR. Nevertheless, the planetary
transits can be distinguished easier from the noise and other intrinsic
fluctuations in the de-noised light curves. These candidates are very
encouraging for further ground or space based observations.

##  5 Conclusions

Transit photometry is a very sensitive method for detecting and characterizing
exoplanets. Photometric data of exoplanetary transits from space based
missions like Kepler/k2 and TESS has considerably expanded our knowledge of
planets beyond our Solar system and the accurate modeling and characterization
of these transits is of great importance.

Singular Spectrum Analysis is a powerful method for decomposing a time series
into its constituting components. Reconstructing a new time series by
selecting different groups of these components based on their weighted
correlation values, can incorporate and separate the behaviour of a specific
trend, periodic signal or different noise factors.

In this work, after discussing and reviewing the principles of basic SSA, we
studied its further applications for the analysis of photometric time series
by simulating different sets of time series incorporating planetary transit
signals with low SNRs bellow the detection criteria ($SNR<10$) . We have shown
that by implementing SSA, we can separate the main variability trends of a
star and the noise factors from its light curve and then use the BLS algorithm
to look for transit-like signals with an increased SNR that could indicate the
existence of a planetary object. Provided that the number of transits in a
light curve is large enough (more than 5 transits per light curve) the
increase in SNR can be as substantial as 30 units.

Additionally, In low SNR transit signals, the BLS algorithm may yield
inaccurate results for the transit depth, duration or mid-transit time which
in turn leads to the inaccurate inference of the astrophysical parameters of
the planet. we have shown that SSA can be used as a noise removal technique in
these circumstances for the accurate modeling of low SNR transiting
exoplanets.

In this paper, we mainly applied our implementation of SSA to photometric time
series from NASA’s TESS mission, but this algorithm can be easily applied to
other current or future ground and space based photometry surveys.

A Python implementation of all the algorithms discussed in this paper can be
found in Fatheddin (2024).

## Acknowledgements

This work was supported by the benefactors of the Swaantje Mondt fund. H. F.
thanks Huub Rottgering, Wouter Schrier and Mozafar Allahyari for their help
and support.

This paper includes data collected by the TESS mission. Funding for the TESS
mission is provided by the NASA’s Science Mission Directorate.

This work has made use of data from the European Space Agency (ESA) mission
Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing
and Analysis Consortium (DPAC,
https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has
been provided by national institutions, in particular the institutions
participating in the Gaia Multilateral Agreement.

This research made use of LightKurve, a Python package for Kepler and TESS
data analysis (Lightkurve Collaboration et al., 2018).

We thank the reviewer for useful comments and suggestions.

## References

  * Agarwal et al. (2020) Agarwal, A., Alomar, A., & Shah, D. 2020, arXiv e-prints, arXiv:2006.13448, doi: 10.48550/arXiv.2006.13448
  * Akeson & Christiansen (2019) Akeson, R., & Christiansen, J. 2019, in American Astronomical Society Meeting Abstracts, Vol. 233, American Astronomical Society Meeting Abstracts #233, 140.09
  * Aoki (1983) Aoki, M. 1983, doi: 10.1007/978-3-642-45565-0
  * Banerjee & Roy (2014) Banerjee, S., & Roy, A. 2014, Linear Algebra and Matrix Analysis for Statistics, 1–559, doi: 10.1201/b17040
  * Bógalo et al. (2017) Bógalo, J., Poncela, P., & Senra, E. 2017.  https://api.semanticscholar.org/CorpusID:55917954
  * Borucki et al. (2010) Borucki, W. J., Koch, D., Basri, G., et al. 2010, Science, 327, 977, doi: 10.1126/science.1185402
  * Broomhead & King (1986) Broomhead, D., & King, G. P. 1986, Physica D: Nonlinear Phenomena, 20, 217, doi: https://doi.org/10.1016/0167-2789(86)90031-X
  * Coussin (2022) Coussin, M. 2022, Journal of International Money and Finance, 120, 102532, doi: https://doi.org/10.1016/j.jimonfin.2021.102532
  * Cox (1974) Cox, J. P. 1974, Reports on Progress in Physics, 37, 563, doi: 10.1088/0034-4885/37/5/001
  * Di Mauro (2016) Di Mauro, M. P. 2016, in Frontier Research in Astrophysics II (FRAPWS2016), 29, doi: 10.22323/1.269.0029
  * Elsner & Tsonis (1996) Elsner, J. B., & Tsonis, A. A. 1996, doi: 10.1007/978-1-4757-2514-8
  * Fatheddin (2024) Fatheddin, H. 2024, AstroFatheddin/SSA_Transits: v2.0.0, Zenodo, doi: 10.5281/zenodo.10938904
  * Fatheddin & Sajadian (2023) Fatheddin, H., & Sajadian, S. 2023, AJ, 166, 140, doi: 10.3847/1538-3881/aced8b
  * Forgan & Rice (2010) Forgan, D., & Rice, K. 2010, MNRAS, 402, 1349, doi: 10.1111/j.1365-2966.2009.15974.x
  * Gaia Collaboration et al. (2023) Gaia Collaboration, Vallenari, A., Brown, A. G. A., et al. 2023, A&A, 674, A1, doi: 10.1051/0004-6361/202243940
  * Gaudi (2000) Gaudi, B. S. 2000, ApJ, 539, L59, doi: 10.1086/312828
  * Gilliland et al. (2011) Gilliland, R. L., Chaplin, W. J., Dunham, E. W., et al. 2011, ApJS, 197, 6, doi: 10.1088/0067-0049/197/1/6
  * Golyandina et al. (2001) Golyandina, N., Nekrutkin, V., & Zhigljavsky, A. 2001, Monographs on Statistics and Applied Probability, 90, doi: 10.1201/9781420035841
  * Golyandina & Zhigljavsky (2013) Golyandina, N., & Zhigljavsky, A. 2013, Singular Spectrum Analysis for Time Series, doi: 10.1007/978-3-642-34913-3
  * Greco et al. (2016) Greco, G., Kondrashov, D., Kobayashi, S., et al. 2016, in Astrophysics and Space Science Proceedings, Vol. 42, The Universe of Digital Sky Surveys, ed. N. R. Napolitano, G. Longo, M. Marconi, M. Paolillo, & E. Iodice, 105, doi: 10.1007/978-3-319-19330-4_16
  * Groth & Ghil (2015) Groth, A., & Ghil, M. 2015, Journal of Climate, 28, 7873, doi: 10.1175/JCLI-D-15-0100.1
  * Hassani (2007) Hassani, H. 2007, University Library of Munich, Germany, MPRA Paper, 5, doi: 10.6339/JDS.2007.05(2).396
  * Heller (2019) Heller, R. 2019, A&A, 623, A137, doi: 10.1051/0004-6361/201834620
  * Howell et al. (2014) Howell, S. B., Sobeck, C., Haas, M., et al. 2014, PASP, 126, 398, doi: 10.1086/676406
  * Kipping (2023) Kipping, D. 2023, MNRAS, 523, 1182, doi: 10.1093/mnras/stad1492
  * Klema & Laub (1980) Klema, V., & Laub, A. 1980, IEEE Transactions on Automatic Control, 25, 164, doi: 10.1109/TAC.1980.1102314
  * Kovács et al. (2002) Kovács, G., Zucker, S., & Mazeh, T. 2002, A&A, 391, 369, doi: 10.1051/0004-6361:20020802
  * Kunimoto et al. (2018) Kunimoto, M., Matthews, J. M., Rowe, J. F., & Hoffman, K. 2018, AJ, 155, 43, doi: 10.3847/1538-3881/aaa005
  * Lightkurve Collaboration et al. (2018) Lightkurve Collaboration, Cardoso, J. V. d. M., Hedges, C., et al. 2018, Lightkurve: Kepler and TESS time series analysis in Python, Astrophysics Source Code Library, record ascl:1812.013.  http://ascl.net/1812.013
  * Majumder & Kanjilal (2019) Majumder, S., & Kanjilal, P. P. 2019, Pure and Applied Geophysics, 176, 3769, doi: 10.1007/s00024-019-02140-4
  * Mao (2012) Mao, S. 2012, Research in Astronomy and Astrophysics, 12, 947, doi: 10.1088/1674-4527/12/8/005
  * Mayor & Queloz (1995) Mayor, M., & Queloz, D. 1995, Nature, 378, 355, doi: 10.1038/378355a0
  * Packard et al. (1980) Packard, N. H., Crutchfield, J. P., Farmer, J. D., & Shaw, R. S. 1980, Phys. Rev. Lett., 45, 712, doi: 10.1103/PhysRevLett.45.712
  * Ricker et al. (2014) Ricker, G. R., Winn, J. N., Vanderspek, R., et al. 2014, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 9143, Space Telescopes and Instrumentation 2014: Optical, Infrared, and Millimeter Wave, ed. J. Oschmann, Jacobus M., M. Clampin, G. G. Fazio, & H. A. MacEwen, 914320, doi: 10.1117/12.2063489
  * Sajadian & Fatheddin (2023) Sajadian, S., & Fatheddin, H. 2023, AJ, 166, 252, doi: 10.3847/1538-3881/ad07d9
  * Savitzky & Golay (1964) Savitzky, A., & Golay, M. J. E. 1964, Analytical Chemistry, 36, 1627, doi: 10.1021/ac60214a047
  * Stapelfeldt et al. (2021) Stapelfeldt, K., Rauer, H., Roberge, A., & Mandell, A. 2021, in ExoFrontiers, 2514-3433 (IOP Publishing), 9–1 to 9–18, doi: 10.1088/2514-3433/abfa8fch9
  * Takens (1981) Takens, F. 1981, in Dynamical Systems and Turbulence, Warwick 1980, ed. D. Rand & L.-S. Young (Berlin, Heidelberg: Springer Berlin Heidelberg), 366–381
  * Thekkeppattu et al. (2023) Thekkeppattu, J. N., Trott, C. M., & McKinley, B. 2023, MNRAS, 520, 6040, doi: 10.1093/mnras/stad522
  * Vautard & Ghil (1989) Vautard, R., & Ghil, M. 1989, Physica D: Nonlinear Phenomena, 35, 395, doi: https://doi.org/10.1016/0167-2789(89)90077-8
  * Vautard et al. (1992) Vautard, R., Yiou, P., & Ghil, M. 1992, Physica D Nonlinear Phenomena, 58, 95, doi: 10.1016/0167-2789(92)90103-T
  * Weinberg & Petersen (2021) Weinberg, M. D., & Petersen, M. S. 2021, MNRAS, 501, 5408, doi: 10.1093/mnras/staa3997
  * Willems et al. (2006) Willems, B., Kolb, U., & Justham, S. 2006, MNRAS, 367, 1103, doi: 10.1111/j.1365-2966.2006.10041.x
  * Wright & Gaudi (2013) Wright, J. T., & Gaudi, B. S. 2013, in Planets, Stars and Stellar Systems. Volume 3: Solar and Stellar Planetary Systems, ed. T. D. Oswalt, L. M. French, & P. Kalas, 489, doi: 10.1007/978-94-007-5606-9_10
  * Yuan-feng & Mei (2005) Yuan-feng, L., & Mei, Z. 2005, Applied Mathematics and Mechanics-english Edition - APPL MATH MECH-ENGL ED, 26, 179, doi: 10.1007/BF02438239
  * Zabalza et al. (2014) Zabalza, J., Ren, J., & Marshall, S. 2014, in 2014 6th Workshop on Hyperspectral Image and Signal Processing: Evolution in Remote Sensing (WHISPERS), 1–4, doi: 10.1109/WHISPERS.2014.8077583