BELLE2-CONF-PH-2022-007

Measurement of decay-time dependent CP violation in $B^{0}\rightarrow
K^{0}_{\scriptscriptstyle S}K^{0}_{\scriptscriptstyle
S}K^{0}_{\scriptscriptstyle S}$ using 2019–2021 Belle II data

The Belle II Collaboration

F. Abudinén, I. Adachi, K. Adamczyk, L. Aggarwal, P. Ahlburg, H. Ahmed, J. K.
Ahn, H. Aihara, N. Akopov, A. Aloisio, F. Ameli, L. Andricek, N. Anh Ky, D. M.
Asner, H. Atmacan, V. Aulchenko, T. Aushev, V. Aushev, T. Aziz, V. Babu, S.
Bacher, H. Bae, S. Baehr, S. Bahinipati, A. M. Bakich, P. Bambade, Sw.
Banerjee, S. Bansal, M. Barrett, G. Batignani, J. Baudot, M. Bauer, A. Baur,
A. Beaubien, A. Beaulieu, J. Becker, P. K. Behera, J. V. Bennett, E. Bernieri,
F. U. Bernlochner, V. Bertacchi, M. Bertemes, E. Bertholet, M. Bessner, S.
Bettarini, V. Bhardwaj, B. Bhuyan, F. Bianchi, T. Bilka, S. Bilokin, D.
Biswas, A. Bobrov, D. Bodrov, A. Bolz, A. Bondar, G. Bonvicini, A. Bozek, M.
Bračko, P. Branchini, N. Braun, R. A. Briere, T. E. Browder, D. N. Brown, A.
Budano, L. Burmistrov, S. Bussino, M. Campajola, L. Cao, G. Casarosa, C.
Cecchi, D. Červenkov, M.-C. Chang, P. Chang, R. Cheaib, P. Cheema, V.
Chekelian, C. Chen, Y. Q. Chen, Y. Q. Chen, Y.-T. Chen, B. G. Cheon, K.
Chilikin, K. Chirapatpimol, H.-E. Cho, K. Cho, S.-J. Cho, S.-K. Choi, S.
Choudhury, D. Cinabro, L. Corona, L. M. Cremaldi, S. Cunliffe, T. Czank, S.
Das, N. Dash, F. Dattola, E. De La Cruz-Burelo, S. A. De La Motte, G. de
Marino, G. De Nardo, M. De Nuccio, G. De Pietro, R. de Sangro, B. Deschamps,
M. Destefanis, S. Dey, A. De Yta-Hernandez, R. Dhamija, A. Di Canto, F. Di
Capua, S. Di Carlo, J. Dingfelder, Z. Doležal, I. Domínguez Jiménez, T. V.
Dong, M. Dorigo, K. Dort, D. Dossett, S. Dreyer, S. Dubey, S. Duell, G.
Dujany, P. Ecker, S. Eidelman, M. Eliachevitch, D. Epifanov, P. Feichtinger,
T. Ferber, D. Ferlewicz, T. Fillinger, C. Finck, G. Finocchiaro, P. Fischer,
K. Flood, A. Fodor, F. Forti, A. Frey, M. Friedl, B. G. Fulsom, M. Gabriel, A.
Gabrielli, N. Gabyshev, E. Ganiev, M. Garcia-Hernandez, R. Garg, A. Garmash,
V. Gaur, A. Gaz, U. Gebauer, A. Gellrich, J. Gemmler, T. Geßler, G.
Ghevondyan, G. Giakoustidis, R. Giordano, A. Giri, A. Glazov, B. Gobbo, R.
Godang, P. Goldenzweig, B. Golob, P. Gomis, G. Gong, P. Grace, W. Gradl, S.
Granderath, E. Graziani, D. Greenwald, T. Gu, Y. Guan, K. Gudkova, J.
Guilliams, C. Hadjivasiliou, S. Halder, K. Hara, T. Hara, O. Hartbrich, K.
Hayasaka, H. Hayashii, S. Hazra, C. Hearty, M. T. Hedges, I. Heredia de la
Cruz, M. Hernández Villanueva, A. Hershenhorn, T. Higuchi, E. C. Hill, H.
Hirata, M. Hoek, M. Hohmann, S. Hollitt, T. Hotta, C.-L. Hsu, K. Huang, T.
Humair, T. Iijima, K. Inami, G. Inguglia, N. Ipsita, J. Irakkathil Jabbar, A.
Ishikawa, S. Ito, R. Itoh, M. Iwasaki, Y. Iwasaki, S. Iwata, P. Jackson, W. W.
Jacobs, D. E. Jaffe, E.-J. Jang, M. Jeandron, H. B. Jeon, Q. P. Ji, S. Jia, Y.
Jin, C. Joo, K. K. Joo, H. Junkerkalefeld, I. Kadenko, J. Kahn, H. Kakuno, M.
Kaleta, A. B. Kaliyar, J. Kandra, K. H. Kang, S. Kang, P. Kapusta, R. Karl, G.
Karyan, Y. Kato, H. Kawai, T. Kawasaki, C. Ketter, H. Kichimi, C. Kiesling,
C.-H. Kim, D. Y. Kim, H. J. Kim, K.-H. Kim, K. Kim, S.-H. Kim, Y.-K. Kim, Y.
Kim, T. D. Kimmel, H. Kindo, K. Kinoshita, C. Kleinwort, B. Knysh, P. Kodyš,
T. Koga, S. Kohani, K. Kojima, I. Komarov, T. Konno, A. Korobov, S. Korpar, N.
Kovalchuk, E. Kovalenko, R. Kowalewski, T. M. G. Kraetzschmar, F. Krinner, P.
Križan, R. Kroeger, J. F. Krohn, P. Krokovny, H. Krüger, W. Kuehn, T. Kuhr, J.
Kumar, M. Kumar, R. Kumar, K. Kumara, T. Kumita, T. Kunigo, M. Künzel, S.
Kurz, A. Kuzmin, P. Kvasnička, Y.-J. Kwon, S. Lacaprara, Y.-T. Lai, C. La
Licata, K. Lalwani, T. Lam, L. Lanceri, J. S. Lange, M. Laurenza, K.
Lautenbach, P. J. Laycock, R. Leboucher, F. R. Le Diberder, I.-S. Lee, S. C.
Lee, P. Leitl, D. Levit, P. M. Lewis, C. Li, L. K. Li, S. X. Li, Y. B. Li, J.
Libby, K. Lieret, J. Lin, Z. Liptak, Q. Y. Liu, Z. A. Liu, D. Liventsev, S.
Longo, A. Loos, A. Lozar, P. Lu, T. Lueck, F. Luetticke, T. Luo, C. Lyu, C.
MacQueen, M. Maggiora, R. Maiti, S. Maity, R. Manfredi, E. Manoni, A. Manthei,
S. Marcello, C. Marinas, L. Martel, A. Martini, L. Massaccesi, M. Masuda, T.
Matsuda, K. Matsuoka, D. Matvienko, J. A. McKenna, J. McNeil, F. Meggendorfer,
F. Meier, M. Merola, F. Metzner, M. Milesi, C. Miller, K. Miyabayashi, H.
Miyake, H. Miyata, R. Mizuk, K. Azmi, G. B. Mohanty, N. Molina-Gonzalez, S.
Moneta, H. Moon, T. Moon, J. A. Mora Grimaldo, T. Morii, H.-G. Moser, M.
Mrvar, F. J. Müller, Th. Muller, G. Muroyama, C. Murphy, R. Mussa, I.
Nakamura, K. R. Nakamura, E. Nakano, M. Nakao, H. Nakayama, H. Nakazawa, A.
Narimani Charan, M. Naruki, Z. Natkaniec, A. Natochii, L. Nayak, M. Nayak, G.
Nazaryan, D. Neverov, C. Niebuhr, M. Niiyama, J. Ninkovic, N. K. Nisar, S.
Nishida, K. Nishimura, M. H. A. Nouxman, K. Ogawa, S. Ogawa, S. L. Olsen, Y.
Onishchuk, H. Ono, Y. Onuki, P. Oskin, F. Otani, E. R. Oxford, H. Ozaki, P.
Pakhlov, G. Pakhlova, A. Paladino, T. Pang, A. Panta, E. Paoloni, S. Pardi, K.
Parham, H. Park, S.-H. Park, B. Paschen, A. Passeri, A. Pathak, S. Patra, S.
Paul, T. K. Pedlar, I. Peruzzi, R. Peschke, R. Pestotnik, F. Pham, M. Piccolo,
L. E. Piilonen, G. Pinna Angioni, P. L. M. Podesta-Lerma, T. Podobnik, S.
Pokharel, L. Polat, V. Popov, C. Praz, S. Prell, E. Prencipe, M. T. Prim, M.
V. Purohit, H. Purwar, N. Rad, P. Rados, S. Raiz, A. Ramirez Morales, R.
Rasheed, N. Rauls, M. Reif, S. Reiter, M. Remnev, I. Ripp-Baudot, M. Ritter,
M. Ritzert, G. Rizzo, L. B. Rizzuto, S. H. Robertson, D. Rodríguez Pérez, J.
M. Roney, C. Rosenfeld, A. Rostomyan, N. Rout, M. Rozanska, G. Russo, D.
Sahoo, Y. Sakai, D. A. Sanders, S. Sandilya, A. Sangal, L. Santelj, P.
Sartori, Y. Sato, V. Savinov, B. Scavino, M. Schnepf, M. Schram, H. Schreeck,
J. Schueler, C. Schwanda, A. J. Schwartz, B. Schwenker, M. Schwickardi, Y.
Seino, A. Selce, K. Senyo, I. S. Seong, J. Serrano, M. E. Sevior, C. Sfienti,
V. Shebalin, C. P. Shen, H. Shibuya, T. Shillington, T. Shimasaki, J.-G. Shiu,
B. Shwartz, A. Sibidanov, F. Simon, J. B. Singh, S. Skambraks, J. Skorupa, K.
Smith, R. J. Sobie, A. Soffer, A. Sokolov, Y. Soloviev, E. Solovieva, S.
Spataro, B. Spruck, M. Starič, S. Stefkova, Z. S. Stottler, R. Stroili, J.
Strube, J. Stypula, Y. Sue, R. Sugiura, M. Sumihama, K. Sumisawa, T.
Sumiyoshi, W. Sutcliffe, S. Y. Suzuki, H. Svidras, M. Tabata, M. Takahashi, M.
Takizawa, U. Tamponi, S. Tanaka, K. Tanida, H. Tanigawa, N. Taniguchi, Y. Tao,
P. Taras, F. Tenchini, R. Tiwary, D. Tonelli, E. Torassa, N. Toutounji, K.
Trabelsi, I. Tsaklidis, T. Tsuboyama, N. Tsuzuki, M. Uchida, I. Ueda, S.
Uehara, Y. Uematsu, T. Ueno, T. Uglov, K. Unger, Y. Unno, K. Uno, S. Uno, P.
Urquijo, Y. Ushiroda, Y. V. Usov, S. E. Vahsen, R. van Tonder, G. S. Varner,
K. E. Varvell, A. Vinokurova, L. Vitale, V. Vobbilisetti, V. Vorobyev, A.
Vossen, B. Wach, E. Waheed, H. M. Wakeling, K. Wan, W. Wan Abdullah, B. Wang,
C. H. Wang, E. Wang, M.-Z. Wang, X. L. Wang, A. Warburton, M. Watanabe, S.
Watanuki, J. Webb, S. Wehle, M. Welsch, C. Wessel, J. Wiechczynski, P.
Wieduwilt, H. Windel, E. Won, L. J. Wu, X. P. Xu, B. D. Yabsley, S. Yamada, W.
Yan, S. B. Yang, H. Ye, J. Yelton, J. H. Yin, M. Yonenaga, Y. M. Yook, K.
Yoshihara, T. Yoshinobu, C. Z. Yuan, Y. Yusa, L. Zani, Y. Zhai, J. Z. Zhang,
Y. Zhang, Y. Zhang, Z. Zhang, V. Zhilich, J. Zhou, Q. D. Zhou, X. Y. Zhou, V.
I. Zhukova, V. Zhulanov, and R. Žlebčík

