00footnotetext: Received xxxx

# Amplitude analysis of $D^{0}\to\pi^{+}\pi^{-}\pi^{+}\pi^{-}$ and
$D^{0}\to\pi^{+}\pi^{-}\pi^{0}\pi^{0}$ ††thanks: 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, 11835012, 11935015, 11935016,
11935018, 11961141012, 12025502, 12035009, 12035013, 12061131003, 12105276,
12122509, 12192260, 12192261, 12192262, 12192263, 12192264, 12192265,
12221005, 12225509, 12235017; the Chinese Academy of Sciences (CAS) Large-
Scale Scientific Facility Program; the CAS Center for Excellence in Particle
Physics (CCEPP); Joint Large-Scale Scientific Facility Funds of the NSFC and
CAS under Contract No. U1732263, U1832103, U1832207, U2032111; 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; 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. 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

###### Abstract

Using $e^{+}e^{-}$ annihilation data corresponding to an integrated luminosity
of 2.93 $\rm fb^{-1}$ taken at the center-of-mass energy $\sqrt{s}=3.773$ GeV
with the BESIII detector, a joint amplitude analysis is performed on the
decays $D^{0}\to\pi^{+}\pi^{-}\pi^{+}\pi^{-}$ and
$D^{0}\to\pi^{+}\pi^{-}\pi^{0}\pi^{0}$(non-$\eta$). The fit fractions of
individual components are obtained, and large interferences among the dominant
components of $D^{0}\to a_{1}(1260)\pi$, $D^{0}\to\pi(1300)\pi$,
$D^{0}\to\rho(770)\rho(770)$ and $D^{0}\to 2(\pi\pi)_{S}$ are found in both
channels. With the obtained amplitude model, the $CP$-even fractions of
$D^{0}\to\pi^{+}\pi^{-}\pi^{+}\pi^{-}$ and
$D^{0}\to\pi^{+}\pi^{-}\pi^{0}\pi^{0}$(non-$\eta$) are determined to be
$(75.2\pm 1.1_{\rm stat.}\pm 1.5_{\rm syst.})\%$ and $(68.9\pm 1.5_{\rm
stat.}\pm 2.4_{\rm syst.})\%$, respectively. The branching fractions of
$D^{0}\to\pi^{+}\pi^{-}\pi^{+}\pi^{-}$ and
$D^{0}\to\pi^{+}\pi^{-}\pi^{0}\pi^{0}$(non-$\eta$) are measured to be
$(0.688\pm 0.010_{\rm stat.}\pm 0.010_{\rm syst.})\%$ and $(0.951\pm
0.025_{\rm stat.}\pm 0.021_{\rm syst.})\%$, respectively. The amplitude
analysis provides an important model for binning strategy in the measurements
of the strong phase parameters of $D^{0}\to 4\pi$ when used to determine the
CKM angle $\gamma(\phi_{3})$ via the $B^{-}\to DK^{-}$ decay.

###### keywords:

BESIII, $D^{0}$ meson decays, amplitude analysis, $CP$-even fraction

M. Ablikim1, M. N. Achasov4,b, P. Adlarson75, O. Afedulidis3, X. C. Ai80, R.
Aliberti35, A. Amoroso74A,74C, Q. An71,58, Y. Bai57, O. Bakina36, I.
Balossino29A, Y. Ban46,g, H.-R. Bao63, V. Batozskaya1,44, K. Begzsuren32, N.
Berger35, M. Berlowski44, M. Bertani28A, D. Bettoni29A, F. Bianchi74A,74C, E.
Bianco74A,74C, A. Bortone74A,74C, I. Boyko36, R. A. Briere5, A. Brueggemann68,
H. Cai76, X. Cai1,58, A. Calcaterra28A, G. F. Cao1,63, N. Cao1,63, S. A.
Cetin62A, J. F. Chang1,58, W. L. Chang1,63, G. R. Che43, G. Chelkov36,a, C.
Chen43, C. H. Chen9, Chao Chen55, G. Chen1, H. S. Chen1,63, M. L. Chen1,58,63,
S. J. Chen42, S. L. Chen45, S. M. Chen61, T. Chen1,63, X. R. Chen31,63, X. T.
Chen1,63, Y. B. Chen1,58, Y. Q. Chen34, Z. J. Chen25,h, Z. Y. Chen1,63, S. K.
Choi10A, X. Chu43, G. Cibinetto29A, F. Cossio74C, J. J. Cui50, H. L. Dai1,58,
J. P. Dai78, A. Dbeyssi18, R. E. de Boer3, D. Dedovich36, C. Q. Deng72, Z. Y.
Deng1, A. Denig35, I. Denysenko36, M. Destefanis74A,74C, F. De Mori74A,74C, B.
Ding66,1, X. X. Ding46,g, Y. Ding34, Y. Ding40, J. Dong1,58, L. Y. Dong1,63,
M. Y. Dong1,58,63, X. Dong76, M. C. Du1, S. X. Du80, Z. H. Duan42, P.
Egorov36,a, Y. H. Fan45, J. Fang1,58, J. Fang59, S. S. Fang1,63, W. X. Fang1,
Y. Fang1, Y. Q. Fang1,58, R. Farinelli29A, L. Fava74B,74C, F. Feldbauer3, G.
Felici28A, C. Q. Feng71,58, J. H. Feng59, Y. T. Feng71,58, K. Fischer69, M.
Fritsch3, C. D. Fu1, J. L. Fu63, Y. W. Fu1, H. Gao63, Y. N. Gao46,g, Yang
Gao71,58, S. Garbolino74C, I. Garzia29A,29B, P. T. Ge76, Z. W. Ge42, C.
Geng59, E. M. Gersabeck67, A. Gilman69, K. Goetzen13, L. Gong40, W. X.
Gong1,58, W. Gradl35, S. Gramigna29A,29B, M. Greco74A,74C, M. H. Gu1,58, Y. T.
Gu15, C. Y. Guan1,63, Z. L. Guan22, A. Q. Guo31,63, L. B. Guo41, M. J. Guo50,
R. P. Guo49, Y. P. Guo12,f, A. Guskov36,a, J. Gutierrez27, K. L. Han63, T. T.
Han1, X. Q. Hao19, F. A. Harris65, K. K. He55, K. L. He1,63, F. H. Heinsius3,
C. H. Heinz35, Y. K. Heng1,58,63, C. Herold60, T. Holtmann3, P. C. Hong12,f,
G. Y. Hou1,63, X. T. Hou1,63, Y. R. Hou63, Z. L. Hou1, B. Y. Hu59, H. M.
Hu1,63, J. F. Hu56,i, T. Hu1,58,63, Y. Hu1, G. S. Huang71,58, K. X. Huang59,
L. Q. Huang31,63, X. T. Huang50, Y. P. Huang1, T. Hussain73, F. Hölzken3, N
Hüsken27,35, N. in der Wiesche68, M. Irshad71,58, J. Jackson27, S. Janchiv32,
J. H. Jeong10A, Q. Ji1, Q. P. Ji19, W. Ji1,63, X. B. Ji1,63, X. L. Ji1,58, Y.
Y. Ji50, X. Q. Jia50, Z. K. Jia71,58, D. Jiang1,63, H. B. Jiang76, P. C.
Jiang46,g, S. S. Jiang39, T. J. Jiang16, X. S. Jiang1,58,63, Y. Jiang63, J. B.
Jiao50, J. K. Jiao34, Z. Jiao23, S. Jin42, Y. Jin66, M. Q. Jing1,63, X. M.
Jing63, T. Johansson75, S. Kabana33, N. Kalantar-Nayestanaki64, X. L. Kang9,
X. S. Kang40, M. Kavatsyuk64, B. C. Ke80, V. Khachatryan27, A. Khoukaz68, R.
Kiuchi1, O. B. Kolcu62A, B. Kopf3, M. Kuessner3, X. Kui1,63, A. Kupsc44,75, W.
Kühn37, J. J. Lane67, P. Larin18, L. Lavezzi74A,74C, T. T. Lei71,58, Z. H.
Lei71,58, H. Leithoff35, M. Lellmann35, T. Lenz35, C. Li47, C. Li43, C. H.
Li39, Cheng Li71,58, D. M. Li80, F. Li1,58, G. Li1, H. Li71,58, H. B. Li1,63,
H. J. Li19, H. N. Li56,i, Hui Li43, J. R. Li61, J. S. Li59, Ke Li1, L. J
Li1,63, L. K. Li1, Lei Li48, M. H. Li43, P. R. Li38,k, Q. M. Li1,63, Q. X.
Li50, R. Li17,31, S. X. Li12, T. Li50, W. D. Li1,63, W. G. Li1, X. Li1,63, X.
H. Li71,58, X. L. Li50, Xiaoyu Li1,63, Y. G. Li46,g, Z. J. Li59, Z. X. Li15,
C. Liang42, H. Liang71,58, H. Liang1,63, Y. F. Liang54, Y. T. Liang31,63, G.
R. Liao14, L. Z. Liao50, Y. P. Liao1,63, J. Libby26, A. Limphirat60, D. X.
Lin31,63, T. Lin1, B. J. Liu1, B. X. Liu76, C. Liu34, C. X. Liu1, F. H. Liu53,
Fang Liu1, Feng Liu6, G. M. Liu56,i, H. Liu38,j,k, H. B. Liu15, H. M. Liu1,63,
Huanhuan Liu1, Huihui Liu21, J. B. Liu71,58, J. Y. Liu1,63, K. Liu38,j,k, K.
Y. Liu40, Ke Liu22, L. Liu71,58, L. C. Liu43, Lu Liu43, M. H. Liu12,f, P. L.
Liu1, Q. Liu63, S. B. Liu71,58, T. Liu12,f, W. K. Liu43, W. M. Liu71,58, X.
Liu38,j,k, X. Liu39, Y. Liu38,j,k, Y. Liu80, Y. B. Liu43, Z. A. Liu1,58,63, Z.
D. Liu9, Z. Q. Liu50, X. C. Lou1,58,63, F. X. Lu59, H. J. Lu23, J. G. Lu1,58,
X. L. Lu1, Y. Lu7, Y. P. Lu1,58, Z. H. Lu1,63, C. L. Luo41, M. X. Luo79, T.
Luo12,f, X. L. Luo1,58, X. R. Lyu63, Y. F. Lyu43, F. C. Ma40, H. Ma78, H. L.
Ma1, J. L. Ma1,63, L. L. Ma50, M. M. Ma1,63, Q. M. Ma1, R. Q. Ma1,63, X. T.
Ma1,63, X. Y. Ma1,58, Y. Ma46,g, Y. M. Ma31, F. E. Maas18, M. Maggiora74A,74C,
S. Malde69, A. Mangoni28B, Y. J. Mao46,g, Z. P. Mao1, S. Marcello74A,74C, Z.
X. Meng66, J. G. Messchendorp13,64, G. Mezzadri29A, H. Miao1,63, T. J. Min42,
R. E. Mitchell27, X. H. Mo1,58,63, B. Moses27, N. Yu. Muchnoi4,b, J.
Muskalla35, Y. Nefedov36, F. Nerling18,d, I. B. Nikolaev4,b, Z. Ning1,58, S.
Nisar11,l, Q. L. Niu38,j,k, W. D. Niu55, Y. Niu 50, S. L. Olsen63, Q.
Ouyang1,58,63, S. Pacetti28B,28C, X. Pan55, Y. Pan57, A. Pathak34, P.
Patteri28A, Y. P. Pei71,58, M. Pelizaeus3, H. P. Peng71,58, Y. Y. Peng38,j,k,
K. Peters13,d, J. L. Ping41, R. G. Ping1,63, S. Plura35, V. Prasad33, F. Z.
Qi1, H. Qi71,58, H. R. Qi61, M. Qi42, T. Y. Qi12,f, S. Qian1,58, W. B. Qian63,
C. F. Qiao63, J. J. Qin72, L. Q. Qin14, X. S. Qin50, Z. H. Qin1,58, J. F.
Qiu1, S. Q. Qu61, Z. H. Qu72, C. F. Redmer35, K. J. Ren39, A. Rivetti74C, M.
Rolo74C, G. Rong1,63, Ch. Rosner18, S. N. Ruan43, N. Salone44, A.
Sarantsev36,c, Y. Schelhaas35, K. Schoenning75, M. Scodeggio29A, K. Y.
Shan12,f, W. Shan24, X. Y. Shan71,58, J. F. Shangguan55, L. G. Shao1,63, M.
Shao71,58, C. P. Shen12,f, H. F. Shen1,8, W. H. Shen63, X. Y. Shen1,63, B. A.
Shi63, H. C. Shi71,58, J. L. Shi12, J. Y. Shi1, Q. Q. Shi55, R. S. Shi1,63, S.
Y. Shi72, X. Shi1,58, X. D. Shi71,58, J. J. Song19, T. Z. Song59, W. M.
Song34,1, Y. J. Song12, Y. X. Song46,g,m, S. Sosio74A,74C, S. Spataro74A,74C,
F. Stieler35, Y. J. Su63, G. B. Sun76, G. X. Sun1, H. Sun63, H. K. Sun1, J. F.
Sun19, K. Sun61, L. Sun76, S. S. Sun1,63, T. Sun51,e, W. Y. Sun34, Y. Sun9, Y.
J. Sun71,58, Y. Z. Sun1, Z. Q. Sun1,63, Z. T. Sun50, C. J. Tang54, G. Y.
Tang1, J. Tang59, Y. A. Tang76, L. Y. Tao72, Q. T. Tao25,h, M. Tat69, J. X.
Teng71,58, V. Thoren75, W. H. Tian59, Y. Tian31,63, Z. F. Tian76, I. Uman62B,
Y. Wan55, S. J. Wang 50, B. Wang1, B. L. Wang63, Bo Wang71,58, D. Y. Wang46,g,
F. Wang72, H. J. Wang38,j,k, J. P. Wang 50, K. Wang1,58, L. L. Wang1, M.
Wang50, Meng Wang1,63, N. Y. Wang63, S. Wang38,j,k, S. Wang12,f, T. Wang12,f,
T. J. Wang43, W. Wang59, W. Wang72, W. P. Wang71,58, X. Wang46,g, X. F.
Wang38,j,k, X. J. Wang39, X. L. Wang12,f, X. N. Wang1, Y. Wang61, Y. D.
Wang45, Y. F. Wang1,58,63, Y. L. Wang19, Y. N. Wang45, Y. Q. Wang1, Yaqian
Wang17, Yi Wang61, Z. Wang1,58, Z. L. Wang72, Z. Y. Wang1,63, Ziyi Wang63, D.
Wei70, D. H. Wei14, F. Weidner68, S. P. Wen1, Y. R. Wen39, U. Wiedner3, G.
Wilkinson69, M. Wolke75, L. Wollenberg3, C. Wu39, J. F. Wu1,8, L. H. Wu1, L.
J. Wu1,63, X. Wu12,f, X. H. Wu34, Y. Wu71, Y. H. Wu55, Y. J. Wu31, Z. Wu1,58,
L. Xia71,58, X. M. Xian39, B. H. Xiang1,63, T. Xiang46,g, D. Xiao38,j,k, G. Y.
Xiao42, S. Y. Xiao1, Y. L. Xiao12,f, Z. J. Xiao41, C. Xie42, X. H. Xie46,g, Y.
Xie50, Y. G. Xie1,58, Y. H. Xie6, Z. P. Xie71,58, T. Y. Xing1,63, C. F.
Xu1,63, C. J. Xu59, G. F. Xu1, H. Y. Xu66, Q. J. Xu16, Q. N. Xu30, W. Xu1, W.
L. Xu66, X. P. Xu55, Y. C. Xu77, Z. P. Xu42, Z. S. Xu63, F. Yan12,f, L.
Yan12,f, W. B. Yan71,58, W. C. Yan80, X. Q. Yan1, H. J. Yang51,e, H. L.
Yang34, H. X. Yang1, Tao Yang1, Y. Yang12,f, Y. F. Yang43, Y. X. Yang1,63,
Yifan Yang1,63, Z. W. Yang38,j,k, Z. P. Yao50, M. Ye1,58, M. H. Ye8, J. H.
Yin1, Z. Y. You59, B. X. Yu1,58,63, C. X. Yu43, G. Yu1,63, J. S. Yu25,h, T.
Yu72, X. D. Yu46,g, C. Z. Yuan1,63, J. Yuan34, L. Yuan2, S. C. Yuan1, Y.
Yuan1,63, Z. Y. Yuan59, C. X. Yue39, A. A. Zafar73, F. R. Zeng50, S. H.
Zeng72, X. Zeng12,f, Y. Zeng25,h, Y. J. Zeng59, Y. J. Zeng1,63, X. Y. Zhai34,
Y. C. Zhai50, Y. H. Zhan59, A. Q. Zhang1,63, B. L. Zhang1,63, B. X. Zhang1, D.
H. Zhang43, G. Y. Zhang19, H. Zhang71, H. C. Zhang1,58,63, H. H. Zhang34, H.
H. Zhang59, H. Q. Zhang1,58,63, H. Y. Zhang1,58, J. Zhang59, J. Zhang80, J. J.
Zhang52, J. L. Zhang20, J. Q. Zhang41, J. W. Zhang1,58,63, J. X. Zhang38,j,k,
J. Y. Zhang1, J. Z. Zhang1,63, Jianyu Zhang63, L. M. Zhang61, Lei Zhang42, P.
Zhang1,63, Q. Y. Zhang39,80, Shuihan Zhang1,63, Shulei Zhang25,h, X. D.
Zhang45, X. M. Zhang1, X. Y. Zhang50, Y. Zhang72, Y. T. Zhang80, Y. H.
Zhang1,58, Y. M. Zhang39, Yan Zhang71,58, Yao Zhang1, Z. D. Zhang1, Z. H.
Zhang1, Z. L. Zhang34, Z. Y. Zhang43, Z. Y. Zhang76, G. Zhao1, J. Y. Zhao1,63,
J. Z. Zhao1,58, Lei Zhao71,58, Ling Zhao1, M. G. Zhao43, R. P. Zhao63, S. J.
Zhao80, Y. B. Zhao1,58, Y. X. Zhao31,63, Z. G. Zhao71,58, A. Zhemchugov36,a,
B. Zheng72, J. P. Zheng1,58, W. J. Zheng1,63, Y. H. Zheng63, B. Zhong41, X.
Zhong59, H. Zhou50, J. Y. Zhou34, L. P. Zhou1,63, X. Zhou76, X. K. Zhou6, X.
R. Zhou71,58, X. Y. Zhou39, Y. Z. Zhou12,f, J. Zhu43, K. Zhu1, K. J.
Zhu1,58,63, L. Zhu34, L. X. Zhu63, S. H. Zhu70, S. Q. Zhu42, T. J. Zhu12,f, W.
J. Zhu12,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 Bochum Ruhr-University, D-44780 Bochum, Germany
4 Budker Institute of Nuclear Physics SB RAS (BINP), Novosibirsk 630090,
Russia
5 Carnegie Mellon University, Pittsburgh, Pennsylvania 15213, USA
6 Central China Normal University, Wuhan 430079, People’s Republic of China
7 Central South University, Changsha 410083, People’s Republic of China
8 China Center of Advanced Science and Technology, Beijing 100190, People’s
Republic of China
9 China University of Geosciences, Wuhan 430074, People’s Republic of China
10 Chung-Ang University, Seoul, 06974, Republic of Korea
11 COMSATS University Islamabad, Lahore Campus, Defence Road, Off Raiwind
Road, 54000 Lahore, Pakistan
12 Fudan University, Shanghai 200433, People’s Republic of China
13 GSI Helmholtzcentre for Heavy Ion Research GmbH, D-64291 Darmstadt, Germany
14 Guangxi Normal University, Guilin 541004, People’s Republic of China
15 Guangxi University, Nanning 530004, People’s Republic of China
16 Hangzhou Normal University, Hangzhou 310036, People’s Republic of China
17 Hebei University, Baoding 071002, People’s Republic of China
18 Helmholtz Institute Mainz, Staudinger Weg 18, D-55099 Mainz, Germany
19 Henan Normal University, Xinxiang 453007, People’s Republic of China
20 Henan University, Kaifeng 475004, People’s Republic of China
21 Henan University of Science and Technology, Luoyang 471003, People’s
Republic of China
22 Henan University of Technology, Zhengzhou 450001, People’s Republic of
China
23 Huangshan College, Huangshan 245000, People’s Republic of China
24 Hunan Normal University, Changsha 410081, People’s Republic of China
25 Hunan University, Changsha 410082, People’s Republic of China
26 Indian Institute of Technology Madras, Chennai 600036, India
27 Indiana University, Bloomington, Indiana 47405, USA
28 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
29 INFN Sezione di Ferrara, (A)INFN Sezione di Ferrara, I-44122, Ferrara,
Italy; (B)University of Ferrara, I-44122, Ferrara, Italy
30 Inner Mongolia University, Hohhot 010021, People’s Republic of China
31 Institute of Modern Physics, Lanzhou 730000, People’s Republic of China
32 Institute of Physics and Technology, Peace Avenue 54B, Ulaanbaatar 13330,
Mongolia
33 Instituto de Alta Investigación, Universidad de Tarapacá, Casilla 7D, Arica
1000000, Chile
34 Jilin University, Changchun 130012, People’s Republic of China
35 Johannes Gutenberg University of Mainz, Johann-Joachim-Becher-Weg 45,
D-55099 Mainz, Germany
36 Joint Institute for Nuclear Research, 141980 Dubna, Moscow region, Russia
37 Justus-Liebig-Universitaet Giessen, II. Physikalisches Institut, Heinrich-
Buff-Ring 16, D-35392 Giessen, Germany
38 Lanzhou University, Lanzhou 730000, People’s Republic of China
39 Liaoning Normal University, Dalian 116029, People’s Republic of China
40 Liaoning University, Shenyang 110036, People’s Republic of China
41 Nanjing Normal University, Nanjing 210023, People’s Republic of China
42 Nanjing University, Nanjing 210093, People’s Republic of China
43 Nankai University, Tianjin 300071, People’s Republic of China
44 National Centre for Nuclear Research, Warsaw 02-093, Poland
45 North China Electric Power University, Beijing 102206, People’s Republic of
China
46 Peking University, Beijing 100871, People’s Republic of China
47 Qufu Normal University, Qufu 273165, People’s Republic of China
48 Renmin University of China, Beijing 100872, 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 Yantai University, Yantai 264005, People’s Republic of China
78 Yunnan University, Kunming 650500, People’s Republic of China
79 Zhejiang University, Hangzhou 310027, People’s Republic of China
80 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 MOE 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
m Also at Ecole Polytechnique Federale de Lausanne (EPFL), CH-1015 Lausanne,
Switzerland

