# Precise measurement of the $e^{+}e^{-}\rightarrow
D_{s}^{\ast+}D_{s}^{\ast-}$ cross sections at center-of-mass energies from
threshold to 4.95 GeV

M. Ablikim1, M. N. Achasov5,b, P. Adlarson75, X. C. Ai81, R. Aliberti36, A.
Amoroso74A,74C, M. R. An40, Q. An71,58, Y. Bai57, O. Bakina37, I.
Balossino30A, Y. Ban47,g, V. Batozskaya1,45, K. Begzsuren33, N. Berger36, M.
Berlowski45, M. Bertani29A, D. Bettoni30A, F. Bianchi74A,74C, E.
Bianco74A,74C, A. Bortone74A,74C, I. Boyko37, R. A. Briere6, A. Brueggemann68,
H. Cai76, X. Cai1,58, A. Calcaterra29A, G. F. Cao1,63, N. Cao1,63, S. A.
Cetin62A, J. F. Chang1,58, T. T. Chang77, W. L. Chang1,63, G. R. Che44, G.
Chelkov37,a, C. Chen44, Chao Chen55, G. Chen1, H. S. Chen1,63, M. L.
Chen1,58,63, S. J. Chen43, S. M. Chen61, T. Chen1,63, X. R. Chen32,63, X. T.
Chen1,63, Y. B. Chen1,58, Y. Q. Chen35, Z. J. Chen26,h, W. S. Cheng74C, S. K.
Choi11A, X. Chu44, G. Cibinetto30A, S. C. Coen4, F. Cossio74C, J. J. Cui50, H.
L. Dai1,58, J. P. Dai79, A. Dbeyssi19, R. E. de Boer4, D. Dedovich37, Z. Y.
Deng1, A. Denig36, I. Denysenko37, M. Destefanis74A,74C, F. De Mori74A,74C, B.
Ding66,1, X. X. Ding47,g, Y. Ding41, Y. Ding35, J. Dong1,58, L. Y. Dong1,63,
M. Y. Dong1,58,63, X. Dong76, M. C. Du1, S. X. Du81, Z. H. Duan43, P.
Egorov37,a, Y. L. Fan76, J. Fang1,58, S. S. Fang1,63, W. X. Fang1, Y. Fang1,
R. Farinelli30A, L. Fava74B,74C, F. Feldbauer4, G. Felici29A, C. Q. Feng71,58,
J. H. Feng59, K Fischer69, M. Fritsch4, C. Fritzsch68, C. D. Fu1, J. L. Fu63,
Y. W. Fu1, H. Gao63, Y. N. Gao47,g, Yang Gao71,58, S. Garbolino74C, I.
Garzia30A,30B, P. T. Ge76, Z. W. Ge43, C. Geng59, E. M. Gersabeck67, A
Gilman69, K. Goetzen14, L. Gong41, W. X. Gong1,58, W. Gradl36, S.
Gramigna30A,30B, M. Greco74A,74C, M. H. Gu1,58, Y. T. Gu16, C. Y Guan1,63, Z.
L. Guan23, A. Q. Guo32,63, L. B. Guo42, M. J. Guo50, R. P. Guo49, Y. P.
Guo13,f, A. Guskov37,a, T. T. Han50, W. Y. Han40, X. Q. Hao20, F. A. Harris65,
K. K. He55, K. L. He1,63, F. H H.. Heinsius4, C. H. Heinz36, Y. K.
Heng1,58,63, C. Herold60, T. Holtmann4, P. C. Hong13,f, G. Y. Hou1,63, X. T.
Hou1,63, Y. R. Hou63, Z. L. Hou1, H. M. Hu1,63, J. F. Hu56,i, T. Hu1,58,63, Y.
Hu1, G. S. Huang71,58, K. X. Huang59, L. Q. Huang32,63, X. T. Huang50, Y. P.
Huang1, T. Hussain73, N Hüsken28,36, W. Imoehl28, M. Irshad71,58, J.
Jackson28, S. Jaeger4, S. Janchiv33, J. H. Jeong11A, Q. Ji1, Q. P. Ji20, X. B.
Ji1,63, X. L. Ji1,58, Y. Y. Ji50, X. Q. Jia50, Z. K. Jia71,58, P. C.
Jiang47,g, S. S. Jiang40, T. J. Jiang17, X. S. Jiang1,58,63, Y. Jiang63, J. B.
Jiao50, Z. Jiao24, S. Jin43, Y. Jin66, M. Q. Jing1,63, T. Johansson75, X. K.1,
S. Kabana34, N. Kalantar-Nayestanaki64, X. L. Kang10, X. S. Kang41, R.
Kappert64, M. Kavatsyuk64, B. C. Ke81, A. Khoukaz68, R. Kiuchi1, R. Kliemt14,
O. B. Kolcu62A, B. Kopf4, M. K. Kuessner4, A. Kupsc45,75, W. Kühn38, J. J.
Lane67, P. Larin19, A. Lavania27, L. Lavezzi74A,74C, T. T. Lei71,k, Z. H.
Lei71,58, H. Leithoff36, M. Lellmann36, T. Lenz36, C. Li44, C. Li48, C. H.
Li40, Cheng Li71,58, D. M. Li81, F. Li1,58, G. Li1, H. Li71,58, H. B. Li1,63,
H. J. Li20, H. N. Li56,i, Hui Li44, J. R. Li61, J. S. Li59, J. W. Li50, K. L.
Li20, Ke Li1, L. J Li1,63, L. K. Li1, Lei Li3, M. H. Li44, P. R. Li39,j,k, Q.
X. Li50, S. X. Li13, T. Li50, W. D. Li1,63, W. G. Li1, X. H. Li71,58, X. L.
Li50, Xiaoyu Li1,63, Y. G. Li47,g, Z. J. Li59, Z. X. Li16, C. Liang43, H.
Liang1,63, H. Liang35, H. Liang71,58, Y. F. Liang54, Y. T. Liang32,63, G. R.
Liao15, L. Z. Liao50, Y. P. Liao1,63, J. Libby27, A. Limphirat60, D. X.
Lin32,63, T. Lin1, B. J. Liu1, B. X. Liu76, C. Liu35, C. X. Liu1, F. H. Liu53,
Fang Liu1, Feng Liu7, G. M. Liu56,i, H. Liu39,j,k, H. B. Liu16, H. M. Liu1,63,
Huanhuan Liu1, Huihui Liu22, J. B. Liu71,58, J. L. Liu72, J. Y. Liu1,63, K.
Liu1, K. Y. Liu41, Ke Liu23, L. Liu71,58, L. C. Liu44, Lu Liu44, M. H.
Liu13,f, P. L. Liu1, Q. Liu63, S. B. Liu71,58, T. Liu13,f, W. K. Liu44, W. M.
Liu71,58, X. Liu39,j,k, Y. Liu81, Y. Liu39,j,k, Y. B. Liu44, Z. A. Liu1,58,63,
Z. Q. Liu50, X. C. Lou1,58,63, F. X. Lu59, H. J. Lu24, J. G. Lu1,58, X. L.
Lu1, Y. Lu8, Y. P. Lu1,58, Z. H. Lu1,63, C. L. Luo42, M. X. Luo80, T. Luo13,f,
X. L. Luo1,58, X. R. Lyu63, Y. F. Lyu44, F. C. Ma41, H. L. Ma1, J. L. Ma1,63,
L. L. Ma50, M. M. Ma1,63, Q. M. Ma1, R. Q. Ma1,63, R. T. Ma63, X. Y. Ma1,58,
Y. Ma47,g, Y. M. Ma32, F. E. Maas19, M. Maggiora74A,74C, S. Malde69, Q. A.
Malik73, A. Mangoni29B, Y. J. Mao47,g, Z. P. Mao1, S. Marcello74A,74C, Z. X.
Meng66, J. G. Messchendorp14,64, G. Mezzadri30A, H. Miao1,63, T. J. Min43, R.
E. Mitchell28, X. H. Mo1,58,63, N. Yu. Muchnoi5,b, Y. Nefedov37, F.
Nerling19,d, I. B. Nikolaev5,b, Z. Ning1,58, S. Nisar12,l, Y. Niu 50, S. L.
Olsen63, Q. Ouyang1,58,63, S. Pacetti29B,29C, X. Pan55, Y. Pan57, A. Pathak35,
P. Patteri29A, Y. P. Pei71,58, M. Pelizaeus4, H. P. Peng71,58, K. Peters14,d,
J. L. Ping42, R. G. Ping1,63, S. Plura36, S. Pogodin37, V. Prasad34, F. Z.
Qi1, H. Qi71,58, H. R. Qi61, M. Qi43, T. Y. Qi13,f, S. Qian1,58, W. B. Qian63,
C. F. Qiao63, J. J. Qin72, L. Q. Qin15, X. P. Qin13,f, X. S. Qin50, Z. H.
Qin1,58, J. F. Qiu1, S. Q. Qu61, C. F. Redmer36, K. J. Ren40, A. Rivetti74C,
V. Rodin64, M. Rolo74C, G. Rong1,63, Ch. Rosner19, S. N. Ruan44, N. Salone45,
A. Sarantsev37,c, Y. Schelhaas36, K. Schoenning75, M. Scodeggio30A,30B, K. Y.
Shan13,f, W. Shan25, X. Y. Shan71,58, J. F. Shangguan55, L. G. Shao1,63, M.
Shao71,58, C. P. Shen13,f, H. F. Shen1,63, W. H. Shen63, X. Y. Shen1,63, B. A.
Shi63, H. C. Shi71,58, J. L. Shi13, J. Y. Shi1, Q. Q. Shi55, R. S. Shi1,63, X.
Shi1,58, J. J. Song20, T. Z. Song59, W. M. Song35,1, Y. J. Song13, Y. X.
Song47,g, S. Sosio74A,74C, S. Spataro74A,74C, F. Stieler36, Y. J. Su63, G. B.
Sun76, G. X. Sun1, H. Sun63, H. K. Sun1, J. F. Sun20, K. Sun61, L. Sun76, S.
S. Sun1,63, T. Sun1,63, W. Y. Sun35, Y. Sun10, Y. J. Sun71,58, Y. Z. Sun1, Z.
T. Sun50, Y. X. Tan71,58, C. J. Tang54, G. Y. Tang1, J. Tang59, Y. A. Tang76,
L. Y Tao72, Q. T. Tao26,h, M. Tat69, J. X. Teng71,58, V. Thoren75, W. H.
Tian52, W. H. Tian59, Y. Tian32,63, Z. F. Tian76, I. Uman62B, S. J. Wang 50,
B. Wang1, B. L. Wang63, Bo Wang71,58, C. W. Wang43, D. Y. Wang47,g, F. Wang72,
H. J. Wang39,j,k, H. P. Wang1,63, J. P. Wang 50, K. Wang1,58, L. L. Wang1, M.
Wang50, Meng Wang1,63, S. Wang39,j,k, S. Wang13,f, T. Wang13,f, T. J. Wang44,
W. Wang59, W. Wang72, W. P. Wang71,58, X. Wang47,g, X. F. Wang39,j,k, X. J.
Wang40, X. L. Wang13,f, Y. Wang61, Y. D. Wang46, Y. F. Wang1,58,63, Y. H.
Wang48, Y. N. Wang46, Y. Q. Wang1, Yaqian Wang18,1, Yi Wang61, Z. Wang1,58, Z.
L. Wang72, Z. Y. Wang1,63, Ziyi Wang63, D. Wei70, D. H. Wei15, F. Weidner68,
S. P. Wen1, C. W. Wenzel4, U. W. Wiedner4, G. Wilkinson69, M. Wolke75, L.
Wollenberg4, C. Wu40, J. F. Wu1,63, L. H. Wu1, L. J. Wu1,63, X. Wu13,f, X. H.
Wu35, Y. Wu71, Y. J. Wu32, Z. Wu1,58, L. Xia71,58, X. M. Xian40, T. Xiang47,g,
D. Xiao39,j,k, G. Y. Xiao43, H. Xiao13,f, S. Y. Xiao1, Y. L. Xiao13,f, Z. J.
Xiao42, C. Xie43, X. H. Xie47,g, Y. Xie50, Y. G. Xie1,58, Y. H. Xie7, Z. P.
Xie71,58, T. Y. Xing1,63, C. F. Xu1,63, C. J. Xu59, G. F. Xu1, H. Y. Xu66, Q.
J. Xu17, Q. N. Xu31, W. Xu1,63, W. L. Xu66, X. P. Xu55, Y. C. Xu78, Z. P.
Xu43, Z. S. Xu63, F. Yan13,f, L. Yan13,f, W. B. Yan71,58, W. C. Yan81, X. Q.
Yan1, H. J. Yang51,e, H. L. Yang35, H. X. Yang1, Tao Yang1, Y. Yang13,f, Y. F.
Yang44, Y. X. Yang1,63, Yifan Yang1,63, Z. W. Yang39,j,k, Z. P. Yao50, M.
Ye1,58, M. H. Ye9, J. H. Yin1, Z. Y. You59, B. X. Yu1,58,63, C. X. Yu44, G.
Yu1,63, J. S. Yu26,h, T. Yu72, X. D. Yu47,g, C. Z. Yuan1,63, L. Yuan2, S. C.
Yuan1, X. Q. Yuan1, Y. Yuan1,63, Z. Y. Yuan59, C. X. Yue40, A. A. Zafar73, F.
R. Zeng50, X. Zeng13,f, Y. Zeng26,h, Y. J. Zeng1,63, X. Y. Zhai35, Y. C.
Zhai50, Y. H. Zhan59, A. Q. Zhang1,63, B. L. Zhang1,63, B. X. Zhang1, D. H.
Zhang44, G. Y. Zhang20, H. Zhang71, H. H. Zhang59, H. H. Zhang35, H. Q.
Zhang1,58,63, H. Y. Zhang1,58, J. J. Zhang52, J. L. Zhang21, J. Q. Zhang42, J.
W. Zhang1,58,63, J. X. Zhang39,j,k, J. Y. Zhang1, J. Z. Zhang1,63, Jianyu
Zhang63, Jiawei Zhang1,63, L. M. Zhang61, L. Q. Zhang59, Lei Zhang43, P.
Zhang1, Q. Y. Zhang40,81, Shuihan Zhang1,63, Shulei Zhang26,h, X. D. Zhang46,
X. M. Zhang1, X. Y. Zhang55, X. Y. Zhang50, Y. Zhang69, Y. Zhang72, Y. T.
Zhang81, Y. H. Zhang1,58, Yan Zhang71,58, Yao Zhang1, Z. H. Zhang1, Z. L.
Zhang35, Z. Y. Zhang76, Z. Y. Zhang44, G. Zhao1, J. Zhao40, J. Y. Zhao1,63, J.
Z. Zhao1,58, Lei Zhao71,58, Ling Zhao1, M. G. Zhao44, S. J. Zhao81, Y. B.
Zhao1,58, Y. X. Zhao32,63, Z. G. Zhao71,58, A. Zhemchugov37,a, B. Zheng72, J.
P. Zheng1,58, W. J. Zheng1,63, Y. H. Zheng63, B. Zhong42, X. Zhong59, H.
Zhou50, L. P. Zhou1,63, X. Zhou76, X. K. Zhou7, X. R. Zhou71,58, X. Y. Zhou40,
Y. Z. Zhou13,f, J. Zhu44, K. Zhu1, K. J. Zhu1,58,63, L. Zhu35, L. X. Zhu63, S.
H. Zhu70, S. Q. Zhu43, T. J. Zhu13,f, W. J. Zhu13,f, Y. C. Zhu71,58, Z. A.
Zhu1,63, J. H. Zou1, J. Zu71,58
(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 Budker Institute of Nuclear Physics SB RAS (BINP), Novosibirsk 630090,
Russia
6 Carnegie Mellon University, Pittsburgh, Pennsylvania 15213, USA
7 Central China Normal University, Wuhan 430079, People’s Republic of China
8 Central South University, Changsha 410083, People’s Republic of China
9 China Center of Advanced Science and Technology, Beijing 100190, People’s
Republic of China
10 China University of Geosciences, Wuhan 430074, People’s Republic of China
11 Chung-Ang University, Seoul, 06974, Republic of Korea
12 COMSATS University Islamabad, Lahore Campus, Defence Road, Off Raiwind
Road, 54000 Lahore, Pakistan
13 Fudan University, Shanghai 200433, People’s Republic of China
14 GSI Helmholtzcentre for Heavy Ion Research GmbH, D-64291 Darmstadt, Germany
15 Guangxi Normal University, Guilin 541004, People’s Republic of China
16 Guangxi University, Nanning 530004, People’s Republic of China
17 Hangzhou Normal University, Hangzhou 310036, People’s Republic of China
18 Hebei University, Baoding 071002, People’s Republic of China
19 Helmholtz Institute Mainz, Staudinger Weg 18, D-55099 Mainz, Germany
20 Henan Normal University, Xinxiang 453007, People’s Republic of China
21 Henan University, Kaifeng 475004, People’s Republic of China
22 Henan University of Science and Technology, Luoyang 471003, People’s
Republic of China
23 Henan University of Technology, Zhengzhou 450001, People’s Republic of
China
24 Huangshan College, Huangshan 245000, People’s Republic of China
25 Hunan Normal University, Changsha 410081, People’s Republic of China
26 Hunan University, Changsha 410082, People’s Republic of China
27 Indian Institute of Technology Madras, Chennai 600036, India
28 Indiana University, Bloomington, Indiana 47405, USA
29 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
30 INFN Sezione di Ferrara, (A)INFN Sezione di Ferrara, I-44122, Ferrara,
Italy; (B)University of Ferrara, I-44122, Ferrara, Italy
31 Inner Mongolia University, Hohhot 010021, People’s Republic of China
32 Institute of Modern Physics, Lanzhou 730000, People’s Republic of China
33 Institute of Physics and Technology, Peace Avenue 54B, Ulaanbaatar 13330,
Mongolia
34 Instituto de Alta Investigación, Universidad de Tarapacá, Casilla 7D, Arica
1000000, Chile
35 Jilin University, Changchun 130012, People’s Republic of China
36 Johannes Gutenberg University of Mainz, Johann-Joachim-Becher-Weg 45,
D-55099 Mainz, Germany
37 Joint Institute for Nuclear Research, 141980 Dubna, Moscow region, Russia
38 Justus-Liebig-Universitaet Giessen, II. Physikalisches Institut, Heinrich-
Buff-Ring 16, D-35392 Giessen, Germany
39 Lanzhou University, Lanzhou 730000, People’s Republic of China
40 Liaoning Normal University, Dalian 116029, People’s Republic of China
41 Liaoning University, Shenyang 110036, People’s Republic of China
42 Nanjing Normal University, Nanjing 210023, People’s Republic of China
43 Nanjing University, Nanjing 210093, People’s Republic of China
44 Nankai University, Tianjin 300071, People’s Republic of China
45 National Centre for Nuclear Research, Warsaw 02-093, Poland
46 North China Electric Power University, Beijing 102206, People’s Republic of
China
47 Peking University, Beijing 100871, People’s Republic of China
48 Qufu Normal University, Qufu 273165, People’s Republic of China
49 Shandong Normal University, Jinan 250014, People’s Republic of China
50 Shandong University, Jinan 250100, People’s Republic of China
51 Shanghai Jiao Tong University, Shanghai 200240, People’s Republic of China
52 Shanxi Normal University, Linfen 041004, People’s Republic of China
53 Shanxi University, Taiyuan 030006, People’s Republic of China
54 Sichuan University, Chengdu 610064, People’s Republic of China
55 Soochow University, Suzhou 215006, People’s Republic of China
56 South China Normal University, Guangzhou 510006, People’s Republic of China
57 Southeast University, Nanjing 211100, People’s Republic of China
58 State Key Laboratory of Particle Detection and Electronics, Beijing 100049,
Hefei 230026, People’s Republic of China
59 Sun Yat-Sen University, Guangzhou 510275, People’s Republic of China
60 Suranaree University of Technology, University Avenue 111, Nakhon
Ratchasima 30000, Thailand
61 Tsinghua University, Beijing 100084, People’s Republic of China
62 Turkish Accelerator Center Particle Factory Group, (A)Istinye University,
34010, Istanbul, Turkey; (B)Near East University, Nicosia, North Cyprus,
99138, Mersin 10, Turkey
63 University of Chinese Academy of Sciences, Beijing 100049, People’s
Republic of China
64 University of Groningen, NL-9747 AA Groningen, The Netherlands
65 University of Hawaii, Honolulu, Hawaii 96822, USA
66 University of Jinan, Jinan 250022, People’s Republic of China
67 University of Manchester, Oxford Road, Manchester, M13 9PL, United Kingdom
68 University of Muenster, Wilhelm-Klemm-Strasse 9, 48149 Muenster, Germany
69 University of Oxford, Keble Road, Oxford OX13RH, United Kingdom
70 University of Science and Technology Liaoning, Anshan 114051, People’s
Republic of China
71 University of Science and Technology of China, Hefei 230026, People’s
Republic of China
72 University of South China, Hengyang 421001, People’s Republic of China
73 University of the Punjab, Lahore-54590, Pakistan
74 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
75 Uppsala University, Box 516, SE-75120 Uppsala, Sweden
76 Wuhan University, Wuhan 430072, People’s Republic of China
77 Xinyang Normal University, Xinyang 464000, People’s Republic of China
78 Yantai University, Yantai 264005, People’s Republic of China
79 Yunnan University, Kunming 650500, People’s Republic of China
80 Zhejiang University, Hangzhou 310027, People’s Republic of China
81 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 75270,
Pakistan

###### Abstract

The process $e^{+}e^{-}\rightarrow D_{s}^{\ast+}D_{s}^{\ast-}$ is studied with
a semi-inclusive method using data samples at center-of-mass energies from
threshold to 4.95 GeV collected with the BESIII detector operating at the
Beijing Electron Positron Collider. The Born cross sections of the process are
measured for the first time with high precision in this energy region. Two
resonance structures are observed in the energy-dependent cross sections
around 4.2 and 4.4 GeV. By fitting the cross sections with a coherent sum of
three Breit-Wigner amplitudes and one phase-space amplitude, the two
significant structures are assigned masses of (4186.5$\pm$9.0$\pm$30)
MeV/$c^{2}$ and (4414.5$\pm$3.2$\pm$6.0) MeV/$c^{2}$, widths of
(55$\pm$17$\pm$53) MeV and (122.6$\pm$7.0$\pm$8.2) MeV, where the first errors
are statistical and the second ones are systematic. The inclusion of a third
Breit-Wigner amplitude is necessary to describe a structure around 4.79 GeV.

###### pacs:

Valid PACS appear here

††preprint: APS/123-QED

In $e^{+}e^{-}$ annihilations, several conventional vector charmonium states
are established in the inclusive hadronic cross sections above the open charm
threshold, such as the $\psi(3770)$, $\psi(4040)$, $\psi(4160)$, and
$\psi(4415)$. However, unexpected vector charmonium-like resonance structures
have also been discovered over the past two decades with a charmonium and
light hadrons in the final state. These include $\psi(4260)$, initially
observed in the $e^{+}e^{-}\to$$\pi^{+}\pi^{-}J/\psi$ process [1, 2], and
$\psi(4360)$ and $\psi(4660)$, observed in
$e^{+}e^{-}\to\pi^{+}\pi^{-}\psi(3686)$ [3, 4] by the $B$-factory experiments
BABAR and Belle.

The $\psi(4260)$ state has been searched for in several alternative decay
modes, involving charmonium states (such as $K\bar{K}J/\psi$ [5, 6], $\eta
J/\psi$ [7, 8], $\eta^{\prime}J/\psi$ [9], $\pi\pi h_{c}$ [10, 11],
$\pi^{+}\pi^{-}\psi(3686)$, and $\omega\chi_{cJ}$ [12]) and charmed mesons
(such as $D\bar{D}$ [13, 14], $D^{\ast 0}D^{\ast-}\pi^{+}+c.c.$ [15] and
$\bar{D}D^{*}\pi+c.c.$ [16]). A structure around 4.26 GeV is observed in most
of the modes except in the $D\bar{D}$ mode [13, 14].

In 2017, a dedicated analysis of $e^{+}e^{-}\to\pi^{+}\pi^{-}J/\psi$ performed
by the BESIII experiment revealed that the so-called $\psi(4260)$ state is
actually a combination of two states [17]. The main component has a lower mass
and a much narrower width than the initially measured state, and is consistent
with the structure observed in $e^{+}e^{-}\to\pi^{+}\pi^{-}h_{c}$ [18] and
$\omega\chi_{c0}$ [12], denoted as $\psi(4230)$. Futhermore, the BESIII
experiment also discovered another structure around
$4.32~{}{\rm~{}GeV}/c^{2}$, with a significance greater than $7.6\sigma$ [17].
Interestingly, the mass of the $\psi(4230)$ state lies just at the production
threshold of the $D_{s}^{\ast+}D_{s}^{\ast-}$ pair, which suggests a possible
correlation between the $\psi(4230)$ and this open-charm decay mode.

The measurements of the exclusive cross sections for charmed strange meson
pairs were performed by CLEO-c at center-of-mass energies ($E_{\textrm{CM}}$)
up to 4.26 GeV [19], and by BABAR [20] and Belle [21] with the initial-state
radiation (ISR) method. Since these results are limited by either energy range
(CLEO-c) or statistics (BABAR and Belle), it could not be concluded yet
whether these vector charmonium-like states decay into charmed strange meson
pairs or not.

With the data samples taken by the BESIII experiment, the cross sections of
$e^{+}e^{-}\to D_{s}^{\ast+}D_{s}^{\ast-}$ can be measured from the production
threshold up to 4.95 GeV, with much improved precision over previous
experiments. This measurement allows to investigate the vector-
chamonium(-like) structures, especially the $\psi(4230)$, therefore shedding
light on their nature.

The BESIII detector is described in detail in [22, 23]. The experimental data
used in this analysis were taken at $E_{\mathrm{CM}}$ ranging from $4.226$ GeV
(just above the $D_{s}^{\ast+}D_{s}^{\ast-}$ production threshold) to $4.95$
GeV with 76 energy points [24, 26, 25] corresponding in total to an integrated
luminosity of $15.67~{}\textrm{fb}^{-1}$ [27, 28, 26]. A GEANT4-based [29]
Monte Carlo (MC) package, which includes the geometric description of the
BESIII detector and its response, is used to produce simulated samples which
are used to estimate the backgrounds and to determine the detection
efficiencies and ISR corrections.

The $e^{+}e^{-}\to D_{s}^{\ast+}D_{s}^{\ast-}$ events are generated using
helicity-amplitude (HELAMP) models in EvtGen [30, 31] at all the energy
points, where the HELAMP inputs (relative magnitudes of the helicity
amplitudes) are extracted from the helicity angular distributions of data, the
width of $D_{s}^{\ast\pm}$ is set to zero, and beam energy spread and ISR are
considered with the generator KKMC [32, 33]. Possible background contributions
are estimated with inclusive MC samples generated by KKMC with integrated
luminosities comparable to data, in which the known decay modes are modeled
with EvtGen using branching fractions taken from the Particle Data Group [34],
and the remaining unknown charmonium decays are modeled with lundcharm [35].
Final-state radiation (FSR) from charged particles is incorporated by the
photos package [36].

To increase efficiency, a semi-inclusive method is performed by reconstructing
only $D_{s}^{\ast+}$ or $D_{s}^{\ast-}$ of $e^{+}e^{-}\to
D_{s}^{\ast+}D_{s}^{\ast-}$, with $D_{s}^{\ast\pm}\to\gamma
D_{s}^{\pm}\to\gamma K^{+}K^{-}\pi^{\pm}$. With this method at least one good
photon, one pair of $K^{+}$$K^{-}$, and one $\pi^{\pm}$ are required in the
final state. The selection criteria for charged tracks and photon candidates
are described in Ref. [37]. To select $D_{s}^{\pm}$ candidates, the invariant
mass $m_{KK\pi}$ of $K^{+}K^{-}\pi^{\pm}$ is required to be in the range
$|m_{KK\pi}-m_{D_{s}}|<15$ MeV/$c^{2}$, where $m_{D_{s}}$ is the nominal
$D_{s}^{\pm}$ mass [34]. All the $\gamma D_{s}^{\pm}$ (i.e. $\gamma
K^{+}K^{-}\pi^{\pm}$) combinations are taken as $D_{s}^{\ast\pm}$ candidates
for further selection. The resolution of the missing mass for
$D_{s}^{\ast\pm}$ candidates can be improved by incorporating
$M_{\mathrm{miss}}\equiv m_{\mathrm{miss}}+m_{\gamma KK\pi}-m_{D_{s}^{\ast}}$,
where $m_{\mathrm{miss}}$ and $m_{\gamma KK\pi}$ are the missing mass and the
invariant mass of the $D_{s}^{\ast\pm}\to\gamma K^{+}K^{-}\pi^{\pm}$
candidates, respectively, and $m_{D_{s}^{\ast}}$ is the nominal mass of
$D_{s}^{\ast\pm}$ [34]. With the signal MC events, the resolution of
$M_{\mathrm{miss}}$ is estimated and parameterized as an energy dependent
function $\sigma^{\mathrm{MC}}_{M_{\mathrm{miss}}}(E_{\mathrm{CM}})$. To
select $D_{s}^{\ast\pm}$ from $e^{+}e^{-}\to D_{s}^{\ast+}D_{s}^{\ast-}$
events and suppress background, the modified missing mass of $D_{s}^{\ast\pm}$
candidates is required to be in a window
$|M_{\mathrm{miss}}-m_{D_{s}^{\ast}}|<5\sigma^{\mathrm{MC}}_{M_{\mathrm{miss}}}(E_{\mathrm{CM}})$.
To improve the mass resolution, $M_{\gamma KK\pi}\equiv m_{\gamma
KK\pi}-m_{KK\pi}+m_{D_{s}}$ is calculated for $D_{s}^{\ast\pm}$ candidates.

The yields of $D_{s}^{\ast\pm}$ signals are determined by fitting the
$M_{\gamma KK\pi}$ distributions. By matching the MC truth, correct $\gamma
KK\pi$ combinations from $e^{+}e^{-}\to D_{s}^{\ast+}D_{s}^{\ast-}$ MC can be
figured out, and the others from MC are called random combinations. To
describe the $D_{s}^{\ast\pm}$ signal in data, the $M_{\gamma KK\pi}$ shape of
correct $\gamma KK\pi$ combinations from $e^{+}e^{-}\to
D_{s}^{\ast+}D_{s}^{\ast-}$ MC is used. This MC signal shape is convolved with
a Gaussian function to take into account the possible mass shift $\Delta m$,
and the resolution (and $D_{s}^{\ast}$ width) difference $\Delta\sigma$
between data and MC simulation. The $M_{\gamma KK\pi}$ shape of random $\gamma
KK\pi$ combinations from $e^{+}e^{-}\to D_{s}^{\ast+}D_{s}^{\ast-}$ MC is also
used as one component and its ratio relative to the one for correct
combinations is fixed according to the MC study. To describe the background in
$M_{\gamma KK\pi}$ distribution for data, a second order Chebyshev function
$1+c_{0}x+c_{1}(2x^{2}-1)$ with two coefficients $c_{0}$ and $c_{1}$ is used.
The parameters $\Delta m$, $\Delta\sigma$, $c_{0}$, and $c_{1}$ are floating
initially in the fits. Since there is no $E_{\mathrm{CM}}$ dependence observed
for $\Delta m$, $\Delta\sigma$, and $c_{1}$, they are finally fixed to the
values averaged over $E_{\mathrm{CM}}$ in the fits to determine the nominal
$D_{s}^{\ast\pm}$ signal yields, leaving only $c_{0}$ as floating. As an
example, the fit of $M_{\gamma KK\pi}$ to data at $E_{\mathrm{CM}}=4.29$ GeV
111This is a nominal value and the measured one is $4287.88\pm 0.06\pm 0.34$
MeV [25]. is shown in Fig. 1.

Figure 1: The $M_{\gamma KK\pi}$ distribution for data at 4.29 GeV. The black
curve represents the fit, the red dashed curve the background, the blue dot-
dashed curve the random combinations of $\gamma KK\pi$ from $e^{+}e^{-}\to
D_{s}^{\ast\pm}D_{s}^{\ast\mp}$, and the blue dashed curve the correct $\gamma
KK\pi$ combinations from the signal.

A study of the inclusive MC samples shows that no peaking background is found
from events without $D_{s}^{\ast\pm}$ after applying all the selection
criteria. The ISR-produced $D_{s}^{\pm}D_{s}^{\ast\mp}$ events can contribute
as a peaking background, which is subtracted by normalized MC according to the
cross sections of $e^{+}e^{-}\to D_{s}^{\pm}D_{s}^{\ast\mp}$ 222A measurement
by BESIII to be published in near future. and the luminosities. The other two-
body processes containing $D_{s}^{\ast\pm}$ are $e^{+}e^{-}\to
D_{s}^{\ast\pm}D_{sJ}^{(\ast)\mp}$, but they do not contribute to the peaking
background since their missing mass (i.e. the mass of $D_{sJ}^{(\ast)\mp}$) is
outside of the $D_{s}^{\ast\mp}$ mass region. For the same reason, the three-
body process $e^{+}e^{-}\to D_{s}^{\ast}D^{(\ast)}K$ does not contribute as
peaking backgrounds either. The contribution from the three-body process
$e^{+}e^{-}\to D_{s}^{\ast\pm}D_{s}^{\mp}\pi^{0}$ is expected to be negligible
due to isospin violation.

The Born cross sections $\sigma_{\mathrm{Born}}$ and the dressed cross
sections $\sigma_{\mathrm{dressed}}$ at each energy point are calculated using

$\displaystyle\sigma_{\mathrm{Born}}$
$\displaystyle=\sigma_{\mathrm{dressed}}|1+\Pi|^{2}$ (1)
$\displaystyle=\frac{N_{D_{s}^{\ast}}^{\mathrm{fit}}-N_{D_{s}^{\pm}D_{s}^{\ast\mp}}}{2\mathcal{B}(D_{s}^{\pm}\to
K^{+}K^{-}\pi^{\pm})\epsilon(1+\delta)\frac{1}{|1+\Pi|^{2}}\mathcal{L}_{\mathrm{int}}},$

where $N_{D_{s}^{\ast}}^{\mathrm{fit}}$ is the fitted $D_{s}^{\ast}$ signal
yield, $N_{D_{s}^{\pm}D_{s}^{\ast\mp}}$ is the number of the estimated peaking
background events from the ISR produced $D_{s}^{\pm}D_{s}^{\ast\mp}$,
$\mathcal{L}_{\mathrm{int}}$ is the integrated luminosity, $1+\delta$ is the
ISR correction factor, $\frac{1}{|1+\Pi|^{2}}$ is the correction factor from
the vacuum polarization (VP) [38], $\mathcal{B}(D_{s}^{\pm}\to
K^{+}K^{-}\pi^{\pm})$ is the branching fraction of the decay $D_{s}^{\pm}\to
K^{+}K^{-}\pi^{\pm}$, which is taken from Ref. [34], and $\epsilon$ is the
detection efficiency containing the branching fraction of
$D_{s}^{\ast\pm}\to\gamma D_{s}^{\pm}$ [34]. As $D_{s}^{\ast+}$ and
$D_{s}^{\ast-}$ are not separated in the fitted signal yields, there is a
factor of 2 in the denominator.

The ISR effect also impacts on the detection efficiency in addition to the
correction factor. Since ISR depends on the dressed cross section line shape,
an iterative procedure is used to determine the dressed cross sections [39].
First, signal MC samples are generated at all the energy points with a flat
dressed cross section line shape using KKMC to determine the initial detection
efficiencies, ISR correction factors and subsequently the preliminary measured
dressed cross sections. The measured dressed cross sections are fitted by a
smooth line shape (defined later in Eq. (2)) as a function of
$E_{\textrm{CM}}$. With the fitted line shape, a MC-weighting method [39] is
used to iteratively update the efficiencies, the ISR correction factors, the
measured dressed cross sections, and the fitted line shape. After four
iterations, the fit to the measured dressed cross section converged. By
including the VP correction, the Born cross sections are summarized in the
Supplemental Material [40].

The Born cross sections of $e^{+}e^{-}\to D_{s}^{\ast+}D_{s}^{\ast-}$ are
fitted to estimate the parameters of the three structures with a coherent sum
of three Breit-Wigner (BW) and one two-body phase-space (PHSP) amplitudes
assuming $D_{s}^{\ast+}D_{s}^{\ast-}$ in
P-wave333$e^{+}e^{-}\to\gamma^{\ast}/Y\to D_{s}^{\ast+}D_{s}^{\ast-}$ is a
vector$\to$vector+vector process, $P$-parity conservation requires the orbit
angular momentum number $L$ of the $D_{s}^{\ast+}D_{s}^{\ast-}$ system to be
odd and $L=1$ (P-wave) is the most likely case.

$\displaystyle\left|\sum_{j=1}^{3}a_{j}\cdot e^{i\phi_{j}}\cdot
BW_{j}(E_{\mathrm{CM}})+\frac{a_{0}}{E_{\mathrm{CM}}^{n}}\right|^{2}\beta^{3}(E_{\mathrm{CM}}),$
(2)

where
$BW_{j}(E_{\mathrm{CM}})=\frac{B_{1}(E_{\mathrm{CM}})}{(E^{2}_{\mathrm{CM}}-M_{j}^{2}+iM_{j}\Gamma_{j})\cdot{E_{\mathrm{CM}}}}$
is the $j$-th BW amplitude, the Blatt-Weisskopf function
$B_{1}(E_{\mathrm{CM}})=\sqrt{1/(1+q^{2}(E_{\mathrm{CM}})R^{2})}$ is used as
the P-wave barrier factor [41, 42] with $q(E_{\mathrm{CM}})$ the momentum of
$D_{s}^{\ast\pm}$ and $R=1.6~{}\mathrm{GeV}^{-1}$ [43] the barrier radius,
$a_{j(0)}$ is the coefficient for the BW (the PHSP amplitude),
$\beta(E_{\mathrm{CM}})$ is the velocity of $D_{s}^{\ast\pm}$, and $\phi_{j}$
is the relative phase. In the nominal fit of the Born cross sections, only
statistical uncertainties are considered. We find three sets of fitting
results with comparable goodness of fit, which are shown in Fig. 2 and Table
1. Due to the limited number of data points around 4.79 GeV, the fitted mass
of the third structure varies from 4786 to 4793.3 MeV/$c^{2}$ and the width
from 27.1 to 60 MeV. The statistical significance of the third structure
exceeds 6.1$\sigma$ in all three fitting procedures, indicating that
incorporating this amplitude leads to a more accurate description of the cross
sections compared to using only two BWs. Besides the three BWs, one additional
BW is tried to better describe the data points around 4.55 GeV, but it is not
kept due to its small statistical significance.

(a) Result 1

(b) Result 2

(c) Result 3

Figure 2:  Three fitting results for the measured Born cross sections of $e^{+}e^{-}\rightarrow D_{s}^{\ast+}D_{s}^{\ast-}$. The black dots with error bars are for the measured Born cross sections. In each plot, the black curve represents the fit; the green dashed, blue two-dashed and red long-dashed ones are for the three BWs, respectively, and the pink dot-dashed is for the PHSP contributions.  Table 1: The fitting results of the Born cross sections. | Result 1 | Result 2 | Result 3
---|---|---|---
$M_{1}$ (MeV/$c^{2}$) | 4186.5$\pm$9.0 | 4193.8$\pm$7.5 | 4195.3$\pm$7.5
$\Gamma_{1}$ (MeV) | 55$\pm$17 | 61.2$\pm$9.0 | 61.8$\pm$9.0
$M_{2}$ (MeV/$c^{2}$) | 4414.5$\pm$3.2 | 4412.8$\pm$3.2 | 4411.0$\pm$3.2
$\Gamma_{2}$ (MeV) | 122.6$\pm$7.0 | 120.3$\pm$7.0 | 120.0$\pm$7.0
$M_{3}$ (MeV/$c^{2}$) | 4793.3$\pm$7.5 | 4789.8$\pm$9.0 | 4786$\pm$10
$\Gamma_{3}$ (MeV) | 27.1$\pm$7.0 | 41$\pm$39 | 60$\pm$35

The systematic uncertainties for the measured Born cross sections are
described as follows, some of which are common along all the energy points,
while others are estimated depending on the $E_{\rm{CM}}$ range.

The systematic uncertainties of tracking (particle identification) efficiency
are estimated to be 0.5% (0.5%) per $K^{\pm}$ and 0.2% (0.4%) per $\pi^{\pm}$
with a control sample of $D_{s}^{\pm}\to K^{+}K^{-}\pi^{\pm}$ decays [44];
thus a 1.2% (1.4%) systematic uncertainty is assigned for the tracking
efficiency (particle identification) in the $D_{s}^{\pm}\to
K^{+}K^{-}\pi^{\pm}$ candidates selection. The systematic uncertainty in the
efficiency for photon reconstruction is set conservatively to 1% based on a
study with a sample of $J/\psi\to\rho\pi$ events [45]. From fits of the
invariant mass spectrum of $D_{s}^{\pm}\to K^{+}K^{-}\pi^{\pm}$ candidates and
of the modified missing mass of $D_{s}^{\ast\pm}\to\gamma K^{+}K^{-}\pi^{\pm}$
candidates, the efficiencies for signal in the mass window for both data and
MC samples can be calculated, and the relative differences in efficiency
between data and MC simulation are taken as systematic uncertainties for the
$D_{s}^{\pm}$ mass window and the modified missing mass window, respectively,
which cover the possible resolution difference and the zero width setting for
$D_{s}^{\ast\pm}$ in MC simulation. The maximum difference in the dressed
cross sections between the last two iterations, $0.2\%$, is taken as the
systematic uncertainty for the stability of the iteration results. The
systematic uncertainty of the branching fraction for $D_{s}^{\pm}\to
K^{+}K^{-}\pi^{\pm}$ and $D_{s}^{\ast\pm}\to\gamma D_{s}^{\pm}$ is taken from
Ref. [34]. The integrated luminosities are measured by QED events [46] with a
systematic uncertainty of 1%. The uncertainty from VP correction is 0.1% [38].
The uncertainty of the peaking background subtraction from the $e^{+}e^{-}\to
D_{s}^{\pm}D_{s}^{\ast\mp}$ process is estimated to be $1\%$ mainly from the
uncertainties of cross section measurement. Instead of the HELAMP model, the
PHSP model [30, 31] is also used to generate $e^{+}e^{-}\to
D_{s}^{\ast\pm}D_{s}^{\ast\mp}$ events and the maximum difference in
efficiency, 2.2%, is taken as the systematic uncertainty for the generation
model. We get the total common systematic uncertainty by adding them in
quadrature, which is 4.0%.

Table 2: Relative systematic uncertainties of the measured Born cross sections. As the effects and the data statistics are energy dependent, some systematic uncertainties are estimated in several energy intervals combining some data sets with low luminosities.  Sources | Systematic uncertainties (%)
---|---
Common | 4.0
Energy (interval) (GeV) | 4.226 | 4.228 | 4.233 | 4.233$\sim$4.24 | 4.24$\sim$4.3 | 4.3$\sim$4.4 | 4.4$\sim$4.82 | 4.843 | 4.86$\sim$4.95
Signal yields fitting with fixed parameters | 2.0 | 5.0 | 2.0
Signal yields fitting range | 4.0 | 5.0 | 4.0
Cross section line shape description | 5.0 | 2.0 | 5.0 | 2.0 | 5.0 | 2.0
$E_{\textrm{CM}}$ uncertainty | 24.5 | 20.0 | 9.4 | 0.8
Total | 25.7 | 21.5 | 12.2 | 7.9 | 6.4 | 9.6 | 6.4 | 7.9 | 6.4
Table 3: The systematic uncertainties of the fitted masses and widths. The systematic uncertainties of the fitted masses and widths. Fitting represent the uncertainty from the multiple fitting results. $R$ represent the uncertainty from the barrier radius, $E_{\textrm{CM}}$ the systematic errors from the $E_{\mathrm{CM}}$ uncertainty and the $\sigma_{\textrm{Born}}$ the systematic uncertainties on the measured Born cross sections. Sources | Fitting | $R$ | $E_{\textrm{CM}}$ | $\sigma_{\textrm{Born}}$ | Total
---|---|---|---|---|---
$M_{1}$ (MeV/$c^{2}$) | 8.8 | 2.6 | 28.6 | 4.8 | 30
$\Gamma_{1}$ (MeV) | 6.8 | 1.9 | 51 | 11.8 | 53
$M_{2}$ (MeV/$c^{2}$) | 3.5 | 0.7 | 3.9 | 2.9 | 6.0
$\Gamma_{2}$ (MeV) | 2.6 | 0.3 | 7.7 | 0.9 | 8.2
$M_{3}$ (MeV/$c^{2}$) | 7.3 | 1.0 | 2.4 | 5.1 | 9.3
$\Gamma_{3}$ (MeV) | 32.9 | 1.1 | 5.3 | 3.4 | 34

The shape related parameters $\Delta m$, $\Delta\sigma$, and $c_{1}$ are fixed
to the averaged values in the nominal fit (see Fig. 1 for an example). The
differences of the fitted signal yields, when these parameters are floating,
are within statistical uncertainties. However, to cover these differences
conservatively, the whole energy range is divided into three intervals
$(4.226,4.3)$ GeV, $(4.3,4.4)$ GeV, and $(4.4,4.95)$ GeV with assigned
systematic uncertainties 2%, 5%, and 2%, respectively, due to the signal and
background shapes. The boundaries of the nominal fitting range for $M_{\gamma
KK\pi}$, which is $[2.02,~{}2.20]$ GeV$/c^{2}$, are changed by $10$ MeV to
estimate the corresponding fitting range uncertainties. These are assigned to
be 4%, 5%, and 4%, respectively, in the energy intervals $(4.226,4.3)$ GeV,
$(4.3,4.4)$ GeV, and $(4.4,4.95)$ GeV. Eq. (2) is used to fit the data
iteratively and get the converged dressed cross sections, during that the
cross section line shape is similar to the one in Fig. 2(a). Other two
different line shapes with comparable fitting goodness, that are similar to
the results shown in Fig. 2(b) and 2(c), plus an additional line shape
obtained by the LOWESS (LOcally WEighted Scatterplot Smoothing) [48, 49], are
used to calculate the systematic uncertainties for the line shape description
by repeating the iterations and taking the differences in the results. The
systematic uncertainty of the measured $E_{\mathrm{CM}}$ is found to be less
than 0.8 MeV (0.6 MeV) for $4.226<E_{\mathrm{CM}}<4.6$ GeV
($4.6<E_{\mathrm{CM}}<4.95$ GeV) [24, 26, 25] and it is used to shift all the
energy points to conservatively estimate the impacts on the measured cross
sections. Since the cross section line shape varies dramatically near the
threshold, the $E_{\mathrm{CM}}$ uncertainty could have significant impact on
the ISR correction factors nearby, subsequently affecting the measured cross
sections and the fitting results around the threshold. The MC samples at the
first four $E_{\mathrm{CM}}$ points near the threshold are regenerated with
the shifted energies to update ISR corrections. After the iterations with the
shifted energy points and the updated ISR correction factors near threshold,
the resulting differences in the Born cross sections are taken as the
systematic uncertainties due to the $E_{\mathrm{CM}}$ uncertainty. These
systematic uncertainties are summarized in Table 2 by assuming no correlation
among different sources.

The sum of three BW amplitudes with the one for the phase-space is an
imperfect model, but it is useful to describe the cross section line shape
smoothly for the iterative ISR correction procedure and to obtain masses and
widths of the resonance(-like) structures for reference. The systematic
uncertainties of the fitted masses and widths of the three resonances are
listed in Table 3. There are three fitting results as listed in Table 1 with
comparable goodness of fit. The fitting result 1 is taken as the nominal one,
and the biggest differences in the fitted central values between result 1 and
the other two are taken as the systematic uncertainties from the multiple
fitting results. The parameter $R$ is changed from 1.6 in the nominal fit to 5
GeV-1 (corresponds the scale of the strong interaction which is about 1 fm)
[43, 47], and the variations in the fitting results are taken as the
corresponding systematic uncertainties. The uncertainty from the measured
$E_{\mathrm{CM}}$ has impact not only on the cross sections but also on the
fitting results. After the iterations with the shifted energies and the
updated ISR correction factors, the differences in the fitting results are
considered as the uncertainty from the $E_{\mathrm{CM}}$ uncertainty. The
systematic uncertainties on the measured Born cross sections consist of a
common part and an uncommon part. The common part is the same for all energies
and has no effect on fitted masses and widths of the BW functions. The
uncommon part is included to redo the cross section fitting and the resulting
differences are taken as systematic uncertainties.

In summary, with the world’s largest $e^{+}e^{-}$ scan data sample between
$4.226$ and $4.95$ GeV accumulated by BESIII, the Born cross sections of
$e^{+}e^{-}\to D_{s}^{\ast+}D_{s}^{\ast-}$ are measured precisely. Two
enhancements in the $E_{\rm CM}$ dependent cross sections are observed around
4.2 and 4.45 GeV. The $E_{\rm CM}$ dependent cross section line shape is
modeled with a sum of three BW and one PHSP amplitudes. The fitted mass and
width for the first resonance are (4186.5$\pm$9.0$\pm$30) MeV/$c^{2}$ and
(55$\pm$17$\pm$53) MeV, respectively. These results are consistent with the
$\psi(4160)$ observed in the inclusive cross section of
$e^{+}e^{-}\to\textrm{hadrons}$ [50] and in the dimuon spectrum of $B^{-}\to
K^{-}\mu^{+}\mu^{-}$ [51]. While considering the systematic uncertainties,
these results are also consistent with $\psi(4230)$ observed in the
$\pi^{+}\pi^{-}J/\psi$ mode. If the contribution is from $\psi(4230)$, it
indicates that the $\psi(4230)$ couples more strongly to the
$D_{s}^{\ast+}D_{s}^{\ast-}$ mode than to the modes with charmonium states,
since the cross section of $e^{+}e^{-}\to D_{s}^{\ast+}D_{s}^{\ast-}$ at 4.23
GeV is roughly one order of magnitude higher than that of
$e^{+}e^{-}\to\pi^{+}\pi^{-}J/\psi$. This information is vital to understand
the nature of $\psi(4230)$, especially considering the fact that its mass is
very close to the production threshold of $D_{s}^{\ast+}D_{s}^{\ast-}$. The
fitted mass and width for the second resonance are (4414.5$\pm$3.2$\pm$6.0)
MeV/$c^{2}$ and (122.6$\pm$7.0$\pm$8.2) MeV, respectively. The mass is
consistent with $\psi(4415)$, and the measured width from this work is a bit
higher than the world averaged value of $\psi(4415)$ [34], but still within
three standard deviations. If we assume that this resonance is the
$\psi(4415)$, this is the first time that this state is observed in the
$D_{s}^{\ast+}D_{s}^{\ast-}$ decay mode. An additional third BW amplitude
describes the $E_{\rm CM}$ dependent cross sections around 4.79 GeV better
than just using two BW functions, with statistical significance greater than
$6.1\sigma$, however, more data points in the vicinity are needed for further
clarification [52]. It is out of the scope of this paper, but a unitary
approach based on $K$-matrix formalism to fit the cross-section results of
various exclusive channels is expected to be carried out as a more
comprehensive and robust analysis of the vector-charmonium(-like) structures
(for instance, an analysis for vector bottomonia with this method can be found
in Ref. [53]).

###### 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.
12150004, 11635010, 11735014, 11835012, 11935015, 11935016, 11935018,
11961141012, 12022510, 12025502, 12035009, 12035013, 12061131003, 12192260,
12192261, 12192262, 12192263, 12192264, 12192265, 12221005, 12225509,
12235017, U1732105, U1632106; Program of Science and Technology Development
Plan of Jilin Province of China under Contract No. 20210508047RQ; The
Fundamental Research Funds of Shandong University; the Chinese Academy of
Sciences (CAS) Large-Scale Scientific Facility Program; the CAS Center for
Excellence in Particle Physics (CCEPP); CAS Key Research Program of Frontier
Sciences under Contracts Nos. QYZDJ-SSW-SLH003, QYZDJ-SSW-SLH040; 100 Talents
Program of CAS; The Institute of Nuclear and Particle Physics (INPAC) and
Shanghai Key Laboratory for Particle Physics and Cosmology; ERC under Contract
No. 758462; European Union’s Horizon 2020 research and innovation programme
under Marie Sklodowska-Curie grant agreement under Contract No. 894790; German
Research Foundation DFG under Contracts Nos. 443159800, 455635585,
Collaborative Research Center CRC 1044, FOR5327, GRK 2149; Istituto Nazionale
di Fisica Nucleare, Italy; Ministry of Development of Turkey under Contract
No. DPT2006K-120470; National Research Foundation of Korea under Contract No.
NRF-2022R1A2C1092335; National Science and Technology fund of Mongolia;
National Science Research and Innovation Fund (NSRF) via the Program
Management Unit for Human Resources & Institutional Development, Research and
Innovation of Thailand under Contract No. B16F640076; Polish National Science
Centre under Contract No. 2019/35/O/ST2/02907; The Swedish Research Council;
U. S. Department of Energy under Contract No. DE-FG02-05ER41374.

## References

  * [1] B. Aubert _et al._ (BaBar Collaboration), Phys. Rev. Lett. 95, 142001 (2005).
  * [2] C. Z. Yuan _et al._ (Belle Collaboration), Phys.Rev.Lett. 99, 182004 (2007).
  * [3] X. L. Wang _et al._ (Belle Collaboration), Phys. Rev. Lett. 99, 142002 (2007).
  * [4] J. P. Lees _et al._ (BaBar Collaboration), Phys. Rev. D 89, 111103(R) (2014).
  * [5] C. P. Shen _et al._ (Belle Collaboration), Phys. Rev. D 89, 072015 (2014).
  * [6] M. Ablikim _et al._ (BESIII Collaboration), Phys. Rev. D 97, 071101 (2018).
  * [7] X. L. Wang _et al._ (Belle Collaboration), Phys. Rev. D 87, 051101 (2013).
  * [8] M. Ablikim _et al._ (BESIII Collaboration), Phys. Rev. D 91, 112005 (2015).
  * [9] M. Ablikim _et al._ (BESIII Collaboration), Phys. Rev. D 94, 032009 (2016).
  * [10] M. Ablikim _et al._ (BESIII Collaboration), Phys. Rev. Lett. 111, 242001 (2013).
  * [11] M. Ablikim _et al._ (BESIII Collaboration), Phys. Rev. Lett. 113, 212002 (2014).
  * [12] M. Ablikim _et al._ (BESIII Collaboration), Phys. Rev. Lett. 114, 092003 (2015).
  * [13] B. Aubert _et al._ (BaBar Collaboration), Phys. Rev. D 76, 111105 (2007).
  * [14] G. Pakhlova _et al._ (Belle Collaboration), Phys. Rev. D 77, 011103 (2008).
  * [15] M. Ablikim _et al._ (BESIII Collaboration), Phys. Rev. Lett. 130, 121901 (2023).
  * [16] M. Ablikim _et al._ (BESIII Collaboration), Phys. Rev. Lett. 122, 102002 (2019).
  * [17] M. Ablikim _et al._ (BESIII Collaboration), Phys. Rev. Lett. 118, 092001 (2017).
  * [18] M. Ablikim _et al._ (BESIII Collaboration), Phys. Rev. Lett. 118, 092002 (2017).
  * [19] D. Cronin-Hennessy _et al._ (CLEO Collaboration), Phys. Rev. D 80, 072001 (2009).
  * [20] P. del Amo Sanchez _et al._ (BaBar Collaboration), Phys. Rev. D 82, 052004 (2010).
  * [21] G. Pakhlova _et al._ (Belle Collaboration), Phys. Rev. D 83, 011101(R) (2011).
  * [22] M. Ablikim et al. (BESIII Collaboration), Nucl. Instrum. Meth. A 614, 245 (2010).
  * [23] K. X. Huang, et al. Nucl. Sci. Tech. 33, 142 (2022).
  * [24] M. Ablikim et al. (BESIII Collaboration), Chin. Phys. C 40, 063001 (2016).
  * [25] M. Ablikim et al. (BESIII Collaboration), Chin. Phys. C 45, 103001 (2021).
  * [26] M. Ablikim et al. (BESIII Collaboration), Chin. Phys. C 46, 113003 (2022).
  * [27] M. Ablikim et al. (BESIII Collaboration), Chin. Phys. C 46, 113002 (2022).
  * [28] M. Ablikim et al. (BESIII Collaboration), Chin. Phys. C 41, 063001 (2017).
  * [29] S. Agostinelli et al. (GEANT4 Collaboration), Nucl. Instrum. Meth. A 506, 250 (2003).
  * [30] D. J. Lange, Nucl. Instrum. Meth. A 462, 152155 (2001).
  * [31] R. G. Ping, Chin. Phys. C 32, 599 (2008).
  * [32] S. Jadach, B. F. L. Ward and Z. Was, Comput. Phys. Commun. 130, 260325 (2000).
  * [33] S. Jadach, B. F. L. Ward and Z. Was, Phys. Rev. D 63, 113009 (2001).
  * [34] R. L. Workman et al. (Particle Data Group), Prog. Theor. Exp. Phys. 2022, 083C01 (2022).
  * [35] 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).
  * [36] E. Richter-Was, Phys. Lett. B 303, 163 (1993).
  * [37] M. Ablikim _et al._ (BESIII Collaboration), Phys. Rev. D. 92, 092006 (2015).
  * [38] S. Actis _et al._ Eur. Phys. J. C 66, 585 (2010).
  * [39] W. Y. Sun, T. Liu, M. Q. Jing, L. L. Wang, B. Zhong, W. M. Song, Front. Phys. 16, 64501 (2021).
  * [40] Supplemental Material [link to be added].
  * [41] J. M. Blatt, V. F. Weisskopf, Sicence, 117, 419420 (1953).
  * [42] F. V. Hippel and C. Quigg, Phys. Rev. D. 5, 624638 (1972).
  * [43] L. Y. Alain, L. Jean-Pierre and R. Patrick, Phys. Rev. D 99, 073010 (2019).
  * [44] M. Ablikim _et al._ (BESIII Collaboration), Phys. Rev. D. 104, 012016 (2015).
  * [45] M. Ablikim _et al._ (BESIII Collaboration), Phys. Rev. D 81, 052005 (2010).
  * [46] M. Ablikim _et al._ (BESIII Collaboration), Chin. Phys. C 39, 093001 (2015).
  * [47] M. Mikhasenko _et al._ , Eur. Phys. J. C 78, 727 (2018).
  * [48] William S, Cleveland, Journal of the American Statistical Association 74, 829836 (1979).
  * [49] William S, Cleveland and Susan J. Devlin, Journal of the American Statistical Association 83, 596610 (1988).
  * [50] M. Ablikim _et al._ (BESII Collaboration), Phys. Lett. B 660, 315319 (2008).
  * [51] R. Aaij et al. (LHCb Collaboration), Phys. Rev. Lett. 111, 112003 (2013).
  * [52] M. Ablikim et al. (BESIII Collaboration), Chin. Phys. C 44, 040001 (2020).
  * [53] N. Hüsken, R. E. Mitchell, and E. S. Swanson, Phys. Rev. D 106, 094013 (2022).