We report a measurement of decay-time dependent CP-violating parameters in
$B^{0}\rightarrow K^{0}_{\scriptscriptstyle S}K^{0}_{\scriptscriptstyle
S}K^{0}_{\scriptscriptstyle S}$ decays. We use $(198.0\pm 3.0)\times 10^{6}\
B\kern 1.79993pt\overline{\kern-1.79993ptB}{}$ pairs collected at the
$\Upsilon(4S)$ resonance with the Belle II detector at the SuperKEKB
asymmetric-energy $e^{+}e^{-}$ collider. The observed mixing-induced and
direct CP violation parameters are $\mathcal{S}=-1.86\
_{-0.46}^{+0.91}\leavevmode\nobreak\ {\rm(stat)}\pm 0.09\leavevmode\nobreak\
{\rm(syst)}$ and $\mathcal{A}=-0.22\ _{-0.27}^{+0.30}\leavevmode\nobreak\
{\rm(stat)}\pm 0.04\leavevmode\nobreak\ {\rm(syst)}$, respectively.

##  1 Introduction

In the Standard Model (SM), the charmless three-body decay $B^{0}\rightarrow
K^{0}_{\scriptscriptstyle S}K^{0}_{\scriptscriptstyle
S}K^{0}_{\scriptscriptstyle S}$ is mediated by the $b\rightarrow sq\mkern
1.5mu\overline{\mkern-1.5muq\mkern-0.1mu}\mkern 0.1mu$ quark transition
corresponding to a “penguin” one-loop amplitude represented by the diagram in
Fig. 1 (charge-conjugate decays are implied hereafter unless specified
otherwise). The three-$K^{0}_{\scriptscriptstyle S}$ final state is CP even.
The small branching fraction $(6.0\pm 0.5)\times 10^{-6}$[1] suppressed by the
penguin loop makes this decay sensitive to a possible contribution from non-SM
physics [2].