##  1 INTRODUCTION

Precision measurements of the elements of the Cabibbo-Kobayashi-Maskawa (CKM)
matrix and the test of the unitarity of the CKM triangle [1, 2] are essential
goals in the field of flavour physics. One of the three angles of the
unitarity triangle,
$\gamma(\phi_{3})\equiv\mathrm{arg}(-V_{ud}V_{ub}^{*}/V_{cd}V_{cb}^{*})$, can
be measured with the tree-level decay $B^{\pm}\to DK^{\pm}$ through the
interference between $B^{-}\to D^{0}K^{-}$ $(b\to c\bar{u}s)$ and
$B^{-}\to\bar{D}^{0}K^{-}$ $(b\to u\bar{c}s)$, which is one of the most
important measurements of LHCb and Belle II. Several approaches have been
proposed to measure the $\gamma$ angle via the decay $B^{\pm}\to DK^{\pm}$ [3,
4, 5]. Here, the relative magnitude and phase between the $D^{0}$ and
$\bar{D}^{0}$ decays into the same final states, and the $D^{0}$ decay
parameters ($e.g.$ $CP$-even fraction $F_{+}$) are the critical inputs. With
more data collected by the LHCb and Belle II experiments in the coming years,
the decay parameters of $D^{0}$ ($\bar{D}^{0}$) will become the dominant
source of uncertainty in the $\gamma$ measurement. Therefore, precision
measurements of the $D^{0}$ decay parameters are urgently required to improve
the precision of the $\gamma$ measurement.

The decay $D^{0}\to 4\pi$ is regarded as a sensitive mode to extract the
$\gamma$ angle via $B^{-}\to DK^{-}$. Their $CP$-even fractions ($F_{+}$) and
relative strong phase parameters in the different phase space (PHSP) bins
($c_{i}$/$s_{i}$) serve as the direct inputs in the GLW [3] and GGSZ [5]
methods, respectively. A reliable decay amplitude model of $D^{0}\to 4\pi$ is
critical to precisely extract $F_{+}$ and a model-independent $c_{i}$/$s_{i}$
[6], and to search for $CP$ violation in $D^{0}\to 4\pi$ [7]. Moreover, the
four-body $D^{0}$ hadronic decays provide an excellent platform to study the
two-body decays $D^{0}\to VV$ and $D^{0}\to AP$, where $V$, $A$, and $P$
denote vector, axial-vector, and pseudo-scalar mesons, respectively. These
decays thus enhance the understanding of the decay dynamics of the $D^{0}$
meson [8, 9].

Currently, there are only a limited number of experimental studies of
$D^{0}\to 4\pi$. The FOCUS experiment performed an amplitude analysis of
$D^{0}\to\pi^{+}\pi^{-}\pi^{+}\pi^{-}$ based on $\sim 6000$ candidate events
with a background fraction of $\sim 10\%$ [10]. An amplitude analysis of
$D^{0}\to\pi^{+}\pi^{-}\pi^{+}\pi^{-}$ was also carried out with the CLEO-c
data, which contain $\sim 7000$ candidate events with a background fraction of
$\sim 20\%$ [11]. However, no amplitude analysis of
$D^{0}\to\pi^{+}\pi^{-}\pi^{0}\pi^{0}$ has been performed yet. BESIII has
collected 2.93 $\rm{fb}^{-1}$ of $e^{+}e^{-}$ collision data at the center-of-
mass energy $\sqrt{s}=3.773$ GeV, where the $D\bar{D}$ is product by pair
without any additional hadrons. This data sample provide an ideal environment
for studying $D$ meson decays with the double tag (DT) technique [12, 13]. In
this method, a single tag (ST) candidate requires that only one $D$ meson is
reconstructed via ST mode. A DT candidate requires that both of $D$ and
$\bar{D}$ are reconstructed via signal mode and ST mode, respectively. Based
on this data sample and DT method with three ST modes $\bar{D}^{0}\to
K^{+}\pi^{-}$, $\bar{D}^{0}\to K^{+}\pi^{-}\pi^{0}$ and $\bar{D}^{0}\to
K^{+}\pi^{-}\pi^{+}\pi^{-}$, we report a joint amplitude analysis of
$D^{0}\to\pi^{+}\pi^{-}\pi^{+}\pi^{-}$ and
$D^{0}\to\pi^{+}\pi^{-}\pi^{0}\pi^{0}$(non-$\eta$). Furthermore, we determine
the model-dependent $CP$-even fractions, the absolute branching fractions, and
the fractions of individual components. Throughout this paper, the charge-
conjugated processes are always implied.

##  2 BESIII DETECTOR AND MONTE CARLO SIMULATION