# Region-to-region kernel interpolation of
acoustic transfer function with directional weighting

###### Abstract

A method of interpolating the acoustic transfer function (ATF) between regions
that takes into account both the physical properties of the ATF and the
directionality of region configurations is proposed. Most spatial ATF
interpolation methods are limited to estimation in the region of receivers. A
kernel method for region-to-region ATF interpolation makes it possible to
estimate the ATFs for both source and receiver regions from a discrete set of
ATF measurements. We newly formulate the reproducing kernel Hilbert space and
associated kernel function incorporating directional weight to enhance the
interpolation accuracy. We also investigate hyperparameter optimization
methods for this kernel function. Numerical experiments indicate that the
proposed method outperforms the method without the use of directional
weighting.

Index Terms—  Acoustic transfer function, Helmholtz equation, kernel ridge
regression, directional weighting, hyperparameter optimization.

##  1 Introduction

The acoustic transfer function (ATF) characterizes sound propagation between
two points in an acoustic environment, which is equivalent to the frequency
response from a source to a receiver at these points. The ATF estimation of an
environment has many applications, such as sound reproduction [1], sound field
equalization [2], echo cancellation [3], and speech dereverberation [4].

The spatial interpolation of the ATF is widely studied because of its
applicability. The ATF is generally represented as the frequency response of
an impulse signal under the assumption of a linear time-invariant system.
Thus, one method is to treat the ATF as a rational transfer function with
poles/zeros [5, 6]. Although the pole locations relative to the receiver
position can be predicted, the relation between the locations of poles and/or
zeros and the source position cannot. In addition, this method is dependent on
the room shape and its eigenmodes. Several attempts have been made to
spatially interpolate the ATF based on the sparsity of planewave components
[7, 8]. However, those methods are limited to variable receiver positions
within a receiver region, with a fixed source.