Decay-time dependent CP violation arises from interference between decay
amplitudes with and without mixing, due to an irreducible phase in the
Cabibbo-Kobayashi-Maskawa (CKM) quark-mixing matrix [3]. If one of the neutral
$B$ mesons produced from the $\Upsilon(4S)$ decays into a CP eigenstate,
$f_{{\it CP}}$, at time $t_{{\it CP}}$, and the other $B$ meson turns into a
flavor-distinguishable final state, $f_{\rm tag}$, at time $t_{\rm tag}$, the
time-dependent decay rate is given by [4, 5, 6]

$\displaystyle\mathcal{P}(\Delta t)=\frac{e^{-|\Delta
t|/\tau_{B^{0}}}}{4\tau_{B^{0}}}(1+q[\mathcal{S}\sin(\Delta m_{d}\Delta
t)+\mathcal{A}\cos(\Delta m_{d}\Delta t)]),$ (1)

where $\Delta t\equiv t_{{\it CP}}-t_{\rm tag}$, and the CP-violating
parameters $\mathcal{S}$ and $\mathcal{A}$ are related to mixing-induced and
direct CP violation, respectively. We refer to the $B$ meson decaying into
$f_{{\it CP}}$ as $B_{{\it CP}}$ and to the other $B$ meson as $B_{\rm tag}$.
The flavor $q$ is $+1$ ($-1$) for $B^{0}_{\rm tag}$ ($\overline{B}^{0}_{\rm
tag}$), $\tau_{B^{0}}$ is the $B^{0}$ lifetime, and $\Delta m_{d}$ is the mass
difference between the two mass eigenstates of the $B^{0}$-$\mkern
1.5mu\overline{\mkern-1.5muB\mkern-0.1mu}\mkern 0.1mu^{0}$ system. The SM
predicts that $\mathcal{S}\approx-\sin 2\phi_{1}$ and $\mathcal{A}=0$ in
$B^{0}\rightarrow K^{0}_{\scriptscriptstyle S}K^{0}_{\scriptscriptstyle
S}K^{0}_{\scriptscriptstyle S}$, where $\phi_{1}$ is defined in terms of the
CKM matrix elements as
$\phi_{1}\equiv\arg[-V_{cd}V_{cb}^{*}/V_{td}V_{tb}^{*}]$. The deviation of
$\mathcal{S}$ from $-\sin 2\phi_{1}$ is predicted to be 0.02 with an
uncertainty smaller than 0.01 [7]. The Belle[8] and BaBar [9] experiments
measured these asymmetries with comparable uncertainties (25% on $\mathcal{S}$
and 20% on $\mathcal{A}$), where the uncertainties are dominated by the size
of their $B^{0}\rightarrow K^{0}_{\scriptscriptstyle
S}K^{0}_{\scriptscriptstyle S}K^{0}_{\scriptscriptstyle S}$ samples. The world
average values are $\mathcal{S}=-0.83\pm 0.17$ and $\mathcal{A}=0.15\pm 0.12$
[10]. The search for non-SM physics with a new CP-violating phase requires
additional independent measurements with improved sensitivity.

We report a measurement of $\mathcal{S}$ and $\mathcal{A}$ in
$B^{0}\rightarrow K^{0}_{\scriptscriptstyle S}K^{0}_{\scriptscriptstyle
S}K^{0}_{\scriptscriptstyle S}$ decays using a sample of $(198.0\pm 3.0)\times
10^{6}$ $B\mkern 1.5mu\overline{\mkern-1.5muB\mkern-0.1mu}\mkern 0.1mu$ pairs
collected by the Belle II experiment. We reconstruct $B^{0}\rightarrow
K^{0}_{\scriptscriptstyle S}K^{0}_{\scriptscriptstyle
S}K^{0}_{\scriptscriptstyle S}$ decays with $K^{0}_{\scriptscriptstyle
S}\rightarrow\pi^{+}\pi^{-}$ as $B_{{\it CP}}$ and suppress background using
two boosted decision-tree (BDT) classifiers. We then measure $q$ using the
remaining charged particles in the event and $\Delta t$ from the difference
between the decay positions of $B_{{\it CP}}$ and $B_{\rm tag}$. Finally,
likelihood fits are performed to determine the signal yield and CP-violating
parameters, respectively. We use $B^{+}\rightarrow K^{0}_{\scriptscriptstyle
S}K^{0}_{\scriptscriptstyle S}K^{+}$ decays as a control channel to determine
fit models and to validate the fit procedure.

Figure 1: Leading decay amplitude for the $B^{0}\rightarrow
K^{0}_{\scriptscriptstyle S}K^{0}_{\scriptscriptstyle
S}K^{0}_{\scriptscriptstyle S}$ decays.

##  2 The Belle II detector and data sample

The Belle II experiment operates at the SuperKEKB asymmetric-energy
$e^{+}e^{-}$ collider [11], at KEK in Tsukuba, Japan. The Belle II detector
[12] is designed to reconstruct the final-state particles from $e^{+}e^{-}$
collisions. Detector subsystems are arranged cylindrically around the beam
interaction point. The innermost subsystem is the vertex detector (VXD), which
consists of two inner layers of silicon pixel detectors (PXD) and four outer
layers of double-sided silicon strip detectors. The second layer of the PXD is
currently partially installed and covers one sixth of the design acceptance.
Compared to Belle, the larger volume of the VXD is expected to provide more
acceptance for $K^{0}_{\scriptscriptstyle S}$ vertexing. The main tracking
device is the central drift chamber, located outside of the vertex detector.
Outside of the drift chamber, two types of Cherenkov ring-imaging detectors
provide hadron identification. The central volume is covered by a time-of-
propagation detector, which uses quartz bars as Cherenkov radiator, and the
forward endcap is covered by an aerogel radiator detector. Outside of this is
the electromagnetic calorimeter, which is made up of CsI(Tl) crystals. A 1.5 T
axial magnetic field is generated by a superconducting solenoid surrounding
the calorimeter. The outermost subsystem is the $K^{0}_{\scriptscriptstyle L}$
and muon detector, which consists of iron plates interspersed with resistive
plate chambers and scintillators. Muons and $K^{0}_{\scriptscriptstyle L}$’s
are detected through their interactions with the material. The iron plates
also serve as the return yoke for the magnetic field.

We use a data set collected at the $\Upsilon(4S)$ resonance in 2019–2021. The
data set corresponds to an integrated luminosity of 189.3$\mbox{\,fb}^{-1}$
and contains $(198.0\pm 3.0)\times 10^{6}\ B\kern
1.79993pt\overline{\kern-1.79993ptB}{}$ pairs. We use two types of Monte Carlo
(MC) simulation samples to optimize selection criteria and to train BDTs used
for event reconstruction. In one sample $B^{0}\kern
1.79993pt\overline{\kern-1.79993ptB}{}^{0}$ pairs are generated, with one of
the $B$ mesons decaying to the signal mode. The signal MC sample is also used
to determine models for the likelihood fit to the $\Delta t$ and $q$
distributions. The other sample includes $e^{+}e^{-}\rightarrow q\overline{q}\
(q=u,d,s,c),B^{0}\kern 1.79993pt\overline{\kern-1.79993ptB}{}^{0},\ {\rm and}\
B^{+}B^{-}$ events that proceed with realistic hadronization and decay
processes. We use the EvtGen [13] package to simulate hadron decays and KKMC
[14] with Pythia [15] for $q\overline{q}$. The detector response is simulated
by Geant4 [16]. Both data and MC samples are analyzed with the Belle II
analysis software framework [17].