The BESIII detector [14] records symmetric $e^{+}e^{-}$ collisions provided by
the BEPCII storage ring [15] in the center-of-mass energy range from 2.0 to
4.95 GeV, with a peak luminosity of $1\times
10^{33}\;\text{cm}^{-2}\text{s}^{-1}$ achieved at $\sqrt{s}=3.77\;\text{GeV}$.
BESIII has collected large data samples in this energy region [16]. The
cylindrical core of the $\rm{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 ${\rm d}E/{\rm d}x$ 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 that in the end cap
region is 110 ps.

Simulated data samples produced with a geant4-based [17] Monte Carlo (MC)
package, which includes the geometric description [18] 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 [19, 20]. The inclusive MC sample including the production of $D\bar{D}$
pairs (including quantum coherence for the neutral $D$ channels), the
non-$D\bar{D}$ decays of the $\psi(3770)$, the ISR production of the $J/\psi$
and $\psi(3686)$ states, and the continuum processes incorporated in kkmc are
generated to estimate the background and ST efficiencies. All particle decays
are modeled with evtgen [21, 22] using branching fractions either taken from
the Particle Data Group (PDG) [23], when available, or otherwise estimated
with lundcharm [24, 25]. Final state radiation from charged final state
particles is incorporated using the photos package [26]. The PHSP signal MC
samples of $D^{0}\to\pi^{+}\pi^{-}\pi^{+}\pi^{-}$ and
$D^{0}\to\pi^{+}\pi^{-}\pi^{0}\pi^{0}$ are generated in PHSP uniformly and are
used to calculate the normalization factor of probability density function
(PDF) in the amplitude analysis. The signal MC samples of
$D^{0}\to\pi^{+}\pi^{-}\pi^{+}\pi^{-}$ and
$D^{0}\to\pi^{+}\pi^{-}\pi^{0}\pi^{0}$(non-$\eta$) are generated according to
the amplitude analysis results and are used to estimate the DT efficiencies.

##  3 EVENT SELECTION

Charged tracks detected in the MDC are required to be within a polar angle
($\theta$) range of $|\rm{cos\theta}|<0.93$, where $\theta$ is defined with
respect to the $z$ axis, which is the symmetry axis of the MDC. The distance
of closest approach of these charged tracks to the interaction point must be
less than 10 cm along the $z$ axis, and less than 1 cm in the transverse
plane. Particle identification (PID) for charged tracks combines measurements
of the specific ionization energy loss in the MDC (d$E$/d$x$) and the flight
time in the TOF to form likelihoods $\mathcal{L}(h)~{}(h=K,\pi)$ for each
hadron $h$ hypothesis. The charged kaons and pions are identified by comparing
the likelihoods for the kaon and pion hypotheses,
$\mathcal{L}(K)>\mathcal{L}(\pi)$ and $\mathcal{L}(\pi)>\mathcal{L}(K)$,
respectively.

Photon candidates are identified using showers in the EMC. The deposited
energy of each shower must be $>25$ MeV in the barrel region
($|\cos\theta|<0.80$) and $>50$ MeV in the end cap region
($0.86<|\cos\theta|<0.92$). To exclude showers that originate from charged
tracks, the angle subtended by the EMC shower and the position of the closest
charged track at the EMC must be $>10^{\circ}$ as measured from the
interaction point. To suppress electronic noises and reject showers unrelated
to the event, the difference between the EMC time and the event start time is
required to be within [0, 700] ns. The $\pi^{0}$ candidates are reconstructed
from pairs of photon candidates with invariant mass being in the interval
(0.115, 0.150) GeV/$c^{2}$. To improve the momentum resolution, a kinematic
fit constraining the two-photon invariant mass to the known $\pi^{0}$ mass
[23] is performed, and the four-momenta updated by this kinematic fit are used
in the following analysis.

The signal candidates of $D^{0}\to\pi^{+}\pi^{-}\pi^{+}\pi^{-}$ and
$D^{0}\to\pi^{+}\pi^{-}\pi^{0}\pi^{0}$ are selected with the DT method [12,
13]. First, the ST $\bar{D}^{0}$ mesons are reconstructed with the three
hadronic decay modes $\bar{D}^{0}\to K^{+}\pi^{-}$, $\bar{D}^{0}\to
K^{+}\pi^{-}\pi^{0}$, and $\bar{D}^{0}\to K^{+}\pi^{-}\pi^{+}\pi^{-}$. Two
kinematic variables, the energy difference with respect to the beam energy
$\Delta E$ and the beam energy constrained mass $M_{\rm{bc}}$ are defined as

$\displaystyle\Delta E=$ $\displaystyle E_{\bar{D}^{0}}-E_{\rm{beam}},$ (1)
$\displaystyle M_{\rm{bc}}=$
$\displaystyle\sqrt{E_{\rm{beam}}^{2}-\vec{p}_{\bar{D}^{0}}^{2}}~{},$ (2)

where $E_{\rm{beam}}$ is the beam energy, and $E_{\bar{D}^{0}}$ and
$\vec{p}_{\bar{D}^{0}}$ are the energy and momentum of the ST $\bar{D}^{0}$
candidate in the $e^{+}e^{-}$ center-of-mass frame. For multiple $\bar{D}^{0}$
candidates in each ST mode, only the one with the smallest $|\Delta E|$ is
kept for further analysis. To reject the backgrounds from cosmic rays and
Bhabha events in the ST mode $\bar{D}^{0}\to K^{+}\pi^{-}$, the requirements
described in Ref. [27] are applied. To reject the peaking background
$\bar{D}^{0}\to K^{+}K_{S}^{0}\pi^{-}$ in the ST mode $\bar{D}^{0}\to
K^{+}\pi^{-}\pi^{+}\pi^{-}$, the events with $|M(\pi^{+}\pi^{-})-0.4976|<0.03$
GeV/$c^{2}$ are vetoed. To further reject combinatorial backgrounds, the
$\bar{D}^{0}$ candidates are required to have $\Delta E$ within a given
interval defined in Table 1, about 3 times the resolution for each ST mode. In
the sub-sample containing ST candidates, the signal candidates of
$D^{0}\to\pi^{+}\pi^{-}\pi^{+}\pi^{-}$ and
$D^{0}\to\pi^{+}\pi^{-}\pi^{0}\pi^{0}$ are reconstructed with the $\pi^{\pm}$
and $\pi^{0}$ candidates which have not been used in the ST side (namely DT
thereafter). Similar kinematic variables $\Delta E$ and $M_{\rm{bc}}$ are
formed for signal $D^{0}$ candidates, and the corresponding $\Delta E$
requirements are listed in Table 1. For multiple signal $D^{0}$ candidates,
only the one giving the smallest $|\Delta E|$ is kept.

Table 1: The $\Delta E$ requirements for different decay modes. Decay mode |  $\Delta E$ (GeV)
---|---
$\bar{D}^{0}\to K^{+}\pi^{-}$ | $(-0.027,0.026)$
$\bar{D}^{0}\to K^{+}\pi^{-}\pi^{0}$ | $(-0.057,0.043)$
$\bar{D}^{0}\to K^{+}\pi^{-}\pi^{+}\pi^{-}$ | $(-0.020,0.018)$
$D^{0}\to\pi^{+}\pi^{-}\pi^{+}\pi^{-}$ | $(-0.032,0.028)$
$D^{0}\to\pi^{+}\pi^{-}\pi^{0}\pi^{0}$ | $(-0.066,0.041)$

To improve the purity of signal candidates in the amplitude analysis, some
further selection criteria are applied. The studies based on the inclusive MC
sample indicate that most backgrounds from the continuum process
$e^{+}e^{-}\to q\bar{q}$ include a $K_{S}^{0}$ in the final state. Therefore,
common and secondary vertex fits are performed on all $\pi^{+}\pi^{-}$ pairs
in the event to reconstruct the $K_{S}^{0}$ candidate. The candidate events
are rejected if there exists any $K_{S}^{0}$ candidate with $\pi^{+}\pi^{-}$
invariant mass within the interval $|M(\pi^{+}\pi^{-})-0.4976|<0.03~{}{\rm
GeV}/c^{2}$ with the decay length greater than twice its resolution. The MC
studies also show that most backgrounds from $e^{+}e^{-}\to D^{+}D^{-}$
contain the decay $D^{-}\to K^{+}\pi^{-}\pi^{-}$ due to its large branching
fraction and a similar topology as the signal. The candidate events are
rejected if there exists any $K^{+}\pi^{-}\pi^{-}$ combinations with
$1.863<M_{\rm{bc}}(K^{+}\pi^{-}\pi^{-})<1.878~{}{\rm GeV}/c^{2}$ and $|\Delta
E(K^{+}\pi^{-}\pi^{-})|<0.03~{}\rm{GeV}$. There is also the background from
the process $e^{+}e^{-}\to D^{0}\bar{D}^{0}$. Events with any
$\pi^{+}\pi^{-}\pi^{0}$ combinations satisfying
$1.859<M_{\rm{bc}}(\pi^{+}\pi^{-}\pi^{0})<1.873~{}{\rm GeV}/c^{2}$ and
$-0.057<\Delta E(\pi^{+}\pi^{-}\pi^{0})<0.043~{}\rm{GeV}$ are rejected to
eliminate the background of $D^{0}\to\pi^{+}\pi^{-}\pi^{0}$ in the signal
process $D^{0}\to\pi^{+}\pi^{-}\pi^{0}\pi^{0}$. All of the above backgrounds
do not form peaks in the $M_{\rm bc}$ distribution of the signal side. There
are also some backgrounds which have the same final states as those in the
signal mode and produce peaks in the $M_{\rm bc}$ distributions of both the ST
and signal sides. To reject the peaking backgrounds $D^{0}\to
K_{S}^{0}(\rightarrow\pi^{+}\pi^{-})\pi^{+}\pi^{-}$ and $D^{0}\to
K^{-}(\rightarrow\pi^{+}\pi^{-}\pi^{-})\pi^{+}$ in the signal process
$D^{0}\to\pi^{+}\pi^{-}\pi^{+}\pi^{-}$, events with
$|M(\pi^{+}\pi^{-})-0.4976|<0.03~{}{\rm GeV}/c^{2}$ or
$M(\pi^{+}\pi^{-}\pi^{-})<0.51~{}{\rm GeV}/c^{2}$ are vetoed. To reject the
peaking backgrounds $D^{0}\to
K_{S}^{0}(\rightarrow\pi^{0}\pi^{0})\pi^{+}\pi^{-}$, $D^{0}\to
K_{S}^{0}(\rightarrow\pi^{+}\pi^{-})\pi^{0}\pi^{0}$ and $D^{0}\to
K^{-}(\rightarrow\pi^{-}\pi^{0})\pi^{+}\pi^{0}$ in the signal process
$D^{0}\to\pi^{+}\pi^{-}\pi^{0}\pi^{0}$, events with
$0.4376<M(\pi^{0}\pi^{0})<0.5276~{}{\rm GeV}/c^{2}$ or
$|M(\pi^{+}\pi^{-})-0.4976|<0.03~{}{\rm GeV}/c^{2}$ or
$0.4677<M(\pi^{-}\pi^{0})<0.5067~{}{\rm GeV}/c^{2}$ are vetoed. Since the
interference between $D^{0}\to\pi^{0}\eta(\rightarrow\pi^{+}\pi^{-}\pi^{0})$
and $D^{0}\to\pi^{+}\pi^{-}\pi^{0}\pi^{0}$ is negligible, events with
$M(\pi^{+}\pi^{-}\pi^{0})<0.57~{}{\rm GeV}/c^{2}$ are also vetoed in the
signal process $D^{0}\to\pi^{+}\pi^{-}\pi^{0}\pi^{0}$. Moreover, we reject
events if the energy of the photon from $\pi^{0}$ decays is below 50 MeV, as
this indicates an incorrect choice of a soft photon.

##  4 BACKGROUND STUDY AND SIGNAL EXTRACTION

With all the above selection criteria, the two-dimensional (2D) distributions
of $M_{\rm bc}^{\rm sig}$ versus $M_{\rm bc}^{\rm tag}$ of the survived events
are shown in the left columns of Figs. 2 and 3, where $M_{\rm bc}^{\rm sig}$
and $M_{\rm bc}^{\rm tag}$ are for the signal and tag sides, respectively.
Typically, the signal events are accumulated around the intersection of
$M_{\rm bc}^{\rm sig}$ and $M_{\rm bc}^{\rm tag}$ at the $D^{0}$ nominal mass.
The background, which arises from miscombinations both in the signal and tag
sides (namely BKGI thereafter), is located at the diagonal band. There are
also the backgrounds from $e^{+}e^{-}\to D^{0}\bar{D}^{0}$ but with the wrong
reconstruction of the signal side (namely BKGII) or tag side (namely BKGIII),
distributed as the vertical and horizontal bands, respectively. Detailed MC
studies indicate that most of the backgrounds from $e^{+}e^{-}\to
D^{0}\bar{D}^{0}$ do not form individual peaks in the distribution of $M_{\rm
bc}^{\rm sig}$ or $M_{\rm bc}^{\rm tag}$. However, there are backgrounds from
$e^{+}e^{-}\to D^{0}\bar{D}^{0}$ with $\bar{D}^{0}\to K^{+}\pi^{-}\pi^{0}$
(namely BKGIV) or $D^{0}\to\pi^{+}\pi^{-}\pi^{0}\pi^{0}$ (namely BKGV) with
the wrong reconstruction of $\pi^{0}$, which are exactly the signal but with
the wrong reconstruction of $\pi^{0}$ and produce a relatively broad peak in
the $M_{\rm bc}^{\rm tag}$ or $M_{\rm bc}^{\rm sig}$ distribution,
respectively. There are also the backgrounds of $D^{0}\to
K_{S}^{0}\pi^{+}\pi^{-}$, $D^{0}\to K_{S}^{0}\pi^{0}\pi^{0}$, $D^{0}\to
K^{-}\pi^{+}\pi^{0}$ and $D^{0}\to K_{S}^{0}\eta^{\prime}$ (namely BKGVI),
which have analogous final states as the signal with the specific decay modes
of $K_{S}^{0}\to\pi^{+}\pi^{-}$, $K_{S}^{0}\to\pi^{0}\pi^{0}$,
$K^{-}\to\pi^{-}\pi^{0}$ and $\eta^{\prime}\to\gamma\pi^{+}\pi^{-}$. These
backgrounds are not directly distinguished from the signal in the $M_{\rm
bc}^{\rm sig}$ and $M_{\rm bc}^{\rm tag}$ distributions and are estimated by
MC simulation. The corresponding yields are summarized in Table 2.

Table 2: The estimated numbers of peaking background events. The uncertainties include the statistical uncertainties of the estimated background yields and the uncertainties of the quoted branching fractions of different background processes. Signal mode | $\pi^{+}\pi^{-}\pi^{+}\pi^{-}$ | $\pi^{+}\pi^{-}\pi^{0}\pi^{0}$
---|---|---
Tag mode | $K^{+}\pi^{-}$ | $K^{+}\pi^{-}\pi^{0}$ | $K^{+}\pi^{-}\pi^{+}\pi^{-}$ | $K^{+}\pi^{-}$ | $K^{+}\pi^{-}\pi^{0}$ | $K^{+}\pi^{-}\pi^{+}\pi^{-}$
$N(K_{S}^{0}\pi^{+}\pi^{-})$ | $24.7\pm 2.3$ | $42.4\pm 4.3$ | $22.9\pm 2.6$ | $12.6\pm 3.1$ | $27.7\pm 6.2$ | $29.2\pm 6.2$
$N(K_{S}^{0}\pi^{0}\pi^{0})$ | - | - | - | $1.7\pm 0.4$ | $4.3\pm 0.9$ | $2.4\pm 0.6$
$N(K^{-}\pi^{+}\pi^{0})$ | - | - | - | $12.6\pm 1.2$ | $23.1\pm 2.2$ | $24.0\pm 1.8$
$N(K_{S}^{0}\eta^{\prime})$ | - | - | - | $3.0\pm 1.2$ | $3.9\pm 1.7$ | $2.8\pm 1.2$

To extract the ST and DT yields, unbinned maximum likelihood fits are
performed on the $M_{\rm bc}^{\rm tag}$ distribution of the remaining ST
candidates and the 2D distributions of $M_{\rm bc}^{\rm sig}$ versus $M_{\rm
bc}^{\rm tag}$ of the remaining DT events in individual tag modes. In the fit
to the $M_{\rm bc}^{\rm tag}$ distribution, the signal shape is described by
the MC simulated shape of the truth-matched ST events, and the background
shape is described by an ARGUS function with the cut-off parameter fixed at
$1.8865$ GeV/$c^{2}$ [28]. For the $\bar{D}^{0}\to K^{+}\pi^{-}\pi^{0}$ tag
mode, an additional peaking background is considered which represents the
wrong reconstruction of $\pi^{0}$ as discussed above. This peaking background
shape is described by the MC-simulated shape of the truth-matched ST events
convolved with a bifurcated Gaussian function with fixed parameters that are
obtained from the fit to the corresponding simulated background events. The
results of the fits to the $M_{\rm bc}^{\rm tag}$ distributions are shown in
Fig. 1.

In the fit to the 2D distribution of $M_{\rm bc}^{\rm sig}$ (labelled as
$y_{1}$) versus $M_{\rm bc}^{\rm tag}$ (labelled as $y_{2}$), the signal is
described by

$\displaystyle S(y_{1},y_{2})\otimes
G(y_{1};\mu_{y_{1}},\sigma_{y_{1}})\otimes
G(y_{2};\mu_{y_{2}},\sigma_{y_{2}}),$

where $S(y_{1},y_{2})$ is the signal MC shape derived from the truth-matched
events, and $G(\mu,\sigma)$ is the Gaussian function describing the resolution
difference between data and MC simulation by fixing the parameters
$\mu_{y_{1(2)}}$ and $\sigma_{y_{1(2)}}$ to those from the fits on the one-
dimensional $M_{\rm bc}$ distributions.

The BKGI is described by

$\displaystyle T(y_{1}-y_{2};\mu,\sigma(y_{1}+y_{2}),n)\times
A(y_{1};m_{y_{1}},z^{\prime}_{y_{1}},\rho^{\prime}_{y_{1}})$
$\displaystyle\times
A(y_{2};m_{y_{2}},z^{\prime}_{y_{2}},\rho^{\prime}_{y_{2}}),$

where $T$ is the student’s function defined as
$T(y;\mu,\sigma,n)=\frac{\Gamma(n/2+0.5)}{\sigma\sqrt{n\pi}\Gamma(n/2)}[1+\frac{1}{2}(\frac{y-\mu}{\sigma})^{2}]^{-\frac{n+1}{2}}$,
$\sigma(y_{1}+y_{2})=\sigma_{0}+\sigma_{1}(y_{1}+y_{2}-m_{y_{1}}-m_{y_{2}})$
and $A$ is the ARGUS function defined as
$A(y;m,z,\rho)=y(1-\frac{y^{2}}{m^{2}})^{\rho}e^{z(1-\frac{y^{2}}{m^{2}})}$.
In the fit, only the cut-off parameter $m_{y_{1(2)}}$ is fixed to
$1.8865~{}{\rm GeV}/c^{2}$, the other parameters are free to vary.

The BKGII and BKGIII are described by

$\displaystyle[S(y_{2(1)})\otimes
G(y_{2(1)};\mu_{y_{2(1)}},\sigma_{y_{2(1)}})]$ $\displaystyle\times
A(y_{1(2)};m_{y_{1(2)}},z_{y_{1(2)}},\rho_{y_{1(2)}}),$

where $A$ is the ARGUS function as described above, $S(y_{2(1)})$ is the
projection of $S(y_{1},y_{2})$ on $y_{2(1)}$ and $G$ is the same Gaussian
function as described above. In the fit, the parameter $m_{y_{1(2)}}$ is fixed
to $1.8865~{}{\rm GeV}/c^{2}$, $\rho_{y_{1(2)}}$ is fixed to the values
obtained from the fits to the inclusive MC sample, and $z_{y_{1(2)}}$ is a
free parameter.

The BKGIV is described by

$\displaystyle B_{1}(y_{1},y_{2})\otimes
G(y_{1};\mu_{y_{1}},\sigma_{y_{1}})\otimes
G(y_{2};\mu_{y_{2}},\sigma_{y_{2}}),$

where $B_{1}(y_{1},y_{2})$ is the truth-matched signal MC shape convolved with
a bifurcated Gaussian function in $y_{2}$ and $G$ is the same Gaussian
function as described above. In the fit, the parameters of the bifurcated
Gaussian function are fixed to the values obtained from the fit to the $M_{\rm
bc}^{\rm tag}$ distribution of the corresponding simulated background events.

The BKGV in $D^{0}\to\pi^{+}\pi^{-}\pi^{0}\pi^{0}$ is described by

$\displaystyle B_{2}(y_{1},y_{2})\otimes
G(y_{1};\mu_{y_{1}},\sigma_{y_{1}})\otimes
G(y_{2};\mu_{y_{2}},\sigma_{y_{2}}),$

where $B_{2}(y_{1},y_{2})$ is the truth-matched signal MC shape convolved with
a bifurcated Gaussian function in $y_{1}$ and $G$ is the same Gaussian
function as described above. In the fit, the parameters of the bifurcated
Gaussian function are fixed to the values obtained from the fit to the $M_{\rm
bc}^{\rm sig}$ distribution of the corresponding simulated background events.

The BKGVI is described by

$\displaystyle B_{3}(y_{1},y_{2})\otimes
G(y_{1};\mu_{y_{1}},\sigma_{y_{1}})\otimes
G(y_{2};\mu_{y_{2}},\sigma_{y_{2}}),$

where $B_{3}(y_{1},y_{2})$ is the truth-matched MC shape convolved with a
bifurcated Gaussian function in $y_{1}$ and $G$ are the Gaussian functions. In
the fit, the parameters of the bifurcated Gaussian functions are fixed to the
value obtained from the fit to the $M_{\rm bc}^{\rm sig}$ distribution of the
corresponding simulated background events from the inclusive MC sample, and
the yields are fixed to those summarized in Table 2.

The projections of the 2D fits on $M_{\rm bc}^{\rm tag}$ and $M_{\rm bc}^{\rm
sig}$ are shown in Figs. 2 and 3, respectively. The obtained ST and DT yields
are summarized in Table 3.

\begin{overpic}[width=138.76157pt]{fig/ST/Kpi.pdf}
\end{overpic}\begin{overpic}[width=138.76157pt]{fig/ST/Kpipi0.pdf}
\end{overpic}\begin{overpic}[width=138.76157pt]{fig/ST/K3pi.pdf} \end{overpic}

Figure 1: The $M_{\rm bc}^{\rm tag}$ distributions of the ST candidates for
$\bar{D}^{0}\to K^{+}\pi^{-}$ (left), $\bar{D}^{0}\to K^{+}\pi^{-}\pi^{0}$
(middle) and $\bar{D}^{0}\to K^{+}\pi^{-}\pi^{+}\pi^{-}$ (right). The dots
with error bars are data, the blue solid curves are the total fit results, and
the red and blue dashed curves are the signal and background, respectively.

\begin{overpic}[width=138.76157pt]{fig/PipPimPipPim/DTfit/Kpi2D.pdf}
\end{overpic}\begin{overpic}[width=138.76157pt]{fig/PipPimPipPim/DTfit/KpiST.pdf}
\end{overpic}\begin{overpic}[width=138.76157pt]{fig/PipPimPipPim/DTfit/KpiDT.pdf}
\end{overpic}\begin{overpic}[width=138.76157pt]{fig/PipPimPipPim/DTfit/Kpipi02D.pdf}
\end{overpic}\begin{overpic}[width=138.76157pt]{fig/PipPimPipPim/DTfit/Kpipi0ST.pdf}
\end{overpic}\begin{overpic}[width=138.76157pt]{fig/PipPimPipPim/DTfit/Kpipi0DT.pdf}
\end{overpic}\begin{overpic}[width=138.76157pt]{fig/PipPimPipPim/DTfit/K3pi2D.pdf}
\end{overpic}\begin{overpic}[width=138.76157pt]{fig/PipPimPipPim/DTfit/K3piST.pdf}
\end{overpic}\begin{overpic}[width=138.76157pt]{fig/PipPimPipPim/DTfit/K3piDT.pdf}
\end{overpic}

Figure 2: The 2D distributions (left) of $M_{\rm bc}^{\rm tag}$ versus $M_{\rm
bc}^{\rm sig}$ and the projections on $M_{\rm bc}^{\rm tag}$ (middle) and
$M_{\rm bc}^{\rm sig}$ (right) of the 2D fits on the DT events tagged by
$\bar{D}^{0}\to K^{+}\pi^{-}$ (first row), $\bar{D}^{0}\to
K^{+}\pi^{-}\pi^{0}$ (second row), and $\bar{D}^{0}\to
K^{+}\pi^{-}\pi^{+}\pi^{-}$ (third row) in the
$D^{0}\to\pi^{+}\pi^{-}\pi^{+}\pi^{-}$ decay. The dots with error bars are
data. The blue solid curves are the total fit results, and the red dashed
curves show the signal. The green, blue, orange, and magenta dashed curves are
BKGI, BKGIII, BKGVI, and the combination of BKGII and BKGV, respectively.

\begin{overpic}[width=138.76157pt]{fig/PipPim2Pi0/DTfit/Kpi2D.pdf}
\end{overpic}\begin{overpic}[width=138.76157pt]{fig/PipPim2Pi0/DTfit/KpiST.pdf}
\end{overpic}\begin{overpic}[width=138.76157pt]{fig/PipPim2Pi0/DTfit/KpiDT.pdf}
\end{overpic}\begin{overpic}[width=138.76157pt]{fig/PipPim2Pi0/DTfit/Kpipi02D.pdf}
\end{overpic}\begin{overpic}[width=138.76157pt]{fig/PipPim2Pi0/DTfit/Kpipi0ST.pdf}
\end{overpic}\begin{overpic}[width=138.76157pt]{fig/PipPim2Pi0/DTfit/Kpipi0DT.pdf}
\end{overpic}\begin{overpic}[width=138.76157pt]{fig/PipPim2Pi0/DTfit/K3pi2D.pdf}
\end{overpic}\begin{overpic}[width=138.76157pt]{fig/PipPim2Pi0/DTfit/K3piST.pdf}
\end{overpic}\begin{overpic}[width=138.76157pt]{fig/PipPim2Pi0/DTfit/K3piDT.pdf}
\end{overpic}

Figure 3: The 2D distributions (left) of $M_{\rm bc}^{\rm tag}$ versus $M_{\rm bc}^{\rm sig}$ and the projections on $M_{\rm bc}^{\rm tag}$ (middle) and $M_{\rm bc}^{\rm sig}$ (right) of the 2D fits on the DT events tagged by $\bar{D}^{0}\to K^{+}\pi^{-}$ (first row), $\bar{D}^{0}\to K^{+}\pi^{-}\pi^{0}$ (second row) and $\bar{D}^{0}\to K^{+}\pi^{-}\pi^{+}\pi^{-}$ (third row) in the $D^{0}\to\pi^{+}\pi^{-}\pi^{0}\pi^{0}$ decay. The dots with error bars are data. The blue solid curves are the total fit results, and the red dashed curves show the signal. The green dashed curves are the combination of BKGI and BKGIV, the magenta dashed curves are the combination of BKGII and BKGV, the blue and orange dashed curves are BKGIII and BKGVI, respectively. Table 3: The ST and DT yields, efficiencies, and quantum correlated correction factors for three tag modes. Tag mode | $\bar{D}^{0}\to K^{+}\pi^{-}$ | $\bar{D}^{0}\to K^{+}\pi^{-}\pi^{0}$ | $\bar{D}^{0}\to K^{+}\pi^{-}\pi^{+}\pi^{-}$
---|---|---|---
$N^{\rm ST}$ | $549586\pm 778$ | $914531\pm 1321$ | $600316\pm 856$
$\epsilon^{\rm ST}$ | $0.6762\pm 0.0004$ | $0.2953\pm 0.0001$ | $0.3365\pm 0.0002$
$N^{\rm DT}_{\pi^{+}\pi^{-}\pi^{+}\pi^{-}}$ | $1719\pm 42$ | $2560\pm 56$ | $1520\pm 43$
$\epsilon^{\rm DT}_{\pi^{+}\pi^{-}\pi^{+}\pi^{-}}$ | $0.2792\pm 0.0006$ | $0.1187\pm 0.0002$ | $0.1221\pm 0.0003$
$N^{\rm DT}_{\pi^{+}\pi^{-}\pi^{0}\pi^{0}}$ (non-$\eta$) | $721\pm 32$ | $917\pm 39$ | $562\pm 29$
$\epsilon^{\rm DT}_{\pi^{+}\pi^{-}\pi^{0}\pi^{0}}$ (non-$\eta$) | $0.0818\pm 0.0003$ | $0.0319\pm 0.0001$ | $0.0334\pm 0.0001$
$\frac{2rR\rm{cos}\delta}{1+r^{2}}$ |  $-0.114\pm 0.004$ [29] |  $-0.067\pm 0.005$ [30] |  $-0.046^{+0.013}_{-0.011}$ [30]

##  5 AMPLITUDE FORMULA

To improve the resolutions of kinematic variables, a one-constraint kinematic
fit with the hypothesis of $D^{0}\to 4\pi$ by constraining the $4\pi$
invariant mass to the known $D^{0}$ mass [23] is performed on the candidate
events, and the updated kinematic variables are used in the amplitude
analysis. Using the GPUPWA framework [31], a joint amplitude analysis is
performed on the candidate events of $D^{0}\to\pi^{+}\pi^{-}\pi^{+}\pi^{-}$
and $D^{0}\to\pi^{+}\pi^{-}\pi^{0}\pi^{0}$(non-$\eta$). The general amplitude
of the $D^{0}\to f$ decay is given by

$\displaystyle A_{f}(p)=\sum_{i}\Lambda_{i}U_{i}(p)~{},$ (3)

where $U_{i}$ is the amplitude of the $i$-th intermediate process with a
complex coupling factor $\Lambda_{i}$ and $p$ is a set of four-momenta of
final states. For the $\bar{D}^{0}$ decay amplitude, assuming $CP$
conservation, the $CP$-conjugate PHSP $\bar{p}$ is defined by interchanging
the charge of the final state and the reversal of three-momenta. Then, the
amplitude of the $\bar{D}^{0}\to\bar{f}$ decay is given as

$\displaystyle\bar{A}_{\bar{f}}(p)=\sum_{i}\Lambda_{i}\bar{U}_{i}(p)=\sum_{i}\Lambda_{i}U_{i}(\bar{p})~{}.$
(4)

Since the $D^{0}\bar{D}^{0}$ pair is produced in the $\psi(3770)$ decay,
quantum correlation between $D^{0}$ and $\bar{D}^{0}$ needs to be considered.
By ignoring the effects of $CP$ violation and $D^{0}$-$\bar{D}^{0}$ mixing,
the observed differential cross section $|M_{f}(p)|^{2}$ of $D^{0}\to f$
together with the tag mode $\bar{D}^{0}\to g$ is

$\displaystyle|A_{f}(p)-r_{\bar{g}}R_{\bar{g}}e^{-i\delta_{\bar{g}}}\bar{A}_{f}(p)|^{2}+r_{\bar{g}}^{2}(1-R_{\bar{g}}^{2})|\bar{A}_{f}(p)|^{2}~{},$
(5)

where
$r_{g}^{2}=\frac{\int|\bar{A}_{g}|^{2}{\mathrm{d}\Phi_{g}}}{\int|{A}_{g}|^{2}{\mathrm{d}\Phi_{g}}}$
and $R_{g}e^{-i\delta_{g}}=\frac{\int
A_{g}^{*}\bar{A}_{g}{\mathrm{d}\Phi_{g}}}{\sqrt{\int|A_{g}|^{2}{\mathrm{d}\Phi_{g}}\int|\bar{A}_{g}|^{2}{\mathrm{d}\Phi_{g}}}}$.
The values of $r$, $R$, $\delta$ for the three tag modes are summarized in
Table 4. In practice, the second term in Eq. (5) is ignored due to the
relatively small value of $r_{\bar{g}}^{2}(1-R_{\bar{g}}^{2})$.

Table 4: The input values of $r$, $R$ and $\delta$ for the three tag modes. Mode | $K^{-}\pi^{+}$ | $K^{-}\pi^{+}\pi^{0}$ | $K^{-}\pi^{+}\pi^{-}\pi^{+}$
---|---|---|---
$r(\%)$ |  $5.86\pm 0.02$ [29] |  $4.41\pm 0.11$ [30] |  $5.50\pm 0.07$ [30]
$R$ | 1 |  $0.79\pm 0.04$ [30] |  $0.44_{-0.10}^{+0.09}$ [30]
$\delta(^{\circ})$ |  $192.1^{+8.6}_{-10.2}$ [29] |  $196\pm+11$ [30] |  $161^{+28}_{-18}$ [30]

The amplitude $U_{i}$ is constructed with the spin factor, Blatt-Weisskopf
barrier factors [32], and propagators of resonances. To construct $U_{i}$ of
the four-body $D^{0}$ decay, the isobar model is applied in which the decay is
factorized into subsequent two-body decay amplitudes [33, 34, 35]. The general
amplitude of the $i$-th intermediate process in a four-body decay is given by

$\displaystyle
U_{i}(p)=S_{i}(p)B_{L_{D}}(p)P_{R_{1}}(p)B_{L_{R_{1}}}(p)P_{R_{2}}(p)B_{L_{R_{2}}}(p)~{},$
(6)

where $S_{i}$ is the spin factor of the $i$-th decay amplitude, $B_{L_{X}}$
($X=D$, $R_{1}$, $R_{2}$) are the Blatt-Weisskopf barrier factors for the $D$
meson and the resonances $R_{1}$ and $R_{2}$, and $P_{R_{1}}$ and $P_{R_{2}}$
are the propagators of $R_{1}$ and $R_{2}$, respectively. The amplitude is
constructed with the exchange symmetry for indistinguishable pions.

The spin factor is constructed with the covariant Zemach (Rarita-Schwinger)
tensor formalism [36, 37, 38, 39] by combining pure-orbital-angular-momentum
covariant tensors $\tilde{t}^{(L)}_{\mu_{1}\cdots\mu_{l}}$ and the momenta of
parent particles together with Minkowski metric $g_{\mu\nu}$ and Levi-Civita
symbol $\epsilon_{\mu\nu\lambda\sigma}$. For a process $a\to bc$, the
covariant tensors $\tilde{t}^{(L)}_{\mu_{1}\cdots\mu_{l}}$ for the final
states of pure orbital angular momentum $L$ are

$\displaystyle\tilde{t}^{(L)}_{\mu_{1}\cdots\mu_{l}}=(-1)^{L}P^{(L)}_{\mu_{1}\cdots\mu_{L}\mu^{\prime}_{1}\cdots\mu^{\prime}_{L}}(p_{a})r^{\mu^{\prime}_{1}}\cdots
r^{\mu^{\prime}_{L}}~{},$ (7)

where $r=p_{b}-p_{c}$ and
$P^{(L)}_{\mu_{1}\cdots\mu_{L}\mu^{\prime}_{1}\cdots\mu^{\prime}_{L}}(p_{a})$
is the spin projection operator of the particle $a$,

$\displaystyle P^{(0)}(p_{a})~{}~{}~{}=$ $\displaystyle 1~{},$ (8)
$\displaystyle P^{(1)}_{\mu\mu^{\prime}}(p_{a})~{}~{}~{}=$ $\displaystyle-
g_{\mu\mu^{\prime}}+\frac{p_{a\mu}p_{a\mu^{\prime}}}{p_{a}^{2}}~{},$ (10)
$\displaystyle P^{(2)}_{\mu\nu\mu^{\prime}\nu^{\prime}}(p_{a})=$
$\displaystyle\frac{1}{2}[P^{(1)}_{\mu\mu^{\prime}}(p_{a})P^{(1)}_{\nu\nu^{\prime}}(p_{a})+P^{(1)}_{\mu\nu^{\prime}}(p_{a})P^{(1)}_{\nu\mu^{\prime}}(p_{a})]$
$\displaystyle-\frac{1}{3}P^{(1)}_{\mu\nu}(p_{a})P^{(1)}_{\mu^{\prime}\nu^{\prime}}(p_{a})~{}.$

Following the isobar model, the spin factors of the four-body decay $D^{0}\to
P_{1}P_{2}P_{3}P_{4}$ are summarized in Table 5.

Table 5: Summary of spin factors in this analysis, where $S$, $P$, $V$, $A$, $T$ and $PT$ denote scalar, pseudo-scalar, vector, axial-vector, tensor, and pseudo-tensor mesons, respectively. $[S]$, $[P]$, and $[D]$ represent orbital angular momenta $L$ = 0, 1, and 2 in the decays, respectively. Decay chain | Spin factor
---|---
$D[S]\to P~{}P_{1},~{}P[S]\to S~{}P_{2},~{}S[S]\to P_{3}~{}P_{4}$ | 1
$D[S]\to P~{}P_{1},~{}P[P]\to V~{}P_{2},~{}V[P]\to P_{3}~{}P_{4}$ | $\tilde{t}_{\mu}(P)\tilde{t}^{\mu}(V)$
$D[S]\to P~{}P_{1},~{}P[D]\to T~{}P_{2},~{}T[D]\to P_{3}~{}P_{4}$ | $\tilde{t}_{\mu\nu}(P)\tilde{t}^{\mu\nu}(T)$
$D[P]\to A~{}P_{1},~{}A[S]\to V~{}P_{2},~{}V[P]\to P_{3}~{}P_{4}$ | $\tilde{t}_{\mu}(D)P^{\mu\nu}(A)\tilde{t}_{\nu}(V)$
$D[P]\to A~{}P_{1},~{}A[D]\to V~{}P_{2},~{}V[P]\to P_{3}~{}P_{4}$ | $\tilde{t}_{\mu}(D)\tilde{t}^{\mu\nu}(A)\tilde{t}_{\nu}(V)$
$D[P]\to A~{}P_{1},~{}A[P]\to S~{}P_{2},~{}S[S]\to P_{3}~{}P_{4}$ | $\tilde{t}_{\mu}(D)\tilde{t}^{\mu}(A)$
$D[P]\to A~{}P_{1},~{}A[P]\to T~{}P_{2},~{}T[D]\to P_{3}~{}P_{4}$ | $\tilde{t}_{\mu}(D)\tilde{t}_{\nu}(A)P^{\mu\nu\rho\sigma}(A)\tilde{t}_{\rho\sigma}(T)$
$D[P]\to V_{1}~{}P_{1},~{}V_{1}[P]\to V_{2}~{}P_{2},~{}V_{2}[P]\to P_{3}~{}P_{4}$ | $\tilde{t}_{\mu}(D)P^{\mu\nu}(V_{1})\epsilon_{\nu\lambda\sigma\alpha}p_{V_{1}}^{\alpha}\tilde{t}^{\lambda}(V_{1})\tilde{t}^{\sigma}(V_{2})$
$D[D]\to PT~{}P_{1},~{}PT[S]\to T~{}P_{2},~{}T[D]\to P_{3}~{}P_{4}$ | $\tilde{t}_{\mu\nu}(D)P^{\mu\nu\alpha\beta}(PT)\tilde{t}_{\alpha\beta}(T)$
$D[D]\to PT~{}P_{1},~{}PT[P]\to V~{}P_{2},~{}V[P]\to P_{3}~{}P_{4}$ | $\tilde{t}_{\mu\nu}(D)P^{\mu\nu\alpha\beta}(PT)\tilde{t}_{\alpha}(PT)\tilde{t}_{\beta}(V)$
$D[D]\to PT~{}P_{1},~{}PT[D]\to S~{}P_{2},~{}S[S]\to P_{3}~{}P_{4}$ | $\tilde{t}_{\mu\nu}(D)\tilde{t}^{\mu\nu}(PT)$
$D[D]\to T~{}P_{1},~{}T[D]\to V~{}P_{2},~{}V[P]\to P_{3}~{}P_{4}$ | $\tilde{t}_{\mu\nu}(D)P^{\mu\nu\alpha\beta}(T)\epsilon_{\alpha\lambda\sigma\rho}p_{T}^{\rho}\tilde{t}_{\beta}^{\lambda}(T)P^{\sigma\gamma}(T)\tilde{t}_{\gamma}(V)$
$D[D]\to T_{1}~{}P_{1},~{}T_{1}[P]\to T_{2}~{}P_{2},~{}T_{2}[D]\to P_{3}~{}P_{4}$ | $\tilde{t}_{\mu\nu}(D)P^{\mu\nu\alpha\beta}(T_{1})\epsilon_{\alpha\lambda\sigma\rho}p_{T_{1}}^{\rho}\tilde{t}^{\lambda}(T_{1})\tilde{t}^{\sigma}_{\beta}(T_{2})$
$D[S]\to S_{1}~{}S_{2},~{}S_{1}[S]\to P_{1}~{}P_{2},~{}S_{2}[S]\to P_{3}~{}P_{4}$ | 1
$D[P]\to V~{}S,~{}V[P]\to P_{1}~{}P_{2},~{}S[S]\to P_{3}~{}P_{4}$ | $\tilde{t}_{\mu}(D)\tilde{t}^{\mu}(V)$
$D[S]\to V_{1}~{}V_{2},~{}V_{1}[P]\to P_{1}~{}P_{2},~{}V_{2}[P]\to P_{3}~{}P_{4}$ | $\tilde{t}_{\mu}(V_{1})\tilde{t}^{\mu}(V_{2})$
$D[P]\to V_{1}~{}V_{2},~{}V_{1}[P]\to P_{1}~{}P_{2},~{}V_{2}[P]\to P_{3}~{}P_{4}$ | $\epsilon_{\mu\nu\alpha\beta}p_{D}^{\beta}\tilde{t}^{\mu}(D)\tilde{t}^{\nu}(V_{1})\tilde{t}^{\alpha}(V_{2})$
$D[D]\to V_{1}~{}V_{2},~{}V_{1}[P]\to P_{1}~{}P_{2},~{}V_{2}[P]\to P_{3}~{}P_{4}$ | $\tilde{t}_{\mu\nu}(D)\tilde{t}^{\mu}(V_{1})\tilde{t}^{\nu}(V_{2})$
$D[D]\to T~{}S,~{}T[D]\to P_{1}~{}P_{2},~{}S[S]\to P_{3}~{}P_{4}$ | $\tilde{t}_{\mu\nu}(D)\tilde{t}^{\mu\nu}(T)$
$D[P]\to T~{}V,~{}T[D]\to P_{1}~{}P_{2},~{}V[P]\to P_{3}~{}P_{4}$ | $\tilde{t}_{\mu}(D)\tilde{t}^{\mu\nu}(T)\tilde{t}_{\nu}(V)$
$D[D]\to T~{}V,~{}T[D]\to P_{1}~{}P_{2},~{}V[P]\to P_{3}~{}P_{4}$ | $\epsilon_{\mu\nu\alpha\beta}p_{D}^{\beta}\tilde{t}^{\mu\rho}(D)\tilde{t}^{\nu}_{\rho}(T)\tilde{t}^{\alpha}(V)$
$D[S]\to T_{1}~{}T_{2},~{}T_{1}[D]\to P_{1}~{}P_{2},~{}T_{2}[D]\to P_{3}~{}P_{4}$ | $\tilde{t}_{\mu\nu}(T_{1})\tilde{t}^{\mu\nu}(T_{2})$
$D[P]\to T_{1}~{}T_{2},~{}T_{1}[D]\to P_{1}~{}P_{2},~{}T_{2}[D]\to P_{3}~{}P_{4}$ | $\epsilon_{\mu\nu\alpha\beta}p_{D}^{\beta}\tilde{t}^{\mu}(D)\tilde{t}^{\nu\rho}(T_{1})\tilde{t}_{\rho}^{\alpha}(T_{2})$
$D[D]\to T_{1}~{}T_{2},~{}T_{1}[D]\to P_{1}~{}P_{2},~{}T_{2}[D]\to P_{3}~{}P_{4}$ | $\tilde{t}_{\mu\nu}(D)\tilde{t}^{\mu\rho}(T_{1})\tilde{t}^{\nu}_{\rho}(T_{2})$

The Blatt-Weisskopf barrier factors $B_{L}(q)$ are derived by assuming a
square well interaction potential as

$\displaystyle B_{L=0}(q)=$ $\displaystyle 1~{},$ (11) $\displaystyle
B_{L=1}(q)=$ $\displaystyle\sqrt{\frac{2}{q^{2}+q_{R}^{2}}}~{},$ (12)
$\displaystyle B_{L=2}(q)=$
$\displaystyle\sqrt{\frac{13}{q^{4}+3q^{2}q_{R}^{2}+9q_{R}^{4}}}~{},$ (13)

where $q$ is the momentum of daughter particle in the rest frame of the mother
particle, $L$ is the orbital angular momentum, and $q_{R}-1/R$ is a hadron
“scale” parameter ($R$ denotes the radius of the centrifugal barrier). In this
analysis, the radius $R$ is taken to be $5.0~{}{\rm{GeV}^{-1}}c$ for $D^{0}$
mesons, and $3.0~{}{\rm{GeV}^{-1}}c$ for other intermediate resonances.

Generally, the propagators of resonances are described by a relativistic
Breit-Wigner function

$\displaystyle P(s)=\frac{1}{m_{0}^{2}-s-im_{0}\Gamma(s)}~{},$ (14)

where $\Gamma(s)$ is the energy dependent width and $m_{0}$ is the nominal
mass of the resonance. For a resonance decaying into two scalar particles
$a\to bc$, the energy dependent partial width $\Gamma_{a\to bc}(s)$ is given
by

$\displaystyle\Gamma_{a\to bc}(s)=\Gamma_{0}^{a\to
bc}\left(\frac{q}{q_{0}}\right)^{2L+1}\left(\frac{m_{0}}{\sqrt{s}}\right)\left(\frac{B_{L}(q)}{B_{L}(q_{0})}\right)^{2}~{},$
(15)

where $\Gamma_{0}^{a\to bc}$ is the nominal width when $s=m_{0}^{2}$ and
$q_{0}$ is the corresponding momentum of the daughter particle in the rest
frame of the mother particle.

In this analysis, the propagators of $\rho(770)$ and $\rho(1450)$ are
described by the Gounaris-Sakurai parametrization [40]. The propagator of
$f_{0}(980)$ in the decay $a_{1}(1420)\to f_{0}(980)\pi$ is described by a
Flatt$\rm{\acute{e}}$ parametrization of $\pi\pi$ and $KK$ coupled channels
with parameters from Ref. [41].

The K-matrix parametrization [42, 43] instead of the Breit-Wigner formula is
adopted for $\pi\pi$ S-wave. The P-vector parametrization of the K-matrix for
the $\pi\pi$ S-wave with five coupled channels $\pi\pi$, $KK$, $\pi\pi\pi\pi$,
$\eta\eta$, $\eta\eta^{\prime}$ and five poles is written as

$\displaystyle F_{\mu}(s)=[I-iK(s)\rho(s)]^{-1}_{\mu\nu}P_{\nu}(s)~{},$ (16)

where $I$ is the identity matrix, $K$ is the K-matrix describing the
scattering process, $\rho$ is the PHSP matrix and $P$ is the initial
production vector (P-vector). The indices $\mu$ and $\nu$ denote the coupled
channels ($\pi\pi$, $KK$, $\pi\pi\pi\pi$, $\eta\eta$, and
$\eta\eta^{\prime}$), and only $F_{\mu}(s)$ of the $\pi\pi$ component is used
as the propagator of the $\pi\pi$ S-wave. The K-matrix is given by

$\displaystyle
K_{\mu\nu}(s)=\left(\sum_{\alpha}\frac{g^{\alpha}_{\mu}g^{\alpha}_{\nu}}{m_{\alpha}^{2}-s}+f_{\mu\nu}^{\rm
scatt}\frac{1-s_{0}^{\rm scatt}}{s-s_{0}^{\rm scatt}}\right)f_{A0}(s)~{},$
(17)

which contains five poles ($\alpha=1$ to $5$). The parameters of the K-matrix
are fixed to those in Ref. [43]. The P-vector is written as

$\displaystyle
P_{\nu}(s)=\left(\sum_{\alpha}\frac{\beta_{\alpha}g^{\alpha}_{\nu}}{m_{\alpha}^{2}-s}+f_{\nu}^{\rm
prod}\frac{1-s_{0}^{\rm prod}}{s-s_{0}^{\rm prod}}\right)~{},$ (18)

where the poles are the same as for the K-matrix. The parameters
$\beta_{\alpha}$, $f_{1\nu}^{\rm prod}$, and $s_{0}^{\rm prod}$ are free in
the fit.

In practice, there are two $\pi\pi$ S-waves in the decay $D^{0}\to
S_{1}S_{2}$, and the corresponding propagator is

$\displaystyle F^{\prime}_{\mu\nu}(s_{1},s_{2})$ $\displaystyle=$
$\displaystyle[I-iK(s_{1})\rho(s_{1})]^{-1}_{\mu\rho}[I-iK(s_{2})\rho(s_{2})]^{-1}_{\nu\sigma}$
(19) $\displaystyle\times P_{\rho\sigma}(s_{1},s_{2})~{},$

where $s_{1}$ and $s_{2}$ are the invariant masses squared of $\pi\pi$,
$P_{\rho\sigma}(s_{1},s_{2})$ is the expansion of the product of two P-vectors
$P_{\rho}(s_{1})$ and $P_{\sigma}(s_{2})$, and the corresponding coefficients
of each term are taken to be independent parameters in the fit. Benefiting
from the exchange symmetry of two $\pi\pi$ S-waves in the amplitude,
$P_{\rho\sigma}(s_{1},s_{2})$ can be written as

$\displaystyle P_{\rho\sigma}(s_{1},s_{2})$ $\displaystyle=$
$\displaystyle\sum_{\alpha,\beta}^{\alpha\leq\beta}a_{\alpha,\beta}\bigg{[}\frac{g^{\alpha}_{\rho}g^{\beta}_{\sigma}}{(m_{\alpha}^{2}-s_{1})(m_{\beta}^{2}-s_{2})}$
(20)
$\displaystyle+\frac{g^{\beta}_{\rho}g^{\alpha}_{\sigma}}{(m_{\beta}^{2}-s_{1})(m_{\alpha}^{2}-s_{2})}\bigg{]}$
$\displaystyle+\sum_{\alpha}b_{\alpha,\rho}\frac{g^{\alpha}_{\sigma}(1-s_{0}^{\rm
prod})}{(s_{1}-s_{0}^{\rm prod})(m_{\alpha}^{2}-s_{2})}$
$\displaystyle+\sum_{\alpha}b_{\alpha,\sigma}\frac{g^{\alpha}_{\rho}(1-s_{0}^{\rm
prod})}{(s_{2}-s_{0}^{\rm prod})(m_{\alpha}^{2}-s_{1})}$
$\displaystyle+c_{[\rho,\sigma]}\frac{(1-s_{0}^{\rm
prod})^{2}}{(s_{1}-s_{0}^{\rm prod})(s_{2}-s_{0}^{\rm prod})},$

where $c_{[\rho,\sigma]}=c_{\rho,\sigma}+c_{\sigma,\rho}$, and
$c_{[\rho,\sigma]}$ with $\rho>\sigma$ is fixed to zero in the fit. Similar to
one $\pi\pi$ S-wave, only ($\pi\pi,\pi\pi$) component of $F^{\prime}_{\mu\nu}$
is used.

The propagators of the resonances $a_{1}(1260)$, $a_{1}(1640)$, $h_{1}(1170)$
and $\pi(1300)$ decaying into $3\pi$ are described by relativistic Breit-
Wigner functions with energy dependent widths with the coupled channel $3\pi$,
where the widths are obtained by integrating the amplitude squared over PHSP

$\displaystyle\Gamma_{a\to bcd}(s)\propto\frac{1}{m_{0}}\int\sum_{\rm
spin}|A_{a\to bcd}|^{2}\mathrm{d}\Phi_{3}~{}~{},$ (21)

based on the amplitudes $A_{R\to 3\pi}$ obtained in this analysis. For other
resonances decaying into $3\pi$, a relativistic Breit-Wigner function with a
constant width is used. The parameters of $a_{1}(1260)$ and $\pi(1300)$ are
determined in the fit and others are fixed to their respective PDG values.
[23] Due to its small contribution, the $a_{1}(1420)$ is described as a
relativistic Breit-Wigner function with a constant width, even though it has
been regarded as the effect of the Triangle Singularity [44]. Only the decay
of $a_{1}(1420)\to f_{0}(980)\pi$ is taken into account, and the corresponding
resonant parameters are taken from Ref. [45].

To take into account the resolution effect for the narrow resonances
$\phi(1020)$ and $\omega(782)$, the shapes of corresponding relativistic
Breit-Wigner functions in the $M(\pi^{+}\pi^{-}\pi^{0})$ distribution are
convolved with a Gaussian function. The parameters of this Gaussian function
are obtained by the fit to the $\eta$ peak in the $M(\pi^{+}\pi^{-}\pi^{0})$
distribution with the Breit-Wigner function of $\eta$ convolved with the same
Gaussian function.

In this analysis, two decay topologies, $D\to R_{1}R_{2}$ and $D\to
R_{1}\pi\to R_{2}\pi\pi$, are considered. Meanwhile, isospin symmetry is
considered by applying the Clebsch-Gordan (CG) coefficients between the
isospin processes. The isospin states for $u$ and $d$ quarks are

$\displaystyle|u\rangle:|\frac{1}{2},\frac{1}{2}\rangle,~{}|\bar{u}\rangle:-|\frac{1}{2},-\frac{1}{2}\rangle,~{}|d\rangle:|\frac{1}{2},-\frac{1}{2}\rangle,~{}|\bar{d}\rangle:|\frac{1}{2},\frac{1}{2}\rangle.$
(22)

Based on the quark components of hadrons, the isospin states of the hadrons
are shown in Table 6. The isospin amplitudes for typical processes are written
as

$\displaystyle A(D^{0}\to a_{1}^{+}\pi^{-})$ $\displaystyle=$ $\displaystyle
A(\rho^{0}_{\pi^{+}\pi^{-}}\pi^{+},\pi^{-})+A(\rho^{+}_{\pi^{+}\pi^{0}}\pi^{0},\pi^{-})-c[2A(f_{\pi^{+}\pi^{-}}\pi^{+},\pi^{-})+A(f_{\pi^{0}\pi^{0}}\pi^{+},\pi^{-})],$
(23) $\displaystyle A(D^{0}\to a_{1}^{0}\pi^{0})~{}$ $\displaystyle=$
$\displaystyle-A(\rho^{+}_{\pi^{+}\pi^{0}}\pi^{-},\pi^{0})+A(\rho^{-}_{\pi^{0}\pi^{-}}\pi^{+},\pi^{0})-c[2A(f_{\pi^{+}\pi^{-}}\pi^{0},\pi^{0})+A(f_{\pi^{0}\pi^{0}}\pi^{0},\pi^{0})],$
(24) $\displaystyle A(D^{0}\to a_{1}^{-}\pi^{+})$ $\displaystyle=$
$\displaystyle-A(\rho^{0}_{\pi^{+}\pi^{-}}\pi^{-},\pi^{+})-A(\rho^{-}_{\pi^{0}\pi^{-}}\pi^{0},\pi^{+})-c[2A(f_{\pi^{+}\pi^{-}}\pi^{-},\pi^{+})+A(f_{\pi^{0}\pi^{0}}\pi^{-},\pi^{+})],$
(25) $\displaystyle A(D^{0}\to h_{1}^{0}\pi^{0})~{}$ $\displaystyle=$
$\displaystyle-A(\rho^{+}_{\pi^{+}\pi^{0}}\pi^{-},\pi^{0})-A(\rho^{-}_{\pi^{0}\pi^{-}}\pi^{+},\pi^{0})+A(\rho^{0}_{\pi^{+}\pi^{-}}\pi^{0},\pi^{0}),$
(26) $\displaystyle A(D^{0}\to ff^{\prime})~{}~{}~{}$ $\displaystyle=$
$\displaystyle
4A(f_{\pi^{+}\pi^{-}},f^{\prime}_{\pi^{+}\pi^{-}})+2[A(f_{\pi^{+}\pi^{-}},f^{\prime}_{\pi^{0}\pi^{0}})+A(f_{\pi^{0}\pi^{0}},f^{\prime}_{\pi^{+}\pi^{-}})]+A(f_{\pi^{0}\pi^{0}},f^{\prime}_{\pi^{0}\pi^{0}}),$
(27) $\displaystyle A(D^{0}\to\rho^{0}f)~{}~{}~{}$ $\displaystyle=$
$\displaystyle
2A(\rho^{0}_{\pi^{+}\pi^{-}},f_{\pi^{+}\pi^{-}})+A(\rho^{0}_{\pi^{+}\pi^{-}},f_{\pi^{0}\pi^{0}}).$
(28)

Here, $a_{1}$ and $\rho$ are isospin vectors, whereas $h_{1}$ and
$f^{(^{\prime})}$ are isospin scalars. The normalization factors are dropped
here and $c$ is the relative difference (including magnitude and phase)
between two isospin processes which is determined in the amplitude analysis.

Table 6: The quark components and isospin states for different hadrons used in the analysis. | $\pi^{+},b^{+},\rho^{+},a^{+}$ | $\pi^{0},b^{0},\rho^{0},a^{0}$ | $\pi^{-},b^{-},\rho^{-},a^{-}$
---|---|---|---
$I=1$ | $|u\bar{d}\rangle:|1,1\rangle$ | $|\frac{u\bar{u}-d\bar{d}}{\sqrt{2}}\rangle:-|1,0\rangle$ | $|d\bar{u}\rangle:-|1,-1\rangle$
| $h,\omega,\phi,f$
$I=0$ | $|(u\bar{u}+d\bar{d})+g(s\bar{s})\rangle:-|0,0\rangle$ |

##  6 AMPLITUDE FIT PROCESS AND RESULT

In the amplitude fit, the sPlot technique [46] is applied to deal with the
effect of background. In the sPlot method, the four-momenta of final states
are chosen as control variables $x$, $M_{\rm bc}^{\rm sig}$ and $M_{\rm
bc}^{\rm tag}$ are chosen as discriminating variables $y$, and the weight for
each candidate event is obtained by an unbinned maximum likelihood fit on the
2D distribution of $M_{\rm bc}^{\rm sig}$ versus $M_{\rm bc}^{\rm tag}$, as
described in Sec. 4.

According to the sPlot technique and based on the above DT fit, the weight
$W_{s}(y)$ of each candidate event is calculated as

$\displaystyle
W_{s}(y)=\frac{\sum_{i}V_{si}F_{i}(y)}{\sum_{i}N_{j}F_{j}(y)+\sum_{k}N^{\prime}_{k}F^{\prime}_{k}(y)},$
(29)

where the index $s$ denotes the signal, $F_{i(j)}(y)$ are the PDFs of the
signal ($i(j)=s$) and backgrounds ($i(j)\neq s$) from BKGI to BGKV with
floating yields, $F^{\prime}_{k}(y)$ are the PDFs of BKGVI with fixed yields,
and $N_{i(j)}$ and $N_{k}^{{}^{\prime}}$ are the corresponding yields for each
component. The inverse of covariance matrix $V_{ij}$ is calculated with

$\displaystyle V_{ij}^{-1}=\sum_{n\in{\rm
Data}}\frac{F_{i}(y_{n})F_{j}(y_{n})}{(\sum_{k}N_{k}F_{k}(y_{n})+\sum_{l}N^{\prime}_{l}F^{\prime}_{l}(y_{n}))^{2}},$
(30)

where the summations for indices $k$ and $l$ run over all the components of
$F_{k}(y)$ and $F^{\prime}_{l}(y)$ as described above, respectively. The
summation for index $n$ runs over all the events in the DT fit.

Based on the sPlot technique, $W_{s}(y)$ defined in Eq. (29) includes the
effects of the backgrounds with the fixed number of events. To take into
account these effects, the coefficients $c_{i}$ are calculated as

$\displaystyle c_{i}=\sum_{j}V_{sj}\nu_{ij},$ (31)

where the summation runs over the signal and backgrounds from BKGI to BKGV
with floating yields in the fit, $i$ represents the different components in
the BKGVI, and $\nu_{ij}$ is

$\displaystyle\nu_{ij}=\sum_{n\in{\rm
Data}}\frac{F^{\prime}_{i}(y_{n})F_{j}(y_{n})}{(\sum_{k}N_{k}F_{k}(y_{n})+\sum_{l}N^{\prime}_{l}F^{\prime}_{l}(y_{n}))^{2}}.$
(32)

With the obtained $W_{s}(y)$ and $c_{i}$, the $x$ distribution of all
candidate events with the weight $W_{s}(y)$ is the real signal together with
the contributions of the peaking background in BKGVI,

$\displaystyle N_{s}P_{s}(x)+\sum_{j}c_{j}N^{\prime}_{j}P^{\prime}_{j}(x),$
(33)

where $P_{s}(x)$ and $P^{\prime}_{j}(x)$ are the $x$ distributions of the
signal and the peaking background $j$ in the BKGVI and $N_{s}$ and
$N^{\prime}_{j}$ are the corresponding yields obtained from the above 2D fit
and MC simulation, respectively.

In practice, the PDF of observing a signal event with the given final
kinematic $p$ is written as

$\displaystyle
P_{s}(p)=\frac{\epsilon_{s}(p)|M_{f}(p)|^{2}\phi_{4}(p)}{\int\epsilon_{s}(p)|M_{f}(p)|^{2}\mathrm{d}\Phi_{4}}~{},$
(34)

where $\epsilon_{s}(p)$ is the signal efficiency, $|M_{f}(p)|^{2}$ is the
differential cross section as discussed in Sec. 5, and $\phi_{4}(p)$ is the
PHSP density. The normalization factor is calculated by MC integration with
the PHSP signal MC sample after event selection,

$\displaystyle\int\epsilon_{s}(p)|M_{f}(p)|^{2}\mathrm{d}\Phi_{4}\propto\frac{1}{N_{\rm
MC}}\sum_{i=1}^{N_{\rm MC}}|M_{f}(p_{i})|^{2}.$ (35)

where the $N_{\rm MC}$ is the number of events of the PHSP signal MC.
According to Eq. (33), the weighted likelihood $\mathrm{ln}L$ is given by

$\displaystyle f\left[\sum_{i\in{\rm
Data}}W_{s}(y_{i})\mathrm{ln}P_{s}(p_{i})-\sum_{j}\sum_{i\in{\rm
BMC}_{j}}\omega_{j}\mathrm{ln}P_{s}(p_{i})\right],$ (36)

where the second term in the bracket is the contribution from the peaking
background, which is estimated with the corresponding simulated background MC.
The normalization factor $\omega_{j}$ is given by

$\displaystyle\omega_{j}=c_{j}N_{{\rm BKG}_{j}}^{\rm Data}/N_{{\rm
BKG}_{j}}^{\rm MC}.$ (37)

Here, $c_{j}$ is obtained from Eq. (31), and $N_{{\rm BKG}_{j}}^{\rm Data}$
and $N_{{\rm BKG}_{j}}^{\rm MC}$ are the background yields in data and
simulated background events, respectively. The factor $f$ in Eq. (36),

$\displaystyle f=\frac{\sum_{i\in{\rm
Data}}W_{s}(y_{i})+\sum_{j}\omega_{j}N^{\rm MC}_{{\rm
BKG}_{j}}}{\sum_{i\in{\rm Data}}W^{2}_{s}(y_{i})+\sum_{j}\omega_{j}^{2}N^{\rm
MC}_{{\rm BKG}_{j}}}~{},$ (38)

is the global factor to correct the statistical bias in the weighted maximum
likelihood fit.

The total likelihood function in this analysis is summed over the two signal
channels and three tag modes,

$\displaystyle\mathrm{ln}L_{\rm total}=\sum_{i\in{\rm
tag}}\left(\mathrm{ln}L_{i}^{\pi^{+}\pi^{-}\pi^{+}\pi^{-}}+\mathrm{ln}L_{i}^{\pi^{+}\pi^{-}\pi^{0}\pi^{0}}\right),$
(39)

and the free parameters are optimized via a maximum likelihood fit using the
MINUIT [47] package.

Generally, all the possible intermediate processes including the resonances
listed in Table 7 and based on $J^{PC}$ conservation are considered in the
fit. The only exception are $D^{0}$ decays, where parity conservation is not
required. Only processes with a significance greater than $5\sigma$ are kept
during the fit unless otherwise noted. Here, the significance of a specific
amplitude is calculated according to the Wilks’s Theorem by comparing the
change of log-likelihood (2$\Delta$ln$L$) to the expected values from the chi-
square distribution ($\chi^{2}_{\Delta N_{\rm para}}$) with the number of
degrees of freedom (NDF) equal to the change of the numbers of fit parameters
($\Delta N_{\rm para}$). If the significance of the amplitudes containing the
isospin vector in the $3\pi$ invariant mass spectrum is greater than
5$\sigma$, the corresponding isospin partners with significance greater than
3$\sigma$ are also kept. For the P-vector of the $\pi\pi$ S-wave, only the
parameters before the 1st, 2nd poles and $\pi\pi$, $KK$ non-resonant terms are
considered, and the others are fixed to zero in the nominal fit. Meanwhile,
only those with significance greater than $3\sigma$ are kept. We assume the
same values of the parameter $s_{0}^{\rm prod}$ in the P-vector of all
$\pi\pi$ S-waves and fix them to $-5~{}{\rm GeV}^{2}/c^{4}$ since they are
insensitive to any choice if $s_{0}^{\rm prod}\leq-5~{}{\rm GeV}^{2}/c^{4}$.

Table 7: Resonances considered in this analysis. | $J^{P}=0^{+}$ | $J^{P}=0^{-}$ | $J^{P}=1^{+}$ | $J^{P}=1^{-}$ | $J^{P}=2^{+}$ | $J^{P}=2^{-}$
---|---|---|---|---|---|---
$\pi\pi$ | $(\pi\pi)_{S}$ |  |  | $\rho(770)$ | $f_{2}(1270)$ |
|  |  | $\rho(1450)$ |  |
$\pi\pi\pi$ |  | $\pi(1300)$ | $a_{1}(1260)$ | $\omega(782)$ | $a_{2}(1320)$ | $\pi_{2}(1670)$
|  | $a_{1}(1420)$ | $\phi(1020)$ |  |
|  | $a_{1}(1640)$ | $\pi_{1}(1400)$ |  |
|  | $h_{1}(1170)$ | $\pi_{1}(1600)$ |  |

To find the optimal solution, the baseline model which contains the processes
from Refs. [10, 11] is built up first. Next, starting from the baseline model,
the significance of each possible process is tested and the most significant
one among those satisfying the significance requirement is added to the
current model. This step is repeated until no additional processes can be
added. After this, the significances of individual processes in the existing
model are tested again and those that do not satisfy the significance
requirement are removed. The above steps are repeated until all the processes
in the model satisfy the significance requirement and no further ones can be
added.

With the nominal amplitude model, the fit fraction (FF) of a specific
amplitude $i$ is calculated as

$\displaystyle{\rm
FF}_{i}=\frac{\int|\Lambda_{i}U_{i}(p)|^{2}\mathrm{d}\Phi}{\int|\sum_{j}\Lambda_{j}U_{j}(p)|^{2}\mathrm{d}\Phi}.$
(40)

The $CP$-even fraction $F_{+}$ is calculated as

$\displaystyle
F_{+}^{f}=\frac{\int|A^{+}_{f}(p)|^{2}\mathrm{d}\Phi}{\int|A^{+}_{f}(p)|^{2}+|A^{-}_{f}(p)|^{2}\mathrm{d}\Phi},$
(41)

where $A^{\pm}_{f}(p)=\frac{1}{\sqrt{2}}[A_{f}(p)\pm\bar{A}_{f}(p)]$ is the
amplitude of the $D^{0}\to f$ decay in a $CP$-even or $CP$-odd state. The FF
results of different amplitudes are summarized in Tables 8, 9 and 10. The
resonant parameters, masses, and widths of $a_{1}(1260)$ and $\pi(1300)$ are
determined by the parameter scans as shown in Fig. 4. The fit results show
large interferences among the dominant intermediate processes $D^{0}\to
a_{1}(1260)\pi$, $D^{0}\to\pi(1300)\pi$, $D^{0}\to\rho(770)\rho(770)$, and
$D^{0}\to 2(\pi\pi)_{S}$. The results for $F_{+}$ obtained in this analysis
are presented in Table 11, which show good agreement with other measurements.
In all above tables, the mean values are obtained based on the output from the
MINUIT fit, while the corresponding statistical uncertainties are estimated by
the bootstrap method [48] based on data, since the weighted maximum likelihood
fit is used, and the statistical uncertainties given by the inverse second
derivative of the negative logarithmic likelihood are no longer asymptotically
correct [49]. In the bootstrap method, the bootstrap samples are generated by
repeatedly resampling the data set with replacement for thousand times, and
both sPlot and amplitude analysis are performed for these samples. The width
of the distribution of estimated parameter values is used as estimator for the
parameter uncertainty. The comparisons between data and MC projection based on
the nominal model for various invariant mass and angle distributions are shown
in Figs. 5 and 6. To determine the quality of the fit, the goodness of fit is
calculated for the amplitude analysis. Here the adaptive bins (with more than
25 events and the sum of weights greater than 0 in each bin) in five-
dimensional PHSP {$M(\pi^{+}\pi^{+})$, $M(\pi^{-}\pi^{-})$, ${\rm
cos}\theta_{+,+}$, ${\rm cos}\theta_{-,-}$, $\phi_{++,--}$} and
{$M(\pi^{+}\pi^{-})$, $M(\pi^{0}\pi^{0})$, ${\rm cos}\theta_{+,-}$, ${\rm
cos}\theta_{0,0}$, $\phi_{+-,00}$} are defined for
$D^{0}\to\pi^{+}\pi^{-}\pi^{+}\pi^{-}$ and
$D^{0}\to\pi^{+}\pi^{-}\pi^{0}\pi^{0}$, respectively, where $\theta$ is the
helicity angle of $\pi$ in the corresponding $\pi\pi$ system, $\phi$ is the
angle between two $\pi\pi$ decay planes in the $D^{0}$ rest frame and the
indices $+$, $-$ and $0$ in subscript represent $\pi^{+}$, $\pi^{-}$ and
$\pi^{0}$, respectively. The $\chi^{2}$ of the amplitude analysis fit is
calculated as

$\displaystyle\chi^{2}=\sum_{i}\frac{(w^{i}_{\rm Data}-w^{i}_{\rm
MC})^{2}}{(w_{\rm Data}^{i})^{2}+(w_{\rm MC}^{i})^{2}},$ (42)

where the summation runs over all the bins in the ten variables, $w^{i}$ and
$(w^{i})^{2}$ are the sum and quadrature sum of weights in the $i$-th bin,
respectively. The goodness of fit is calculated to be $\chi^{2}/{\rm
NDF}=382/396=0.96$, where NDF is equal to $(N_{\rm bin}-N_{\rm para}-6)$.

Table 8: The fit parameters, FFs and significances of individual amplitudes, where the first uncertainties are statistical and the second in FFs are systematic. Amplitude | Magnitude | Phase (rad) | FF (%) | Significance ($\sigma$)
---|---|---|---|---
$\pi^{+}\pi^{-}\pi^{+}\pi^{-}$ | $\pi^{+}\pi^{-}\pi^{0}\pi^{0}$
$a_{1}(1260)^{+}\pi^{-}$ | 100(fixed) | 0(fixed) | $82.2\pm 3.3\pm 16.0$ | $57.4\pm 2.7\pm 7.3$ |  $>$20
$a_{1}(1260)^{-}\pi^{+}$ | $35.3\pm 2.7$ | $\phantom{+}0.23\pm 0.07$ | $10.3\pm 1.5\pm 2.5$ | $7.2\pm 1.1\pm 2.2$ | 15.6
$a_{1}(1260)^{0}\pi^{0}$ | $50.9\pm 3.1$ | $-2.99\pm 0.06$ | - | $32.9\pm 3.2\pm 8.6$ |  $>$20
$a_{1}(1420)^{+}\pi^{-}$ | $19.0\pm 3.6$ | $\phantom{+}2.70\pm 0.18$ | $0.6\pm 0.2\pm 0.2$ | $0.3\pm 0.1\pm 0.1$ | 6.0
$a_{1}(1640)^{+}\pi^{-}$ | $20.1\pm 3.0$ | $-2.07\pm 0.16$ | $1.7\pm 0.5\pm 0.9$ | $1.1\pm 0.3\pm 0.7$ | 7.3
$a_{1}(1640)^{-}\pi^{+}$ | $10.5\pm 2.8$ | $-1.26\pm 0.29$ | $0.5\pm 0.3\pm 0.3$ | $0.3\pm 0.2\pm 0.2$ | 5.2
$a_{2}(1320)^{+}\pi^{-}$ | $0.23\pm 0.07$ | $-2.92\pm 0.30$ | $0.2\pm 0.1\pm 0.1$ | $0.2\pm 0.1\pm 0.1$ | 4.6
$a_{2}(1320)^{-}\pi^{+}$ | $0.30\pm 0.05$ | $-0.47\pm 0.21$ | $0.4\pm 0.1\pm 0.1$ | $0.3\pm 0.1\pm 0.1$ | 6.4
$h_{1}(1170)^{0}\pi^{0}$ | $9.7\pm 2.2$ | $-0.59\pm 0.27$ | - | $1.3\pm 0.6\pm 1.1$ | 6.5
$\pi(1300)^{+}\pi^{-}$ | $76.3\pm 3.6$ | $-2.325\pm 0.044$ | $32.3\pm 2.6\pm 4.5$ | $15.6\pm 1.4\pm 2.6$ |  $>$20
$\pi(1300)^{-}\pi^{+}$ | $65.1\pm 3.4$ | $-2.631\pm 0.045$ | $23.5\pm 2.3\pm 4.1$ | $11.4\pm 1.1\pm 2.2$ | 19.0
$\pi(1300)^{0}\pi^{0}$ | $61.1\pm 3.2$ | $\phantom{+}0.61\pm 0.05$ | - | $23.2\pm 2.8\pm 3.5$ | 16.8
$\pi_{2}(1670)^{0}\pi^{0}$ | $12.2\pm 1.5$ | $-1.11\pm 0.14$ | - | $1.1\pm 0.2\pm 0.4$ | 6.9
$\rho(770)^{0}\rho(770)^{0}$ | - | - | $28.0\pm 1.9\pm 3.2$ | - |  $>$20
$\rho(770)^{0}\rho(770)^{0}[S]$ | $6.1\pm 1.1$ | $-3.10\pm 0.17$ | $1.7\pm 0.6\pm 0.5$ | - | 6.5
$\rho(770)^{0}\rho(770)^{0}[P]$ | $6.17\pm 0.36$ | $\phantom{+}1.62\pm 0.07$ | $9.8\pm 1.0\pm 0.9$ | - | 19.5
$\rho(770)^{0}\rho(770)^{0}[D]$ | $4.54\pm 0.22$ | $-3.06\pm 0.05$ | $23.1\pm 2.1\pm 2.4$ | - |  $>$20
$\rho(770)^{0}\rho(1450)^{0}$ | - | - | $2.5\pm 0.9\pm 1.1$ | - | 8.0
$\rho(770)^{0}\rho(1450)^{0}[P]$ | $13.9\pm 2.5$ | $\phantom{+}0.68\pm 0.20$ | $1.0\pm 0.4\pm 0.5$ | - | 6.4
$\rho(770)^{0}\rho(1450)^{0}[D]$ | $5.6\pm 1.3$ | $\phantom{+}3.08\pm 0.20$ | $1.5\pm 0.9\pm 1.1$ | - | 5.0
$\rho(770)^{+}\rho(770)^{-}$ | - | - | - | $90.9\pm 3.9\pm 7.9$ |  $>$20
$\rho(770)^{+}\rho(770)^{-}[S]$ | $13.7\pm 1.2$ | $\phantom{+}3.03\pm 0.09$ | - | $13.0\pm 2.0\pm 3.5$ | 13.0
$\rho(770)^{+}\rho(770)^{-}[P]$ | $7.10\pm 0.36$ | $-1.69\pm 0.07$ | - | $19.6\pm 1.3\pm 2.0$ |  $>$20
$\rho(770)^{+}\rho(770)^{-}[D]$ | $4.59\pm 0.22$ | $\phantom{+}0.06\pm 0.05$ | - | $36.0\pm 3.0\pm 2.5$ |  $>$20
$\rho(770)^{+}\rho(1450)^{-}[D]$ | $8.1\pm 1.7$ | $-1.01\pm 0.18$ | - | $1.69\pm 0.8\pm 1.7$ | 6.3
$\rho(770)^{0}(\pi\pi)_{S}$ | - | - | $2.7\pm 0.6\pm 1.8$ | $1.0\pm 0.2\pm 0.4$ | 12.8
$\beta_{1}$ | $8.4\pm 3.6$ | $-1.68\pm 0.50$ | - | - | -
$f_{\pi\pi}^{\rm prod}$ | $40.7\pm 5.0$ | $-0.50\pm 0.14$ | - | - | -
$f_{KK}^{\rm prod}$ | $121\pm 25$ | $\phantom{+}1.73\pm 0.23$ | - | - | -
$(\pi^{+}\pi^{-})_{S}(\pi\pi)_{S}$ | - | - | $62.8\pm 4.6\pm 9.6$ | $37.4\pm 3.0\pm 5.1$ |  $>$20
$a_{1,1}$ | $2224\pm 35$ | $-1.044\pm 0.019$ | - | - | -
$a_{1,2}$ | $7287\pm 62$ | $\phantom{+}1.727\pm 0.009$ | - | - | -
$b_{2,\pi\pi}$ | $8816\pm 120$ | $-1.107\pm 0.014$ | - | - | -
$c_{[\pi\pi,\pi\pi]}$ | $2433\pm 96$ | $\phantom{+}1.796\pm 0.043$ | - | - | -
$c_{[\pi\pi,KK]}$ | $5417\pm 477$ | $\phantom{+}2.68\pm 0.10$ | - | - | -
$f_{2}(1270)^{0}(\pi\pi)_{S}$ | - | - | $1.8\pm 0.4\pm 1.2$ | $1.1\pm 0.2\pm 0.8$ | 9.1
$f_{\pi\pi}^{\rm prod}$ | $18.3\pm 1.8$ | $-1.39\pm 0.10$ | - | - | -
$f_{KK}^{\rm prod}$ | $56\pm 10$ | $\phantom{+}2.29\pm 0.20$ | - | - | -
$\omega(782)\pi^{0}$ | $1.58\pm 0.30$ | $-0.50\pm 0.44$ | - | $0.9\pm 0.4\pm 0.2$ | 6.1
$\phi(1020)\pi^{0}$ | $0.44\pm 0.06$ | $\phantom{+}2.51\pm 0.41$ | - | $1.5\pm 0.4\pm 0.3$ | 7.4
Table 9: The fit parameters, FFs and significance for the three-body decays of
$a_{1}(1260)$, $a_{1}(1420)$, $a_{1}(1640)$, $a_{2}(1320)$, $h_{1}(1170)$,
$\pi(1300)$ and $\pi_{2}(1670)$. The first uncertainties are statistical and
the second in FFs are systematic.

Amplitude Magnitude Phase (rad) Relative FF (%) Significance ($\sigma$)
$\pi^{+}\pi^{-}\pi^{+}\pi^{-}$ $\pi^{+}\pi^{-}\pi^{0}\pi^{0}$ charge $=\pm 1$
charge $=\pm 1$ charge $=0$ $a_{1}(1260)\to\rho(770)\pi[S]$ 1(fixed) 0(fixed)
$79.7\pm 2.2\pm 3.5$ $81.2\pm 2.0\pm 3.2$ $78.9\pm 2.2\pm 3.8$ $>$20
$a_{1}(1260)\to\rho(770)\pi[D]$ $0.060\pm 0.009$ $-0.01\pm 0.17$ $1.3\pm
0.3\pm 0.8$ $1.2\pm 0.3\pm 0.7$ $1.3\pm 0.3\pm 0.8$ 7.5 $a_{1}(1260)\to
f_{2}(1270)\pi[P]$ $0.311\pm 0.033$ $-1.58\pm 0.11$ $1.8\pm 0.4\pm 0.6$
$0.8\pm 0.2\pm 0.2$ $1.2\pm 0.2\pm 0.4$ 10.6
$a_{1}(1260)\to(\pi^{+}\pi^{-})_{S}\pi[P]$ - - $5.4\pm 0.6\pm 1.0$ $3.1\pm
0.4\pm 0.8$ $3.9\pm 0.5\pm 0.9$ 15.8 $\beta_{1}$ $0.83\pm 0.13$ $-2.18\pm
0.14$ - - - - $f_{\pi\pi}^{\rm prod}$ $2.47\pm 0.16$ $\phantom{+}0.34\pm 0.08$
- - - - $f_{KK}^{\rm prod}$ $6.59\pm 0.84$ $\phantom{+}1.76\pm 0.12$ - - - -
$\pi(1300)\to\rho(770)\pi$ 1(fixed) 0(fixed) $53.8\pm 2.9\pm 8.8$ $79.0\pm
2.0\pm 6.8$ $73.6\pm 2.6\pm 8.9$ $>$20 $\pi(1300)\to(\pi^{+}\pi^{-})_{S}\pi$ -
- $51.1\pm 2.9\pm 8.7$ $26.2\pm 2.1\pm 6.5$ $36.7\pm 2.6\pm 8.6$ 14.4
$\beta_{1}$ $5.0\pm 0.5$ $-0.64\pm 0.10$ - - - - $f_{KK}^{\rm prod}$ $43.8\pm
2.8$ $\phantom{+}0.34\pm 0.06$ - - - - $a_{1}(1420)\to f_{0}(980)\pi[P]$
1(fixed) 0(fixed) 100 100 - 6.0 $a_{1}(1640)\to\rho(770)\pi[S]$ 1(fixed)
0(fixed) 100 100 - 9.1 $a_{2}(1320)\to\rho(770)\pi[D]$ 1(fixed) 0(fixed) 100
100 - 7.3 $h_{1}(1170)\to\rho(770)\pi[S]$ 1(fixed) 0(fixed) - - 100 6.5
$\pi_{2}(1670)\to f_{2}(1270)\pi[S]$ 1(fixed) 0(fixed) - - 100 6.9

Table 10: Masses and widths of $a_{1}(1260)$ and $\pi(1300)$ obtained via parameter scans and presented in the PDG [23]. | This work | PDG
---|---|---
| ${\rm Mass}~{}({\rm GeV}/c^{2})$ | ${\rm Width}~{}({\rm GeV})$ | ${\rm Mass}~{}({\rm GeV}/c^{2})$ | ${\rm Width}~{}({\rm GeV})$
$a_{1}(1260)$ | $1.193\pm 0.005\pm 0.023$ | $0.487\pm 0.009\pm 0.041$ | $1.230\pm 0.040$ | $0.250-0.600$
$\pi(1300)$ | $1.534\pm 0.011\pm 0.022$ | $0.610\pm 0.030\pm 0.093$ | $1.300\pm 0.100$ | $0.200-0.600$

\begin{overpic}[width=106.23477pt]{fig/resscan/Ma11260.pdf}
\end{overpic}\begin{overpic}[width=106.23477pt]{fig/resscan/Wa11260.pdf}
\end{overpic}\begin{overpic}[width=106.23477pt]{fig/resscan/Mpi1300.pdf}
\end{overpic}\begin{overpic}[width=106.23477pt]{fig/resscan/Wpi1300.pdf}
\end{overpic}

Figure 4: Likelihood scans for masses and widths of $a_{1}(1260)$ and $\pi(1300)$. Table 11: The $CP$-even fractions obtained in this work and comparisons with the CLEO-c and prior BESIII measurements. | $F_{+}^{\pi^{+}\pi^{-}\pi^{+}\pi^{-}}$ |  $F_{+}^{\pi^{+}\pi^{-}\pi^{0}\pi^{0}}$(non-$\eta$)
---|---|---
This work (model-dependent) | $(75.2\pm 1.1_{\rm stat.}\pm 1.5_{\rm syst.})\%$ | $(68.9\pm 1.5_{\rm stat.}\pm 2.4_{\rm syst.})\%$
CLEO-c (model-dependent) |  $(72.9\pm 0.9_{\rm stat.}\pm 1.5_{\rm syst.}\pm 1.0_{\rm model})\%$ [11] | -
CLEO-c (model-independent, global) |  $(73.7\pm 2.8)\%$ [50] | -
CLEO-c (model-independent, binned) |  $(76.9\pm 2.1_{\rm stat.}\pm 1.0_{\rm syst.}\pm 0.2_{K_{S}~{}{\rm veto}})\%$ [6] | -
BESIII (model-independent, global) |  $(73.4\pm 1.5_{\rm stat.}\pm 0.8_{\rm syst.})\%$ [51] |  $(68.2\pm 7.7)\%$ [52]

\begin{overpic}[width=138.76157pt]{fig/PipPimPipPim/PWA/MPipPim.pdf}
\end{overpic}\begin{overpic}[width=138.76157pt]{fig/PipPimPipPim/PWA/MPipPipPim.pdf}
\end{overpic}\begin{overpic}[width=138.76157pt]{fig/PipPimPipPim/PWA/MPimPimPip.pdf}
\end{overpic}\begin{overpic}[width=138.76157pt]{fig/PipPimPipPim/PWA/costheta_pm.pdf}
\end{overpic}\begin{overpic}[width=138.76157pt]{fig/PipPimPipPim/PWA/chi_pm_pm.pdf}
\end{overpic}\begin{overpic}[width=138.76157pt]{fig/PipPimPipPim/PWA/leg.pdf}
\end{overpic}

Figure 5: The mass and angular distributions for
$D^{0}\to\pi^{+}\pi^{-}\pi^{+}\pi^{-}$, where $\theta_{a}^{ab}$ is the
helicity angle of $a$ in the $ab$ system, $\phi_{ab,cd}$ is the angle between
the decay planes of $ab$ and $cd$ systems in the $D^{0}$ rest frame.

\begin{overpic}[width=138.76157pt]{fig/PipPim2Pi0/PWA/MPipPimPi0.pdf}
\end{overpic}\begin{overpic}[width=138.76157pt]{fig/PipPim2Pi0/PWA/MPipPi01Pi02.pdf}
\end{overpic}\begin{overpic}[width=138.76157pt]{fig/PipPim2Pi0/PWA/MPimPi01Pi02.pdf}
\end{overpic}\begin{overpic}[width=138.76157pt]{fig/PipPim2Pi0/PWA/MPipPim.pdf}
\end{overpic}\begin{overpic}[width=138.76157pt]{fig/PipPim2Pi0/PWA/MPi01Pi02.pdf}
\end{overpic}\begin{overpic}[width=138.76157pt]{fig/PipPim2Pi0/PWA/MPipPi0.pdf}
\end{overpic}\begin{overpic}[width=138.76157pt]{fig/PipPim2Pi0/PWA/MPimPi0.pdf}
\end{overpic}\begin{overpic}[width=138.76157pt]{fig/PipPim2Pi0/PWA/costheta_pm.pdf}
\end{overpic}\begin{overpic}[width=138.76157pt]{fig/PipPim2Pi0/PWA/costheta_z.pdf}
\end{overpic}\begin{overpic}[width=138.76157pt]{fig/PipPim2Pi0/PWA/costheta_p.pdf}
\end{overpic}\begin{overpic}[width=138.76157pt]{fig/PipPim2Pi0/PWA/costheta_m.pdf}
\end{overpic}\begin{overpic}[width=138.76157pt]{fig/PipPim2Pi0/PWA/chi_pm_12.pdf}
\end{overpic}\begin{overpic}[width=138.76157pt]{fig/PipPim2Pi0/PWA/chi_p_m.pdf}
\end{overpic}\begin{overpic}[width=138.76157pt]{fig/PipPim2Pi0/PWA/leg.pdf}
\end{overpic}\begin{overpic}[width=138.76157pt]{fig/PipPim2Pi0/PWA/null.pdf}
\end{overpic}

Figure 6: The mass and angular distributions in
$D^{0}\to\pi^{+}\pi^{-}\pi^{0}\pi^{0}$, where $\theta_{a}^{ab}$ is the
helicity angle of $a$ in the $ab$ system and $\phi_{ab,cd}$ is the angle
between the decay planes of $ab$ and $cd$ systems in the $D^{0}$ rest frame.

##  7 SYSTEMATIC UNCERTAINTY OF THE AMPLITUDE ANALYSIS

The systematic uncertainties of the amplitude analysis come from two aspects.
One is the experimental systematic uncertainty which includes those from
background estimation and detection efficiency over PHSP. Another is the
model-dependent uncertainty which includes those from resonance line shape,
radii of Blatt-Weisskopf barrier factors, quantum correlation parameters,
extra amplitudes, and fit bias. To estimate these uncertainties for each
source, the fit is performed with alternative conditions, and the deviations
from the nominal results are taken as the corresponding uncertainties. Table
12 summarizes the systematic uncertainties on the FFs, resonance parameters,
and $CP$-even fractions in the unit of the statistical uncertainty, where the
total systematic uncertainties are the square roots of the quadrature sums of
the individual contributions. The individual uncertainties are obtained as
follows.

  * •

Background estimation

The sPlot technique is used to estimate the signal weight and background, the
corresponding sources of systematic uncertainty are all the PDFs and the
magnitudes of peaking backgrounds in BKGVI in the 2D fit on the $M_{\rm
bc}^{\rm tag}$ versus $M_{\rm bc}^{\rm sig}$ distribution. To estimate the
corresponding uncertainties, an alternative fit on the $M_{\rm bc}^{\rm tag}$
versus $M_{\rm bc}^{\rm sig}$ distribution is performed by varying the means
and widths of two Gaussian functions by $\pm 1\sigma$ for signal, varying the
fixed parameters of the ARGUS function by $\pm 1\sigma$ for BKGI and BKGII,
changing the Student’s function to the bifurcated Student’s function with
different $n$ and $\sigma$ on the left and right sides of the maximum value
and two ARGUS functions to one ARGUS function on the ($M_{\rm bc}^{\rm
sig}+M_{\rm bc}^{\rm tag}$) dimension for BKGIII, varying the means and widths
of the bifurcated Gaussian functions by $\pm 1\sigma$ for BKGIV and BKGV and
varying the fixed peaking background yields by $\pm 1\sigma$. The output
results are used in the sPlot technique and amplitude analysis.

  * •

Detection efficiency

The systematic uncertainties due to detection efficiencies are from the
$\pi^{\pm}$ tracking/PID and $\pi^{0}$ reconstruction. The corresponding
uncertainties are obtained by weighting the PHSP signal MC sample with a
factor $\epsilon_{\rm Data}/\epsilon_{\rm MC}$ when normalizing the PDF in the
amplitude analysis, where $\epsilon_{\rm Data}$ and $\epsilon_{\rm MC}$ are
the efficiencies of data and MC simulation, respectively. The $\pi^{\pm}$
tracking/PID efficiencies are quoted from Ref. [53], while the efficiency of
the $\pi^{0}$ reconstruction is studied with the control samples of $D^{0}\to
K^{-}\pi^{+}\pi^{0}$ and $D^{0}\to\pi^{+}\pi^{-}\pi^{0}$ versus three tag
modes used in this analysis.

  * •

The parameters and models of resonances

The systematic uncertainties associated with the masses and widths of
resonances in the amplitude analysis are estimated by shifting the
corresponding fixed PDG [23] values or optimized values by $\pm 1\sigma$. The
systematic uncertainties associated with the model of $\pi\pi$ S-wave are
estimated by replacing the K-matrix formula with the sum of three independent
resonances of $f_{0}(500)$, $f_{0}(980)$, and $f_{0}(1370)$, which are
described with the formula in Ref. [54], a Flatt$\rm{\acute{e}}$
parametrization, and a relative Breit-Wigner function with energy dependent
width. The overall uncertainties are the square roots of the quadrature sums
of the individual contributions.

  * •

Radii of Blatt-Weisskopf barrier factors

The systematic uncertainties due to radii of Blatt-Weisskopf barrier factors
are obtained by varying the radii of the $D^{0}$ meson and other resonances by
$\pm 1~{}{\rm{GeV}^{-1}}c$ in the amplitude analysis, and the square root of
the quadrature sum of the individual effects is taken as the uncertainty.

  * •

Quantum correlation correction

The systematic uncertainty due to the quantum correlation correction is
estimated by varying the input quantum correlation parameters by $\pm 1\sigma$
in the amplitude analysis.

  * •

Extra amplitudes

The systematic uncertainties due to the extra amplitudes are estimated with
the alternative fits with the models including an additional amplitude of
$f_{2}(1270)f_{2}(1270)[D]$ or $1^{++}[\rho(770)^{+}\pi^{-}[S]]_{\rm
NR}\pi^{0}$ (“NR” represents non-resonant contribution), which are of
3$\sigma$ to 5$\sigma$ significance based on the nominal model. The square
root of the quadrature sum of the individual effects is taken as the
uncertainty.

  * •

Fit bias

The sPlot technique, MC integration, detection resolution, and many other
potential problems from fit tools may lead to fit bias. To estimate this fit
bias, 100 sets of samples with the same size and signal purity as the data are
generated, where the signal events are generated according to the the nominal
model and the background events are from the inclusive MC sample. The
amplitude fit is performed on these samples, and the average deviations of the
fitted parameters from their input values are taken as uncertainties due to
fit bias.

Table 12:  Systematic uncertainties in the FFs, resonance parameters and
$CP$-even fractions in units of statistical standard deviations. 1: Background
estimation. 2: Detection efficiency. 3: Resonance parameters. 4: Radii of
Blatt-Weisskopf barrier factors. 5: Quantum correlation correction. 6: Extra
amplitudes. 7: Fit bias. For the items with “/”, the quantities before and
after “/” are for $D^{0}\to\pi^{+}\pi^{-}\pi^{+}\pi^{-}$ and
$D^{0}\to\pi^{+}\pi^{-}\pi^{0}\pi^{0}$, respectively.

1 2 3 4 5 6 7 Total FF($D^{0}\to a_{1}(1260)^{+}\pi^{-}$) 0.6/1.1 0.0/0.1
3.8/0.5 3.0/2.1 0.0/0.1 0.2/1.1 0.3/0.1 4.9/2.7 FF($D^{0}\to
a_{1}(1260)^{-}\pi^{+}$) 0.1/0.1 0.1/0.1 0.7/1.5 1.4/1.2 0.1/0.1 0.7/0.7
0.2/0.2 1.7/2.1 FF($D^{0}\to a_{1}(1260)^{0}\pi^{0}$) -/0.4 -/0.1 -/1.8 -/1.7
-/0.0 -/0.9 -/0.2 -/2.6 FF($D^{0}\to a_{1}(1420)^{+}\pi^{-}$) 0.1/0.2 0.0/0.0
0.4/0.7 0.4/0.5 0.0/0.1 0.6/0.5 0.1/0.0 0.8/1.0 FF($D^{0}\to
a_{1}(1640)^{+}\pi^{-}$) 0.7/0.8 0.2/0.2 0.6/0.9 1.1/1.1 0.2/0.2 1.2/1.3
0.2/0.1 1.9/2.1 FF($D^{0}\to a_{1}(1640)^{-}\pi^{+}$) 0.1/0.1 0.1/0.1 0.4/0.2
0.6/0.6 0.2/0.2 0.9/1.0 0.0/0.0 1.2/1.2 FF($D^{0}\to a_{2}(1320)^{+}\pi^{-}$)
0.1/0.2 0.2/0.2 0.5/0.4 0.2/0.1 0.1/0.1 0.3/0.2 0.3/0.3 0.7/0.6 FF($D^{0}\to
a_{2}(1320)^{-}\pi^{+}$) 0.1/0.0 0.0/0.0 0.7/0.3 0.5/0.6 0.1/0.1 0.4/0.5
0.1/0.1 1.0/0.9 FF($D^{0}\to h_{1}(1170)^{0}\pi^{0}$) -/0.7 -/0.1 -/1.1 -/0.7
-/0.0 -/1.0 -/0.2 -/1.8 FF($D^{0}\to\pi(1300)^{+}\pi^{-}$) 0.6/0.9 0.1/0.1
0.9/1.0 1.1/1.2 0.2/0.1 0.6/0.2 0.0/0.0 1.7/1.8
FF($D^{0}\to\pi(1300)^{-}\pi^{+}$) 0.1/0.4 0.2/0.2 1.1/1.4 1.2/1.5 0.1/0.1
0.5/0.1 0.0/0.1 1.8/2.1 FF($D^{0}\to\pi(1300)^{0}\pi^{0}$) -/0.5 -/0.2 -/0.9
-/0.6 -/0.0 -/0.2 -/0.0 -/1.2 FF($D^{0}\to\pi_{2}(1670)^{0}\pi^{0}$) -/0.8
-/0.1 -/0.3 -/0.4 -/0.1 -/1.4 -/0.1 -/1.7
FF($D^{0}\to\rho(770)^{0}\rho(770)^{0}$) 0.3/- 0.0/- 1.5/- 0.2/- 0.1/- 0.6/-
0.0/- 1.7/- FF($D^{0}\to\rho(770)^{0}\rho(770)^{0}[S]$) 0.1/- 0.0/- 0.7/-
0.3/- 0.1/- 0.3/- 0.1/- 0.8/- FF($D^{0}\to\rho(770)^{0}\rho(770)^{0}[P]$)
0.0/- 0.0/- 0.5/- 0.5/- 0.3/- 0.1/- 0.4/- 0.9/-
FF($D^{0}\to\rho(770)^{0}\rho(770)^{0}[D]$) 0.3/- 0.0/- 0.9/- 0.3/- 0.1/-
0.5/- 0.2/- 1.1/- FF($D^{0}\to\rho(770)^{0}\rho(1450)^{0}$) 0.0/- 0.0/- 0.7/-
0.7/- 0.1/- 0.8/- 0.0/- 1.2/- FF($D^{0}\to\rho(770)^{0}\rho(1450)^{0}[P]$)
0.1/- 0.2/- 1.3/- 0.6/- 0.3/- 0.1/- 0.2/- 1.5/-
FF($D^{0}\to\rho(770)^{0}\rho(1450)^{0}[D]$) 0.1/- 0.1/- 0.2/- 0.9/- 0.1/-
0.8/- 0.1/- 1.2/- FF($D^{0}\to\rho(770)^{+}\rho(770)^{-}$) -/0.9 -/0.0 -/1.3
-/0.5 -/0.1 -/1.0 -/0.3 -/2.0 FF($D^{0}\to\rho(770)^{+}\rho(770)^{-}[S]$)
-/0.6 -/0.1 -/1.1 -/0.6 -/0.0 -/1.0 -/0.1 -/1.7
FF($D^{0}\to\rho(770)^{+}\rho(770)^{-}[P]$) -/1.0 -/0.2 -/0.5 -/0.4 -/0.1
-/0.8 -/0.1 -/1.5 FF($D^{0}\to\rho(770)^{+}\rho(770)^{-}[D]$) -/0.1 -/0.1
-/0.4 -/0.3 -/0.1 -/0.6 -/0.2 -/0.8
FF($D^{0}\to\rho(770)^{+}\rho(1450)^{-}[D]$) -/0.7 -/0.1 -/1.9 -/0.8 -/0.1
-/0.7 -/0.4 -/2.3 FF($D^{0}\to\rho(770)^{0}(\pi\pi)_{S}$) 0.5/0.3 0.0/0.0
2.7/1.8 0.1/0.1 0.0/0.0 1.0/0.7 0.2/0.2 3.0/2.0
FF($D^{0}\to(\pi^{+}\pi^{-})_{S}(\pi\pi)_{S}$) 0.1/0.6 0.0/0.0 2.0/1.3 0.6/0.8
0.1/0.0 0.3/0.5 0.1/0.2 2.1/1.7 FF($D^{0}\to f_{2}(1270)^{0}(\pi\pi)_{S}$)
0.0/0.1 0.0/0.0 2.6/2.9 1.1/1.0 0.0/0.0 1.4/1.6 0.1/0.1 3.1/3.5
FF($D^{0}\to\omega(782)\pi^{0}$) -/0.1 -/0.0 -/0.3 -/0.1 -/0.0 -/0.4 -/0.1
-/0.5 FF($D^{0}\to\phi(1020)\pi^{0}$) -/0.2 -/0.1 -/0.3 -/0.2 -/0.1 -/0.5
-/0.4 -/0.8 FF($a_{1}(1260)^{\pm}\to\rho(770)\pi[S]$) 0.8/0.8 0.1/0.1 1.2/1.2
0.4/0.4 0.1/0.1 0.6/0.6 0.1/0.0 1.6/1.6
FF($a_{1}(1260)^{\pm}\to\rho(770)\pi[D]$) 0.1/0.1 0.0/0.0 2.2/2.1 0.5/0.4
0.0/0.0 0.1/0.1 0.2/0.2 2.2/2.2 FF($a_{1}(1260)^{\pm}\to f_{2}(1270)\pi[P]$)
0.7/0.7 0.0/0.0 0.4/0.3 1.2/1.2 0.1/0.1 0.3/0.3 0.1/0.1 1.5/1.5
FF($a_{1}(1260)^{\pm}\to(\pi^{+}\pi^{-})_{S}\pi[P]$) 0.8/0.9 0.1/0.0 1.0/1.6
0.3/0.1 0.2/0.1 0.9/0.8 0.2/0.2 1.6/2.0
FF($a_{1}(1260)^{0}\to\rho(770)\pi[S]$) -/0.7 -/0.1 -/1.4 -/0.3 -/0.1 -/0.6
-/0.0 -/1.7 FF($a_{1}(1260)^{0}\to\rho(770)\pi[D]$) -/0.1 -/0.0 -/2.2 -/0.4
-/0.0 -/0.1 -/0.2 -/2.2 FF($a_{1}(1260)^{0}\to f_{2}(1270)\pi[P]$) -/0.7 -/0.0
-/0.4 -/1.4 -/0.1 -/0.3 -/0.1 -/1.7
FF($a_{1}(1260)^{0}\to(\pi^{+}\pi^{-})_{S}\pi[P]$) -/0.8 -/0.0 -/1.6 -/0.1
-/0.1 -/0.8 -/0.3 -/2.0 FF($\pi(1300)^{\pm}\to\rho(770)\pi$) 0.0/0.2 0.1/0.1
2.0/0.8 2.2/3.3 0.1/0.1 0.4/0.3 0.3/0.3 3.1/3.4
FF($\pi(1300)^{\pm}\to(\pi^{+}\pi^{-})_{S}\pi$) 0.1/0.1 0.1/0.1 1.8/1.2
2.3/2.9 0.1/0.1 0.3/0.2 0.3/0.4 3.0/3.1 FF($\pi(1300)^{0}\to\rho(770)\pi$)
-/0.2 -/0.1 -/1.0 -/3.3 -/0.1 -/0.4 -/0.2 -/3.5
FF($\pi(1300)^{0}\to(\pi^{+}\pi^{-})_{S}\pi$) -/0.0 -/0.1 -/1.9 -/2.6 -/0.1
-/0.1 -/0.4 -/3.3 $M_{a_{1}(1260)}$ 0.5 0.1 4.0 1.2 0.1 1.6 0.0 4.5
$\Gamma_{a_{1}(1260)}$ 1.7 0.2 3.2 2.3 0.4 1.7 0.0 4.6 $M_{\pi(1300)}$ 0.8 0.1
1.1 1.1 0.0 1.0 0.0 2.0 $\Gamma_{\pi(1300)}$ 0.7 0.1 2.9 0.4 0.2 0.5 0.0 3.1
$F^{+}$ 0.7/1.0 0.0/0.2 0.7/0.8 0.3/0.8 0.1/0.1 1.0/0.6 0.1/0.1 1.4/1.6

##  8 MEASUREMENT OF BRANCHING FRACTIONS

The absolute branching fractions of $D^{0}\to\pi^{+}\pi^{-}\pi^{+}\pi^{-}$ and
$D^{0}\to\pi^{+}\pi^{-}\pi^{0}\pi^{0}$(non-$\eta$) are measured with the DT
method. The numbers of ST events ($N_{g}^{\rm ST}$) for the tag mode $g$ and
DT events ($N_{fg}^{\rm DT}$) for the self-conjugated signal model $f$ with
the tag mode $g$ are given by

$\displaystyle N_{g}^{\rm ST}=$ $\displaystyle
2N_{D^{0}\bar{D}^{0}}\mathcal{B}_{g}\epsilon_{g}^{ST}(1+y_{D}^{2})(1+r_{g}^{2}-2r_{g}R_{g}y_{D}\mathrm{cos}\delta_{g}),$
(43) $\displaystyle N_{fg}^{\rm DT}=$ $\displaystyle
2N_{D^{0}\bar{D}^{0}}\mathcal{B}_{f}\mathcal{B}_{g}\epsilon_{fg}^{DT}(1+y_{D}^{2})[1+r_{g}^{2}-2r_{g}R_{g}\mathrm{cos}\delta_{g}(2F_{+}^{f}-1)].$
(44)

Here, $N_{D^{0}\bar{D}^{0}}$ is the total number of $D^{0}\bar{D}^{0}$ pairs
in data, $\mathcal{B}_{f}$ and $\mathcal{B}_{g}$ are the branching fractions
of the signal mode $f$ and tag mode $g$, respectively, $\epsilon_{f}^{\rm ST}$
and $\epsilon_{fg}^{\rm ST}$ are the corresponding ST and DT efficiencies, and
$y_{D}$ is the $D^{0}$-$\bar{D}^{0}$ mixing parameter. By combining Eqs. (43)
and (44) and ignoring the term $2r_{g}R_{g}y_{D}\mathrm{cos}\delta_{g}$ in Eq.
(43), the branching fraction of $D^{0}\to f$ is given by

$\displaystyle\mathcal{B}_{f}=\frac{\sum_{g}N_{fg}^{\rm
DT}}{\sum_{g}N_{g}^{\rm ST}(\epsilon_{fg}^{\rm DT}/\epsilon_{g}^{\rm
ST})\left[1-\frac{2r_{g}R_{g}\mathrm{cos}\delta_{g}}{1+r_{g}^{2}}(2F_{+}^{f}-1)\right]}.$
(45)

As described in Sec. 4, the ST and DT yields are extracted by performing an
unbinned maximum likelihood fit on the $M_{\rm bc}^{\rm tag}$ distribution and
the 2D distribution of $M_{\rm bc}^{\rm tag}$ versus $M_{\rm bc}^{\rm sig}$,
respectively. The corresponding ST and DT efficiencies are obtained with the
similar fit processes on the signal MC sample as described in Sec. 2. The
corresponding ST and DT yields as well as the ST and DT efficiencies are
summarized in Table 3. According to Eq. (45) and the values in Table 3, the
branching fractions of $D^{0}\to\pi^{+}\pi^{-}\pi^{+}\pi^{-}$ and
$D^{0}\to\pi^{+}\pi^{-}\pi^{0}\pi^{0}$(non-$\eta$) are calculated and
summarized in Table 13. According to the FFs obtained in the amplitude
analysis, the branching fractions of the intermediate processes are also
determined. The obtained results are summarized in Table 14.

Table 13: Branching fractions of $D^{0}\to\pi^{+}\pi^{-}\pi^{+}\pi^{-}$ and $D^{0}\to\pi^{+}\pi^{-}\pi^{0}\pi^{0}$(non-$\eta$) from this work compared to PDG [23] values, where the first uncertainties are statistical and the second systematic. | This work | PDG
---|---|---
$D^{0}\to\pi^{+}\pi^{-}\pi^{+}\pi^{-}$ | $(0.688\pm 0.010\pm 0.010)\%$ | $(0.756\pm 0.020)\%$
$D^{0}\to\pi^{+}\pi^{-}\pi^{0}\pi^{0}$(non-$\eta$) | $(0.951\pm 0.025\pm 0.021)\%$ | $(1.005\pm 0.090)\%$
Table 14:  Branching fractions of the intermediate processes with FFs$>$1% in $D^{0}\to\pi^{+}\pi^{-}\pi^{+}\pi^{-}$ and $D^{0}\to\pi^{+}\pi^{-}\pi^{0}\pi^{0}$.  Component | Branching fraction (%)
---|---
$\pi^{+}\pi^{-}\pi^{+}\pi^{-}$ | $\pi^{+}\pi^{-}\pi^{0}\pi^{0}$
$D^{0}\to a_{1}(1260)^{+}\pi^{-}$ | $0.566\pm 0.024\pm 0.110$ | $0.546\pm 0.027\pm 0.070$
$D^{0}\to a_{1}(1260)^{-}\pi^{+}$ | $0.071\pm 0.010\pm 0.017$ | $0.068\pm 0.011\pm 0.021$
$D^{0}\to a_{1}(1260)^{0}\pi^{0}$ | - | $0.313\pm 0.031\pm 0.082$
$D^{0}\to a_{1}(1640)^{+}\pi^{-}$ | $0.012\pm 0.003\pm 0.006$ | $0.010\pm 0.003\pm 0.007$
$D^{0}\to h_{1}(1170)^{0}\pi^{0}$ | - | $0.012\pm 0.006\pm 0.010$
$D^{0}\to\pi(1300)^{+}\pi^{-}$ | $0.222\pm 0.018\pm 0.031$ | $0.148\pm 0.014\pm 0.025$
$D^{0}\to\pi(1300)^{-}\pi^{+}$ | $0.162\pm 0.016\pm 0.028$ | $0.108\pm 0.011\pm 0.021$
$D^{0}\to\pi(1300)^{0}\pi^{0}$ | - | $0.221\pm 0.027\pm 0.034$
$D^{0}\to\pi_{2}(1670)^{0}\pi^{0}$ | - | $0.010\pm 0.002\pm 0.004$
$D^{0}\to\rho(770)^{0}\rho(770)^{0}$ | $0.193\pm 0.013\pm 0.022$ | -
$D^{0}\to\rho(770)^{0}\rho(770)^{0}[S]$ | $0.012\pm 0.004\pm 0.003$ | -
$D^{0}\to\rho(770)^{0}\rho(770)^{0}[P]$ | $0.067\pm 0.007\pm 0.006$ | -
$D^{0}\to\rho(770)^{0}\rho(770)^{0}[D]$ | $0.159\pm 0.015\pm 0.017$ | -
$D^{0}\to\rho(770)^{0}\rho(1450)^{0}$ | $0.017\pm 0.006\pm 0.008$ | -
$D^{0}\to\rho(770)^{0}\rho(1450)^{0}[P]$ | $0.007\pm 0.003\pm 0.003$ | -
$D^{0}\to\rho(770)^{0}\rho(1450)^{0}[D]$ | $0.010\pm 0.006\pm 0.008$ | -
$D^{0}\to\rho(770)^{+}\rho(770)^{-}$ | - | $0.864\pm 0.040\pm 0.077$
$D^{0}\to\rho(770)^{+}\rho(770)^{-}[S]$ | - | $0.124\pm 0.019\pm 0.033$
$D^{0}\to\rho(770)^{+}\rho(770)^{-}[P]$ | - | $0.186\pm 0.013\pm 0.019$
$D^{0}\to\rho(770)^{+}\rho(770)^{-}[D]$ | - | $0.342\pm 0.029\pm 0.025$
$D^{0}\to\rho(770)^{+}\rho(1450)^{-}[D]$ | - | $0.016\pm 0.008\pm 0.016$
$D^{0}\to\rho(770)^{0}(\pi\pi)_{S}$ | $0.019\pm 0.004\pm 0.012$ | $0.010\pm 0.002\pm 0.004$
$D^{0}\to(\pi^{+}\pi^{-})_{S}(\pi\pi)_{S}$ | $0.432\pm 0.032\pm 0.066$ | $0.356\pm 0.029\pm 0.049$
$D^{0}\to f_{2}(1270)^{0}(\pi\pi)_{S}$ | $0.012\pm 0.003\pm 0.008$ | $0.010\pm 0.002\pm 0.008$
$D^{0}\to\omega(782)\pi^{0}$ | - | $0.009\pm 0.004\pm 0.002$
$D^{0}\to\phi(1020)\pi^{0}$ | - | $0.014\pm 0.004\pm 0.003$

Benefiting from the DT method, most of the systematic uncertainties associated
with the ST selection are cancelled. The relative systematic uncertainties in
the branching fraction measurement are summarized in Table 15, where the total
uncertainties are the square roots of the quadrature sums of the individual
ones. Details of the systematic uncertainties in the branching fraction
measurements are described below.

  * •

Tracking efficiency

The tracking efficiency of $\pi^{\pm}$ has been studied in Ref. [53]. The
systematic uncertainties are assigned by re-weighting the data-MC difference
in $\pi^{\pm}$ tracking efficiencies according to the momentum distribution of
$\pi^{\pm}$ in the signal MC sample. They are assigned to be 0.4% and 0.2% for
$D^{0}\to\pi^{+}\pi^{-}\pi^{+}\pi^{-}$ and
$D^{0}\to\pi^{+}\pi^{-}\pi^{0}\pi^{0}$, respectively.

  * •

PID efficiency

The PID efficiency of $\pi^{\pm}$ has been studied in Ref. [53]. The
systematic uncertainties are assigned by re-weighting the data-MC difference
in $\pi^{\pm}$ PID efficiencies according to the momentum distribution of
$\pi^{\pm}$ in the signal MC sample. They are assigned to be 1.2% and 0.6% for
$D^{0}\to\pi^{+}\pi^{-}\pi^{+}\pi^{-}$ and
$D^{0}\to\pi^{+}\pi^{-}\pi^{0}\pi^{0}$, respectively.

  * •

$\pi^{0}$ reconstruction efficiency

The $\pi^{0}$ reconstruction efficiency is studied by the control samples
$\bar{D}^{0}\to K^{+}\pi^{-}$, $K^{+}\pi^{-}\pi^{0}$,
$K^{+}\pi^{-}\pi^{+}\pi^{-}$ versus $D^{0}\to K^{-}\pi^{+}\pi^{0}$,
$\pi^{+}\pi^{-}\pi^{0}$. The corresponding systematic uncertainty for
$D^{0}\to\pi^{+}\pi^{-}\pi^{0}\pi^{0}$ is assigned to be 1.3%, by re-weighting
the data-MC difference in $\pi^{0}$ reconstruction efficiencies according to
the momentum distribution of $\pi^{0}$ in the signal MC sample.

  * •

$\Delta E$ requirement in signal side

To study the systematic uncertainty from the $\Delta E$ requirement of the
signal side, the $\Delta E$ of the signal MC sample is smeared by a Gaussian
function, while the mean and width are obtained by performing a fit to the
$\Delta E$ distribution of data with the signal MC shape convolved with a
Gaussian function. The resultant changes of the efficiencies with respect to
the nominal values, 0.1% for both $D^{0}\to\pi^{+}\pi^{-}\pi^{+}\pi^{-}$ and
$D^{0}\to\pi^{+}\pi^{-}\pi^{0}\pi^{0}$, are taken as the systematic
uncertainties.

  * •

ST fit

The systematic uncertainty of the ST fit is studied by changing fit range,
signal shape and background shape. The fit range is enlarged by $0.005~{}{\rm
GeV}/c^{2}$. The signal shape is changed by replacing the Gaussian resolution
with a Crystal Ball function [55]. The background shape is changed with a
floating cut-off parameter for the ARGUS function. The square root of the
quadrature sum of the relative change of ST yields, which gives 0.5%, is taken
as the systematic uncertainty.

  * •

DT fit

As described in Sec. 7, the systematic uncertainties of the DT fit are
obtained by varying the signal and background shapes, and the yields of
peaking backgrounds in BKGVI in the fit. The square roots of the quadrature
sums of the relative changes of DT yields in the individual changes are 0.6%
and 1.4% for $D^{0}\to\pi^{+}\pi^{-}\pi^{+}\pi^{-}$ and
$D^{0}\to\pi^{+}\pi^{-}\pi^{0}\pi^{0}$, respectively, which are taken as the
systematic uncertainties.

  * •

Amplitude model

The DT efficiencies are obtained by varying the parameters of the amplitude
model within their uncertainties. The resultant standard deviations of DT
efficiencies, which are 0.3% and 0.6% for
$D^{0}\to\pi^{+}\pi^{-}\pi^{+}\pi^{-}$ and
$D^{0}\to\pi^{+}\pi^{-}\pi^{0}\pi^{0}$, respectively, are taken as the
systematic uncertainties.

  * •

Quantum correlation corrections

The uncertainties associated with the input parameters of quantum correlation
($r$, $R$, $\delta$, and $F_{+}$) are 0.3% and 0.4% for
$D^{0}\to\pi^{+}\pi^{-}\pi^{+}\pi^{-}$ and
$D^{0}\to\pi^{+}\pi^{-}\pi^{0}\pi^{0}$, respectively, which are assigned
according to uncertainty propagation.

  * •

MC statistics

The systematic uncertainties related with the limited statistics of the signal
MC samples are 0.1% and 0.2% for $D^{0}\to\pi^{+}\pi^{-}\pi^{+}\pi^{-}$ and
$D^{0}\to\pi^{+}\pi^{-}\pi^{0}\pi^{0}$, respectively.

Table 15: Relative systematic uncertainties (%) of the measured branching fractions. Source | $\mathcal{B}(D^{0}\to\pi^{+}\pi^{-}\pi^{+}\pi^{-})$ | $\mathcal{B}(D^{0}\to\pi^{+}\pi^{-}\pi^{0}\pi^{0})$
---|---|---
$\pi^{\pm}$ tracking | 0.4 | 0.2
$\pi^{\pm}$ PID | 1.2 | 0.6
$\pi^{0}$ reconstruction | - | 1.3
$\Delta E$ cut in signal side | 0.1 | 0.1
ST fit | 0.5 | 0.5
DT fit | 0.6 | 1.4
Amplitude model | 0.3 | 0.6
Quantum correlation correction | 0.3 | 0.4
MC statistics | 0.1 | 0.2
Total | 1.5 | 2.2

##  9 SUMMARY

Using 2.93 $\rm fb^{-1}$ of $e^{+}e^{-}$ collision data taken at
$\sqrt{s}=3.773~{}\rm{GeV}$ with the BESIII detector, a joint amplitude
analysis of $D^{0}\to\pi^{+}\pi^{-}\pi^{+}\pi^{-}$ and
$D^{0}\to\pi^{+}\pi^{-}\pi^{0}\pi^{0}$(non-$\eta$) is performed. Large
interferences between the dominant amplitudes of $D^{0}\to a_{1}(1260)\pi$,
$D^{0}\to\pi(1300)\pi$, $D^{0}\to\rho(770)\rho(770)$ and $D^{0}\to
2(\pi\pi)_{S}$ are observed. Based on the amplitude model, the model dependent
$CP$-even fractions of $D^{0}\to\pi^{+}\pi^{-}\pi^{+}\pi^{-}$ and
$D^{0}\to\pi^{+}\pi^{-}\pi^{0}\pi^{0}$(non-$\eta$) are determined to be
$(75.2\pm 1.1_{\rm stat.}\pm 1.5_{\rm syst.})\%$ and $(68.9\pm 1.5_{\rm
stat.}\pm 2.4_{\rm syst.})\%$, respectively, which are consistent with the
previous measurements carried out by the CLEO [11, 50, 6] and BESIII [51, 52]
collaboration. The branching fractions of
$D^{0}\to\pi^{+}\pi^{-}\pi^{+}\pi^{-}$ and
$D^{0}\to\pi^{+}\pi^{-}\pi^{0}\pi^{0}$(non-$\eta$) are measured to be
$(0.688\pm 0.010_{\rm stat.}\pm 0.010_{\rm syst.})\%$ and $(0.951\pm
0.025_{\rm stat.}\pm 0.021_{\rm syst.})\%$, respectively, where the former one
is $3\sigma$ lower than the PDG [23] value. This $3\sigma$ deviation may be
due to the differences in the amplitude models which affect the global
reconstruction efficiency. These results provide essential information for
search of $CP$ violation [7] and measurements of the binned strong phase
parameter [6], which are important inputs in the $\gamma$ measurement via the
$B^{-}\to DK^{-}$ decay.

The BESIII Collaboration thanks the staff of BEPCII, the IHEP computing center
and the supercomputing center of the University of Science and Technology of
China (USTC) for their strong support.

## References

  * [1] Cabibbo, N.,  Phys. Rev. Lett. 10 (1963) 531.
  * [2] Kobayashi, M. and Maskawa, T.,  Prog. Theor. Phys. 49 (1973) 652.
  * [3] Gronau, M. and Wyler, D.,  Phys. Lett. B 265 (1991) 172.
  * [4] Atwood, D., Dunietz, I., and Soni, A.,  Phys. Rev. Lett. 78 (1997) 3257.
  * [5] Giri, A., Grossman, Y., Soffer, A., and Zupan, J.,  Phys. Rev. D 68 (2003) 054018.
  * [6] Harnew, S., Naik, P., Prouve, C., Rademacker, J., and Asner, D.,  JHEP 01 (2018) 144.
  * [7] Aaij, R. et al.,  Phys. Lett. B 769 (2017) 345.
  * [8] Cheng, H.-Y.,  Phys. Rev. D 67 (2003) 094007.
  * [9] Cheng, H.-Y. and Chiang, C.-W.,  Phys. Rev. D 81 (2010) 074021.
  * [10] Link, J. M. et al.,  Phys. Rev. D 75 (2007) 052003.
  * [11] d’Argent, P. et al.,  JHEP 05 (2017) 143.
  * [12] Baltrusaitis, R. M. et al.,  Phys. Rev. Lett. 56 (1986) 2140.
  * [13] Li, H.-B. and Lyu, X.-R.,  Natl. Sci. Rev. 8 (2021) nwab181.
  * [14] Ablikim, M. et al.,  Nucl. Instrum. Meth. A 614 (2010) 345.
  * [15] Yu, C. et al.,  BEPCII Performance and Beam Dynamics Studies on Luminosity,  in 7th International Particle Accelerator Conference, page TUYA01, 2016.
  * [16] Ablikim, M. et al.,  Chin. Phys. C 44 (2020) 040001.
  * [17] Agostinelli, S. et al.,  Nucl. Instrum. Meth. A 506 (2003) 250.
  * [18] Huang, K.-X. et al.,  Nucl. Sci. Tech. 33 (2022) 142.
  * [19] Jadach, S., Ward, B. F. L., and Was, Z.,  Phys. Rev. D 63 (2001) 113009.
  * [20] Jadach, S., Ward, B. F. L., and Was, Z.,  Comput. Phys. Commun. 130 (2000) 260.
  * [21] Lange, D. J.,  Nucl. Instrum. Meth. A 462 (2001) 152.
  * [22] Ping, R.-G.,  Chin. Phys. C 32 (2008) 599.
  * [23] Workman, R. L. et al.,  PTEP 2022 (2022) 083C01.
  * [24] Chen, J. C., Huang, G. S., Qi, X. R., Zhang, D. H., and Zhu, Y. S.,  Phys. Rev. D 62 (2000) 034003.
  * [25] Yang, R.-L., Ping, R.-G., and Chen, H.,  Chin. Phys. Lett. 31 (2014) 061301.
  * [26] Richter-Was, E.,  Phys. Lett. B 303 (1993) 163.
  * [27] Ablikim, M. et al.,  Phys. Lett. B 734 (2014) 227.
  * [28] Albrecht, H. et al.,  Phys. Lett. B 241 (1990) 278.
  * [29] Amhis, Y. S. et al.,  Eur. Phys. J. C 81 (2021) 226.
  * [30] Ablikim, M. et al.,  JHEP 05 (2021) 164.
  * [31] Berger, N., Liu, B., and Wang, J.,  J. Phys. Conf. Ser. 219 (2010) 042031.
  * [32] Von Hippel, F. and Quigg, C.,  Phys. Rev. D 5 (1972) 624.
  * [33] Mandelstam, S., Paton, J. E., Peierls, R. F., and Sarker, A. Q.,  Annals Phys. 18 (1962) 198.
  * [34] Herndon, D., Soding, P., and Cashmore, R. J.,  Phys. Rev. D 11 (1975) 3165.
  * [35] Brehm, J. J.,  Annals Phys. 108 (1977) 454.
  * [36] Rarita, W. and Schwinger, J.,  Phys. Rev. 60 (1941) 61.
  * [37] Zemach, C.,  Phys. Rev. 140 (1965) B97.
  * [38] Chung, S. U.,  Phys. Rev. D 57 (1998) 431.
  * [39] Zou, B. S. and Bugg, D. V.,  Eur. Phys. J. A 16 (2003) 537.
  * [40] Gounaris, G. J. and Sakurai, J. J.,  Phys. Rev. Lett. 21 (1968) 244.
  * [41] Ablikim, M. et al.,  Phys. Lett. B 607 (2005) 243.
  * [42] Anisovich, V. V. and Sarantsev, A. V.,  Eur. Phys. J. A 16 (2003) 229.
  * [43] Aubert, B. et al.,  Phys. Rev. D 78 (2008) 034023.
  * [44] Alexeev, G. D. et al.,  Phys. Rev. Lett. 127 (2021) 082501.
  * [45] Aghasyan, M. et al.,  Phys. Rev. D 98 (2018) 092003.
  * [46] Pivk, M. and Le Diberder, F. R.,  Nucl. Instrum. Meth. A 555 (2005) 356.
  * [47] James, F. and Roos, M.,  Comput. Phys. Commun. 10 (1975) 343.
  * [48] Efron, B.,  The Annals of Statistics 7 (1979) 1 .
  * [49] Langenbruch, C.,  Eur. Phys. J. C 82 (2022) 393.
  * [50] Malde, S. et al.,  Phys. Lett. B 747 (2015) 9.
  * [51] Ablikim, M. et al.,  Phys. Rev. D 106 (2022) 092004.
  * [52] Ablikim, M. et al.,  Phys. Rev. D 106 (2022) 092005.
  * [53] Ablikim, M. et al.,  Phys. Rev. D 97 (2018) 072004.
  * [54] Bugg, D. V.,  J. Phys. G 34 (2007) 151.
  * [55] Skwarnicki, T.,  A study of the radiative CASCADE transitions between the Upsilon-Prime and Upsilon resonances,  PhD thesis, Cracow, INP, 1986.