Sound field estimation or reconstruction is a problem similar to the spatial
interpolation of the ATF, and is aimed at estimating the continuous pressure
distribution from a discrete set of microphones. Many sound field estimation
methods are based on series expansions of finite-dimensional basis functions
using planewaves and spherical wavefunctions [9, 10]. An alternative approach
is the kernel method, where the solution space of the homogeneous Helmholtz
equation is defined as the reproducing kernel Hilbert space (RKHS), and the
estimate is obtained by kernel ridge regression [11, 12]. This infinite-
dimensional analysis of the sound field makes it possible to estimate it
without truncating the expansion order and has been used in various
applications, such as spatial audio reproduction [13, 14], active noise
cancellation [15], and sensor placement optimization [16].

In a prior work of the authors [17], an ATF interpolation method for both
variable source and receiver positions, i.e., region-to-region ATF
interpolation, based on the kernel method was proposed. The definition of the
RKHS is based on acoustic properties of ATFs: the constraint of satisfying the
homogeneous Helmholtz equation for the reverberant component and acoustic
reciprocity. This method significantly outperformed the method based on a
finite-dimensional expansion of spherical wavefunctions proposed in [18].

In this paper, we extend the kernel method of region-to-region ATF
interpolation to incorporate directional weight. The directional weight has
recently been introduced in the kernel method for sound field estimation to
take prior information on source directions into consideration [19, 20]. We
introduce the directional weight for the region-to-region ATF interpolation to
incorporate knowledge of the configuration of source and receiver regions.
Since the use of directional weight requires the determination of several
hyperparameters of the reproducing kernel functions, we also investigate
hyperparameter optimization methods. Numerical experiments are conducted to
compare the proposed method with the method without directional weighting in
order to evaluate the effect of directionality on the estimations.