##  3 Event reconstruction

The $\Upsilon(4S)$ is produced at the collision point with a Lorentz boost
(${\beta\gamma}$) of 0.287, and subsequently decays to $B$ and $\mkern
1.5mu\overline{\mkern-1.5muB\mkern-0.1mu}\mkern 0.1mu$ mesons, which are
nearly at rest in the center-of-mass (CM) frame. Therefore, the $B$ meson
pairs propagate in the laboratory nearly along the boost direction with a
known boost factor, which enables us to approximate the decay-time difference
between them as $\Delta t=(\ell_{{\it CP}}-\ell_{\rm tag})/\beta\gamma$. Here,
$\ell_{{\it CP}({\rm tag})}$ is the decay position of $B_{{\it CP}({\rm
tag})}$ which is projected onto the boost axis.

Pairs of oppositely charged particles with dipion mass between
457.6${\mathrm{\,Me\kern-1.00006ptV\\!/}c^{2}}$ and
537.6${\mathrm{\,Me\kern-1.00006ptV\\!/}c^{2}}$ are used to reconstruct
$K^{0}_{\scriptscriptstyle S}\rightarrow\pi^{+}\pi^{-}$ candidates. The
$K^{0}_{\scriptscriptstyle S}$ properties are obtained from a kinematic fit of
the $\pi^{+}$ and $\pi^{-}$ trajectories. To reduce the combinatorial
background from incorrectly reconstructed (“fake”) $K^{0}_{\scriptscriptstyle
S}$ candidates, a discriminant variable is formed from a BDT classifier with
22 input variables that include kinematic quantities, particle identification
information, and the number of hits in the VXD associated to the $\pi^{\pm}$
tracks, which are referred to as VXD hits. The most discriminating variables
are the angle between the directions of $K^{0}_{\scriptscriptstyle S}$
momentum and the decay position seen from the IP in the laboratory frame, and
the flight length of $K^{0}_{\scriptscriptstyle S}$ normalized by its
uncertainty. A selection based on this BDT discriminant accepts $95\%$ of
correctly reconstructed (“true”) $K^{0}_{\scriptscriptstyle S}$ mesons and
$0.59\%$ of fake $K^{0}_{\scriptscriptstyle S}$ mesons. Using a fit to the
dipion-mass distribution in data, we confirm that the signal efficiency is
consistent with MC. Although the selection slightly biases $\tau_{B^{0}}$, it
does not significantly affect $\mathcal{S}$ and $\mathcal{A}$. We take this
effect into account as a source of systematic uncertainty.

We reconstruct $B$ candidates by combining three $K^{0}_{\scriptscriptstyle
S}$ candidates. We select $B$ candidates using the invariant mass
$M(K^{0}_{\scriptscriptstyle S}K^{0}_{\scriptscriptstyle
S}K^{0}_{\scriptscriptstyle S})$ and beam-energy-constrained mass $M_{\rm
bc}\equiv\sqrt{E_{\rm beam}^{2}-|\vec{p}_{B}|^{2}}$, where $E_{\rm beam}$ and
$\vec{p}_{B}$ are the beam energy and the momentum of $B$ meson in the CM
frame. The difference between the beam energy ($\sqrt{s}/2$) and the
reconstructed $B$ energy in the CM frame, $\Delta E\equiv E_{\rm beam}-E_{B}$,
is not used because of its correlation with $M_{\rm bc}$. We retain the $B$
candidates satisfying $5.2<M_{\rm
bc}<5.29{\mathrm{\,Ge\kern-1.00006ptV\\!/}c^{2}}$ and
$5.08<M(K^{0}_{\scriptscriptstyle S}K^{0}_{\scriptscriptstyle
S}K^{0}_{\scriptscriptstyle S})<5.48{\mathrm{\,Ge\kern-1.00006ptV\\!/}c^{2}}$,
but exclude the candidates satisfying $5.265<M_{\rm
bc}<5.29{\mathrm{\,Ge\kern-1.00006ptV\\!/}c^{2}}$ and
$5.08<M(K^{0}_{\scriptscriptstyle S}K^{0}_{\scriptscriptstyle
S}K^{0}_{\scriptscriptstyle S})<5.181{\mathrm{\,Ge\kern-1.00006ptV\\!/}c^{2}}$
to avoid contamination by background due to $B^{0(+)}\rightarrow
K^{0}_{\scriptscriptstyle S}K^{0}_{\scriptscriptstyle S}K^{*0(+)}$.

The dominant source of background arises from continuum $e^{+}e^{-}\rightarrow
q\mkern 1.5mu\overline{\mkern-1.5muq\mkern-0.1mu}\mkern 0.1mu$ events. We
suppress the continuum background using another BDT classifier
$\mathcal{O}^{\prime}_{CS}$ with the following input variables related to
event topology: the cosine of the angle between the thrust axes of $B_{{\it
CP}}$ and $B_{\rm tag}$; the magnitude of the thrust of $B_{\rm tag}$; the sum
of the transverse momenta of the particles in the event; missing mass squared;
and modified Fox-Wolfram moments [18]. Here, the thrust axis of a $B$ meson is
a unit vector $\vec{t}$ that maximizes the thrust magnitude
$T\equiv(\sum_{i}\left|\vec{t}\cdot\vec{p}_{i}\right|)/(\sum_{i}\left|\vec{p}_{i}\right|)$,
where $\vec{p}_{i}$ is the momentum of the $B$ meson $i$th daughter. The
selection on $\mathcal{O}^{\prime}_{CS}$ rejects 49% of background and retains
98% of signal. The selection criteria for $K^{0}_{\scriptscriptstyle S}$
candidates are determined by maximizing a figure of merit, $N_{\rm
sig}/\sqrt{N_{\rm sig}+N_{\rm bkg}}$, where $N_{\rm sig}$ and $N_{\rm bkg}$
are the yields of simulated signal and background events, respectively,
satisfying $5.27<M_{\rm bc}<5.29{\mathrm{\,Ge\kern-1.00006ptV\\!/}c^{2}}$,
$5.18<M(K^{0}_{\scriptscriptstyle S}K^{0}_{\scriptscriptstyle
S}K^{0}_{\scriptscriptstyle S})<5.38{\mathrm{\,Ge\kern-1.00006ptV\\!/}c^{2}}$,
and $\mathcal{O}^{\prime}_{CS}>0.5$.

