Analysis of the decay D 0 → K 0 S K + K −
BESIII Collaboration, M. Ablikim, M. N. Achasov, P. Adlarson, S. Ahmed, M. Albrecht, M. Alekseev, A. Amoroso, F. F. An, Q. An, Y. Bai, O. Bakina, R. Baldini Ferroli, I. Balossino, Y. Ban, K. Begzsuren, J. V. Bennett, N. Berger, M. Bertani, D. Bettoni, F. Bianchi, J. Biernat, J. Bloms, I. Boyko, R. A. Briere, H. Cai, X. Cai, A. Calcaterra, G. F. Cao, N. Cao, S. A. Cetin, J. Chai, J. F. Chang, W. L. Chang, G. Chelkov, D. Y. Chen, G. Chen, H. S. Chen, J. Chen, M. L. Chen, S. J. Chen, Y. B. Chen, W. Cheng, G. Cibinetto, F. Cossio, X. F. Cui, H. L. Dai, J. P. Dai, X. C. Dai, A. Dbeyssi, D. Dedovich, Z. Y. Deng, A. Denig, I. Denysenko, M. Destefanis, F. De Mori, Y. Ding, C. Dong, J. Dong, L. Y. Dong, M. Y. Dong, Z. L. Dou, S. X. Du, J. Z. Fan, J. Fang, S. S. Fang, Y. Fang, R. Farinelli, L. Fava, F. Feldbauer, G. Felici, C. Q. Feng, M. Fritsch, C. D. Fu, Y. Fu, Q. Gao, X. L. Gao, Y. Gao, Y. Gao, Y. G. Gao, B. Garillon, I. Garzia, E. M. Gersabeck, A. Gilman, K. Goetzen, L. Gong, W. X. Gong, W. Gradl, M. Greco, L. M. Gu, M. H. Gu, S. Gu, Y. T. Gu, A. Q. Guo, L. B. Guo, R. P. Guo, Y. P. Guo, A. Guskov, S. Han, X. Q. Hao, et al. (386 additional authors not shown)
AAnalysis of the decay D → K S K + K − (Dated: June 5, 2020)M. Ablikim , M. N. Achasov ,d , P. Adlarson , S. Ahmed , M. Albrecht , M. Alekseev A, C , A. Amoroso A, C , F. F. An ,Q. An , , Y. Bai , O. Bakina , R. Baldini Ferroli A , I. Balossino A , Y. Ban ,l , K. Begzsuren , J. V. Bennett ,N. Berger , M. Bertani A , D. Bettoni A , F. Bianchi A, C , J Biernat , J. Bloms , I. Boyko , R. A. Briere , H. Cai ,X. Cai , , A. Calcaterra A , G. F. Cao , , N. Cao , , S. A. Cetin B , J. Chai C , J. F. Chang , , W. L. Chang , ,G. Chelkov ,b,c , D. Y. Chen , G. Chen , H. S. Chen , , J. Chen , M. L. Chen , , S. J. Chen , Y. B. Chen , ,W. Cheng C , G. Cibinetto A , F. Cossio C , X. F. Cui , H. L. Dai , , J. P. Dai ,h , X. C. Dai , , A. Dbeyssi ,D. Dedovich , Z. Y. Deng , A. Denig , I. Denysenko , M. Destefanis A, C , F. De Mori A, C , Y. Ding , C. Dong ,J. Dong , , L. Y. Dong , , M. Y. Dong , , , Z. L. Dou , S. X. Du , J. Z. Fan , J. Fang , , S. S. Fang , , Y. Fang ,R. Farinelli A, B , L. Fava B, C , F. Feldbauer , G. Felici A , C. Q. Feng , , M. Fritsch , C. D. Fu , Y. Fu , Q. Gao ,X. L. Gao , , Y. Gao , Y. Gao , Y. G. Gao , B. Garillon , I. Garzia A , E. M. Gersabeck , A. Gilman , K. Goetzen ,L. Gong , W. X. Gong , , W. Gradl , M. Greco A, C , L. M. Gu , M. H. Gu , , S. Gu , Y. T. Gu , A. Q. Guo ,L. B. Guo , R. P. Guo , Y. P. Guo , A. Guskov , S. Han , X. Q. Hao , F. A. Harris , K. L. He , , F. H. Heinsius ,T. Held , Y. K. Heng , , , M. Himmelreich ,g , Y. R. Hou , Z. L. Hou , H. M. Hu , , J. F. Hu ,h , T. Hu , , , Y. Hu ,G. S. Huang , , J. S. Huang , X. T. Huang , X. Z. Huang , N. Huesken , T. Hussain , W. Ikegami Andersson ,W. Imoehl , M. Irshad , , Q. Ji , Q. P. Ji , X. B. Ji , , X. L. Ji , , H. L. Jiang , X. S. Jiang , , , X. Y. Jiang ,J. B. Jiao , Z. Jiao , D. P. Jin , , , S. Jin , Y. Jin , T. Johansson , N. Kalantar-Nayestanaki , X. S. Kang ,R. Kappert , M. Kavatsyuk , B. C. Ke , I. K. Keshk , A. Khoukaz , P. Kiese , R. Kiuchi , R. Kliemt , L. Koch ,O. B. Kolcu B,f , B. Kopf , M. Kuemmel , M. Kuessner , A. Kupsc , M. Kurth , M. G. Kurth , , W. Kühn , J. S. Lange ,P. Larin , L. Lavezzi C , H. Leithoff , T. Lenz , C. Li , C. H. Li , Cheng Li , , D. M. Li , F. Li , , G. Li , H. B. Li , ,H. J. Li ,j , J. C. Li , Ke Li , L. K. Li , Lei Li , P. L. Li , , P. R. Li , W. D. Li , , W. G. Li , X. H. Li , , X. L. Li ,X. N. Li , , Z. B. Li , Z. Y. Li , H. Liang , , H. Liang , , Y. F. Liang , Y. T. Liang , G. R. Liao , L. Z. Liao , ,J. Libby , C. X. Lin , D. X. Lin , Y. J. Lin , B. Liu ,h , B. J. Liu , C. X. Liu , D. Liu , , D. Y. Liu ,h , F. H. Liu ,Fang Liu , Feng Liu , H. B. Liu , H. M. Liu , , Huanhuan Liu , Huihui Liu , J. B. Liu , , J. Y. Liu , , K. Liu ,K. Y. Liu , Ke Liu , L. Y. Liu , Q. Liu , S. B. Liu , , T. Liu , , X. Liu , X. Y. Liu , , Y. B. Liu , Z. A. Liu , , ,Zhiqing Liu , Y. F. Long ,l , X. C. Lou , , , H. J. Lu , J. D. Lu , , J. G. Lu , , Y. Lu , Y. P. Lu , , C. L. Luo ,M. X. Luo , P. W. Luo , T. Luo ,j , X. L. Luo , , S. Lusso C , X. R. Lyu , F. C. Ma , H. L. Ma , L. L. Ma ,M. M. Ma , , Q. M. Ma , X. N. Ma , X. X. Ma , , X. Y. Ma , , Y. M. Ma , F. E. Maas , M. Maggiora A, C ,S. Maldaner , S. Malde , Q. A. Malik , A. Mangoni B , Y. J. Mao ,l , Z. P. Mao , S. Marcello A, C , Z. X. Meng ,J. G. Messchendorp , G. Mezzadri A , J. Min , , T. J. Min , R. E. Mitchell , X. H. Mo , , , Y. J. Mo , C. MoralesMorales , N. Yu. Muchnoi ,d , H. Muramatsu , A. Mustafa , S. Nakhoul ,g , Y. Nefedov , F. Nerling ,g , I. B. Nikolaev ,d ,Z. Ning , , S. Nisar ,k , S. L. Niu , , S. L. Olsen , Q. Ouyang , , , S. Pacetti B , Y. Pan , , M. Papenbrock ,P. Patteri A , M. Pelizaeus , H. P. Peng , , K. Peters ,g , J. Pettersson , J. L. Ping , R. G. Ping , , A. Pitka , R. Poling ,V. Prasad , , M. Qi , S. Qian , , C. F. Qiao , X. P. Qin , X. S. Qin , Z. H. Qin , , J. F. Qiu , S. Q. Qu ,K. H. Rashid ,i , K. Ravindran , C. F. Redmer , M. Richter , A. Rivetti C , V. Rodin , M. Rolo C , G. Rong , ,Ch. Rosner , M. Rump , A. Sarantsev ,e , M. Savrié B , Y. Schelhaas , K. Schoenning , W. Shan , X. Y. Shan , ,M. Shao , , C. P. Shen , P. X. Shen , X. Y. Shen , , H. Y. Sheng , X. Shi , , X. D Shi , , J. J. Song , Q. Q. Song , ,X. Y. Song , S. Sosio A, C , C. Sowa , S. Spataro A, C , F. F. Sui , G. X. Sun , J. F. Sun , L. Sun , S. S. Sun , ,X. H. Sun , Y. J. Sun , , Y. K Sun , , Y. Z. Sun , Z. J. Sun , , Z. T. Sun , Y. T Tan , , C. J. Tang , G. Y. Tang ,X. Tang , V. Thoren , B. Tsednee , I. Uman D , B. Wang , B. L. Wang , C. W. Wang , D. Y. Wang ,l , K. Wang , ,L. L. Wang , L. S. Wang , M. Wang , M. Z. Wang ,l , Meng Wang , , P. L. Wang , R. M. Wang , W. P. Wang , ,X. Wang ,l , X. F. Wang , X. L. Wang ,j , Y. Wang , Y. Wang , , Y. F. Wang , , , Y. Q. Wang , Z. Wang , ,Z. G. Wang , , Z. Y. Wang , Z. Y. Wang , Zongyuan Wang , , T. Weber , D. H. Wei , P. Weidenkaff , H. W. Wen ,S. P. Wen , U. Wiedner , G. Wilkinson , M. Wolke , L. H. Wu , L. J. Wu , , Z. Wu , , L. Xia , , Y. Xia , S. Y. Xiao ,Y. J. Xiao , , Z. J. Xiao , Y. G. Xie , , Y. H. Xie , T. Y. Xing , , X. A. Xiong , , Q. L. Xiu , , G. F. Xu , J. J. Xu ,L. Xu , Q. J. Xu , W. Xu , , X. P. Xu , F. Yan , L. Yan A, C , W. B. Yan , , W. C. Yan , Y. H. Yan , H. J. Yang ,h ,H. X. Yang , L. Yang , R. X. Yang , , S. L. Yang , , Y. H. Yang , Y. X. Yang , Yifan Yang , , Z. Q. Yang , M. Ye , ,M. H. Ye , J. H. Yin , Z. Y. You , B. X. Yu , , , C. X. Yu , J. S. Yu , T. Yu , C. Z. Yuan , , X. Q. Yuan ,l , Y. Yuan ,C. X. Yue , A. Yuncu B,a , A. A. Zafar , Y. Zeng , B. X. Zhang , B. Y. Zhang , , C. C. Zhang , D. H. Zhang ,H. H. Zhang , H. Y. Zhang , , J. Zhang , , J. L. Zhang , J. Q. Zhang , J. W. Zhang , , , J. Y. Zhang , J. Z. Zhang , ,K. Zhang , , L. Zhang , L. Zhang , S. F. Zhang , T. J. Zhang ,h , X. Y. Zhang , Y. Zhang , , Y. H. Zhang , ,Y. T. Zhang , , Yang Zhang , Yao Zhang , Yi Zhang ,j , Yu Zhang , Z. H. Zhang , Z. P. Zhang , Z. Y. Zhang , G. Zhao ,J. Zhao , J. W. Zhao , , J. Y. Zhao , , J. Z. Zhao , , Lei Zhao , , Ling Zhao , M. G. Zhao , Q. Zhao , S. J. Zhao ,T. C. Zhao , Y. B. Zhao , , Z. G. Zhao , , A. Zhemchugov ,b , B. Zheng , J. P. Zheng , , Y. Zheng ,l , Y. H. Zheng ,B. Zhong , L. Zhou , , L. P. Zhou , , Q. Zhou , , X. Zhou , X. K. Zhou , X. R. Zhou , , Xiaoyu Zhou , Xu Zhou ,A. N. Zhu , , J. Zhu , J. Zhu , K. Zhu , K. J. Zhu , , , S. H. Zhu , W. J. Zhu , X. L. Zhu , Y. C. Zhu , ,Y. S. Zhu , , Z. A. Zhu , , J. Zhuang , , B. S. Zou , J. H. Zou (BESIII Collaboration) a r X i v : . [ h e p - e x ] J un Institute of High Energy Physics, Beijing 100049, People’s Republic of China Beihang University, Beijing 100191, People’s Republic of China Beijing Institute of Petrochemical Technology, Beijing 102617, People’s Republic of China Bochum Ruhr-University, D-44780 Bochum, Germany Carnegie Mellon University, Pittsburgh, Pennsylvania 15213, USA Central China Normal University, Wuhan 430079, People’s Republic of China China Center of Advanced Science and Technology, Beijing 100190, People’s Republic of China COMSATS University Islamabad, Lahore Campus, Defence Road, Off Raiwind Road, 54000 Lahore, Pakistan Fudan University, Shanghai 200443, People’s Republic of China G.I. Budker Institute of Nuclear Physics SB RAS (BINP), Novosibirsk 630090, Russia GSI Helmholtzcentre for Heavy Ion Research GmbH, D-64291 Darmstadt, Germany Guangxi Normal University, Guilin 541004, People’s Republic of China Guangxi University, Nanning 530004, People’s Republic of China Hangzhou Normal University, Hangzhou 310036, People’s Republic of China Helmholtz Institute Mainz, Johann-Joachim-Becher-Weg 45, D-55099 Mainz, Germany Henan Normal University, Xinxiang 453007, People’s Republic of China Henan University of Science and Technology, Luoyang 471003, People’s Republic of China Huangshan College, Huangshan 245000, People’s Republic of China Hunan Normal University, Changsha 410081, People’s Republic of China Hunan University, Changsha 410082, People’s Republic of China Indian Institute of Technology Madras, Chennai 600036, India Indiana University, Bloomington, Indiana 47405, USA (A)INFN Laboratori Nazionali di Frascati, I-00044, Frascati, Italy; (B)INFN and University of Perugia, I-06100, Perugia,Italy (A)INFN Sezione di Ferrara, I-44122, Ferrara, Italy; (B)University of Ferrara, I-44122, Ferrara, Italy Institute of Physics and Technology, Peace Ave. 54B, Ulaanbaatar 13330, Mongolia Johannes Gutenberg University of Mainz, Johann-Joachim-Becher-Weg 45, D-55099 Mainz, Germany Joint Institute for Nuclear Research, 141980 Dubna, Moscow region, Russia Justus-Liebig-Universitaet Giessen, II. Physikalisches Institut, Heinrich-Buff-Ring 16, D-35392 Giessen, Germany KVI-CART, University of Groningen, NL-9747 AA Groningen, The Netherlands Lanzhou University, Lanzhou 730000, People’s Republic of China Liaoning Normal University, Dalian 116029, People’s Republic of China Liaoning University, Shenyang 110036, People’s Republic of China Nanjing Normal University, Nanjing 210023, People’s Republic of China Nanjing University, Nanjing 210093, People’s Republic of China Nankai University, Tianjin 300071, People’s Republic of China Peking University, Beijing 100871, People’s Republic of China Shandong Normal University, Jinan 250014, People’s Republic of China Shandong University, Jinan 250100, People’s Republic of China Shanghai Jiao Tong University, Shanghai 200240, People’s Republic of China Shanxi University, Taiyuan 030006, People’s Republic of China Sichuan University, Chengdu 610064, People’s Republic of China Soochow University, Suzhou 215006, People’s Republic of China Southeast University, Nanjing 211100, People’s Republic of China State Key Laboratory of Particle Detection and Electronics, Beijing 100049, Hefei 230026, People’s Republic of China Sun Yat-Sen University, Guangzhou 510275, People’s Republic of China Tsinghua University, Beijing 100084, People’s Republic of China (A)Ankara University, 06100 Tandogan, Ankara, Turkey; (B)Istanbul Bilgi University, 34060 Eyup, Istanbul, Turkey;(C)Uludag University, 16059 Bursa, Turkey; (D)Near East University, Nicosia, North Cyprus, Mersin 10, Turkey University of Chinese Academy of Sciences, Beijing 100049, People’s Republic of China University of Hawaii, Honolulu, Hawaii 96822, USA University of Jinan, Jinan 250022, People’s Republic of China University of Manchester, Oxford Road, Manchester, M13 9PL, United Kingdom University of Minnesota, Minneapolis, Minnesota 55455, USA University of Muenster, Wilhelm-Klemm-Str. 9, 48149 Muenster, Germany University of Oxford, Keble Rd, Oxford, UK OX13RH University of Science and Technology Liaoning, Anshan 114051, People’s Republic of China University of Science and Technology of China, Hefei 230026, People’s Republic of China University of South China, Hengyang 421001, People’s Republic of China University of the Punjab, Lahore-54590, Pakistan (A)University of Turin, I-10125, Turin, Italy; (B)University of Eastern Piedmont, I-15121, Alessandria, Italy; (C)INFN, I-10125, Turin, Italy Uppsala University, Box 516, SE-75120 Uppsala, Sweden Wuhan University, Wuhan 430072, People’s Republic of China Xinyang Normal University, Xinyang 464000, People’s Republic of China Zhejiang University, Hangzhou 310027, People’s Republic of China Zhengzhou University, Zhengzhou 450001, People’s Republic of China a Also at Bogazici University, 34342 Istanbul, Turkey b Also at the Moscow Institute of Physics and Technology, Moscow 141700, Russia c Also at the Functional Electronics Laboratory, Tomsk State University, Tomsk, 634050, Russia d Also at the Novosibirsk State University, Novosibirsk, 630090, Russia e Also at the NRC "Kurchatov Institute", PNPI, 188300, Gatchina, Russia f Also at Istanbul Arel University, 34295 Istanbul, Turkey g Also at Goethe University Frankfurt, 60323 Frankfurt am Main, Germany h Also at Key Laboratory for Particle Physics, Astrophysics and Cosmology, Ministry of Education; Shanghai Key Laboratoryfor Particle Physics and Cosmology; Institute of Nuclear and Particle Physics, Shanghai 200240, People’s Republic of China i Also at Government College Women University, Sialkot - 51310. Punjab, Pakistan. j Also at Key Laboratory of Nuclear Physics and Ion-beam Application (MOE) and Institute of Modern Physics, FudanUniversity, Shanghai 200443, People’s Republic of China k Also at Harvard University, Department of Physics, Cambridge, MA, 02138, USA l Also at State Key Laboratory of Nuclear Physics and Technology, Peking University, Beijing 100871, People’s Republic ofChina
Using a data sample of 2 .
93 fb − of e + e − collisions collected at √ s = 3 .
773 GeV in the BESIIIexperiment, we perform an analysis of the decay D → K S K + K − . The Dalitz plot is analyzedusing 1856 ±
45 flavor-tagged signal decays. We find that the Dalitz plot is well described by aset of six resonances: a (980) , a (980) + , φ (1020), a (1320) + , a (1320) − and a (1450) − . Theirmagnitudes, phases and fit fractions are determined as well as the coupling of a (980) to KK , g KK = (3 . ± . ± . D → K S K + K − ismeasured using 11 660 ±
118 untagged signal decays to be (4 . ± . ± . × − .Both measurements are limited by their systematic uncertainties. I. INTRODUCTION
The decay D → K S K + K − is a self-conjugate channel with a resonant substructure containing CP eigenstatesas well as non- CP eigenstates. An accurate measurementof the decay and its substructure has implications forvarious fields. The substructure of D → K S K + K − isdominated by the KK S-wave which can be studied in analmost background-free environment. In particular, lightscalar mesons are of interest since their spectrum is notfree of doubt [1, p. 658ff.]. The branching fractions of theresonant substructure and the total branching fractionare inputs to a better theoretical understanding of D - D mixing. Furthermore, the strong phase difference betweenthe decays D → K S K + K − and D → K S K + K − can bedetermined from the amplitude model. This phase is aninput to a measurement of the angle γ of the CKM uni-tarity triangle using the decay of B − → D K − with D → K S K + K − [2]. A model-independent determination ofthe strong phase will be presented in a separate paper [3]. Charge conjugation is implied throughout this work, except whereexplicitly noted otherwise. Where beneficial we abbreviate the final state K S K + K − by 3 K . The most recent analysis of D → K S K + K − was per-formed by the B A B AR experiment [4]. Using 12 500 flavor-tagged D decays, the total branching fraction was mea-sured relative to the decay D → K S π + π − . The currentvalue given by the Particle Data Group (PDG) [1] is de-rived from that measurement. Furthermore, a Dalitz plotanalysis was performed and it was found that the resonantsubstructure is well described by a set of four resonances: a (980) , a (980) + , φ (1020) and f (1370).In this work, the decay D → K S K + K − is analyzedusing a data sample of e + e − collisions corresponding toan integrated luminosity of 2 .
93 fb − [5] collected withthe BESIII detector at √ s = 3 .
773 GeV. At this energy,the produced ψ (3770) decays predominantly to D D and D + D − . The pair of neutral D mesons is produced ina quantum entangled state. The flavor of one mesoncan be inferred from the decay of the other meson if itis reconstructed in a flavor-specific decay channel. Thesubsample in which both D ’s are fully reconstructed isreferred to as the ‘tagged sample’. The sample in whichonly the reconstruction of the signal decay is required isdenoted as the ‘untagged sample’.The paper is structured as follows: we introduce the de-tector and the Monte-Carlo (MC) simulation in Section IIfollowed by the description of the event selection in Sec-tion III. The Dalitz plot analysis is presented in Section IVand the branching fraction measurement in Section V.Finally, we summarize our results in Section VI. II. DETECTOR AND DATA SETS
The BESIII detector records symmetric e + e − collisionswith high luminosity provided by the BEPCII storagering [6]. The center-of-mass energy ranges from 2 GeVto 4 . E/ d x .The momentum of a charged track with a transversemomentum of 1 GeV / c is measured with a resolutionof 0 . E/ d x resolution for electrons fromBhabha scattering is 6 %. A plastic scintillator time-of-flight (TOF) system provides a time resolution of 68 ps(110 ps) in the barrel (end cap) part and is used for particleidentification. The electromagnetic calorimeter (EMC)measures the energy of electromagnetic showers with aresolution of better than 2 . Geant4 [9] which includes the geometric descriptionof the BESIII detector and the detector response. Thesimulation includes the beam energy spread and initialstate radiation (ISR) in the e + e − annihilations. Theinclusive MC samples are simulated using KKMC [10]and consist of the production of DD pairs, non- DD decays of the ψ (3770), the ISR production of the J/ψ andΨ(3686) states, and continuum processes. The knowndecay modes are modelled with
EvtGen [11, 12] usingbranching fractions taken from the PDG [13], and theremaining unknown decays from the charmonium stateswith
LundCharm [14, 15]. The final state radiations(FSR) from charged final state particles are incorporatedwith the
Photos package [16]. All samples correspondto a luminosity of 5-10 times that of the data sample. Inaddition, the Dalitz plot analysis requires a large sample The current record is 1 . × cm − s − at √ s = 3 .
773 GeV of phase space distributed events of the signal channel.In this case, the decay ψ (3770) → D D with D → K S K + K − and D → (tag) is simulated. ‘Tag’ refers to anumber of flavor-specific channels that are used for flavortagging (see Table I). The measurement of the branchingfraction requires an accurate decay model for the signaldecay, and the result of the Dalitz plot analysis is usedto generate appropriate signal events which are used tosubstitute signal events in the inclusive MC sample. III. DATA PREPARATIONA. Event selection
Charged tracks are reconstructed from hits in the MDCand their momenta are determined from the track curva-ture. We require that each track has a point-of-closestapproach to the interaction point of 10 cm along the beamline ( V z ) and 1 cm perpendicular to it ( V r ). Furthermore,we require that the reconstructed polar angle is withinthe acceptance of the MDC of | cos θ | < .
93. The particlespecies is determined from d E/ d x measured by the MDCand TOF information. The combined χ ( H ) for a particlehypothesis H is given by χ ( H ) = χ E/ d x ( H ) + χ T OF ( H ) . (1)Using the corresponding number of degrees of freedom,a probability P H is calculated, and we require that allkaon and pion candidates satisfy P K > P π and P π > P K ,respectively.Two pions with opposite charge are combined to forma K S candidate. The requirement on V r is removed andthe requirement on V z is loosened to 20 cm for thesetracks. No particle identification requirement is imposed.The reconstructed K S invariant mass is denoted by m ks .The common vertex and the signed K S candidate flightdistance are estimated using a secondary vertex fit. The χ of the secondary vertex fit is required to be smaller than100. To suppress combinatorial background and decaysof the type K + K − π + π − , we require that K S candidateshave a ratio of flight distance over its uncertainty largerthan 2 for the untagged sample and larger than 0 for thetagged sample. The flavor tag final states include π ’sand η ’s which are reconstructed from their decays to apair of photons [17].We combine a K S candidate and two oppositely chargedkaon candidates to form a D signal candidate. Dueto the charm threshold decay kinematics, a convenientvariable to discriminate signal from background is thebeam-constrained mass m bc c = E − | ~p D | c . (2)The combined 3-momenta of the daughter tracks in therest frame of the ψ (3770) is denoted by ~p D and the beam TABLE I. Yields of the tagged sample for various tag channels.Yields include background candidates.Flavor tag Signal yield K − π + K − π + π K − π + π π K − π + π + π − K − π + π + π − π K − π + η energy by E beam . A kinematic fit with the nominal D mass as a constraint is applied, and candidates are re-quired to have a χ smaller than 20. In the untaggedsample, only the decay of one D meson to K S K + K − isreconstructed and the one with the smallest differencebetween the reconstructed energy of the candidate andhalf the center-of-mass energy of the e + e − beams is cho-sen. In the tagged sample, both D and D decays arereconstructed. One D meson is reconstructed using decaychannels specific to the flavor of the decaying meson (seeTable I). Usually multiple tag-signal candidate combina-tions are found, and we select the best combination ofone tag and one signal candidate via the average beam-constrained mass m bc closest to the nominal D mass.Finally, m bc and m ks are required to be within the axisboundaries of Fig. 1.Tagged and untagged samples contain 1935 and 13 209candidates, respectively. The contributions to the taggedsample from individual tag channels are listed in Table I.The signal yields are determined using a two-dimensionalfit to m bc and m ks . B. Signal and background
The analysis requires accurate determination of thesignal yields of tagged and untagged samples. The back-ground that passes our selection can be categorized ac-cording to its distribution in the plane of m ks versus m bc :• Non- K S background: The final state of the sig-nal decay is K + K − ( π + π − ) K S . Events that do notcontain the intermediate decay of a K S show up asa peak in m bc and a flat distribution in m ks .• Combinatorial background:
Events which donot contain the correct final state. These eventscome from qq production and from misreconstructed D decays. The m bc distribution has a phase spacecomponent and a wide peak component. The π + π − - - ) [[GeV]] E (D D E n t r i e s [ . ] Dtag_deltaE_data
Entries 2025Mean 0.005742 - Std Dev 0.02403 - tag decay mode E n t r i e s [ . ] Dtag_decayMode_data
Entries 2025Mean 1.77Std Dev 1.382 ] [GeV/c bc m ] [ G e V / c k s m E n t r i e s FIG. 1. Distribution of m ks versus m bc of the untagged sample. pair can, but does not have to originate from a K S decay. In case of an intermediate K S decay theevents show up as a band along m ks . Otherwisethey are broadly distributed.The distribution of m ks versus m bc of the untagged sam-ple is shown in Fig. 1. Since the signal peaks in thecenter of the m bc versus m ks plane, the signal and bothbackground components can be distinguished in data bya fit procedure.We use an unbinned extended two-dimensional maxi-mum likelihood fit to determine the signal yields. Sincethe variables m ks and m bc are almost independent thetwo-dimensional probability density function (PDF) isconstructed as a product of the one-dimensional PDFs.The signal components S bc ( m bc ) and S ks ( m ks ) are bothmodeled by a Crystal Ball function [18] with two-sidedpower law tails S ( x ) = ( n L | α L | ) n L e − | αL | ( n L | α L | − | α L | − x ) − n L e − x ( n R | α R | ) n R e − | αR | ( n R | α R | − | α R | + x ) − n R , (3)with x = ( m − µ ) /σ . The three function componentsare defined for the lower tail x < α L , the central part α L < x < α R and the higher tail x > α R . S bc and S ks have separate shape parameters µ, σ, n R , n L , α R and α L .The combinatorial background model for m bc consistsof an ARGUS phase space shape A ( m bc ) [19] and a Gaus-sian G c ( m bc ). The m ks shape is described by a polynomialof first order P (1) and a peaking component modeled bythe m ks shape of the signal B bcc ( m bc ) = f bcc A ( m bc ) + (1 − f bcc ) G c ( m bc ) (4) B ksc ( m ks ) = f ksc P (1) c ( m ks ) + (1 − f ksc ) S ks ( m ks ) . The non- K S background model has the same shape asthe signal in m bc and a polynomial of first order P (1) in m ks B bck ( m bc ) = S bc ( m bc ) (5) B ksk ( m ks ) = P (1) k ( m ks ) . The untagged sample additionally contains non- K S back-ground candidates which were reconstructed from tracksof the ’tag’ decay. We model these candidates by ad-ditional terms A ( m bc ) and S ks ( m ks ) in B bck ( m bc ) and B ksk ( m ks ), respectively.Shape parameters are determined using a simultaneousfit to MC samples representing signal and backgroundcomponents. The shape parameters, except the param-eter σ of the K S signal peak, are fixed afterwards. Thecomplete PDF is given by F ( m bc , m ks ) = N s S bc ( m bc ) S ks ( m ks )+ N c B bcc ( m bc ) B ksc ( m ks ) (6)+ N k B bck ( m bc ) B ksk ( m ks ) . The yields N s , N c and N k are determined by an extendedmaximum likelihood fit to the data sample. This proce-dure is used for the determination of the untagged signalyield for the branching fraction measurement. The taggedsignal fraction for the Dalitz plot analysis is determinedby a maximum likelihood fit with fixed sample yield. Thefit to the untagged sample is shown in Fig. 2 and yields11 660 ±
118 signal candidates. For the tagged sample werequire additionally to the description in Section III Athat all candidates are within a box in the m ks versus m bc plane. We choose a box of ± σ around the peak maximum in m bc and m ks , respectively.The signal yield determined by the fit within this boxis 1856 ±
45 tagged signal candidates with a purity of96 .
37 %.
IV. DALITZ PLOT ANALYSIS
The Dalitz plot of the decay D → K S K + K − afterreconstruction and selection is shown in Fig. 3. Trackmomenta are updated according to the kinematic fit,described in Section III A. The distribution is dominatedby the φ (1020) and the KK S -wave which in turn isusually described by the charged and neutral a (980)resonances. The a (980) mass is below the KK thresholdand therefore only its high-mass tail is visible. From thedistribution along the K S K + invariant mass the vectornature of the φ (1020) can be observed.In the following, the free parameters of the Dalitz am-plitude model are denoted with β and the Dalitz variableswith ξ . The latter can be either two invariant masses (asused in Fig. 3) or an invariant mass and the correspondinghelicity angle (see Eq. (24)). The Dalitz plot analysis isperformed using the ComPWA framework [20, 21]. ] [GeV/c BC m020040060080010001200140016001800 E v e n t s / ( . ) DataBackgroundModel ] [GeV/c bc m - ] s D e v i a t i o n [ ] [GeV/c KS m200400600800100012001400 E v e n t s / ( . ) ] [GeV/c ks m - ] s D e v i a t i o n [ FIG. 2. Projections of the untagged data sample and thefit model. Below the deviation between fit model and datasample is shown in units of its uncertainty.
A. Background
The Dalitz plot of background candidates is studiedusing an MC background sample as well as data and MCsideband samples. Sideband sample events are requiredto lie outside a box region of ± m bc and m ks , and within 1 . < m bc < . and 0 . < m ks < .
528 GeV/c . The comparison of theMC samples with the data sample shows a good agree-ment (Fig. 4), and the MC background sample is usedin the following to fix the shape of a phenomenologicalbackground model. A peaking component originates fromthe decay D → ( KK ) φ ( ππ ) ρ . The φ (1020) is generatedunpolarized; the effects of any spin alignment are expectedto be negligible. We describe it by a Breit-Wigner modelwith mass and width parameters of the φ (1020) and spinzero. An additional Breit-Wigner component with free ] /c [GeV - K + K m ] / c [ G e V + K S K m E n t r i e s FIG. 3. Dalitz plot of the decay D → K S K + K − from thetagged sample. The phase space boundary is indicated. ] /c [GeV - K + K m E n t r i e s / . Data sidebandMC backgroundBackground model
FIG. 4. Projection of the MC background sample (full dottedpoints) and the background model (line) on m KK , both scaledto the expected number of background candidates. The datasideband sample (open triangles) is scaled to approximate theMC background and is also re-binned due to low statistics.Note that zero entries are not plotted. parameters improves the fit quality close to the φ (1020)peak. Combinatorial background is described by a phasespace component. The background model shows goodagreement with sideband data, as shown in Fig. 4. Wefind χ / ndf = 116 /
99 for the comparison of the modeland sideband data.
B. Quantum entangled D D decays We analyze D mesons produced in the reaction e + e − → D D via a ψ (3770) as intermediate state. Incontrast to an isolated D decay, this has implications forthe decay rate since fundamental conservation laws holdfor the combined D D decay amplitude and not just forthe decay amplitude of one D . We mention especiallythe conservation of charge-parity ( CP ) which we assumeto be strictly conserved in the D system. We follow thephase convention CP (cid:12)(cid:12) D (cid:11) = − (cid:12)(cid:12) D (cid:11) . (7)The decay is mediated by the decay operator H . Inthe following, the transition amplitude (cid:10) j (cid:12)(cid:12) H (cid:12)(cid:12) D (cid:11) of anisolated D decay to the final state j is denoted by A j .From CP conservation it follows A j = (cid:10) j (cid:12)(cid:12) H (cid:12)(cid:12) D (cid:11) = − (cid:10) (cid:12)(cid:12) H (cid:12)(cid:12) D (cid:11) = −A A j = (cid:10) j (cid:12)(cid:12) H (cid:12)(cid:12) D (cid:11) = − (cid:10) (cid:12)(cid:12) H (cid:12)(cid:12) D (cid:11) = −A . (8)In the case that j is a CP eigenstate we have j = and we include the CP eigenvalue η of the final state j : A j = − η A j . We describe the amplitude ratio of D and D to the same final state j using its magnitude r j andphase δ j λ j = A j A j = − r j e − iδ j . (9)In general, the amplitudes A j depend on the phase spaceposition: A j is constant only for two-body decays. Wedenote those final states with r j ≤ j and theircharge-conjugates with r > .The combined wave function of D D is anti-symmetricdue to the negative parity of the e + e − reaction. Thematrix element of the decay of D and D to the finalstates i and j at decay times t and t , respectively, isgiven by M ij ( t , t ) = 1 √ (cid:20) (cid:10) i (cid:12)(cid:12) H (cid:12)(cid:12) D ( t ) (cid:11) (cid:10) j (cid:12)(cid:12) H (cid:12)(cid:12) D ( t ) (cid:11) − (cid:10) i (cid:12)(cid:12) H (cid:12)(cid:12) D ( t ) (cid:11) (cid:10) j (cid:12)(cid:12) H (cid:12)(cid:12) D ( t ) (cid:11) (cid:21) . (10)The BESIII experiment does not give access to the D decay time and therefore we are only interested in the time-integrated transition matrix element. The integration ofEq. (10) over the D decay time difference yields |M ij | = Z ∞−∞ |M ij ( | t − t | ) | d ( | t − t | ) ≈ (cid:12)(cid:12) A j A i − A j A i (cid:12)(cid:12) . (11) TABLE II. Hadronic parameters of tag channels. The mixingparameters and the hadronic parameters for the final state Kπ are obtained by a global fit [23]. The measurement for Kππ and K − π + π − π + are from [24]. For the other tag channels nomeasurements exist today.Parameter Value x (0 . ± .
14) % y (0 . ± .
07) % (cid:0) r KπD (cid:1) (0 . ± . δ Kπ (9 . ± . r Kππ D (4 . ± .
12) % δ Kππ (19 ± r K3 πD (5 . ± .
12) % δ K π ( − ± We choose the normalization such that |M ij | and |A i | are branching fractions when integrated over the phasespace. The D mixing parameters x and y are of O (10 − ) [1, p. 691ff] and thus are neglected in secondorder in the previous expression.Using Eq. (9) we can write Eq. (11) as |M ij | ≈ (cid:12)(cid:12) A j A i − A j A i (cid:12)(cid:12) = (cid:12)(cid:12) A i (cid:12)(cid:12) (cid:12)(cid:12) A j λ i − A j (cid:12)(cid:12) . (12)For our case we get |M tag , K | = (cid:12)(cid:12) A tag (cid:12)(cid:12) (cid:12)(cid:12) A K λ tag − A K (cid:12)(cid:12) . (13)The decay amplitudes of D → K S K + K − and D → K S K + K − are connected via [22] A K ( m K K + ,m K K − )= A K ( m K K − , m K K + ) (14)if CP is conserved. The amplitude of the flavor tag decay A tag does not depend on the phase space position of thesignal decay and is therefore constant in Eq. (13). Theratio of D to D amplitude of the tag decay is denotedby λ tag = − r D e − iδ D . (15)We use experimental input for the magnitude and phaseof λ tag . Since these parameters have not yet been mea-sured for each tag channel separately, we set them tocommon values for all tag channels. As nominal value wechoose the experimental average for the final state K − π + .Experimental results are summarized in Table II.A detailed derivation of the decay amplitude of quan-tum entangled D mesons is given in [1, 25]. To sum-marize, the result of an analysis of the D → K S K + K − Dalitz plot needs to be the amplitude model of an iso-lated D decay in order to be comparable with other D production reactions (e. g. D ∗ → D π ). Therefore, theeffect of the production mechanism needs to be consideredin the amplitude model. For many D final states thiseffect can be neglected, but it needs to be consideredin this analysis since K S K + K − is a self-conjugate finalstate. The effect of the quantum entanglement on themeasurement of the branching fraction is discussed inSection V A. C. Resonance model
The signal decay amplitude A K is parameterized in theisobar model using the helicity formalism. The dynamicparts are described by a Breit-Wigner formula and, inthe case of the KK S -wave, by a Flatté description. Weconsider an intermediate resonance R produced in theinitial state i and decaying to the final state f with finalstate particles a and b . The third (spectator) particlein the three-body decay is denoted by c . The resonancehas the angular momentum J and its parameterizationdepends on the center-of-mass energy squared s of thefinal state particles a and b .The Breit-Wigner description suggested by the PDG [1]is R J ( s ) = − g D → R g R → f m R − s + i √ s Γ( s ) , (16)with the mass dependent widthΓ( s ) = Γ R (cid:18) q ( s ) q ( m R ) (cid:19) J +1 (cid:18) m R √ s (cid:19) F J ( z ) F J ( z R ) . (17)The center-of-mass daughter momentum q ( s ) is imag-inary below threshold. Therefore, we derive it from ananalytic continuation of the phase-space factor: iρ = − ˆ ρπ log (cid:12)(cid:12)(cid:12) ρ − ˆ ρ (cid:12)(cid:12)(cid:12) , s < − ρπ arctan ρ , < s < s th − ˆ ρπ log (cid:12)(cid:12)(cid:12) ρ − ˆ ρ (cid:12)(cid:12)(cid:12) + i ˆ ρ, s th < s . (18)with ˆ ρ ( s ) = 116 π q | ˆ q ( s ) | √ s , (19)and ˆ q ( s ) = ( s − ( m a + m b ) )( s − ( m a − m b ) )4 s . (20)Eq. (19) is input to Eq. (18) and is in turn used to calculate q ( s ) from the resulting ρ . This is the parameterizationsuggested by the PDG [1, Section 48.2.3].The coupling constants for the production and decay g i in Eq. (16) are related to the partial width Γ i via g R → f = 1 q J ( s R ) F J ( z R ) s m R Γ R → f ρ ( s ) . (21)This relation holds for narrow and isolated resonances.The Breit-Wigner model assumes a point-like object.The effect of an extended resonance object is taken intoaccount via the angular momentum barrier factors F J ( z ).The Blatt-Weisskopf barrier factors [26] are widely used: F ( z ) = 1 F ( z ) = 2 zz + 1 (22) F ( z ) = 13 z ( z − + 9 z , where z ( s ) = q ( s ) R . Experience shows that the influ-ence of the resonance radius R is rather small. We use R = 1 . − .The KK S -wave involves the charged and neutral a (980). The a (980) couples strongly to the channel KK as well as to the channel ηπ . We describe both by acoupled channel formula, the so-called Flatté formula [27].We use the parameterization from Ref. [4]: R J ch ( s ) = − g D → R g KK m R − s + i ( g KK ρ KK + g ηπ ρ ηπ ) . (23)The coupling constants are denoted by g i and the phasespace factor ρ i is given in Eq. (18). The a (980) couplesalso to the channel K K . Consequently, we add a thirdchannel to the denominator of Eq. (23) with the samecoupling as the K + K − channel. Other than this difference,the a (980) and a (980) + use the same values for theircoupling constants.The decay D → K S K + K − is a decay of a pseudo-scalar particle to three final state pseudo-scalar particles.The angular distribution of an intermediate resonance istherefore given by the Legendre polynomials P J (cos θ R ),which depend on the helicity anglecos θ R = − m ac − m a − m c − E ∗ R E ∗ c q ∗ R q ∗ c . (24)Starred quantities are measured in the rest frame of theresonance.The total decay amplitude A j to a final state j is thecoherent sum over the individual resonances A j ( ξ , β ) = X J √ J +14 π × X i c i n i R Ji ( s i ) P J (cos θ i ) . (25) The prefactor originates from the normalization of theLegendre polynomials. The Dalitz plot variables are de-noted by ξ and consist of invariant mass s i and helicityangle θ i for each subsystem.The magnitude and phase are given by the complexcoefficient c i , and all resonances are normalized with thefactors n − i = Z d ξ R Ji ( ξ , β ) . (26)We insert A ( ξ , β ) into Eq. (13) to obtain the decay am-plitude including the effect of the quantum entanglementof D and D . D. Likelihood function
The probability function is given by L ( ξ , β ) = f · |M ( ξ , β ) | R |M ( ξ , β ) | (cid:15) ( ξ ) d ξ + (1 − f ) · | B ( ξ ) | . (27) B ( ξ ) denotes the Dalitz plot background model. The effi-ciency function is denoted by (cid:15) ( ξ ) and the signal purityby f = (96 . ± . (cid:15) ( ξ ). The likelihood function is evaluatedfor each event, and the logarithm of its products can bewritten as − log L ( β ) = − N X ev log L ( ξ ev , β ) , (28)where N is the size of the data sample. The interestingphysics parameters are the fit fractions which are definedas f i = | c i | R d ξ n i | R i ( ξ , β ) | R d ξ |M ( ξ , β ) | , (29)where c i is the magnitude of resonance R i . The integral R d ξ n i | R i ( ξ ) | is equal to one, due to our choice for theresonance normalization (Eq. (26)). A precise calculationof the statistical uncertainty of the fit fractions requiresthe propagation of the full covariance matrix through theintegration which is achieved using an MC approach. E. Model selection
The PDG [1] lists 12 resonances that could potentiallycontribute to the decay D → K S K + K − . Due to the0 TABLE III. Overview of resonances that could appear asintermediate states. Mass and width are the Breit-Wignerparameters. In case that these are channel dependent wequote the parameters of the KK final state. The parametersof f (980) and a (980) are weighted averages of previous mea-surements of Refs. [28–30] and Refs. [31–36], respectively. Theother values are from the PDG [13]. Resonance I G ( J PC ) Mass [MeV / c ] Coupling f (980) 0 + (0 ++ ) 971 ± g KK = (3 . ± .
05) GeV g ππ = (1 . ± .
1) GeV a (980) 1 − (0 ++ ) 994 +6 − g ηπ = (2 . ± .
04) GeV φ (1020) 0 − (1 −− ) 1019 . ± .
019 Γ = (4 . ± . f (1270) 0 + (2 ++ ) 1275 . ± . . +2 . − . )MeV a (1320) 1 − (2 ++ ) 1318 . ± . KK = (109 . ± .
4) MeV f (1370) [37] 0 + (0 ++ ) 1440 ± KK = (121 ±
15) MeV a (1450) 1 − (0 ++ ) 1474 ±
19 Γ = (256 ±
13) MeV limited size of the data sample we consider only resonancethat were found to decay to KK by previous experiments.An overview of known resonances is given in Table III.The choice of resonances that are included in the modelis a common problem in amplitude analysis. Generally,increasing the complexity of the model improves the fitquality. In the usual approach, a minimum statisticalsignificance or a minimum fit fraction (or both) is requiredfor each resonance. Those requirements are somewhatarbitrary parameters and furthermore, the order in whichresonances are added (or removed) from the model canlead to different sets of resonances. We apply a moreabstract method for resonance selection.Balancing a model between fit quality and model com-plexity is a common problem in the area of machine learn-ing and in statistics in general. One approach to solvesuch a problem is the so-called Least Absolute Shrinkageand Selection Operator (LASSO) method [38]. The basicidea is to penalize undesired behavior of the objectivefunction. In the original approach the objective functionis a least square model and the penalty function is thesum of the absolute values of the free parameters. Inthe context of particle physics this approach is describedin [39]. In our case the objective function is the logarithmof the likelihood, and the undesired behavior is a largesum over all fit fractions (which indicates strong interfer-ences). Therefore, we use the sum of the square-root offit fractions as penalty function. We modify Eq. (28) − log L λ P ( β ) = − log L ( β ) (30)+ λ P X i vuut R d ξ N i R i R ∗ i R d ξ P m,n N m N n R m R ∗ n . The square-root is used since it favors the suppression ofsmall contributions, in contrast to, for example, the sumof the fit fractions which would favor solutions with equalvalues. P l penalty scale - - - - - - - A k a i k e I n f o r m a t i o n C r i t er i o n = 9.00 minP l FIG. 5. Distribution of
AIC λ P c versus the penalty scale λ P .A minimum is found at λ minP = 9 .
0. For each value of λ P several fits with different starting parameters are performed(gray circles). Only the result with the lowest value (red dots)is considered. The parameter λ P regularizes the model complexity. Alarge value suppresses the sum of fit fractions and smallvalues allow for larger interference terms. It is a nui-sance parameter and we have to determine its optimalvalue. Again, this is a common problem in statisticsand a possible solution is the use of so-called informationcriteria. Information criteria are mathematical formula-tions of the ‘principle of parsimony’ [40]. This means inhypothesis testing that we prefer the model with fewerparameters over a more complicated model, given thesame goodness-of-fit. The criteria suggested by [39] arethe Akaike information criteria (AIC) [41] and Bayesianinformation criteria (BIC) [40]. We choose a slightly mod-ified version of the AIC which takes the size of the datasample into account:AIC λ P c = − L λ P + 2 r λ P + 2 r ( r + 1) N − r − . (31)The number of events in data is denoted by N and thecoefficient r is related to the complexity of the model.We follow the suggestion in Ref. [39] and use the numberof resonances as parameter r . We consider only reso-nances with a fit fraction larger than a minimum value.A minimum value of 2 × − gives a stable result.A scan for different values of λ P over a wide range isperformed to map out the minima of AIC λ P c . The scanis shown in Fig. 5 and a minimum at λ minP = 9 . λ P . After a set of resonances is selected the penaltyterm is removed from the likelihood.We obtain a model consisting of a (980) , a (980) + , φ (1020), a (1320) + , a (1320) − and a (1450) − . Themodel contains the a (1320) − and a (1450) − which areexpected to be doubly Cabibbo-suppressed with regard totheir positively charged partner. We suspect that thosecontributions appear as an artefact of an imperfect modeldescription in some phase space regions. We thereforedecide to quote additionally a ‘basic’ model of a (980) , a (980) + and φ (1020).1 F. Goodness-of-fit
The distribution of data events across the Dalitz plotis not uniform, and in a larger area of the phase spacealmost no events are observed (see Fig. 3). Therefore,we apply a goodness-of-fit test which is more suitable forthis situation than the widely used χ test. We choose apoint-to-point dissimilarity method [42] which providesan unbinned goodness-of-fit test. The test variable Φ isdefined asΦ = 12 Z d ξ Z d ξ [ ρ m ( ξ ) − ρ n ( ξ )] (cid:2) ρ m ( ξ ) − ρ n ( ξ ) (cid:3) R ( (cid:12)(cid:12) ξ − ξ (cid:12)(cid:12) ) . (32)We use the Euclidean metric to calculate the distance be-tween two points in phase space. For the general distancefunction R ( (cid:12)(cid:12) ξ − ξ (cid:12)(cid:12) ) we choose a Gaussian function with awidth which is proportional to the amplitude value. Theunderlying (in general unknown) PDFs of two samplesare denoted ρ m and ρ n . Therefore, we estimate ρ m and ρ n by MC integration using samples from each PDFΦ = 1 N ( N + 1) N X j>i R ( | n i − n j | ) − N M
M,N X j,i R ( | n i − m i | ) (33)+ 1 M ( M + 1) M X j>i R ( | m i − m j | ) . Elements of both samples are denoted by m i and n i andthe total sample size by M and N , respectively. Weidentify one sample with our data sample, and the secondone is generated using the final amplitude model. Tocalculate a probability that a certain model fits the data,the distribution of the test variable is needed. This cannot be analytically derived and we simulate it using anMC approach; the result is given below. G. Systematics
Systematic uncertainties on the Dalitz plot amplitudemodel arise from various sources: background description,amplitude model, inaccuracies of the MC simulation, ex-ternal parameters and the fit procedure. For each sourceof uncertainty, we rerun the fit with a different configura-tion and add the deviations from the nominal amplitudemodel in quadrature. An overview of the systematicuncertainties is given in Table IV.
1. Background
Uncertainties from the background treatment comefrom the background model as well as from the uncertaintyon the signal purity. The fit quality of the backgroundmodel is good as illustrated in Fig. 4, and we do not assignan uncertainty due to our choice of the model but weuse different samples to determine the shape parameters.The nominal sample is the MC background sample andwe additionally test MC and data sideband samples. Thedifference of the fit result in comparison to the nominalmodel is taken as systematic uncertainty. Note that thecontribution to the φ (1020) peak is different for the signaland sideband region. Thus, the systematic uncertaintyis a conservative assumption. The effect on the Dalitzplot analysis of the uncertainty on the signal purity isestimated by varying the signal purity by two times itsstatistical uncertainty to larger and smaller values.
2. Amplitude model
A source of uncertainty of the amplitude model is theresonance radius that is used in the barrier factors. Wevary it in steps of 1 GeV − from 0 GeV − to 5 GeV − .Our nominal value is 1 . − .The quantum entanglement of D D is included in theDalitz amplitude model. We use external measurementsof the magnitude and phase of λ tag (Eq. (15)). Theexperimental averages for r D and δ D from the final state K − π + are used as nominal values. The influence on theresult is studied using the value for r D from K − π + π andtwice the K − π + nominal value. The phase δ D is set tozero, twice the K − π + nominal value and to the measuredvalue of K − π + π − π + .
3. Monte-Carlo simulation
The efficiency correction of the data sample is obtainedfrom MC simulation. Differences between data and MCsimulation in track reconstruction and particle identifi-cation can influence the result. Especially, regions withlow momentum K ± tracks are prone to inaccuracies. Wecorrect for these differences using momentum dependentcorrection factors obtained from hadronic D decays. Wetest the influence of the tracking correction by rerun-ning the fit without correction. The influence is foundto be negligible for the Dalitz plot analysis and thus nosystematic uncertainty is assigned.Another effect comes from different momentum reso-lutions in data and MC simulation. The φ (1020) hasa width that is of the same order as the mass resolu-tion. We study the influence of the mass resolution byrerunning the minimization with a free width parameter2 TABLE IV. Overview of uncertainties for the Dalitz plot amplitude model. We list the fit parameters and their statistical andsystematic uncertainties. The fit parameters and fit fractions are corrected for their fitting biases and are denoted by ‘correctedvalue’. Systematic uncertainties are given in units of the statistical uncertainty of the parameter ¯ σ (we use the average value ofthe asymmetric uncertainties). Parameter g KK a (980) a (980) + φ (1020) a (1320) + a (1320) − a (1450) − [GeV] FF [%] | c | φ [rad] FF [%] | c | φ [rad] FF [%] | c | φ [rad] FF [%] | c | φ [rad] FF [%] | c | φ [rad] FF [%]Fit value 3.80 93 0.59 2.96 33 0.71 1.67 47 0.13 − − − − − σ σ Background 0.25 0.71 0.49 0.18 0.53 0.49 0.17 0.29 0.38 0.34 0.27 0.11 0.15 0.17 0.22 0.53 0.21Amplitude model 0.16 0.27 0.26 0.17 0.33 0.18 0.07 0.09 0.15 0.04 0.10 0.04 0.08 0.12 0.14 0.28 0.14Quantum correlation 0.04 1.21 0.56 0.24 0.45 1.14 1.65 1.91 0.13 0.32 0.31 0.19 0.21 0.12 0.21 0.52 0.36External parameters 1.44 0.49 0.15 0.17 0.48 0.36 1.96 0.27 0.09 1.27 0.19 0.56 1.17 0.59 0.31 0.34 0.55Fitting procedure 0.07 0.17 0.20 0.05 0.13 0.24 0.05 0.23 0.12 0.11 0.08 0.10 0.17 0.09 0.29 0.30 0.41Sys. uncertainty 1.46 1.50 0.79 0.34 0.86 1.31 2.56 1.97 0.43 1.35 0.46 0.61 1.21 0.64 0.53 0.87 0.80 which approximates a resolution difference. The parame-ter changes from 4 .
266 MeV to (5 . ± .
3) MeV. We addthe deviation from the nominal model to the systematicuncertainty. We keep the parameter fixed in the nominalfit.
4. External parameters
External parameters are listed in Tables II and III. Weshift each parameter by its uncertainty to smaller andlarger values and rerun the minimization. The deviationfrom the nominal model is taken as systematic uncertainty.The influence of the a (980) coupling to ηπ is estimated byrerunning the minimization with both couplings as free pa-rameters. We obtain a value of g ηπ = (2 . ± .
16) GeV.
5. Fit procedure
We validate that the analysis routine is bias free andthat the fit routine provides a correct estimate of thestatistical uncertainty. We use our nominal fit result togenerate a signal MC sample. This sample passes detec-tor simulation and reconstruction as well as the eventselection procedure. Then, we add the expected amountof background from MC simulation and rerun the mini-mization procedure. We calculate the difference betweenthe parameters of the nominal model and the fit result inunits of the statistical uncertainty of the parameter. Theprocedure is repeated with 200 statistically independentsamples.We find that the error estimate is correct but smallbiases for some parameters are present, especially for the parameters of the a (1450) − . We correct each fitparameter for its bias and add half of the correction tothe systematic uncertainty.Furthermore, we check that no better minimum existsin the parameter space. We do so by rerunning theminimization with start values chosen randomly acrossthe whole parameter space. From 200 fits, no fit with avalid minimum exhibits a smaller negative logarithm ofthe likelihood value than the nominal fit. H. Results
We find that the Dalitz plot is well described by amodel with six resonances: a (980) , a (980) + , φ (1020), a (1320) + , a (1320) − and a (1450) − . The Dalitz plotprojections and the fit model are shown in Fig. 6 and thefit parameters are listed in Table V. The magnitude andphase of the a (980) are fixed to 1 and 0, respectively,as a reference.The projections of the model and the data sample showan excellent fit quality. The probability of the goodness-of-fit test of the model is (71 ±
3) %. The fit fractionsfor the interference terms are listed in Table VI. Thelargest interference is a destructive interference betweenthe neutral and charged a (980). The total fit fraction ofthe interference terms sums up to 106 %.Furthermore, we study a model with a reduced set ofresonances that only includes a (980) , a (980) + , and φ (1020). The results are listed in Table VII. The prob-ability of the goodness-of-fit test of the reduced modelis (68 ±
3) %. The nominal amplitude model is used be-low in the branching fraction measurement to obtain the3
TABLE V. Result from the Dalitz plot analysis. The first uncertainty is statistical followed by systematic uncertainty. Thecoupling constant a (980) → KK is determined to be g KK = (3 . ± . ± . a (1320) + , a (1320) − and a (1320) + the upper limits and the central values (CV) of the fit fractions are quoted as well as their combined significance.Final state Magnitude Phase [rad] Fit fraction [%] Sign.[ σ ] a (980) K S ± ± > a (980) + K − . +0 . − . ± .
09 2 . +0 . − . ± .
06 34 ± ± > φ (1020) K S . +0 . − . ± .
08 1 . ± . ± .
19 48 ± ± > a (1320) + K − . ± . ± . − . +0 . − . ± . < . . . . . . a (1320) − K + . ± . ± . − . ± . ± . < . . a (1450) − K + . +0 . − . ± .
04 0 . ± . ± . < . . ± D and D amplitudes are omitted. a (980) + φ (1020) a (1320) + a (1320) − a (1450) − a (980) − .
83 0 . − .
40 0 . − . a (980) + − .
93 0 . − .
62 8 . φ (1020) − .
59 0 .
10 1 . a (1320) + − . − . a (1320) − − . TABLE VII. Result from the Dalitz plot analysis using amodel with resonant contributions from a (980) , a (980) + and φ (1020). The first uncertainty is statistical followed bysystematic uncertainty. The coupling constant a (980) → KK is determined to be g KK = (3 . ± . ± . Final state Magnitude Phase [rad] Fit fraction [%] a (980) K S ± ± a (980) + K − . ± . ± . − . +0 . − . ± .
19 36 ± ± φ (1020) K S . ± . ± .
04 1 . ± . ± .
20 48 ± ± ± signal efficiency. V. BRANCHING FRACTION MEASUREMENT
The branching fraction of D → K S K + K − is measuredusing the untagged sample. The branching fraction isgiven by B K = N K N D D · f QC · (cid:15) K · B K S → π + π − . (34)Here, the signal yield is denoted by N K , which iscorrected for the efficiency of reconstruction and se- lection (cid:15) K . We correct for the branching fraction ofthe K S reconstruction mode using B ( K S → π + π − ) =(69 . ± .
05) % [13]. Our branching fraction result isnormalized to the number of D D decays in the datasample [17]: N D D = (10 597 ± ± × . (35)The branching fraction in an untagged D D sample islinked to the branching fraction of an isolated D decayvia the correction factor f QC which is derived in thefollowing section. A. Quantum entanglement
We consider a pair of D mesons, of which one mesondecays to the signal final state and the other to an arbi-trary final state. In the following, i and j are differentfinal states. The branching fraction is given by B jX = |M jX | = X i (cid:16) |M ji | + |M jı | (cid:17) = X i ( B ji + B jı ) . (36)We sum over all possible final states of one D meson.As mentioned before we use a normalization in whichthe phase space integral over the norm of an amplitudecorresponds to a branching fraction. Using Eq. (9) wefind that B jX = X i B j B i (cid:20) h r i i + h r j i + h r i i h r j i− h √ r i cos δ i i (cid:10) √ r j cos δ j (cid:11) (cid:21) . (37)Here, h·i denotes phase space averaged values. The branch-ing fractions of isolated D decay sum up to one X i (cid:0) B i + B i (cid:1) = X i (cid:0) B i + B i h r i i (cid:1) = 1 (38)4 ] /c [GeV - K + K m E n t r i e s / . /ndf = 71.52/99 c ] /c [GeV - K + K m - ] s D e v [ ] /c [GeV + K S0 K m E n t r i e s / . /ndf = 92.15/99 c DataBkgModel ] /c [GeV + K S0 K m - ] s D e v [ ] /c [GeV - K S0 K m E n t r i e s / . /ndf = 93.97/99 c ] /c [GeV - K S0 K m - ] s D e v [ FIG. 6. Dalitz plot projections of data sample (full dots)and amplitude model (blue line). Below each projection thedeviation between model and data sample is shown in unitsof its uncertainty. The inset in the first plot shows a zoom onthe
KK S -wave contribution. and the mixing parameter y can be expressed as [25] y = 2 X i B i h√ r i cos δ i i . (39)Thus, Eq. (37) gives B jX = B j (cid:2) h r j i − (cid:10) √ r j cos δ j (cid:11) y (cid:3) . (40)The correction factor that links the branching fractionof a quantum entangled D D pair B jX to the branchingfraction of an isolated D decay B j is then given by2 f QC = 1 + h r j i − y (cid:10) √ r j cos δ j (cid:11) . (41) The quantities r j and δ j depend on the phase space po-sition, and we use the phase space averaged values forthe calculation of f QC . From the Dalitz amplitude modela value of f QC = 1 . ± .
015 is obtained for the finalstate j = K S K + K − . The statistical uncertainty of theDalitz amplitude model is propagated to f QC via an MCapproach. The limited statistics of the tagged samplecause a rather large uncertainty on f QC of 1 .
45 %.
B. Systematic uncertainties
Systematic uncertainties on the branching fraction mea-surement arise from several sources. An overview is givenin Table VIII.Deviations between data and MC simulation can leadto different resolutions in specific variables, thus leadingto different efficiencies for the selection criteria. Mostselection variables are already included in the uncertaintyon track reconstruction (see below). For the remainingrequirement on the χ of the D vertex fit we find an un-certainty of 0 . K S mass resolution and therefore the K S width isa free parameter in the fit. The D mass resolution isconsistent between data and MC simulation.The branching fraction measurement requires the to-tal efficiency for reconstruction and selection which issensitive to the substructure of the decay. We use theDalitz plot model to generate signal events which we usefor efficiency determination. Since the fit quality of theDalitz model is excellent we do not assign an additionaluncertainty.The signal yield is determined using models for sig-nal and background. The model shape is determinedusing MC simulation and discrepancies between data andsimulation can therefore lead to a bias in the yield deter-mination. We use the covariance matrix of the fit anda multi-dimensional Gaussian to generate sets of shapeparameters and recalculate signal and background yieldsusing these sets of parameters. We find that the system-atic uncertainty is less than 0 . K S reconstructionefficiency is studied using J/ψ → K ∗± K ∓ and J/ψ → φK S K ∓ π ± control samples [43]. We assign an uncertaintyof 1 . DD decays[44]. We assign 1 % uncertainty per charged kaontrack.The systematic uncertainty, and also the total uncer-tainty of the measurement is dominated by the contribu-tions due to track reconstruction and particle identifica-tion. In total the systematic uncertainty on the branchingfraction measurement is 3 . TABLE VIII. Overview of systematic uncertainties.Systematic uncertainties [%]Quantum entanglement 1 . . . E ffi c i e n c y K S reconstruction 1 . K ± tracking 2 . K ± particle identification 2 . . E x t . Number of D D decays 1 . B ( K S → π + π − ) 0 . . C. Result
The signal yield N K is determined by a two-dimensional fitting procedure. The projections to m bc and m ks of the data sample and the fit model are shownin Fig. 2. We obtain a signal yield of 11 660 ±
118 events.Using the inclusive MC sample we find an efficiency forreconstruction and selection of (cid:15) K = (17 . ± .
04) %where the uncertainty is due to limited MC statistics.According to Eq. (34) the branching fraction of D → K S K + K − is B ( D → K S K + K − ) = (42)(4 . ± . ± . × − . The relative statistical and systematic uncertainties are1 . .
67 %, respectively. The total uncertainty is3 .
81 %.
VI. CONCLUSION
In summary, we investigate the decay D → K S K + K − using 2 .
93 fb − of e + e − collisions collected at √ s =3 .
773 GeV recorded with the BESIII experiment. Weanalyse the Dalitz plot and measure its branching frac-tion.The K S K + K − Dalitz plot is described using an isobaramplitude model. We select the optimal set of resonancesusing a ‘penalty term’ method and find that the Dalitzplot is well described using an amplitude model withsix resonances of a (980) , a (980) + , φ (1020), a (1320) + , a (1320) − and a (1450) − . The largest contribution tothe total intensity comes from the a (980) that, togetherwith its charged partner, describes the KK threshold.Both resonances show a strong interference which leadsto a sum of fit fractions of the Dalitz amplitude model of(176 ±
20) %. The f (980) could appear as an intermediate resonance,but our strategy for resonance selection does not favora model that includes the f (980). With respect to thenominal model its significance is 1 . σ . The a (980)couples strongly to the channel KK as well as to thechannel πη . We measure its coupling to KK to be g KK = (3 . ± . ± . D mesons in each eventare reconstructed. Therefore the sample size is limitedand statistical and systematic uncertainties are of thesame order. The result is influenced by the quantumentanglement of D and D with respect to the measure-ments of isolated D decays. We include this effect inour amplitude models in order to quote parameters of anisolated D decay. The magnitude and phase of the ratioof D and D amplitudes of the tag decays are necessaryto describe this effect. Since those are not measured forall tag channels we use the parameters of D → K + π − for all channels. The effect of this substitution on thefinal result is included in the systematic uncertainties.The model includes the a (1320) − and the a (1450) − which are expected to be doubly Cabibbo-suppressed and,therefore, should have a significantly smaller fit fractionthan their positively charged partners. The fact that wedo not see this could be a hint that those contributionsare artefacts of an imperfect model description in partsof the phase space. Those states and the a (1320) + havea combined statistical significance of 5 . σ . Each of theirisospin partners a (1320) , a (1450) and a (1450) + havea statistical significance of 2 . σ or less with respect to thenominal model and are not included by our method forresonance selection. Because of the small fit fractions of a (1320) + , a (1320) − and a (1450) − we decide to reportupper limits. Due to these problems of the model, weadditionally quote a model built from the ‘visible’ resonantstates a (980) , a (980) + and φ (1020). The result is givenin Table VII. In comparison with the result from B A B AR [4]we use a different set of resonances which leads to strongerinterference terms.We measure the branching fraction of the decay D → K S K + K − to be (4 . ± . ± . × − .This is the first absolute measurement. We use the Dalitzamplitude model to accurately describe the signal decayin simulation and also to obtain the quantum entangle-ment correction factor. The measurement is in goodagreement with previous measurements and we are ableto reduce the uncertainty significantly. The measurementis systematically limited. VII. ACKNOWLEDGMENTS
The BESIII collaboration thanks the staff of BEPCIIand the IHEP computing center for their strong sup-6port. This work is supported in part by National KeyBasic Research Program of China under Contract No.2015CB856700; National Natural Science Foundation ofChina (NSFC) under Contracts Nos. 11625523, 11635010,11735014, 11822506, 11835012; the Chinese Academy ofSciences (CAS) Large-Scale Scientific Facility Program;Joint Large-Scale Scientific Facility Funds of the NSFCand CAS under Contracts Nos. U1532257, U1532258,U1732263, U1832207; CAS Key Research Program ofFrontier Sciences under Contracts Nos. QYZDJ-SSW-SLH003, QYZDJ-SSW-SLH040; 100 Talents Program ofCAS; INPAC and Shanghai Key Laboratory for ParticlePhysics and Cosmology; ERC under Contract No. 758462;German Research Foundation DFG under Contracts Nos. Collaborative Research Center CRC 1044, FOR 2359;Istituto Nazionale di Fisica Nucleare, Italy; Koninkli-jke Nederlandse Akademie van Wetenschappen (KNAW)under Contract No. 530-4CDP03; Ministry of Develop-ment of Turkey under Contract No. DPT2006K-120470;National Science and Technology fund; STFC (UnitedKingdom); The Knut and Alice Wallenberg Foundation(Sweden) under Contract No. 2016.0157; The Royal So-ciety, UK under Contracts Nos. DH140054, DH160214;The Swedish Research Council; U. S. Department of En-ergy under Contracts Nos. DE-FG02-05ER41374, DE-SC-0010118, DE-SC-0012069; University of Groningen (RuG)and the Helmholtzzentrum fuer SchwerionenforschungGmbH (GSI), Darmstadt. [1] M. Tanabashi et al. (Particle Data Group), Phys. Rev.
D98 , 030001 (2018).[2] B. Aubert et al. (BaBar Collaboration), Phys. Rev.
D78 ,034023 (2008).[3] M. Ablikim et al. (BESIII Collaboration), “Improvedmodel-independent determination of the strong-phasedifference between D and D to K S,L K + K − decays,”(2020), to be submitted to PRD.[4] B. Aubert et al. (BaBar Collaboration), Phys. Rev. D72 ,052008 (2005).[5] M. Ablikim et al. (BESIII Collaboration), Chin. Phys.
C37 , 123001 (2013); Phys. Lett.
B753 , 629 (2016).[6] C. Yu et al. , in
Proceedings, 7th International Particle Ac-celerator Conference (IPAC 2016): Busan, Korea (2016).[7] M. Ablikim et al. (BESIII Collaboration), Chin. Phys.
C44 , 040001 (2020).[8] M. Ablikim et al. (BESIII Collaboration), Nucl. Instrum.Meth.
A614 , 345 (2010).[9] S. Agostinelli et al. (GEANT4), Nucl. Instrum. Meth.
A506 , 250 (2003).[10] S. Jadach, B. F. L. Ward, and Z. Was, Phys. Rev.
D63 ,113009 (2001), arXiv:hep-ph/0006359 [hep-ph].[11] D. J. Lange, Nucl. Instrum. Meth.
A462 , 152 (2001).[12] R.-G. Ping, Chin. Phys.
C32 , 599 (2008).[13] C. Patrignani et al. (Particle Data Group), Chin. Phys.
C40 , 100001 (2016).[14] J. C. Chen, G. S. Huang, X. R. Qi, D. H. Zhang, andY. S. Zhu, Phys. Rev.
D62 , 034003 (2000).[15] R.-L. Yang, R.-G. Ping, and H. Chen, Chin. Phys. Lett. , 061301 (2014).[16] E. Richter-Was, Phys. Lett. B303 , 163 (1993).[17] M. Ablikim et al. (BESIII Collaboration), Chin. Phys.
C42 , 083001 (2018).[18] J. E. Gaiser,
Charmonium Spectroscopy From RadiativeDecays of the
J/ψ and ψ , Ph.D. thesis, SLAC (1982).[19] H. Albrecht et al. (ARGUS), Phys. Lett. B229 , 304(1989).[20] M. Michel et al. , J. Phys. Conf. Ser. , 022025 (2014).[21] M. Fritsch, K. Goetzen, W. Gradl, M. Michel, K. Peters,S. Pflüger, A. Pitka, P. Weidenkaff, and J. Zhang, “Com-pwa: The common partial wave analysis framework,” (2019).[22] J. Libby et al. (CLEO Collaboration), Phys. Rev.
D82 ,112006 (2010).[23] Y. S. Amhis et al. (HFLAV), (2019), arXiv:1909.12524[hep-ex].[24] T. Evans, S. Harnew, J. Libby, S. Malde, J. Rademacker,and G. Wilkinson, Phys. Lett.
B757 , 520 (2016), [Erra-tum: Phys. Lett.B765,402(2017)].[25] D. M. Asner and W. M. Sun, Phys. Rev.
D73 , 034024(2006), [Erratum: Phys. Rev.D77,019901(2008)].[26] J. M. Blatt and V. F. Weisskopf,
Theoretical NuclearPhysics (Springer, New York, 1979).[27] S. M. Flatte, Phys. Lett. , 224 (1976).[28] F. Ambrosino et al. (KLOE Collaboration), Eur. Phys. J.
C49 , 473 (2007).[29] R. Garcia-Martin, R. Kaminski, J. R. Pelaez, andJ. Ruiz de Elvira, Phys. Rev. Lett. , 072001 (2011).[30] M. Ablikim et al. (BES Collaboration), Phys. Lett.
B607 ,243 (2005).[31] S. Teige et al. (E852 Collaboration), Phys. Rev.
D59 ,012001 (1999).[32] A. Abele et al. , Phys. Rev.
D57 , 3860 (1998).[33] D. V. Bugg, Phys. Rev.
D78 , 074023 (2008).[34] F. Ambrosino et al. (KLOE Collaboration), Phys. Lett.
B681 , 5 (2009).[35] S. B. Athar et al. (CLEO Collaboration), Phys. Rev.
D75 ,032002 (2007).[36] G. S. Adams et al. (CLEO Collaboration), Phys. Rev.