##  2 Problem statement and preliminaries

Given a space $\Omega\subset\mathbb{R}^{3}$ with stationary acoustic
properties, our objective is to estimate the ATF
$h\colon\Omega_{\mathrm{R}}\times\Omega_{\mathrm{S}}\to\mathbb{C}$ between any
source/receiver pair of positions
$\mathbf{r}\in\Omega_{\mathrm{R}}\subset\Omega$ and
$\mathbf{s}\in\Omega_{\mathrm{S}}\subset\Omega$, where $\Omega_{\mathrm{S}}$
and $\Omega_{\mathrm{R}}$ are the source and receiver regions, respectively
(see Fig. 1).

###  2.1 Preliminaries

We assume that any ATF $h$ can be separated into two components: a known
direct component $h_{\mathrm{D}}$ given as the free-field Green’s function
$G_{0}$ and an unknown reverberant component $h_{\mathrm{R}}$, as shown in
[18, 17]. $h_{\mathrm{R}}$ satisfies the homogeneous Helmholtz equation [21].
These assumptions are represented as

$\displaystyle
h(\mathbf{r}|\mathbf{s},k)=h_{\mathrm{D}}(\mathbf{r}|\mathbf{s},k)+h_{\mathrm{R}}(\mathbf{r}|\mathbf{s},k)$
(1) $\displaystyle
h_{\mathrm{D}}(\mathbf{r}|\mathbf{s},k)=G_{0}(\mathbf{r}|\mathbf{s},k)=\frac{e^{\mathrm{i}k\|\mathbf{r}-\mathbf{s}\|}}{4\pi\|\mathbf{r}-\mathbf{s}\|}$
(2)
$\displaystyle(\nabla^{2}_{\mathbf{r}}+k^{2})h_{\mathrm{R}}(\mathbf{r}|\mathbf{s},k)=(\nabla^{2}_{\mathbf{s}}+k^{2})h_{\mathrm{R}}(\mathbf{r}|\mathbf{s},k)=0,$
(3)