In addition to the non-resonant decay amplitude, quasi-two-body decays
$B^{0}\rightarrow X(\rightarrow K^{0}_{\scriptscriptstyle
S}K^{0}_{\scriptscriptstyle S})K^{0}_{\scriptscriptstyle S}$ due to
$b\rightarrow s$ and $b\rightarrow c$ transitions contribute to
$B^{0}\rightarrow K^{0}_{\scriptscriptstyle S}K^{0}_{\scriptscriptstyle
S}K^{0}_{\scriptscriptstyle S}$ decays. Since we regard the decays via
$b\rightarrow s$ as signal, we veto the $b\rightarrow c$ contribution to
measure the CP violation in a pure $b\rightarrow s$ process. We only expect a
significant $b\rightarrow c$ contribution from
$\chi_{c0}K^{0}_{\scriptscriptstyle S}$. The branching ratio of
$B^{0}\rightarrow\chi_{c0}(\rightarrow K^{0}_{\scriptscriptstyle
S}K^{0}_{\scriptscriptstyle S})K^{0}_{\scriptscriptstyle S}$ is around 5% of
the signal branching ratio. We reject $B^{0}$ candidates if the invariant mass
of any combination of two $K^{0}_{\scriptscriptstyle S}$ candidates is in the
range $3379<M(K^{0}_{\scriptscriptstyle S}K^{0}_{\scriptscriptstyle
S})<3447{\mathrm{\,Me\kern-1.00006ptV\\!/}c^{2}}$. This requirement rejects
90% of the background from $B^{0}\rightarrow\chi_{c0}K^{0}_{\scriptscriptstyle
S}$ and 7.5% of signal.

##  4 Measurement of $B$-meson flavor and decay-time difference

We use a BDT-based algorithm to identify the $B_{\rm tag}$ flavor [19]. It
uses 13 BDTs, each of which extracts a specific signature of $b\rightarrow
c\rightarrow s$ cascade decays from the particle identification and kinematic
variables of particles not belonging to $B_{{\it CP}}$. The outputs from the
BDTs are combined by a higher-level BDT to return the value of $q$ (defined
earlier) and a tagging quality variable $r$. The variable $r$ varies from zero
for no tagging information to one for unambiguous flavor assignment. The
probability density function (PDF) for signal events is represented as a
modified version of Eq. (1) by the probability to misidentify the flavor, $w$,
and its difference between $B^{0}$ and $\mkern
1.5mu\overline{\mkern-1.5muB\mkern-0.1mu}\mkern 0.1mu^{0}$, $\Delta w$,

$\displaystyle\mathcal{P}^{\rm TD}_{\rm sig}(\Delta t,q)$
$\displaystyle=\frac{e^{-|\Delta
t|/\tau_{B^{0}}}}{4\tau_{B^{0}}}\left(\mathstrut{}^{\mathstrut}_{\mathstrut}1-q\Delta
w+(1-2w)q[\mathcal{S}\sin(\Delta m_{d}\Delta t)+\mathcal{A}\cos(\Delta
m_{d}\Delta t)]\right).$ (2)

The events are classified into seven independent $r$ intervals. For each of
these intervals, $w$ and $\Delta w$ are determined using flavor-specific
$B$-meson decays with large branching fractions [19].

To measure $\Delta t$, we reconstruct the decay vertices of $B_{{\it CP}}$ and
$B_{\rm tag}$ using information about the beam interaction point (IP), which
is modeled by a three-dimensional Gaussian distribution. The $B_{{\it CP}}$
vertex position is reconstructed from $K^{0}_{\scriptscriptstyle S}$ daughter
tracks and the reconstructed $B_{{\it CP}}$ trajectory that originates from
the IP and points toward the reconstructed $B_{{\it CP}}$ momentum. Often
$K^{0}_{\scriptscriptstyle S}$ mesons decay outside of the VXD volume
resulting in less well-measured decay positions. Hence, the $B_{{\it CP}}$
vertex resolution largely depends on the number of $K^{0}_{\scriptscriptstyle
S}$ mesons that have associated VXD hits. In the signal MC, the fractions of
events where zero, one, two, and three $K^{0}_{\scriptscriptstyle S}$
candidates have VXD hits are 0.4%, 8.0%, 37.7%, and 54.0%, respectively. When
only one $K^{0}_{\scriptscriptstyle S}$ has VXD hits, the $B_{{\it CP}}$
trajectory helps to significantly improve the $B_{{\it CP}}$ vertex
resolution, reducing the average vertex position uncertainty from around
270${\,\mu\rm m}$ to 130${\,\mu\rm m}$.

We use the tracks that do not belong to $B_{{\it CP}}$ to reconstruct the
$B_{\rm tag}$ vertex, excluding the tracks without an associated PXD hit and
those that, combined in opposite charge pairs, yield the
$K^{0}_{\scriptscriptstyle S}$ mass. Similarly to the $B_{{\it CP}}$
vertexing, we reconstruct the $B_{\rm tag}$ trajectory using the IP
information and $B_{\rm tag}$ momentum, which is calculated as the difference
between the beam momentum and the $B_{{\it CP}}$ momentum [20]. The $B_{\rm
tag}$ trajectory is included in the vertex fit to improve the vertex
resolution and reconstruction efficiency. We use the $\chi^{2}$ per degree of
freedom of the vertex fit $\chi^{2}/N$ and the vertex position uncertainty
$\sigma_{\ell}$ as indicators of the quality of the $\Delta t$ measurement.
The number $N_{\rm tracks}$ of tracks in the fit, typically six, determines
$N$ as $N=2N_{\rm tracks}-1$. We require $\chi^{2}/N<100$ and
$\sigma_{\ell}<500{\,\mu\rm m}$ for $B_{\rm tag}$.

The multiplicity of $B_{{\it CP}}$ candidates in a selected event is 1.06 on
average. For events with multiple candidates, we choose the one with the
smallest $B_{{\it CP}}$ vertex fit $\chi^{2}$, which is not correlated with
the true $\Delta t$.

We apply the following selection criteria, related to the $B_{{\it CP}}$
vertex and $\Delta t$, to the remaining candidates: one or more
$K^{0}_{\scriptscriptstyle S}$ from $B_{{\it CP}}$ is associated with VXD
hits; if a $K^{0}_{\scriptscriptstyle S}$ daughter track from $B_{{\it CP}}$
is associated with a hit in the first PXD layer (layer-1 hit), its partner
from the same $K^{0}_{\scriptscriptstyle S}$ should also have a layer-1 hit;
otherwise it is likely the layer-1 hit belongs to a particle that was
incorrectly associated to $B_{{\it CP}}$, degrading the $B_{{\it CP}}$ vertex
resolution; the $\chi^{2}$ probability of the vertex fit should be larger than
0.001 for $B_{{\it CP}}$; $\sigma_{\ell}<500{\,\mu\rm m}$ for $B_{{\it CP}}$;
and $-30<\Delta t<30{\rm\,ps}$. We call the events passing these criteria
time-differential (TD) events and the others time-integrated (TI) events. We
use the $B_{\rm tag}$ flavor information for TI events, but not the $\Delta t$
information. Therefore, the PDF in Eq. (2) is integrated over $\Delta t$ for
TI events,

$\displaystyle\mathcal{P}_{\rm sig}^{\rm TI}(q)=\frac{1}{2}\left(1-q\Delta
w+(1-2w)q\mathcal{A}\frac{1}{1+\Delta m_{d}^{2}\tau_{B^{0}}^{2}}\right).$ (3)

##  5 Determination of signal yield