where $\mathbf{r}|\mathbf{s}$ is the source/receiver pair of positions,
$\mathrm{i}$ is the imaginary unit, $\nabla^{2}_{\mathbf{r}}$ and
$\nabla^{2}_{\mathbf{s}}$ are the Laplacian operators on the coordinates of
$\mathbf{r}$ and $\mathbf{s}$, respectively, $k=2\pi\breve{f}/c$ is the
wavenumber, $\breve{f}$ is the frequency, and $c$ is the speed of sound.
Hereafter, $k$ in the argument of the ATFs is omitted for notational
simplicity.

###  2.2 Region-to-region ATF interpolation problem

We distribute a set of $M$ receivers at points
$\\{\mathbf{r}_{m}\\}_{m=1}^{M}$ and $L$ sources at points
$\\{\mathbf{s}_{l}\\}_{l=1}^{L}$ to obtain a total of $N$ ($=LM$) ATF
measurement values. We collectively denote
$\mathbf{q}_{n}=\mathbf{r}_{m}|\mathbf{s}_{l}$ for the position pairs with
index $n=m+(l-1)M$ ($\in\\{1,\ldots,N\\}$). A set of $N$ ATFs is given, and
then the direct component is removed from them to obtain the measurement
vector $\mathbf{y}=[y_{1},\ldots,y_{N}]^{\mathsf{T}}$ corresponding to each
$\mathbf{q}_{n}$. We define our optimization problem as

$\displaystyle\hat{h}_{\mathrm{R}}=\mathop{\rm
arg~{}min}\limits_{f\in\mathscr{H}}\mathcal{J}(f)$
$\displaystyle\mathcal{J}(f)\coloneqq\sum_{n=1}^{N}|y_{n}-f(\mathbf{q}_{n})|^{2}+\lambda\|f\|_{\mathscr{H}}^{2},\
f\in\mathscr{H},$ (4)

where $\lambda>0$ is the regularization constant and $\mathscr{H}$ is the
feature space to which the interpolation of the reverberant component belongs.
This interpolation function is then added to the direct component to obtain
$\hat{h}$:

$\hat{h}(\mathbf{r}|\mathbf{s})=h_{\mathrm{D}}(\mathbf{r}|\mathbf{s})+\hat{h}_{\mathrm{R}}(\mathbf{r}|\mathbf{s}).$
(5)

Fig. 1: Schematic diagram of a region-to-region ATF interpolation problem.

###  2.3 Kernel ridge regression

We consider the space $\mathscr{H}$ to be a reproducing kernel Hilbert space
(RKHS). That is, it is a Hilbert space
$(\mathscr{H},\langle\cdot,\cdot\rangle_{\mathscr{H}})$ that admits a
reproducing kernel function $\kappa$. The reproducing kernel is a bivariate
function that satisfies:

$\langle\kappa(\cdot,\mathbf{r}|\mathbf{s}),f\rangle_{\mathscr{H}}=f(\mathbf{r}|\mathbf{s}),\
\forall\mathbf{r}\in\Omega_{\mathrm{R}},\
\forall\mathbf{s}\in\Omega_{\mathrm{S}}.$ (6)

In that case, the interpolation function $\hat{h}_{\mathrm{R}}$ is given by
kernel ridge regression:

$\hat{h}_{\mathrm{R}}(\mathbf{r}|\mathbf{s})=\bm{\kappa}(\mathbf{r}|\mathbf{s})(\mathbf{K}+\lambda\mathbf{I})^{-1}\mathbf{y},$
(7)

where
$\bm{\kappa}(\mathbf{r}|\mathbf{s})=[\kappa(\mathbf{r}|\mathbf{s},\mathbf{q}_{1}),\dots,\kappa(\mathbf{r}|\mathbf{s},\mathbf{q}_{N})]$
is the kernel function vector, $\mathbf{I}$ is the identity matrix, and
$\mathbf{K}$ is the Gram matrix defined as

$\mathbf{K}=\begin{bmatrix}\kappa(\mathbf{q}_{1},\mathbf{q}_{1})&\kappa(\mathbf{q}_{1},\mathbf{q}_{2})&\dots&\kappa(\mathbf{q}_{1},\mathbf{q}_{N})\\\
\vdots&\vdots&\ddots&\vdots\\\
\kappa(\mathbf{q}_{N},\mathbf{q}_{1})&\kappa(\mathbf{q}_{N},\mathbf{q}_{2})&\dots&\kappa(\mathbf{q}_{N},\mathbf{q}_{N})\end{bmatrix}.$
(8)