We extract the signal yields for TD and TI events separately from three-
dimensional likelihood fits to the unbinned distributions of $M_{\rm bc}$,
$M(K^{0}_{\scriptscriptstyle S}K^{0}_{\scriptscriptstyle
S}K^{0}_{\scriptscriptstyle S})$, and $\mathcal{O}^{\prime}_{CS}$. The
likelihood function accounts for two sample components, signal and background.
For the signal component, the $M_{\rm bc}$ distribution is modeled with a
Gaussian function, the $M(K^{0}_{\scriptscriptstyle
S}K^{0}_{\scriptscriptstyle S}K^{0}_{\scriptscriptstyle S})$ distribution with
the sum of two Gaussian functions, and the $\mathcal{O}^{\prime}_{CS}$
distribution with an asymmetric Gaussian function. We determine the parameters
for the signal shapes with fits to the signal MC sample. We use different
parameter sets for the $M(K^{0}_{\scriptscriptstyle
S}K^{0}_{\scriptscriptstyle S}K^{0}_{\scriptscriptstyle S})$ distribution in
TD and TI events because the latter have a broader distribution due to poorly
reconstructed candidates. For the background component, the $M_{\rm bc}$
distribution is modeled with an ARGUS function [21], the
$M(K^{0}_{\scriptscriptstyle S}K^{0}_{\scriptscriptstyle
S}K^{0}_{\scriptscriptstyle S})$ distribution with a linear function, and the
$\mathcal{O}^{\prime}_{CS}$ distribution with the sum of a Gaussian function
and an asymmetric Gaussian function.

We use $B^{+}\rightarrow K^{0}_{\scriptscriptstyle S}K^{0}_{\scriptscriptstyle
S}K^{+}$ decays to determine the background parameters because their kinematic
properties are similar to those of the signal decay. The background PDF shapes
are confirmed to be consistent between the two decay modes using MC samples.
The fit for $B^{+}\rightarrow K^{0}_{\scriptscriptstyle
S}K^{0}_{\scriptscriptstyle S}K^{+}$ gives both the corresponding yield, used
as a validation of our fitting procedures, and the background PDF parameters
for the $B^{0}\rightarrow K^{0}_{\scriptscriptstyle
S}K^{0}_{\scriptscriptstyle S}K^{0}_{\scriptscriptstyle S}$ fit.

Figure 2 shows the results of fits to $B^{0}\rightarrow
K^{0}_{\scriptscriptstyle S}K^{0}_{\scriptscriptstyle
S}K^{0}_{\scriptscriptstyle S}$ decays, separated into TD and TI samples. We
define the signal region as $5.271<M_{\rm
bc}<5.288{\mathrm{\,Ge\kern-1.00006ptV\\!/}c^{2}}$,
$5.181<M(K^{0}_{\scriptscriptstyle S}K^{0}_{\scriptscriptstyle
S}K^{0}_{\scriptscriptstyle
S})<5.366{\mathrm{\,Ge\kern-1.00006ptV\\!/}c^{2}}$, and
$-3.945<\mathcal{O}^{\prime}_{CS}<5.807$ so that each range retains 99.73% of
signal TD events. The signal yield is $53\pm 8$ events and the purity is 54%
in the signal region for TD events, and $48_{-7}^{+8}$ events and 45% for TI
events.

2

Figure 2: Distributions of (top) $M_{\rm bc}$, (middle)
$M(K^{0}_{\scriptscriptstyle S}K^{0}_{\scriptscriptstyle
S}K^{0}_{\scriptscriptstyle S})$, and (bottom) $\mathcal{O}^{\prime}_{CS}$ for
(left) TD candidates and (right) TI candidates with fit projections overlaid.
The $M_{\rm bc}$ distributions are restricted to events in the
$M(K^{0}_{\scriptscriptstyle S}K^{0}_{\scriptscriptstyle
S}K^{0}_{\scriptscriptstyle S})$ signal region. The
$M(K^{0}_{\scriptscriptstyle S}K^{0}_{\scriptscriptstyle
S}K^{0}_{\scriptscriptstyle S})$ and $\mathcal{O}^{\prime}_{CS}$ distributions
are restricted to events in the $M_{\rm bc}$ signal region.

##  6 Determination of CP-violating parameters

We determine the CP-violating parameters $\mathcal{S}$ and $\mathcal{A}$ by an
unbinned maximum-likelihood fit to the $\Delta t$ and $q$ distributions
combining TD and TI events restricted to the signal region. The contribution
to the likelihood function from the $j$th TD event is

$\displaystyle\mathcal{P}^{\rm TD}_{j}(\Delta t_{j},q_{j})=f^{\rm sig}_{j}\int
d(\Delta t^{\prime})R(\Delta t_{j}-\Delta t^{\prime})\mathcal{P}_{\rm
sig}(\Delta t^{\prime},q_{j})+(1-f^{\rm sig}_{j})\mathcal{P}_{\rm bkg}(\Delta
t_{j}),$ (4)

where $R(\delta\Delta t)$ is the response function of the $\Delta t$
measurement, $f_{j}^{\rm sig}$ is the signal fraction of the $j$th event, and
$\mathcal{P}_{\rm bkg}$ is the $\Delta t$ distribution of background events.
The response function consists of three components: detector resolutions for
$B_{{\it CP}}$ and $B_{\rm tag}$, bias due to secondary tracks from charmed
intermediate states for $B_{\rm tag}$, and a correction to the boost factors
due to their small, but non-zero, CM momentum. The parameters for the response
function are fixed to the values obtained from a fit to signal MC events. We
calculate $f_{j}^{\rm sig}$ from the signal-extraction fit including the
tagging probability $r$ and the cosine of the polar angle of the $B_{{\it
CP}}$ momentum in the CM frame. The additional variables are introduced to
avoid fit biases due to implicitly considering equal distributions that differ
across sample components [22]. The $r$ distribution for the background
component is obtained from the $M_{\rm
bc}<5.265{\mathrm{\,Ge\kern-1.00006ptV\\!/}c^{2}}$ sideband of the control
channel. The distribution $\mathcal{P}_{\rm bkg}$ is modeled in the same way
as the response function component for the detector resolution. Its parameters
are determined by a fit to the data sideband, $M_{\rm
bc}<5.265{\mathrm{\,Ge\kern-1.00006ptV\\!/}c^{2}}$. For TI events, we use the
likelihood function of Eq. (4) integrated over $\Delta t$,

$\displaystyle\mathcal{P}^{\rm TI}_{j}(q_{j})=f^{\rm sig}_{j}\mathcal{P}^{\rm
TI}_{\rm sig}(q_{j})+\frac{1-f^{\rm sig}_{j}}{2}.$ (5)

To validate the analysis procedure, we reconstruct $B^{+}\rightarrow
K^{0}_{\scriptscriptstyle S}K^{0}_{\scriptscriptstyle S}K^{+}$ decays without
using the position information of $K^{+}$ in the vertex fit and extract
$\mathcal{S}$ while fixing $\mathcal{A}$ at zero. The result
$\mathcal{S}=0.37_{-0.33}^{+0.31}$ is consistent with no CP violation and thus
supports the robustness of our analysis procedure. Only TD events are used in
the fit for the control channel.