##  3 Region-to-region ATF interpolation with directional weighting

We previously formulated an RKHS for solving (4) based on the spherical
wavefunction expansion in a previous study [17]. We now formulate an RKHS
based on the planewave expansion incorporating directional weight to take into
consideration the behavior of the ATF with respect to the relative positions
of $\Omega_{\mathrm{R}}$ and $\Omega_{\mathrm{S}}$.

###  3.1 Feature space definition

Since the reverberant component $h_{\mathrm{R}}$ does not include any sources
in $\Omega_{\mathrm{S}}$ and $\Omega_{\mathrm{R}}$, $h_{\mathrm{R}}$ can be
represented by a planewave expansion, i.e., Herglotz wavefunction [22], as

$\displaystyle
h_{\mathrm{R}}(\mathbf{r}|\mathbf{s})=\mathcal{I}\left(\tilde{h}_{\mathrm{R}};\mathbf{r}|\mathbf{s}\right),$
(9)
$\displaystyle\mathcal{I}\left(f;\mathbf{r}|\mathbf{s}\right)\coloneqq\int_{\mathbb{S}^{2}\times\mathbb{S}^{2}}e^{\mathrm{i}k(\hat{\mathbf{r}}\cdot\mathbf{r}+\hat{\mathbf{s}}\cdot\mathbf{s})}f(\hat{\mathbf{r}},\hat{\mathbf{s}})\mathrm{d}\hat{\mathbf{r}}\mathrm{d}\hat{\mathbf{s}},$
(10)

where $\mathcal{I}$ is the operator for the planewave expansion, and
$\mathbb{S}^{2}$ is the set of vectors in $\mathbb{R}^{3}$ of unit norm,
representing directions. Therefore,
$\tilde{h}_{\mathrm{R}}(\hat{\mathbf{r}},\hat{\mathbf{s}})$ is the complex
amplitude of the planewave component of $h_{\mathrm{R}}$ from the source
direction $\hat{\mathbf{s}}$ to the receiver direction $\hat{\mathbf{r}}$. The
reciprocity of $h_{\mathrm{R}}$ imposes the condition that
$\tilde{h}_{\mathrm{R}}(\hat{\mathbf{r}},\hat{\mathbf{s}})=\tilde{h}_{\mathrm{R}}(\hat{\mathbf{s}},\hat{\mathbf{r}})$.
Thus, our inner-product space
$(\mathscr{H},\langle\cdot,\cdot\rangle_{\mathscr{H}})$ is defined as

$\begin{split}\mathscr{H}=&\left\\{h_{\mathrm{R}}=\mathcal{I}\left(\tilde{h}_{\mathrm{R}};\mathbf{r}|\mathbf{s}\right)\colon\tilde{h}_{\mathrm{R}}\in
L^{2}(W,\mathbb{S}^{2}\times\mathbb{S}^{2}),\right.\\\ &\left.\hskip
70.0pt\tilde{h}_{\mathrm{R}}(\hat{\mathbf{r}},\hat{\mathbf{s}})=\tilde{h}_{\mathrm{R}}(\hat{\mathbf{s}},\hat{\mathbf{r}})\
\forall\hat{\mathbf{r}},\hat{\mathbf{s}}\in\mathbb{S}^{2}\vphantom{\hat{h}}\right\\}\end{split}$
(11) $\langle
f,g\rangle_{\mathscr{H}}=\int_{\mathbb{S}^{2}\times\mathbb{S}^{2}}\frac{\overline{\tilde{f}(\hat{\mathbf{r}},\hat{\mathbf{s}})}\tilde{g}(\hat{\mathbf{r}},\hat{\mathbf{s}})}{W(\hat{\mathbf{r}},\hat{\mathbf{s}})}\mathrm{d}\hat{\mathbf{r}}\mathrm{d}\hat{\mathbf{s}},\
\forall f,g\in\mathscr{H},$ (12)

where $\overline{\ \cdot\ }$ is the complex conjugate,
$W\colon\mathbb{S}^{2}\times\mathbb{S}^{2}\to\mathbb{R}_{+}$ is a directional
weighting function, and $L^{2}(W,\mathbb{S}^{2}\times\mathbb{S}^{2})$ is the
space of functions of bounded square integral for $W$. Under these conditions,
$(\mathscr{H},\langle\cdot,\cdot\rangle_{\mathscr{H}})$ inherits the
completeness of $L^{2}$ spaces [23] and as such is a Hilbert space. We can
show that it is an RKHS by showing that the following function is its
reproducing kernel:

To show that this function is the reproducing kernel, we only have to show
that it satisfies the projection property for any given $f\in\mathscr{H}$:

$\displaystyle\langle\kappa(\cdot,\mathbf{r}|\mathbf{s}),f\rangle_{\mathscr{H}}$
$\displaystyle=\int_{\mathbb{S}^{2}\times\mathbb{S}^{2}}\frac{e^{\mathrm{i}k(\hat{\mathbf{r}}\cdot\mathbf{r}+\hat{\mathbf{s}}\cdot\mathbf{s})}\tilde{f}+e^{\mathrm{i}k(\hat{\mathbf{r}}\cdot\mathbf{s}+\hat{\mathbf{s}}\cdot\mathbf{r})}\tilde{f}}{2}\mathrm{d}\hat{\mathbf{r}}\mathrm{d}\hat{\mathbf{s}}$
$\displaystyle=\frac{1}{2}\left(\mathcal{I}\left(\tilde{f};\mathbf{r}|\mathbf{s}\right)+\mathcal{I}\left(\tilde{f};\mathbf{s}|\mathbf{r}\right)\right)$
$\displaystyle=\frac{f(\mathbf{r}|\mathbf{s})+f(\mathbf{s}|\mathbf{r})}{2}.$
(13)

Since $f$ is reciprocal, we have

$\langle\kappa(\cdot,\mathbf{r}|\mathbf{s}),f\rangle_{\mathscr{H}}=f(\mathbf{r}|\mathbf{s}).$
(14)

We also consider that the directional weighting function $W$ is separable for
$\hat{\mathbf{r}}$ and $\hat{\mathbf{s}}$, that is,

$W(\hat{\mathbf{r}},\hat{\mathbf{s}})=w(\hat{\mathbf{r}})w(\hat{\mathbf{s}}),$
(15)

where $w\colon\mathbb{S}^{2}\to\mathbb{R}_{+}$ is a directional weighting
function in a single direction.

###  3.2 Proposed directional weighting function

Since the direct component is removed from the measurements to obtain the
reverberant component, the directional weighting function for the region-to-
region ATF interpolation should have minimal gain in the direction connecting
the centers of both regions, which is defined as $\hat{\mathbf{v}}_{0}$.
Therefore, we define the directional weighting function $w$ as

$w(\hat{\mathbf{v}})=\frac{1}{4\pi}\left(1+\gamma^{2}-\frac{\cosh(\beta\hat{\mathbf{v}}\cdot\hat{\mathbf{v}}_{0})}{\cosh(\beta)}\right),\
\hat{\mathbf{v}}\in\mathbb{S}^{2}.$ (16)

The hyperparameters $\beta$ and $\gamma$ characterize $\Omega$ acoustically.
Since the direct components are expected to be weaker, the early reflections
are expected to be significantly more influential. An increase in $\beta$
makes the weight more selective in regards to the lateral components, while
$\gamma$ sets the minimum gain baseline for all directions.

To better illustrate this, an example of the directional weighting function is
plotted in Fig. 2, which shows that the weight is minimal in the chosen
direction $\hat{\mathbf{v}}_{0}$ but maximal in the lateral directions.
$\beta$ controls the width of the cavity passing through
$\hat{\mathbf{v}}_{0}$, while $\gamma$ controls the depth.

###  3.3 Relation to prior work using uniform weight

Fig. 2: Directional weighting function $w$ when
$\hat{\mathbf{v}}_{0}=[-1/\sqrt{2},-1/\sqrt{2},0]^{\mathsf{T}}$,
$\gamma=0.01$, and $\beta=100$. The gain is given by the distance from the
origin to the surface, quantified according to the bar on the right.

By setting $\gamma=1$ and $\beta=0$, the directional weighting function
becomes uniform, i.e., $w=1/4\pi$. The reproducing kernel of this uniform
weight has the closed form

$\begin{split}\kappa(\mathbf{r}|\mathbf{s},\mathbf{r}^{\prime}|\mathbf{s}^{\prime})=\frac{1}{2}\left(j_{0}(k\|\mathbf{r}-\mathbf{r}^{\prime}\|)j_{0}(k\|\mathbf{s}-\mathbf{s}^{\prime}\|)\right.\\\
\ \ \ \
\left.+j_{0}(k\|\mathbf{s}-\mathbf{r}^{\prime}\|)j_{0}(k\|\mathbf{r}-\mathbf{s}^{\prime}\|)\right),\end{split}$
(17)

where $j_{0}$ is the $0$th-order spherical Bessel function of the first kind.
This kernel function is identical to that used in [17], making both
estimations equivalent.

##  4 Hyperparameter optimization for directional weight

The introduction of hyperparameters in (16) also necessitates a criterion for
choosing them. Although the profile of the reverberant component is
understood, the exact balance between the direct and lateral component gains
is not known outright. There are several methods of hyperparameter
optimization used for kernel ridge regression [24, 25]. For this application,
we employ leave-one-out cross-validation (LOO) because of its simplicity and
nearly unbiased nature [26].

We begin the computation of LOO by selecting a data point and measurement
point pair $(\mathbf{q}_{n},y_{n})$ to remove from the data set. We then
derive the desired model using the remaining $N-1$ elements of the data set
and compute the error in estimating $y_{n}$ with this partial model. The value
of LOO is the average of all errors calculated by repeating this process
exhaustively:

$\mathrm{LOO(\mathbf{y},\ell)}=\frac{1}{N}\sum_{n=1}^{N}\ell\left(\hat{f}_{n}(\mathbf{q}_{n})-y_{n}\right),$
(18)

where $\ell$ is the chosen loss function, and $\hat{f}_{n}$ is the model
derived when considering all data pairs except $(\mathbf{q}_{n},y_{n})$. We
consider two types of loss function $\ell$: square error (SQE) and Tukey loss.

$\displaystyle\mathrm{SQE}(z)=|z|^{2}$ (19)
$\displaystyle\mathrm{Tukey}(z)=\begin{cases}\displaystyle\frac{\sigma^{2}}{6}\left(1-\left(1-\frac{|z|^{2}}{\sigma^{2}}\right)^{3}\right),\
|z|\leq\sigma\\\ \displaystyle\frac{\sigma^{2}}{6},\ |z|>\sigma\end{cases},$
(20)

where $z$ is a complex variable, and $\sigma$ is the selectivity parameter.
Tukey loss is more selective than SQE, which means that these loss functions
have similar behavior near $0$, but Tukey loss has slower growth for higher
discrepancies. This selectivity makes Tukey loss more resilient to outliers.

Both SQE and Tukey loss are differentiable with respect to the
hyperparameters, meaning that LOO can be optimized using gradient descent
methods. We applied the improved robust back-propagation algorithm introduced
in [27].

##  5 Numerical simulations

Fig. 3: Experimental setup utilized for the simulations. Red diamonds
represent sources and blue circles represent receivers.

The proposed method using directional weight was evaluated using 3D acoustic
simulations based on the image source method [28] by comparing it with the
method using uniform weight. The room used for the evaluation was shoebox-
shaped and had $3.2~{}\mathrm{m}\times 4.0~{}\mathrm{m}\times
2.7~{}\mathrm{m}$ dimensions. The reflection coefficients of the walls were
set so that the reverberation time $T_{60}$ was $0.45~{}\mathrm{s}$. The
source and receiver regions were spheres of $0.2~{}\mathrm{m}$ radius whose
centers were $\mathbf{s}_{0}=[0.35,0.43,0.29]^{\mathsf{T}}~{}\mathrm{m}$ for
$\Omega_{\mathrm{S}}$ and $\mathbf{r}_{0}=-\mathbf{s}_{0}$ for
$\Omega_{\mathrm{R}}$. The origin of the coordinate system was set at the
center of the room. We used
$\hat{\mathbf{v}}_{0}=(\mathbf{r}_{0}-\mathbf{s}_{0})/\|\mathbf{r}_{0}-\mathbf{s}_{0}\|$
for the directional kernel.

The measurement points were given for a total of $L=M=41$ points distributed
on two spherical layers. The point distribution was given by the spherical
$t$-design [29]: $t=4$ for the outer layer and $t=3$ for the inner layer. We
also added Gaussian noise so that the signal-to-noise ratio was
$20~{}\mathrm{dB}$. The parameter $\sigma$ in the Tukey loss (20) was set to
$0.4$, which was determined so that the probability of the error falling
within four standard deviations was above $99\%$. The regularization parameter
$\lambda$ in (7) was set to $10^{-2}$.

The evaluation measure was defined as the normalized mean square error (NMSE),

$\mathrm{NMSE}=10\log_{10}\left(\frac{\sum_{n}\left|\hat{h}(\mathbf{q}^{\prime}_{n})-h(\mathbf{q}^{\prime}_{n})\right|^{2}}{\sum_{n}|h(\mathbf{q}^{\prime}_{n})|^{2}}\right),$
(21)

where $\mathbf{q}^{\prime}_{n}$ denote the $n$th evaluation source/receiver
pair. We set $9025$ source/receiver pairs of evaluation.

The NMSE with respect to the frequency ranging from $100$ to
$1150~{}\mathrm{Hz}$ is shown in Fig. 4. We also show the estimated ATF in
$\Omega_{\mathrm{R}}$ generated by a single source point at
$[0.35,0.43,0.29]^{\mathsf{T}}~{}\mathrm{m}$ for $950~{}\mathrm{Hz}$ in Fig.
5. The normalized squared error distribution of this plot is also shown in
Fig. 6, which is defined as

$\mathrm{NSE}(\mathbf{r})=10\log_{10}\left(\frac{\left|h(\mathbf{r}|\mathbf{s}_{0})-\hat{h}(\mathbf{r}|\mathbf{s}_{0})\right|^{2}}{|h(\mathbf{r}|\mathbf{s}_{0})|^{2}}\right).$
(22)

The proposed method using directionally weighted kernels achieved a lower NMSE
than the uniform counterpart for every frequency. For the SNR of
$20~{}\mathrm{dB}$, the more selective Tukey loss also significantly
outperformed the simpler SQE used in LOO. These results indicate that the
directional weight is effective for the region-to-region ATF interpolation
problem. In addition, LOO with Tukey loss is more useful for hyperparameter
optimization than the simple SQE loss.

Fig. 4: NMSE performance of the methods compared.

(a) True (b) Uniform

(c) Directional - SQE (d) Directional - Tukey

Fig. 5: Distributions of true and estimated ATFs in $\Omega_{\mathrm{R}}$ from
the center of $\Omega_{\mathrm{S}}$ at
$[0.35,0.43,0.29]^{\mathsf{T}}~{}\mathrm{m}$ for $950~{}\mathrm{Hz}$. The
black circle indicates the bounds of $\Omega_{\mathrm{R}}$.

(a) Uniform

(b) Directional - SQE (c) Directional - Tukey

Fig. 6: Distributions of NSEs in $\Omega_{\mathrm{R}}$.

##  6 Conclusion

We proposed a kernel interpolation method for region-to-region ATF
interpolation with directional weighting. The reproducing kernel Hilbert space
and associated reproducing kernel function are formulated on the basis of
planewave decomposition, having properties of the Helmholtz equation
constraint, acoustic reciprocity, and directionality. The spatial
interpolation of the ATF is achieved by kernel ridge regression using this
kernel function. Hyperparameters included in the kernel function are optimized
by leave-one-out cross-validation. In the numerical experiments, the proposed
method achieved highly accurate interpolation compared with the method using
uniform weight. Furthermore, robustness to the noise was significantly
improved by using Tukey loss in LOO.

##  7 Acknowledgements

This work was supported by JSPS KAKENHI Grant Number JP19H01116 and JST PRESTO
Grant Number JPMJPR18J4.

## References

  * [1] F. Borra, I. D. Gebru, and D. Markovic,  “Soundfield reconstruction in reverberant environments using higher-order microphones and impulse response measurements,”  in Proc. IEEE Int. Conf. Acoust., Speech, Signal Process. (ICASSP), 2019, pp. 281–285.
  * [2] R. Mazur, F. Katzberg, and A. Mertins,  “Robust room equalization using sparse sound-field reconstruction,”  in Proc. IEEE Int. Conf. Acoust., Speech, Signal Process. (ICASSP), 2019, pp. 4230–4234.
  * [3] P. Gil-Cacho, T. van Waterschoot, M. Moonen, and S. Jensen,  “Multi-microphone acoustic echo cancellation using multi-channel warped linear prediction of common acoustical poles,”  in Proc. European Signal Process. Conf. (EUSIPCO), 2010, vol. 2010, pp. 2121–2125.
  * [4] N. Mohanan, R. Velmurugan, and P. Rao,  “Speech dereverberation using NMF with regularized room impulse response,”  in Proc. IEEE Int. Conf. Acoust., Speech, Signal Process. (ICASSP), 2017, pp. 4955–4959.
  * [5] Y. Haneda, S. Makino, Y. Kaneda, and N. Koizumi,  “ARMA modeling of a room transfer function at low frequencies,”  J. Acoust. Soc. Jpn. (E), vol. 15, pp. 353–355, 1994.
  * [6] Y. Haneda, Y. Kaneda, and N. Kitawaki,  “Common-acoustical-pole and residue model and its application to spatial interpolation and extrapolation of a room transfer function,”  IEEE Trans. Speech Audio Process., vol. 7, no. 6, pp. 709–717, 1999.
  * [7] R. Mignot, G. Chardon, and L. Daudet,  “Low frequency interpolation of room impulse responses using compressed sensing,”  IEEE/ACM Trans. Audio, Speech, Lang. Process., vol. 22, no. 1, pp. 205–216, 2014.
  * [8] N. Antonello, E. De Sena, M. Moonen, P. A. Naylor, and T. van Waterschoot,  “Room impulse response interpolation using a sparse spatio-temporal representation of the sound field,”  IEEE/ACM Trans. Audio, Speech, Lang. Process., vol. 25, no. 10, pp. 1929–1941, 2017.
  * [9] M. A. Poletti,  “Three-dimensional surround sound systems based on spherical harmonics,”  J. Audio Eng. Soc., vol. 53, no. 11, pp. 1004–1025, 2005.
  * [10] S. Koyama, K. Furuya, Y. Hiwasaki, and Y. Haneda,  “Analytical approach to wave field reconstruction filtering in spatio-temporal frequency domain,”  IEEE Trans. Audio, Speech, Lang. Process., vol. 21, no. 4, pp. 685–696, 2013.
  * [11] N. Ueno, S. Koyama, and H. Saruwatari,  “Kernel ridge regression with constraint of Helmholtz equation for sound field interpolation,”  in Proc. Int. Workshop Acoust. Signal Enhancement (IWAENC), 2018, pp. 436–440.
  * [12] N. Ueno, S. Koyama, and H. Saruwatari,  “Sound field recording using distributed microphones based on harmonic analysis of infinite order,”  IEEE Signal Process. Lett., vol. 25, no. 1, pp. 135–139, 2018\.
  * [13] N. Iijima, S. Koyama, and H. Saruwatari,  “Binaural rendering from microphone array signals of arbitrary geometry,”  J. Acoust. Soc. Am., 2021,  (in press).
  * [14] S. Koyama, K. Kimura, and N. Ueno,  “Sound field reproduction with weighted mode matching and infinite-dimensional harmonic analysis: An experimental evaluation,”  in Proc. International Conference on Immersive and 3D Audio (I3DA), 2021.
  * [15] S. Koyama, J. Brunnström, H. Ito, N. Ueno, and H. Saruwatari,  “Spatial active noise control based on kernel interpolation of sound field,”  IEEE/ACM Trans. Audio, Speech, Lang. Process., 2021,  (in press).
  * [16] K. Ariga, T. Nishida, S. Koyama, N. Ueno, and H. Saruwatari,  “Mutual-information-based sensor placement for spatial sound field recording,”  in Proc. IEEE Int. Conf. Acoust., Speech, Signal Process. (ICASSP), 2020, pp. 166–170.
  * [17] J. G. C. Ribeiro, N. Ueno, S. Koyama, and H. Saruwatari,  “Kernel interpolation of acoustic transfer function between regions considering reciprocity,”  in IEEE Sensor Array Multichannel Signal Process. Workshop (SAM), 2020, pp. 1–5.
  * [18] P. N. Samarasinghe, T. D. Abhayapala, M. A. Poletti, and T. Betlehem,  “An efficient parameterization of the room transfer function,”  IEEE/ACM Trans. Audio, Speech, Lang. Process., vol. 23, no. 12, pp. 2217–2227, 2015.
  * [19] H. Ito, S. Koyama, N. Ueno, and H. Saruwatari,  “Spatial active noise control based on kernel interpolation with directional weighting,”  in Proc. IEEE Int. Conf. Acoust., Speech, Signal Process. (ICASSP), 2020, pp. 8399–8403.
  * [20] N. Ueno, S. Koyama, and H. Saruwatari,  “Directionally weighted wave field estimation exploiting prior information on source direction,”  IEEE Trans. Signal Process., vol. 69, pp. 2383–2395, 2021.
  * [21] E. G. Williams,  Fourier Acoustics,  Academic Press, London, 1999.
  * [22] D. Colton and P. Monk,  Herglotz Wave Functions in Inverse Electromagnetic Scattering Theory, pp. 367–394,  Springer Berlin Heidelberg, Berlin, Heidelberg, 2003.
  * [23] W. Rudin,  Functional Analysis,  McGraw-Hill, New York City, 1991.
  * [24] R. Horiuchi, S. Koyama, J. G. C. Ribeiro, N. Ueno, and H. Saruwatari,  “Kernel learning for sound field estimation with l1 and l2 regularizations,”  in Proc. IEEE Int. Workshop Appl. Signal Process. Audio Acoust. (WASPAA), 2021 (to appear).
  * [25] D. Caviedes Nozal, N. Riis, F. Heuchel, J. Brunskog, P. Gerstoft, and E. Fernandez-Grande,  “Gaussian processes for sound field reconstruction,”  J. Acoust. Soc. Am., vol. 149, pp. 1107–1119, 2021.
  * [26] A. Elisseeff and M. Pontil,  “Leave-one-out error and stability of learning algorithms with applications,”  NATO Sci. Ser., III: Comput. Systems Sci., vol. 190, pp. 111–130, 2003.
  * [27] A. Anastasiadis, G. Magoulas, and M. Vrahatis,  “An efficient improvement of the Rprop algorithm,”  in Proc. Int. Workshop Artificial Neural Netw. Pattern Recognit. (ANNPR), 2003.
  * [28] J. B. Allen and D. A. Berkley,  “Image method for efficiently simulating small-room acoustics,”  J. Acoust. Soc. Am., vol. 65, no. 4, pp. 943–950, 1979.
  * [29] X. Chen and R. Womersley,  “Spherical t-design with d=(t+1)`^`2 points,”  (Accessed on Oct. 26, 2018).  [Online].  Available: http://www.polyu.edu.hk/ama/staff/xjchen/sphdesigns.html.