Using $B^{0}\rightarrow K^{0}_{\scriptscriptstyle S}K^{0}_{\scriptscriptstyle
S}K^{0}_{\scriptscriptstyle S}$ decays, we obtain
$\mathcal{S}=-1.86_{-0.52}^{+0.60}$ and $\mathcal{A}=-0.22_{-0.21}^{+0.22}$,
where the uncertainties are statistical. The uncertainties are known to be
underestimated by the fit due to the small sample size. We reevaluate them
using a parametric bootstrap method, in which we generate simplified simulated
experiments obtained by sampling the likelihood, with the most probable
$\mathcal{S}$ and $\mathcal{A}$ within the physical region
$\mathcal{S}^{2}+\mathcal{A}^{2}\leq 1$ as input parameters. We obtain the
distribution of $\mathcal{S}$ and $\mathcal{A}$ from the simplified simulated
experiments and define the statistical uncertainty using 16 and 84 percentiles
of the distribution. The estimated uncertainties are ${}_{-0.46}^{+0.91}$ for
$\mathcal{S}$ and ${}_{-0.27}^{+0.30}$ for $\mathcal{A}$. Figure 3 shows the
signal component of the $\Delta t$ distribution separated for $q=\pm 1$ using
an ${}_{s}\mathcal{P}lot$ technique [23] and the asymmetry of the
distribution. The asymmetry is defined as $\frac{N_{+}(\Delta t)-N_{-}(\Delta
t)}{N_{+}(\Delta t)+N_{-}(\Delta t)}$, where $N_{\pm}(\Delta t)$ represent the
number of entries with $q=\pm 1$ in the corresponding $\Delta t$ bin. The
plots show only TD events.

Figure 3: (top) $\Delta t$ distribution of the signal component extracted with
the ${}_{s}\mathcal{P}lot$ technique and (bottom) its asymmetry for (left)
$B^{+}\rightarrow K^{0}_{\scriptscriptstyle S}K^{0}_{\scriptscriptstyle
S}K^{+}$ candidates and (right) $B^{0}\rightarrow K^{0}_{\scriptscriptstyle
S}K^{0}_{\scriptscriptstyle S}K^{0}_{\scriptscriptstyle S}$ candidates. Only
TD events are shown. In the $\Delta t$ distribution graphs, the red solid
curve and filled circles represent the fit result and data for $q=+1$, while
the blue dashed curve and open circles represent the fit result and data for
$q=-1$, respectively. The asymmetry is defined as $\frac{N_{+}(\Delta
t)-N_{-}(\Delta t)}{N_{+}(\Delta t)+N_{-}(\Delta t)}$, where $N_{\pm}(\Delta
t)$ represent the number of entries with $q=\pm 1$ in the corresponding
$\Delta t$ bin. In the asymmetry graphs, the points represent data and the
solid curve represents the result of the fit.

##  7 Systematic uncertainties

We consider various sources of systematic uncertainties and summarize them in
Table 1. To evaluate the systematic uncertainty on $\mathcal{S}$ and
$\mathcal{A}$ related to assumptions made on parameters of the fit model, we
repeat the fit on data using alternative values of the parameters randomly
sampled based on auxiliary knowledge. This approach is used for $w$ and
$\Delta w$ (referred to as flavor tagging in the table), the parameters
describing the resolution function, $\tau_{B^{0}}$ and $\Delta m_{d}$ (physics
parameters), the parameters for the $M_{\rm bc}$, $M(K^{0}_{\scriptscriptstyle
S}K^{0}_{\scriptscriptstyle S}K^{0}_{\scriptscriptstyle S})$, and
$\mathcal{O}^{\prime}_{CS}$ shapes (signal fraction), and the parameters for
the background $\Delta t$ shape. The widths of the resulting distributions of
$\mathcal{S}$ and $\mathcal{A}$ are taken as contributions to the systematic
uncertainty. We use the world-average values and uncertainties of
$\tau_{B^{0}}$ and $\Delta m_{d}$ [1] and modify the $\tau_{B^{0}}$
uncertainty considering the bias due to the $K^{0}_{\scriptscriptstyle S}$
selection.

Table 1: Systematic uncertainties Source |  $\delta\mathcal{S}$ | $\delta\mathcal{A}$
---|---|---
Vertex reconstruction | 0.025 | 0.022
Flavor tagging | 0.079 | 0.030
Resolution function | 0.012 | 0.006
Physics parameters | 0.008 | 0.000
Fit bias | 0.003 | 0.002
Signal fraction | 0.011 | 0.007
Background $\Delta t$ shape | 0.011 | 0.001
Detector misalignment | 0.002 | 0.004
Resolution model | 0.001 | 0.003
Tag-side interference | 0.014 | 0.015
Total | 0.087 | 0.042

The systematic uncertainty due to the vertex reconstruction is determined by
varying the parameters for the IP profile and boost vector, track requirements
for the $B_{\rm tag}$ vertex reconstruction, criteria to select TD events, and
correction of helix parameter uncertainties for vertexing. We use MC samples
simulated with a misaligned detector geometry to evaluate the misalignment
effect. The systematic uncertainty on the resolution model is determined by
analyzing a MC sample with alternative models of $\Delta t$ response
functions. Correlations are observed between $M(K^{0}_{\scriptscriptstyle
S}K^{0}_{\scriptscriptstyle S}K^{0}_{\scriptscriptstyle S})$ and $\chi^{2}/N$
for $B_{{\it CP}}$, and between $\mathcal{O}^{\prime}_{CS}$ and flavor tag.
For the systematic uncertainty due to the fit bias, two sets of simplified
simulated experiments are generated with and without these correlations; the
fits for $\mathcal{S}$ and $\mathcal{A}$ are performed ignoring these
correlations. We take the difference between the mean value of $\mathcal{S}$
and $\mathcal{A}$ for the two sets as a systematic uncertainty. For tag-side
interference [24], simplified simulated experiments are generated with and
without tag-side interference and the difference is taken as a systematic
uncertainty. The systematic uncertainty is dominated by that of the flavor
tagging performance owing to the limited size of the calibration sample.

##  8 Conclusion

In summary, we report a measurement of decay-time dependent CP violation in
$B^{0}\rightarrow K^{0}_{\scriptscriptstyle S}K^{0}_{\scriptscriptstyle
S}K^{0}_{\scriptscriptstyle S}$ decays using a data set corresponding to
$(198.0\pm 3.0)\times 10^{6}$ $B\kern 1.79993pt\overline{\kern-1.79993ptB}{}$
pairs collected with the Belle II experiment. The measured CP-violating
parameters are

$\displaystyle\mathcal{S}$ $\displaystyle=-1.86\
_{-0.46}^{+0.91}\leavevmode\nobreak\ {\rm(stat)}\pm 0.09\leavevmode\nobreak\
{\rm(syst)}\ {\rm and}$ (6) $\displaystyle\mathcal{A}$ $\displaystyle=-0.22\
_{-0.27}^{+0.30}\leavevmode\nobreak\ {\rm(stat)}\pm 0.04\leavevmode\nobreak\
{\rm(syst)}.$ (7)

Figure 4 shows the confidence regions based on likelihood-ratio ordering,
where both the statistical and systematic uncertainties are taken into account
[25]. Here, we constrain $\mathcal{S}$ and $\mathcal{A}$ within the physical
boundary, $\mathcal{S}^{2}+\mathcal{A}^{2}\leq 1$. The filled and open black
circles indicate the most probable values for $\mathcal{S}$ and $\mathcal{A}$
in the physical region and the SM prediction $(\mathcal{S},\mathcal{A})=(-\sin
2\phi_{1},0)$ based on measurements in $B^{0}\rightarrow(c\overline{c})K^{0}$
decays, respectively [10]. The results are consistent with the latest
measurements at Belle and BaBar.

Figure 4:  The red, orange, and green contours represent the 68.27%, 95.45%,
and 99.73% confidence regions for $\mathcal{S}$ and $\mathcal{A}$ given the
physical constraint $\mathcal{S}^{2}+\mathcal{A}^{2}\leq 1$. The filled and
open black circles indicate the most probable values for $\mathcal{S}$ and
$\mathcal{A}$ in the physical region and the SM prediction
$(\mathcal{S},\mathcal{A})=(-\sin 2\phi_{1},0)$ based on measurements in
$B^{0}\rightarrow(c\overline{c})K^{0}$ decays, respectively [10].

## Acknowledgements

We thank the SuperKEKB group for the excellent operation of the accelerator,
the KEK cryogenics group for the efficient operation of the solenoid, and the
KEK computer group for on-site computing support.

## References

  * [1] P. A. Zyla et al., (Particle Data Group), Review of Particle Physics, PTEP 2020 (2020) 083C01.
  * [2] Y. Grossman and M. P. Worah, CP asymmetries in B decays with new physics in decay amplitudes, Phys. Lett. B 395 (1997) 241, arXiv:hep-ph/9612269.
  * [3] M. Kobayashi and T. Maskawa, CP Violation in the Renormalizable Theory of Weak Interaction, Prog. Theor. Phys. 49 (1973) 652.
  * [4] A. B. Carter and A. I. Sanda, CP Violation in Cascade Decays of B Mesons, Phys. Rev. Lett. 45 (1980) 952.
  * [5] A. B. Carter and A. I. Sanda, CP Violation in B Meson Decays, Phys. Rev. D 23 (1981) 1567.
  * [6] I. I. Y. Bigi and A. I. Sanda, Notes on the Observability of CP Violations in B Decays, Nucl. Phys. B 193 (1981) 85.
  * [7] H. Y. Cheng, C. K. Chua, and A. Soni, CP-violating asymmetries in $B^{0}$ decays to $K^{+}K^{-}K^{0}_{S(L)}$ and $K^{0}_{S}K^{0}_{S}K^{0}_{S(L)}$, Phys. Rev. D 72 (2005) 094003, arXiv:hep-ph/0506268.
  * [8] K. H. Kang et al., (Belle collaboration), Measurement of time-dependent $CP$ violation parameters in $B^{0}\rightarrow K_{S}^{0}K_{S}^{0}K_{S}^{0}$ decays at Belle, Phys. Rev. D 103 (2021) 032003, arXiv:2011.00793 [hep-ex].
  * [9] J. P. Lees et al., (BaBar collaboration), Amplitude analysis and measurement of the time-dependent CP asymmetry of $B^{0}\rightarrow K_{S}^{0}K_{S}^{0}K_{S}^{0}$ decays, Phys. Rev. D 85 (2012) 054023, arXiv:1111.3636 [hep-ex].
  * [10] Y. S. Amhis et al., (HFLAV), Averages of $b$-hadron, $c$-hadron, and $\tau$-lepton properties as of 2018, Eur. Phys. J. C81 (2021) 226, arXiv:1909.12524 [hep-ex]. updated results and plots available at https://hflav.web.cern.ch/.
  * [11] K. Akai, K. Furukawa, and H. Koiso, (SuperKEKB), SuperKEKB Collider, Nucl. Instrum. Meth. A907 (2018) 188, arXiv:1809.01958 [physics.acc-ph].
  * [12] T. Abe et al., (Belle II Collaboration), Belle II Technical Design Report, arXiv:1011.0352 [physics.ins-det].
  * [13] D. J. Lange, The EvtGen particle decay simulation package, Nucl. Instrum. Meth. A462 (2001) 152.
  * [14] S. Jadach, B. F. L. Ward, and Z. Was, The Precision Monte Carlo event generator K K for two fermion final states in $e^{+}e^{-}$ collisions, Comput. Phys. Commun. 130 (2000) 260, arXiv:hep-ph/9912214 [hep-ph].
  * [15] T. Sjöstrand, S. Ask, J. R. Christiansen, R. Corke, N. Desai, P. Ilten, S. Mrenna, S. Prestel, C. O. Rasmussen, and P. Z. Skands, An Introduction to PYTHIA 8.2, Comput. Phys. Commun. 191 (2015) 159, arXiv:1410.3012 [hep-ph].
  * [16] S. Agostinelli et al., (GEANT4), GEANT4: A Simulation toolkit, Nucl. Instrum. Meth. A506 (2003) 250.
  * [17] T. Kuhr, C. Pulvermacher, M. Ritter, T. Hauth, and N. Braun, (Belle-II Framework Software Group), The Belle II Core Software, Comput. Softw. Big Sci. 3 (2019) 1, arXiv:1809.04299 [physics.comp-ph].
  * [18] S. H. Lee et al., (Belle collaboration), Evidence for $B^{0}\rightarrow\pi^{0}\pi^{0}$, Phys. Rev. Lett. 91 (2003) 261801, arXiv:hep-ex/0308040.
  * [19] F. Abudinén et al., (Belle II collaboration), B-flavor tagging at Belle II, Eur. Phys. J. C 82 (2022) 283, arXiv:2110.00790 [hep-ex].
  * [20] S. Dey and A. Soffer, Beam-Constrained Vertexing for $B$ Physics at the Belle II Experiment, Springer Proc. Phys. 248 (2020) 411.
  * [21] H. Albrecht et al., (ARGUS collaboration), Search for Hadronic $b\rightarrow u$ Decays, Phys. Lett. B 241 (1990) 278.
  * [22] G. Punzi, Comments on likelihood fits with variable resolution, eConf C030908 (2003) WELT002, arXiv:physics/0401045.
  * [23] M. Pivk and F. R. Le Diberder, SPlot: A Statistical tool to unfold data distributions, Nucl. Instrum. Meth. A 555 (2005) 356, arXiv:physics/0402083.
  * [24] O. Long, M. Baak, R. N. Cahn, and D. P. Kirkby, Impact of tag side interference on time dependent CP asymmetry measurements using coherent $B^{0}\kern 1.79993pt\overline{\kern-1.79993ptB}{}^{0}$ pairs, Phys. Rev. D 68 (2003) 034010, arXiv:hep-ex/0303030.
  * [25] G. J. Feldman and R. D. Cousins, A Unified approach to the classical statistical analysis of small signals, Phys. Rev. D 57 (1998) 3873, arXiv:physics/9711021.