The Belle Collaboration

#
Measurement of time-dependent $C\\!P$ violation parameters in $B^{0}\to
K^{0}_{S}K^{0}_{S}K^{0}_{S}$ decays at Belle

K. H. Kang  Kyungpook National University, Daegu 41566  H. Park  Kyungpook
National University, Daegu 41566  T. Higuchi  Kavli Institute for the Physics
and Mathematics of the Universe (WPI), University of Tokyo, Kashiwa 277-8583
K. Miyabayashi  Nara Women’s University, Nara 630-8506  K. Sumisawa  High
Energy Accelerator Research Organization (KEK), Tsukuba 305-0801  SOKENDAI
(The Graduate University for Advanced Studies), Hayama 240-0193  I. Adachi
High Energy Accelerator Research Organization (KEK), Tsukuba 305-0801
SOKENDAI (The Graduate University for Advanced Studies), Hayama 240-0193  J.
K. Ahn  Korea University, Seoul 02841  H. Aihara  Department of Physics,
University of Tokyo, Tokyo 113-0033  S. Al Said  Department of Physics,
Faculty of Science, University of Tabuk, Tabuk 71451  Department of Physics,
Faculty of Science, King Abdulaziz University, Jeddah 21589  D. M. Asner
Brookhaven National Laboratory, Upton, New York 11973  V. Aulchenko  Budker
Institute of Nuclear Physics SB RAS, Novosibirsk 630090  Novosibirsk State
University, Novosibirsk 630090  T. Aushev  Higher School of Economics (HSE),
Moscow 101000  R. Ayad  Department of Physics, Faculty of Science, University
of Tabuk, Tabuk 71451  V. Babu  Deutsches Elektronen–Synchrotron, 22607
Hamburg  S. Bahinipati  Indian Institute of Technology Bhubaneswar, Satya
Nagar 751007  A. M. Bakich  School of Physics, University of Sydney, New South
Wales 2006  P. Behera  Indian Institute of Technology Madras, Chennai 600036
C. Beleño  II. Physikalisches Institut, Georg-August-Universität Göttingen,
37073 Göttingen  J. Bennett  University of Mississippi, University,
Mississippi 38677  V. Bhardwaj  Indian Institute of Science Education and
Research Mohali, SAS Nagar, 140306  T. Bilka  Faculty of Mathematics and
Physics, Charles University, 121 16 Prague  J. Biswal  J. Stefan Institute,
1000 Ljubljana  G. Bonvicini  Wayne State University, Detroit, Michigan 48202
A. Bozek  H. Niewodniczanski Institute of Nuclear Physics, Krakow 31-342  M.
Bračko  University of Maribor, 2000 Maribor  J. Stefan Institute, 1000
Ljubljana  T. E. Browder  University of Hawaii, Honolulu, Hawaii 96822  M.
Campajola  INFN - Sezione di Napoli, 80126 Napoli  Università di Napoli
Federico II, 80055 Napoli  L. Cao  University of Bonn, 53115 Bonn  D.
Červenkov  Faculty of Mathematics and Physics, Charles University, 121 16
Prague  M.-C. Chang  Department of Physics, Fu Jen Catholic University, Taipei
24205  V. Chekelian  Max-Planck-Institut für Physik, 80805 München  A. Chen
National Central University, Chung-li 32054  B. G. Cheon  Department of
Physics and Institute of Natural Sciences, Hanyang University, Seoul 04763  K.
Chilikin  P.N. Lebedev Physical Institute of the Russian Academy of Sciences,
Moscow 119991  K. Cho  Korea Institute of Science and Technology Information,
Daejeon 34141  S.-K. Choi  Gyeongsang National University, Jinju 52828  Y.
Choi  Sungkyunkwan University, Suwon 16419  D. Cinabro  Wayne State
University, Detroit, Michigan 48202  S. Cunliffe  Deutsches
Elektronen–Synchrotron, 22607 Hamburg  N. Dash  Indian Institute of Technology
Bhubaneswar, Satya Nagar 751007  G. De Nardo  INFN - Sezione di Napoli, 80126
Napoli  Università di Napoli Federico II, 80055 Napoli  Z. Doležal  Faculty of
Mathematics and Physics, Charles University, 121 16 Prague  T. V. Dong  Key
Laboratory of Nuclear Physics and Ion-beam Application (MOE) and Institute of
Modern Physics, Fudan University, Shanghai 200443  S. Eidelman  Budker
Institute of Nuclear Physics SB RAS, Novosibirsk 630090  Novosibirsk State
University, Novosibirsk 630090  P.N. Lebedev Physical Institute of the Russian
Academy of Sciences, Moscow 119991  J. E. Fast  Pacific Northwest National
Laboratory, Richland, Washington 99352  T. Ferber  Deutsches
Elektronen–Synchrotron, 22607 Hamburg  B. G. Fulsom  Pacific Northwest
National Laboratory, Richland, Washington 99352  R. Garg  Panjab University,
Chandigarh 160014  V. Gaur  Virginia Polytechnic Institute and State
University, Blacksburg, Virginia 24061  N. Gabyshev  Budker Institute of
Nuclear Physics SB RAS, Novosibirsk 630090  Novosibirsk State University,
Novosibirsk 630090  A. Garmash  Budker Institute of Nuclear Physics SB RAS,
Novosibirsk 630090  Novosibirsk State University, Novosibirsk 630090  A. Giri
Indian Institute of Technology Hyderabad, Telangana 502285  P. Goldenzweig
Institut für Experimentelle Teilchenphysik, Karlsruher Institut für
Technologie, 76131 Karlsruhe  B. Golob  Faculty of Mathematics and Physics,
University of Ljubljana, 1000 Ljubljana  J. Stefan Institute, 1000 Ljubljana
D. Greenwald  Department of Physics, Technische Universität München, 85748
Garching  Y. Guan  University of Cincinnati, Cincinnati, Ohio 45221  O.
Hartbrich  University of Hawaii, Honolulu, Hawaii 96822  K. Hayasaka  Niigata
University, Niigata 950-2181  H. Hayashii  Nara Women’s University, Nara
630-8506  M. Hernandez Villanueva  University of Mississippi, University,
Mississippi 38677  W.-S. Hou  Department of Physics, National Taiwan
University, Taipei 10617  C.-L. Hsu  School of Physics, University of Sydney,
New South Wales 2006  K. Inami  Graduate School of Science, Nagoya University,
Nagoya 464-8602  G. Inguglia  Institute of High Energy Physics, Vienna 1050
A. Ishikawa  High Energy Accelerator Research Organization (KEK), Tsukuba
305-0801  SOKENDAI (The Graduate University for Advanced Studies), Hayama
240-0193  M. Iwasaki  Osaka City University, Osaka 558-8585  Y. Iwasaki  High
Energy Accelerator Research Organization (KEK), Tsukuba 305-0801  W. W. Jacobs
Indiana University, Bloomington, Indiana 47408  E.-J. Jang  Gyeongsang
National University, Jinju 52828  H. B. Jeon  Kyungpook National University,
Daegu 41566  S. Jia  Beihang University, Beijing 100191  Y. Jin  Department of
Physics, University of Tokyo, Tokyo 113-0033  K. K. Joo  Chonnam National
University, Gwangju 61186  A. B. Kaliyar  Tata Institute of Fundamental
Research, Mumbai 400005  G. Karyan  Deutsches Elektronen–Synchrotron, 22607
Hamburg  T. Kawasaki  Kitasato University, Sagamihara 252-0373  H. Kichimi
High Energy Accelerator Research Organization (KEK), Tsukuba 305-0801  C.
Kiesling  Max-Planck-Institut für Physik, 80805 München  C. H. Kim  Department
of Physics and Institute of Natural Sciences, Hanyang University, Seoul 04763
D. Y. Kim  Soongsil University, Seoul 06978  K.-H. Kim  Yonsei University,
Seoul 03722  S. H. Kim  Department of Physics and Institute of Natural
Sciences, Hanyang University, Seoul 04763  Y.-K. Kim  Yonsei University, Seoul
03722  K. Kinoshita  University of Cincinnati, Cincinnati, Ohio 45221  P.
Kodyš  Faculty of Mathematics and Physics, Charles University, 121 16 Prague
S. Korpar  University of Maribor, 2000 Maribor  J. Stefan Institute, 1000
Ljubljana  D. Kotchetkov  University of Hawaii, Honolulu, Hawaii 96822  P.
Križan  Faculty of Mathematics and Physics, University of Ljubljana, 1000
Ljubljana  J. Stefan Institute, 1000 Ljubljana  R. Kroeger  University of
Mississippi, University, Mississippi 38677  P. Krokovny  Budker Institute of
Nuclear Physics SB RAS, Novosibirsk 630090  Novosibirsk State University,
Novosibirsk 630090  R. Kulasiri  Kennesaw State University, Kennesaw, Georgia
30144  R. Kumar  Punjab Agricultural University, Ludhiana 141004  A. Kuzmin
Budker Institute of Nuclear Physics SB RAS, Novosibirsk 630090  Novosibirsk
State University, Novosibirsk 630090  Y.-J. Kwon  Yonsei University, Seoul
03722  J. S. Lange  Justus-Liebig-Universität Gießen, 35392 Gießen  I. S. Lee
Department of Physics and Institute of Natural Sciences, Hanyang University,
Seoul 04763  S. C. Lee  Kyungpook National University, Daegu 41566  J. Li
Kyungpook National University, Daegu 41566  L. K. Li  Institute of High Energy
Physics, Chinese Academy of Sciences, Beijing 100049  Y. B. Li  Peking
University, Beijing 100871  L. Li Gioi  Max-Planck-Institut für Physik, 80805
München  J. Libby  Indian Institute of Technology Madras, Chennai 600036  K.
Lieret  Ludwig Maximilians University, 80539 Munich  D. Liventsev  Virginia
Polytechnic Institute and State University, Blacksburg, Virginia 24061  High
Energy Accelerator Research Organization (KEK), Tsukuba 305-0801  J.
MacNaughton  University of Miyazaki, Miyazaki 889-2192  C. MacQueen  School of
Physics, University of Melbourne, Victoria 3010  M. Masuda  Earthquake
Research Institute, University of Tokyo, Tokyo 113-0032  T. Matsuda
University of Miyazaki, Miyazaki 889-2192  D. Matvienko  Budker Institute of
Nuclear Physics SB RAS, Novosibirsk 630090  Novosibirsk State University,
Novosibirsk 630090  P.N. Lebedev Physical Institute of the Russian Academy of
Sciences, Moscow 119991  M. Merola  INFN - Sezione di Napoli, 80126 Napoli
Università di Napoli Federico II, 80055 Napoli  H. Miyata  Niigata University,
Niigata 950-2181  R. Mizuk  P.N. Lebedev Physical Institute of the Russian
Academy of Sciences, Moscow 119991  Higher School of Economics (HSE), Moscow
101000  G. B. Mohanty  Tata Institute of Fundamental Research, Mumbai 400005
T. J. Moon  Seoul National University, Seoul 08826  T. Mori  Graduate School
of Science, Nagoya University, Nagoya 464-8602  T. Morii  Kavli Institute for
the Physics and Mathematics of the Universe (WPI), University of Tokyo,
Kashiwa 277-8583  M. Mrvar  Institute of High Energy Physics, Vienna 1050  R.
Mussa  INFN - Sezione di Torino, 10125 Torino  M. Nakao  High Energy
Accelerator Research Organization (KEK), Tsukuba 305-0801  SOKENDAI (The
Graduate University for Advanced Studies), Hayama 240-0193  H. Nakazawa
Department of Physics, National Taiwan University, Taipei 10617  Z. Natkaniec
H. Niewodniczanski Institute of Nuclear Physics, Krakow 31-342  M. Nayak
School of Physics and Astronomy, Tel Aviv University, Tel Aviv 69978  N. K.
Nisar  University of Pittsburgh, Pittsburgh, Pennsylvania 15260  S. Nishida
High Energy Accelerator Research Organization (KEK), Tsukuba 305-0801
SOKENDAI (The Graduate University for Advanced Studies), Hayama 240-0193  K.
Nishimura  University of Hawaii, Honolulu, Hawaii 96822  K. Ogawa  Niigata
University, Niigata 950-2181  S. Ogawa  Toho University, Funabashi 274-8510
H. Ono  Nippon Dental University, Niigata 951-8580  Niigata University,
Niigata 950-2181  Y. Onuki  Department of Physics, University of Tokyo, Tokyo
113-0033  P. Oskin  P.N. Lebedev Physical Institute of the Russian Academy of
Sciences, Moscow 119991  P. Pakhlov  P.N. Lebedev Physical Institute of the
Russian Academy of Sciences, Moscow 119991  Moscow Physical Engineering
Institute, Moscow 115409  G. Pakhlova  Higher School of Economics (HSE),
Moscow 101000  P.N. Lebedev Physical Institute of the Russian Academy of
Sciences, Moscow 119991  S. Pardi  INFN - Sezione di Napoli, 80126 Napoli
S.-H. Park  Yonsei University, Seoul 03722  S. Patra  Indian Institute of
Science Education and Research Mohali, SAS Nagar, 140306  S. Paul  Department
of Physics, Technische Universität München, 85748 Garching  T. K. Pedlar
Luther College, Decorah, Iowa 52101  R. Pestotnik  J. Stefan Institute, 1000
Ljubljana  L. E. Piilonen  Virginia Polytechnic Institute and State
University, Blacksburg, Virginia 24061  T. Podobnik  Faculty of Mathematics
and Physics, University of Ljubljana, 1000 Ljubljana  J. Stefan Institute,
1000 Ljubljana  V. Popov  Higher School of Economics (HSE), Moscow 101000  E.
Prencipe  Forschungszentrum Jülich, 52425 Jülich  M. T. Prim  Institut für
Experimentelle Teilchenphysik, Karlsruher Institut für Technologie, 76131
Karlsruhe  M. V. Purohit  Okinawa Institute of Science and Technology, Okinawa
904-0495  M. Ritter  Ludwig Maximilians University, 80539 Munich  M. Röhrken
Deutsches Elektronen–Synchrotron, 22607 Hamburg  A. Rostomyan  Deutsches
Elektronen–Synchrotron, 22607 Hamburg  N. Rout  Indian Institute of Technology
Madras, Chennai 600036  M. Rozanska  H. Niewodniczanski Institute of Nuclear
Physics, Krakow 31-342  G. Russo  Università di Napoli Federico II, 80055
Napoli  D. Sahoo  Tata Institute of Fundamental Research, Mumbai 400005  Y.
Sakai  High Energy Accelerator Research Organization (KEK), Tsukuba 305-0801
SOKENDAI (The Graduate University for Advanced Studies), Hayama 240-0193  S.
Sandilya  University of Cincinnati, Cincinnati, Ohio 45221  A. Sangal
University of Cincinnati, Cincinnati, Ohio 45221  L. Santelj  Faculty of
Mathematics and Physics, University of Ljubljana, 1000 Ljubljana  J. Stefan
Institute, 1000 Ljubljana  T. Sanuki  Department of Physics, Tohoku
University, Sendai 980-8578  V. Savinov  University of Pittsburgh, Pittsburgh,
Pennsylvania 15260  G. Schnell  University of the Basque Country UPV/EHU,
48080 Bilbao  IKERBASQUE, Basque Foundation for Science, 48013 Bilbao  J.
Schueler  University of Hawaii, Honolulu, Hawaii 96822  C. Schwanda  Institute
of High Energy Physics, Vienna 1050  A. J. Schwartz  University of Cincinnati,
Cincinnati, Ohio 45221  Y. Seino  Niigata University, Niigata 950-2181  K.
Senyo  Yamagata University, Yamagata 990-8560  M. E. Sevior  School of
Physics, University of Melbourne, Victoria 3010  M. Shapkin  Institute for
High Energy Physics, Protvino 142281  V. Shebalin  University of Hawaii,
Honolulu, Hawaii 96822  J.-G. Shiu  Department of Physics, National Taiwan
University, Taipei 10617  E. Solovieva  P.N. Lebedev Physical Institute of the
Russian Academy of Sciences, Moscow 119991  S. Stanič  University of Nova
Gorica, 5000 Nova Gorica  M. Starič  J. Stefan Institute, 1000 Ljubljana  Z.
S. Stottler  Virginia Polytechnic Institute and State University, Blacksburg,
Virginia 24061  T. Sumiyoshi  Tokyo Metropolitan University, Tokyo 192-0397
W. Sutcliffe  University of Bonn, 53115 Bonn  M. Takizawa  Showa
Pharmaceutical University, Tokyo 194-8543  J-PARC Branch, KEK Theory Center,
High Energy Accelerator Research Organization (KEK), Tsukuba 305-0801
Theoretical Research Division, Nishina Center, RIKEN, Saitama 351-0198  U.
Tamponi  INFN - Sezione di Torino, 10125 Torino  K. Tanida  Advanced Science
Research Center, Japan Atomic Energy Agency, Naka 319-1195  F. Tenchini
Deutsches Elektronen–Synchrotron, 22607 Hamburg  K. Trabelsi  LAL, Univ.
Paris-Sud, CNRS/IN2P3, Université Paris-Saclay, Orsay 91898  M. Uchida  Tokyo
Institute of Technology, Tokyo 152-8550  T. Uglov  P.N. Lebedev Physical
Institute of the Russian Academy of Sciences, Moscow 119991  Higher School of
Economics (HSE), Moscow 101000  Y. Unno  Department of Physics and Institute
of Natural Sciences, Hanyang University, Seoul 04763  S. Uno  High Energy
Accelerator Research Organization (KEK), Tsukuba 305-0801  SOKENDAI (The
Graduate University for Advanced Studies), Hayama 240-0193  P. Urquijo  School
of Physics, University of Melbourne, Victoria 3010  Y. Ushiroda  High Energy
Accelerator Research Organization (KEK), Tsukuba 305-0801  SOKENDAI (The
Graduate University for Advanced Studies), Hayama 240-0193  G. Varner
University of Hawaii, Honolulu, Hawaii 96822  A. Vinokurova  Budker Institute
of Nuclear Physics SB RAS, Novosibirsk 630090  Novosibirsk State University,
Novosibirsk 630090  V. Vorobyev  Budker Institute of Nuclear Physics SB RAS,
Novosibirsk 630090  Novosibirsk State University, Novosibirsk 630090  P.N.
Lebedev Physical Institute of the Russian Academy of Sciences, Moscow 119991
C. H. Wang  National United University, Miao Li 36003  E. Wang  University of
Pittsburgh, Pittsburgh, Pennsylvania 15260  M.-Z. Wang  Department of Physics,
National Taiwan University, Taipei 10617  P. Wang  Institute of High Energy
Physics, Chinese Academy of Sciences, Beijing 100049  X. L. Wang  Key
Laboratory of Nuclear Physics and Ion-beam Application (MOE) and Institute of
Modern Physics, Fudan University, Shanghai 200443  S. Watanuki  Department of
Physics, Tohoku University, Sendai 980-8578  E. Won  Korea University, Seoul
02841  X. Xu  Soochow University, Suzhou 215006  B. D. Yabsley  School of
Physics, University of Sydney, New South Wales 2006  W. Yan  University of
Science and Technology of China, Hefei 230026  S. B. Yang  Korea University,
Seoul 02841  H. Ye  Deutsches Elektronen–Synchrotron, 22607 Hamburg  J. H. Yin
Institute of High Energy Physics, Chinese Academy of Sciences, Beijing 100049
C. Z. Yuan  Institute of High Energy Physics, Chinese Academy of Sciences,
Beijing 100049  Y. Yusa  Niigata University, Niigata 950-2181  Z. P. Zhang
University of Science and Technology of China, Hefei 230026  V. Zhilich
Budker Institute of Nuclear Physics SB RAS, Novosibirsk 630090  Novosibirsk
State University, Novosibirsk 630090  V. Zhukova  P.N. Lebedev Physical
Institute of the Russian Academy of Sciences, Moscow 119991  V. Zhulanov
Budker Institute of Nuclear Physics SB RAS, Novosibirsk 630090  Novosibirsk
State University, Novosibirsk 630090

###### Abstract

We measure the time-dependent $C\\!P$ violation parameters in $B^{0}\to
K^{0}_{S}K^{0}_{S}K^{0}_{S}$ decays using $772\times 10^{6}~{}B\mkern
1.5mu\overline{\mkern-1.5muB\mkern-0.1mu}\mkern 0.1mu$ pairs collected at the
$\Upsilon(4S)$ resonance with the Belle detector at the KEKB asymmetric-energy
$e^{+}e^{-}$ collider. The obtained mixing-induced and direct $C\\!P$
asymmetries are $-0.71\pm 0.23~{}{\rm(stat)}\pm 0.05~{}{\rm(syst)}$ and
$0.12\pm 0.16~{}{\rm(stat)}\pm 0.05~{}{\rm(syst)}$, respectively. These values
are consistent with the Standard Model predictions. The significance of
$C\\!P$ violation differs from zero by 2.5 standard deviations.

###### pacs:

13.20.He, 14.40.Nd

††preprint:  Belle preprint 2020-17 KEK preprint 2020-34

In the Standard Model (SM), $C\\!P$ violation in the quark sector is described
by an irreducible phase in the Kobayashi-Maskawa (KM) mechanism KM . The
charmless three-body decay $B^{0}\to K^{0}_{S}K^{0}_{S}K^{0}_{S}$ is mediated
by the $b\to sq\mkern 1.5mu\overline{\mkern-1.5muq\mkern-0.1mu}\mkern 0.1mu$
quark transition, which is prohibited in the lowest-order SM interaction.
Instead, this $C\\!P$-even decay occurs via a “penguin” amplitude, as shown in
Fig. 1. Deviations from the SM expectations for $C\\!P$-violating parameters
provide sensitivity to new physics np1 .

Time-dependent $C\\!P$ violation can be caused by interference between the
decay and mixing amplitudes. When one of the neutral $B$ mesons produced from
the $\Upsilon(4S)$ decays into a $C\\!P$ eigenstate, $f_{C\\!P}$, at time
$t_{C\\!P}$, and the other into a flavor-distinguishable final state, $f_{\rm
tag}$, at time $t_{\rm tag}$, the time-dependent decay rate is given by Sanda

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

where $\Delta t\equiv t_{C\\!P}-t_{\rm tag}$, measured in the center-of-mass
(CM) frame, and the $C\\!P$-violating parameters $\mathcal{S}$ and
$\mathcal{A}$ are related to mixing-induced and direct $C\\!P$ violation,
respectively. Here the flavor $q$ is $+1$ ($-1$) for $\mkern
1.5mu\overline{\mkern-1.5muB\mkern-0.1mu}\mkern 0.1mu^{0}$ ($B^{0}$),
$\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}=-\sin 2\phi_{1}$ and $\mathcal{A}=0$ in $B^{0}\to
K^{0}_{S}K^{0}_{S}K^{0}_{S}$, where
$\phi_{1}\equiv\arg[-V_{cd}V_{cb}^{*}/V_{td}V_{tb}^{*}]$ phi1notation .

Previous measurements of $\mathcal{S}$ at Belle and BaBar have yielded values
of $-0.30\pm 0.32~{}{\rm(stat)}\pm 0.08~{}{\rm(syst)}$ using $535\times
10^{6}$ $B\mkern 1.5mu\overline{\mkern-1.5muB\mkern-0.1mu}\mkern 0.1mu$ pairs,
and $-0.94^{+0.24}_{-0.21}~{}{\rm(stat)}\pm 0.06~{}{\rm(syst)}$ using
$468\times 10^{6}$ $B\mkern 1.5mu\overline{\mkern-1.5muB\mkern-0.1mu}\mkern
0.1mu$ pairs, respectively b03ksbelle ; b03ksbabar . To search for physics
beyond the SM containing a new $C\\!P$-violating phase, we measure
$\mathcal{S}$ and $\mathcal{A}$ in $B^{0}\to K^{0}_{S}K^{0}_{S}K^{0}_{S}$
decays with the final Belle data set of $772\times 10^{6}$ $B\mkern
1.5mu\overline{\mkern-1.5muB\mkern-0.1mu}\mkern 0.1mu$ pairs.

Figure 1:  Penguin amplitude for the $B^{0}\to K^{0}_{S}K^{0}_{S}K^{0}_{S}$
decays.

The Belle detector is a large-solid-angle magnetic spectrometer that consists
of a silicon vertex detector (SVD), a 50-layer central drift chamber (CDC), an
array of aerogel threshold Cherenkov counters (ACC), a barrel-like arrangement
of time-of-flight scintillation counters (TOF), and an electromagnetic
calorimeter comprised of CsI(Tl) crystals (ECL). These detector components are
located inside a superconducting solenoid coil that provides a 1.5 T magnetic
field. An iron flux-return located outside the magnetic coil is instrumented
to detect $K_{L}^{0}$ mesons and identify muons (KLM). The detector is
described in detail elsewhere Belle . Two inner detector configurations were
used. A 2.0 cm radius beampipe with a double-wall beryllium structure and a
three-layer SVD were used for the first sample of $152\times 10^{6}~{}B\mkern
1.5mu\overline{\mkern-1.5muB\mkern-0.1mu}\mkern 0.1mu$ pairs, while a 1.5 cm
radius beampipe, a four-layer SVD and a small-inner-cell CDC were used to
record the remaining $620\times 10^{6}~{}B\mkern
1.5mu\overline{\mkern-1.5muB\mkern-0.1mu}\mkern 0.1mu$ pairs svd2 . The latter
data sample has been reprocessed with improved software, which incorporates an
improved vertex reconstruction pbf ; ccbar .

The $\Upsilon(4S)$ is produced at the KEKB asymmetric-energy $e^{+}e^{-}$
collider KEKB  with a Lorentz boost (${\beta\gamma}$) of 0.425; it
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 CM frame. The Lorentz boost introduces a sufficient
distance between the $B$ and $\mkern
1.5mu\overline{\mkern-1.5muB\mkern-0.1mu}\mkern 0.1mu$ decay vertices to be
measurable nearly along with the $z$ axis, which is antiparallel to the
$e^{+}$ beam direction. The distance is related to $\Delta
t\approx(z_{C\\!P}-z_{\rm tag})/c\beta\gamma$, where $z_{C\\!P}$ and $z_{\rm
tag}$ are the coordinates of the decay positions of $f_{C\\!P}$ and $f_{\rm
tag}$, respectively. To avoid the large backgrounds accompanying $\gamma$ and
$\pi^{0}$ detection, we reconstruct the $K_{S}^{0}$ only through its decay to
two charged pions. The event selection and measurement of $C\\!P$ violation
parameters are optimized using Monte Carlo (MC) events. The MC events are
generated by EVTGEN Evtgen , and the detector response is modeled using GEANT3
GEANT3 . We simulate the $B$-meson decay to three $K_{S}^{0}$ as uniformly
distributed in the available phase space.

The $K_{S}^{0}$ is selected from charged pion pairs using a neural network
(NN) nisks ; NB  with 13 inputs: the $K_{S}^{0}$ momentum in the lab frame
($>$ 0.06 GeV/$c$); the distance between the pion tracks in the $z$-direction
($<$ 20 cm); the flight length in the $x$-$y$ plane; the angle between the
$K_{S}^{0}$ momentum and the vertex displacement vector; the angle between the
$K_{S}^{0}$ momentum and the pion momentum; for each daughter pion, the
distance of closest approach to the interaction point (IP); the existence of
the SVD hits; and the number of axial- and stereo-wire hits in the CDC. The
mass ranges allowed are 0.474 GeV/$c^{2}$ $<M(\pi^{+}\pi^{-})<$ 0.522
GeV/$c^{2}$ when only one pion hits the SVD, and 0.478 GeV/$c^{2}$
$<M(\pi^{+}\pi^{-})<$ 0.517 GeV/$c^{2}$ otherwise.

To identify the signal $B$-decay, we use the energy difference $\Delta E\equiv
E_{\rm beam}-E_{B}$ and the beam-energy-constrained mass ${M_{\rm
bc}\equiv\sqrt{(E_{\rm beam})^{2}-|\vec{p}_{B}|^{2}c^{2}}/c^{2}}$, where
$E_{\rm beam}$ is the beam energy, and $E_{B}$ and $\vec{p}_{B}$ are the
energy and momentum, respectively, of the $B^{0}$ candidate. All quantities
are evaluated in the CM frame. The $B^{0}$ candidates are required to lie in
the region of $M_{\rm bc}>5.2$ GeV$/c^{2}$ and $|\Delta E|<0.2$ GeV.

We find that seven percent of the events have more than one $B^{0}$ candidate.
When there are multiple $B^{0}$ candidates in an event, we choose the one with
the smallest $\chi^{2}$ as defined by
$\chi^{2}=\sum_{i=1}^{3}[(M_{i}(\pi^{+}\pi^{-})-m_{K^{0}_{S}})/\sigma_{i}]^{2}$,
where $M_{i}(\pi^{+}\pi^{-})$ and $\sigma_{i}$ are the invariant mass and mass
resolution for the $i$-th $K^{0}_{S}$, respectively, and $m_{K^{0}_{S}}$ is
the nominal $K^{0}_{S}$ mass PDG .

The dominant source of background is continuum $e^{+}e^{-}\to q\mkern
1.5mu\overline{\mkern-1.5muq\mkern-0.1mu}\mkern 0.1mu~{}(q=u,d,s,c)$ events.
To suppress this background, we use another NN with the following inputs: the
cosine of the polar angle of the $B^{0}$-candidate flight direction in the CM
frame ($\cos\theta_{B}$); the cosine of the angle between the thrust axis of
the $B^{0}$ candidate and that of the rest of the event ($\cos\theta_{T}$);
and a likelihood ratio obtained from modified Fox-Wolfram moments KSFW . The
NN outputs ($\mathcal{O}_{\rm NN}$) range between $+1$ and $-1$, where
$\mathcal{O}_{\rm NN}$ close to $+1$ ($-1$) indicates a signal-like
(background-like) event. The $\mathcal{O}_{\rm NN}$ criterion is obtained by
maximizing a figure-of-merit (FOM = $N_{\rm sig}/\sqrt{N_{\rm sig}+N_{q\mkern
1.5mu\overline{\mkern-1.5muq\mkern-0.1mu}\mkern 0.1mu}}$), where $N_{\rm sig}$
and $N_{q\mkern 1.5mu\overline{\mkern-1.5muq\mkern-0.1mu}\mkern 0.1mu}$ are
the number of signal and continuum MC events. The FOM is maximal at the value
of $\mathcal{O}_{\rm NN}$ = 0.08. The signal region is defined as $M_{\rm
bc}>5.27$ GeV$/c^{2}$ and $|\Delta E|<0.1$ GeV. The $\mathcal{O}_{\rm NN}$
requirement retains 81% of the signal and reduces continuum background by a
factor of 10 in the signal region.

The $\mathcal{O}_{\rm NN}$ is transformed by using the formula

${\mathcal{O}^{\prime}_{\rm NN}}=\log(\frac{{\mathcal{O}_{\rm
NN}~{}-~{}\mathcal{O}_{\rm NN,min}}}{{\mathcal{O}_{\rm
NN,max}~{}-~{}\mathcal{O}_{\rm NN}}}).$ (2)

The value of $\mathcal{O}_{\rm NN,min}$ is selected to be 0.08, thus
maximizing the FOM. The value of $\mathcal{O}_{\rm NN,max}$ is set to 0.99,
the highest observed value of $\mathcal{O}_{\rm NN}$.

Among the $B^{0}\to K^{0}_{S}K^{0}_{S}K^{0}_{S}$ decays, there are quasi-two-
body intermediate states, where both $b\to c$ and $b\to s$ transitions
contribute. The former contaminates the measurement of $C\\!P$ violation in
$b\to s$ transitions. Among possible $b\to c$ transition-induced $B$ decays,
we expect significant contributions solely from
$B^{0}\to\chi_{c0}(\chi_{c0}\to K^{0}_{S}K^{0}_{S})K^{0}_{S}$ decays; this
contribution is estimated to be $16.3\pm 3.1$ events. We use the invariant
mass, $M_{K^{0}_{S}K^{0}_{S}}$, to veto $B^{0}\to\chi_{c0}K^{0}_{S}$ decays:
we reject the $B^{0}$ candidate if any $K^{0}_{S}$ pair among its decay
products has an invariant mass within $\pm 2\sigma$ of the nominal mass, where
$\sigma$ is the reconstructed $\chi_{c0}$ mass resolution from simulation. The
veto range is 3.388 GeV/$c^{2}$ $<$ $M_{K^{0}_{S}K^{0}_{S}}$ $<$ 3.444
GeV/$c^{2}$, and this removes 83% of the $B^{0}\to\chi_{c0}K^{0}_{S}$ decays.

To identify the $B$ meson flavor, a flavor tagging algorithm TaggingNIM  is
used that utilizes inclusive properties of particles not associated with the
signal decay. This algorithm returns the value of $q$ (defined earlier) and a
tagging quality variable $r$. The latter varies from $r=0$ for no tagging
information to $r=1$ for unambiguous flavor assignment. The probability
density function (PDF) for signal events modifies Eq. (1) by taking the wrong-
tag fraction, $w$, and its difference between $B^{0}$ and $\mkern
1.5mu\overline{\mkern-1.5muB\mkern-0.1mu}\mkern 0.1mu^{0}$, $\Delta w$, into
account:

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

The events are categorized into seven $r$ bins. For each of these bins, $w$
and $\Delta w$ are determined by high statistics flavor-specific $B$ meson
decays ccbar .

The parameter $\Delta t$ in Eq. (3) is determined through vertex
reconstruction for the signal $B^{0}$ ($B_{C\\!P}$) and the accompanying
$B^{0}$ ($B_{\rm tag}$). For reconstruction of the $B_{C\\!P}$ vertex, we use
those $K_{S}^{0}$ trajectories in which both daughter pions have SVD hits in
the $z$ direction. $B_{C\\!P}$ can have up to three $K^{0}_{S}$ that satisfy
this requirement. According to the signal MC, 14% of events do not have any
$K^{0}_{S}$ producing sufficient SVD hits. The IP information is incorporated
as a virtual straight track along the $z$ axis, called the “IP tube”, to
provide a constraint for kinematical fits to reconstruct the $B$ decay vertex.
The $B_{C\\!P}$ vertex is obtained from the available $K^{0}_{S}$ trajectories
and the IP tube. The IP tube is also used to reconstruct the $B_{\rm tag}$
vertex using the charged tracks not assigned to $B_{C\\!P}$, as described in
more detail in Ref. vertexres . Because of this treatment, we can reconstruct
the $B_{C\\!P}$ ($B_{\rm tag}$) vertex with only one $K^{0}_{S}$ trajectory
(charged track).

Events with poorly reconstructed vertices are rejected by requiring:

  1. (1)

$|\Delta t|<70$ ps;

  2. (2)

a vertex quality for each of $B_{C\\!P}$ and $B_{\rm tag}$ of less than 50,
when the vertex is constrained by multiple tracks; and

  3. (3)

uncertainties on the $z$ position of the vertices for both $B_{C\\!P}$ and
$B_{\rm tag}$ of less than 0.2 mm when the vertex is constrained by multiple
tracks, and less than 0.5 mm when the $B_{C\\!P}$ ($B_{\rm tag}$) vertex is
constrained by a single $K^{0}_{S}$ trajectory (single track).

The vertex quality is $\chi^{2}$ per degree of freedom, where $\chi^{2}$ is
obtained from the vertex fit without the IP tube constraint. After the poorly
reconstructed events are discarded, the remaining data events amount to 73% of
the total number of events in the signal region.

We determine the signal yield by performing an extended unbinned maximum-
likelihood fit. The fit is done in the $\Delta E$-$M_{\rm
bc}$-$\mathcal{O}^{\prime}_{\rm NN}$ three-dimensional space, with each PDF
expressed as the product of the one-dimensional PDFs. The signal PDFs are
modeled as a double Gaussian, a Gaussian, and an asymmetric Gaussian,
respectively, and the background PDFs are modeled as a first-order polynomial,
an ARGUS function ARGUS , and an asymmetric Gaussian, respectively. The
parameters of the signal PDF are fixed according to fits to the signal MC; the
parameters of the background PDF are left free in the fit. The fit result is
shown in Fig. 2. In the signal region, the signal yield is 258 $\pm$ 17 events
and the purity is 74%, where $-4.72<$ $\mathcal{O}^{\prime}_{\rm NN}$ $<7.24$
is further required for the signal region.

Figure 2:  Result of signal-extraction fit to $\Delta E$ (top) in the
redefined signal region of $M_{\rm bc}$ and $\mathcal{O}^{\prime}_{\rm NN}$,
$M_{\rm bc}$ (middle) in the redefined signal region of $\Delta E$ and
$\mathcal{O}^{\prime}_{\rm NN}$, and the $\mathcal{O}^{\prime}_{\rm NN}$
(bottom) in the redefined signal region of $\Delta E$ and $M_{\rm bc}$
distributions. The red, blue, and green dashed lines represent the total,
signal, and background PDFs. The points with error bars represent data.

We determine the $C\\!P$-violating parameters $\mathcal{A}$ and $\mathcal{S}$
by performing a second unbinned maximum-likelihood fit. The contribution to
the likelihood function from the $j$th event is

$\begin{split}\mathcal{P}_{j}&=(1-f_{\rm ol})\times\left[f^{\rm
sig}_{j}\left\\{\int d(\Delta t^{\prime})R(\Delta t_{j}-\Delta
t^{\prime})\times\right.\right.\\\ &\left.\left.\mathcal{P}_{\rm sig}(\Delta
t^{\prime},q_{j})\vphantom{\int}\right\\}+(1-f^{\rm sig}_{j})\mathcal{P}_{\rm
bkg}(\Delta t_{j})\right]\\\ &+f_{\rm ol}\mathcal{P}_{\rm ol}(\Delta
t_{j}),\end{split}$ (4)

where $R(\Delta t)$ is a resolution function. The resolution function consists
of three components: detector resolutions for $B_{C\\!P}$ and $B_{\rm tag}$,
nonprimary track effects for $B_{\rm tag}$, and a kinematical approximation
due to the difference in the lab momentum of $B_{C\\!P}$ and $B_{\rm tag}$
owing to non-zero CM momentum vertexres . In Eq. (4), $f_{\rm ol}$ and
$\mathcal{P}_{\rm ol}$ are the fraction and PDF, respectively, of outlier
events for a very long tail shape in the $\Delta t$ distribution, and
$f_{j}^{\rm sig}$ is the signal fraction obtained from the signal extraction.
$\mathcal{P}_{\rm bkg}$ is the $\Delta t$ distribution of background events.
This PDF is the sum of a $\delta$-function and an exponential function, both
convolved with a background resolution function. The parameters of these
functions are determined by fitting events in a data sideband region, defined
as 0.15 GeV $<|\Delta E|<$ 0.20 GeV or $M_{\rm bc}<5.26$ GeV$/c^{2}$. The
unbinned maximum-likelihood fit in the signal region is used to determine the
$C\\!P$ violation parameters, where the world average values are used for
$\tau_{B^{0}}$ and $\Delta m_{d}$ PDG . The measured $\mathcal{S}$ and
$\mathcal{A}$ are $-0.71\pm 0.23$ and $0.12\pm 0.16$, respectively, where the
uncertainties are statistical. The background-subtracted $\Delta t$ and
asymmetry distributions are shown in Fig. 3 ccbar .

Figure 3:  Background-subtracted $\Delta t$ distribution (top) and asymmetry
distribution (bottom) obtained from data. In the $\Delta t$ distribution
graph, the red solid line and open circles represent the fitted curve and data
for $\mkern 1.5mu\overline{\mkern-1.5muB\mkern-0.1mu}\mkern 0.1mu^{0}$, while
the blue dashed line and filled circles represent the fitted curve and data
for $B^{0}$, respectively. In the asymmetry distribution graph, the points
represent data and the solid line represents the fitted curve.

The estimated systematic uncertainties are summarized in Table 1. The
systematic uncertainty is calculated by varying the fixed parameters in the
fit for $\mathcal{S}$ and $\mathcal{A}$, and all uncertainties are summed in
quadrature. For inputs obtained from data and MC, we use $1\sigma$ and
$2\sigma$ variations, respectively.

Table 1: Systematic uncertainties Source | $\delta\mathcal{S}$ | $\delta\mathcal{A}$
---|---|---
Vertex reconstruction | 0.031 | 0.038
Flavor tagging | 0.002 | 0.004
Resolution function | 0.016 | 0.014
Physics parameters | 0.004 | 0.001
Fit bias | 0.012 | 0.009
Signal fraction | 0.024 | 0.021
Background $\Delta t$ shape | 0.016 | 0.001
SVD misalignment | 0.004 | 0.005
$\Delta z$ bias | 0.002 | 0.004
Tag-side interference | 0.001 | 0.008
Total | 0.047 | 0.047

The systematic uncertainty on the vertex reconstruction is determined by
varying the $x$-$y$ plane smearing parameter for the IP profile, charged track
requirements for the $B_{\rm tag}$ vertex reconstruction, criteria to discard
poorly reconstructed vertices for measurement of $C\\!P$ violation, and
correction of helix parameter errors for vertexing. The systematic
uncertainties due to the parameters $w$ and $\Delta w$ are estimated by
varying these parameters by their uncertainties. We vary each resolution
function parameter by its uncertainty. For physics parameters, we calculate
differences in $\mathcal{S}$ and $\mathcal{A}$ by varying world average values
of $\tau_{B^{0}}$ and $\Delta m_{B^{0}}$. For the systematic uncertainty on
the fit bias, we measure $C\\!P$ violation parameters using the signal MC
events, mixed signal MC with continuum toy MC events, and mixed signal toy MC
with continuum toy MC events at the ratio expected from data. The larger of
the difference between the input $\mathcal{S}$ ($\mathcal{A}$) and the output
$\mathcal{S}$ ($\mathcal{A}$) and the statistical error on the output
$\mathcal{S}$ ($\mathcal{A}$) is considered as the uncertainty due to the fit
bias. The systematic uncertainties due to the signal fraction and background
$\Delta t$ shape are obtained by varying the parameter $f^{\rm sig}_{j}$ and
$\mathcal{P}_{\rm bkg}$ in Eq. (4), respectively. For possible SVD
misalignment, $\Delta z$ bias, and tag-side interference, we quote the
systematic uncertainties obtained from a study of a large control sample of
$B\to(c\mkern 1.5mu\overline{\mkern-1.5muc\mkern-0.1mu}\mkern 0.1mu)K_{S}^{0}$
decays ccbar .

The significance, taking both statistical and systematic uncertainties into
account, is calculated using a two-dimensional Feldman-Cousins approach FD .
The significance of $C\\!P$ violation is determined to be 2.5$\sigma$ away
from $(0,0)$ as depicted in Fig. 4, which shows the two-dimensional confidence
contour in the $\mathcal{S}$ and $\mathcal{A}$ plane.

Figure 4:  Confidence contour for $\mathcal{S}$ and $\mathcal{A}$. The
contours represent 1$\sigma$, 2$\sigma$, and 3$\sigma$ from inner to outer.
The point with error bars represents the measured values of $\mathcal{S}$ and
$\mathcal{A}$. The circle is the physical region defined by
$\mathcal{S}^{2}+\mathcal{A}^{2}\leq 1$.

In summary, we have studied time-dependent $C\\!P$ violation in $B^{0}\to
K^{0}_{S}K^{0}_{S}K^{0}_{S}$ decays using the final data set containing
$772\times 10^{6}$ $B\mkern 1.5mu\overline{\mkern-1.5muB\mkern-0.1mu}\mkern
0.1mu$ collected at the Belle experiment. The NN methods for $K^{0}_{S}$
selection and background suppression and improved vertex reconstruction, along
with increased statistics, result in a more precise measurement than the
previous Belle one. The measured values of $\mathcal{S}$ and $\mathcal{A}$ are

$\displaystyle\mathcal{S}$ $\displaystyle=$ $\displaystyle-0.71\pm
0.23\rm{(stat)}\pm 0.05\rm{(syst)},$ $\displaystyle\mathcal{A}$
$\displaystyle=$ $\displaystyle 0.12\pm 0.16\rm{(stat)}\pm 0.05\rm{(syst)}.$

The results are consistent with the world average value of $-\sin 2\phi_{1}$
($-0.70$) HFLV  as well as with the SM prediction. These result supersede our
previous measurements in Ref. b03ksbelle . We thank the KEKB group for the
excellent operation of the accelerator; the KEK cryogenics group for the
efficient operation of the solenoid; and the KEK computer group, and the
Pacific Northwest National Laboratory (PNNL) Environmental Molecular Sciences
Laboratory (EMSL) computing group for strong computing support; and the
National Institute of Informatics, and Science Information NETwork 5 (SINET5)
for valuable network support. We acknowledge support from the Ministry of
Education, Culture, Sports, Science, and Technology (MEXT) of Japan, the Japan
Society for the Promotion of Science (JSPS), and the Tau-Lepton Physics
Research Center of Nagoya University; the Australian Research Council
including grants DP180102629, DP170102389, DP170102204, DP150103061,
FT130100303; Austrian Science Fund (FWF); the National Natural Science
Foundation of China under Contracts No. 11435013, No. 11475187, No. 11521505,
No. 11575017, No. 11675166, No. 11705209; Key Research Program of Frontier
Sciences, Chinese Academy of Sciences (CAS), Grant No. QYZDJ-SSW-SLH011; the
CAS Center for Excellence in Particle Physics (CCEPP); the Shanghai Pujiang
Program under Grant No. 18PJ1401000; the Ministry of Education, Youth and
Sports of the Czech Republic under Contract No. LTT17020; the Carl Zeiss
Foundation, the Deutsche Forschungsgemeinschaft, the Excellence Cluster
Universe, and the VolkswagenStiftung; the Department of Science and Technology
of India; the Istituto Nazionale di Fisica Nucleare of Italy; National
Research Foundation (NRF) of Korea Grant Nos. 2016R1D1A1B01010135,
2016R1D1A1B02012900, 2018R1A2B3003643, 2018R1A6A1A06024970,
2018R1D1A1B07047294, 2019K1A3A7A09033840, 2019R1I1A3A01058933,
2019R1A6A3A01096585; Radiation Science Research Institute, Foreign Large-size
Research Facility Application Supporting project, the Global Science
Experimental Data Hub Center of the Korea Institute of Science and Technology
Information and KREONET/GLORIAD; the Polish Ministry of Science and Higher
Education and the National Science Center; the Ministry of Science and Higher
Education of the Russian Federation, Agreement 14.W03.31.0026; University of
Tabuk research grants S-1440-0321, S-0256-1438, and S-0280-1439 (Saudi
Arabia); the Slovenian Research Agency; Ikerbasque, Basque Foundation for
Science, Spain; the Swiss National Science Foundation; the Ministry of
Education and the Ministry of Science and Technology of Taiwan; and the United
States Department of Energy and the National Science Foundation.

## References

  * (1) M. Kobayashi and T. Maskawa, Prog. Theor. Phys. 49, 652 (1973).
  * (2) Y. Grossman and M. P. Worah, Phys. Lett. B 395, 241 (1997).
  * (3) A. B. Carter and A. I. Sanda, Phys. Rev. Lett. 45, 952 (1980); A. B. Carter and A. I. Sanda, Phys. Rev. D 23, 1567 (1981); I. I. Bigi and A. I. Sanda, Nucl. Phys. 193, 85 (1981).
  * (4) Another naming convention, $\beta(=\phi_{1})$, is also used in the literature.
  * (5) K.-F. Chen et al. (Belle Collaboration), Phys. Rev. Lett. 98, 031802 (2007).
  * (6) J.P. Lees et al. (BaBar Collaboration), Phys. Rev. D 85, 054023 (2012).
  * (7) A. Abashian et al. (Belle Collaboration), Nucl. Instrum. Methods Phys. Res., Sect. A 479, 117 (2002); also, see the detector section in J. Brodzicka et al., Prog. Theor. Exp. Phys. 2012, 04D001 (2012).
  * (8) Z. Natkaniec et al. (Belle SVD2 Group), Nucl. Instrum. Methods Phys. Res., Sect. A 560, 1 (2006).
  * (9) Ed. A.J. Bevan, B. Golob, Th. Mannel, S. Prell, and B.D. Yabsley, Eur. Phys. J. C74 (2014) 3026, SLAC-PUB-15968, KEK Preprint 2014-3, p. 40.
  * (10) I. Adachi et al. (Belle Collaboration), Phys. Rev. Lett. 108, 171802 (2012).
  * (11) S. Kurokawa and E. Kikutani, Nucl. Instrum. Methods Phys. Res., Sect. A 499, 1 (2003), and other papers included in this volume; T. Abe et al., Prog. Theor. Exp. Phys. 2013, 03A001 (2013), and following articles up to 03A011.
  * (12) D. J. Lange et al., Nucl. Instrum. Methods Phys. Res., Sect. A 462, 152 (2001).
  * (13) R. Brun et al., CERN DD/EE/84-1 (1984).
  * (14) H. Nakano, Ph.D Thesis, Tohoku University (2014) Chapter 4, unpublished, https://tohoku.repo.nii.ac.jp/?action=pages_view_main&active_action=repository_view_main_item_detail&item_id=70563&item_no=1&page_id=33&block_id=38.
  * (15) M. Feindt and U. Kerzel, Nucl. Instrum. Methods Phys. Res., Sect. A 559, 190 (2006).
  * (16) P. Zyla et al. (Particle Data Group), Prog. Theor. Exp. Phys. 2020, 083C01 (2020).
  * (17) S. H. Lee et al. (Belle Collaboration), Phys. Rev. Lett. 91, 261801 (2003).
  * (18) H. Kakuno et al., Nucl. Instrum. Methods Phys. Res., Sect. A 533, 516 (2004).
  * (19) H. Tajima et al., Nucl. Instrum. Methods Phys. Res., Sect. A 533, 370 (2004).
  * (20) H. Albrecht et al. (ARGUS Collaboration), Phys. Lett. B 241, 278 (1990).
  * (21) G. J. Feldman and R. D. Cousins, Phys. Rev. D 57, 3873 (1998).
  * (22) Y. S. Amhis et al. (HFLAV Group), arXiv:1909.12524 (2019); and updates at https://hflav.web.cern.ch/.


# Teaching a Robot to Walk
Using Reinforcement Learning

Jack Dibachi and Jacob Azoulay
Stanford University
AA228: Decision Making under Uncertainty
<EMAIL_ADDRESS>—<EMAIL_ADDRESS>

###### Abstract

> Classical control techniques such as PID and LQR have been used effectively
> in maintaining a system state, but these techniques become more difficult to
> implement when the model dynamics increase in complexity and sensitivity.
> For adaptive robotic locomotion tasks with several degrees of freedom, this
> task becomes infeasible with classical control techniques. Instead,
> reinforcement learning can train optimal walking policies with ease. We
> apply deep Q-learning and augmented random search (ARS) to teach a simulated
> two-dimensional bipedal robot how to walk using the OpenAI Gym
> `BipedalWalker-v3` environment. Deep Q-learning did not yield a high reward
> policy, often prematurely converging to suboptimal local maxima likely due
> to the coarsely discretized action space. ARS, however, resulted in a better
> trained robot, and produced an optimal policy which officially “solves” the
> `BipedalWalker-v3` problem. Various naive policies—including a random
> policy, a manually encoded inch forward policy, and a stay still policy—were
> used as benchmarks to evaluate the proficiency of the learning algorithm
> results.

## Problem Statement

Creating feedback controllers capable of maintaining stability in a bipedal
robot is a complicated task, particularly because of model complexity and
sensitivity. Even when controllers can be made, they are extremely sensitive
to variations in the environment, such as irregularities in the walking
surface or unexpected obstacles (?). Moreover, mitigating state uncertainty by
including range sensors and contact sensors requires the control designer to
account for any number of possible environments. In short, creating a multi-
input-multi-output feedback controller to address this control task is
intractable.

Instead, the robot’s dynamics can be modeled and tested in numerous simulation
trials, and use reinforcement learning to determine the optimal policy mapping
the robot’s state to the action required to walk forward. In the OpenAI Gym’s
`BipedalWalker-v3` environment, the reward system is defined by incrementally
gaining 300 points if the finish line is reached, losing 100 points if robot
falls over, and gradually losing points based on the amount of control effort
exerted from the motors. Solving the problem is defined by collecting over 300
reward points for 100 consecutive episodes. We use OpenAI Gym’s
`BipedalWalker-v3` for the robot dynamics and environment, and an augmented
random search (ARS) strategy for determining the optimal policy.

## Other Work

Robot locomotion has been treated as a learning task more frequently in the
past few years. ? apply a model-based approach to a simple five-link bipedal
robot using Poincare maps. The policy, however, is only determined using the
robot’s intrinsic state, such as the joint angles, positions, and velocities,
rather than including measurements such as foot contacts and range
measurements. This approach, then, cannot be used in environments with surface
irregularities (`BipedalWalker-v3`) or obstacles (`BipedalWalkerHardcore-v3`).
? effectively implement reinforcement learning to teach the Cassie robot to
walk, though their approach uses a gait library as input to a deep
reinforcement learning network (?). While effective, their objective is to
enable robust online planning for environmental uncertainty on physical
robotic platforms. ? outline an ARS strategy, validated in a number of MuJoCo
locomotion environments via OpenAI Gym with promising results. ARS is a
particularly useful learning strategy for locomotion tasks, especially
`BipedalWalker-v3`, as the policy space is filled with locally optimal
locomotion strategies (?). Other model-free learning strategies require many
more episodes to solve these kinds of locomotion tasks, if they solve them at
all. For this reason, we use ARS to train the bipedal walker.

## Approach

OpenAI Gym’s `BipedalWalker-v3` environment provides a model of a five-link
bipedal robot, depicted in Figure 1. The robot state is a vector with 24
elements: $\theta,\dot{x},\dot{y},\omega$ of the hull center of mass (white),
$\theta,\omega$ of each joint (two green, two orange), contacts with the
ground (red), and 10 lidar scans of the ground (red line). The action space,
then, is a four element vector representing the torque command to the motor at
each of the four joints (?).

Figure 1: The robot in the BipedalWalker-v3 environment. Red flag indicates
starting position.

This robot has to navigate through the environment, which is a randomly seeded
ground with an uneven surface. The lidar scans can help detect this
unevenness, though they are most useful in the hardcore version of the
environment where obstacles and stairs are present. To find the optimal
policy, which maps the robot states to the robot actions yielding the highest
reward, we try both deep Q-learning and ARS.

### Deep Q-Learning

Traditional Q-learning techniques incrementally estimate the action value
function $Q(s,a)$, using information at each step to update the Q-function.
This populates a lookup table containing action values for every possible
state-action pair. This works well for Markov decision process (MDP) problems
with discrete state and action spaces, however is insufficient for problems
with large continuous state and action spaces. Fortunately, the Q-learning
algorithm can be adapted to problems with continuous spaces by using a
parametric approximation of the action value function $Q_{\theta}(s,a)$ such
as a neural network with the weights and biases of the network serving as the
function hyperparameters. The input layer consists of one node for each state
variable and the output layer consists of nodes representing utilities for
each discretized action.

Deep Q-learning relies on the minimization of the loss between the approximate
action value function $Q_{\theta}(s,a)$ and the true optimal action value
function $Q^{*}(s,a)$ denoted as

$l(\theta)=\frac{1}{2}\operatorname{\mathbb{E}}_{(s,a)\sim\pi^{*}}[(Q^{*}(s,a)-Q_{\theta}(s,a))^{2}]$

Applying gradient descent to minimize this loss and using samples to
approximate the expectation results in the update rule

$\theta\leftarrow\theta+\alpha(Q^{*}(s,a)-Q_{\theta}(s,a))\nabla_{\theta}Q_{\theta}(s,a)$

where $\alpha$ is the learning rate which determines the size of each update
step. Lastly, because the true optimal policy $Q^{*}(s,a)$ is unknown, it can
be approximated, yielding the finalized update rule:

$\theta\leftarrow\theta+\alpha(r+\gamma\max_{a^{\prime}}Q_{\theta}(s^{\prime},a^{\prime})-Q_{\theta}(s,a))\nabla_{\theta}Q_{\theta}(s,a)$

To implement this algorithm for the bipedal walker, the continuous action
space is first discretized into a predetermined number of bins. The continuous
action space consists of four action variables (one for each joint), each
taking a value between $-1$ and $1$. These were discretized and constrained to
only take one of three values $-1,0,$ and $1$, which results in $3^{4}=81$
possible actions represented as a vector of four numbers.

The neural network representation has 24 input nodes, each representing one of
the 24 state variables observed by the agent. The output layer has 81 nodes
(one per possible action) which correspond to the approximated utilities of
taking an action given the observed state. The output layer uses a linear
activation function in order to produce unbounded linear utility outputs.
There are two hidden layers, each with 55 nodes using $tanh$ activation
functions. A visual representation of the neural network structure is depicted
in Figure 2.

Figure 2: Q-learning neural network structure. Algorithm 1 Deep Q-Learning
Algorithm

Initialize hyperparameters: $\alpha$, $\gamma$, action bin size

Initialize neural net Q with weights $\theta$

for batch in training duration do

$\epsilon\leftarrow$ decay($\epsilon$)

$Q\\_\leftarrow Q$

for episode in batch size do

$s\leftarrow$ reset environment

while episode is not $done$ do

$a=\begin{cases}\operatorname*{arg\,max}_{a^{\prime}}Q(s,a^{\prime}),&\mbox{with
prob }1-\epsilon\\\ \mbox{random action},&\mbox{with prob
}\epsilon\end{cases}$

$s^{\prime},r$ = observed new state, reward

$Q_{target}=r+(1-done)\gamma\max_{a^{\prime}}Q\\_(s^{\prime},a^{\prime})$

$Q_{output}=Q(s,a)$

$\theta$ updated with gradient descent:

$d\theta\leftarrow
d\theta+\frac{\partial(Q\\_target-Q\\_output)^{2}}{\partial\theta}$

$s\leftarrow s^{\prime}$

end while

end for

end for

To find the optimal walking policy, we encourage exploration of the action
space using $\epsilon$-greedy exploration, where the robot takes some random
action with probability $\epsilon$. The amount of exploration needed decreases
as the robot improves its gait, so we decay the $\epsilon$ term over time.
Furthermore, to minimize unstable learning, training episodes are run in
batches, with each batch using a static copy of the current estimation of the
action value function (represented as a neural network) when computing
$r+\gamma\max_{a^{\prime}}Q_{\theta}(s^{\prime},a^{\prime})$.

The neural network representations and gradient descent loss minimization are
implemented using the Python module `pytorch` (?).

### Augmented Random Search (ARS)

While Q-learning required the discretization of the action space, ARS allows
for both a continuous state space and action space, often resulting in
efficient optimal policies in fewer training iterations (?). ARS explores the
policy parameter space, as opposed to the action and state space, by taking a
set of parameter samples with zero mean Gaussian noise and executing rollouts
for each of those samples.

The parameters of the model $\theta$ take the form of an $n{\times}m$ matrix
where $n$ is the number of action variables and $m$ is the number of state
variables. A state observation is matrix multiplied to produce a desired
action to take. In ARS, several parameter noise matrices $\delta$ are sampled,
resulting in $\theta+v\delta$ and $\theta-v\delta$ pairs where $\theta$ is the
current parameterization, $v$ is a step size scalar that determines how much
noise is introduced to the model, and $\delta$ is a random noise matrix of the
same dimensions as $\theta$. If in total there are $h$ number of random noise
matrices $\delta$ then there are $2h$ matrices, each of which are used to
conduct rollouts. These noisy parameterizations all lie within a hypersphere
of a specified radius around the current parameterization $\theta$ (?).

Rewards for each rollout are collected, and their standard deviation $\sigma$
is computed and used to normalize the final update step. The $h$ $v\delta$
matrices are then sorted in descending order of maximum reward using the
criteria $\max(r_{\theta+v\delta_{i}},r_{\theta-v\delta_{i}})$ for each
parameterization pair $i$, where $r_{\theta+v\delta_{i}}$ denotes the reward
received from the rollout using a model parameterization $\theta+v\delta_{i}$,
and $r_{\theta-v\delta_{i}}$ denotes the reward received from the rollout
using a model parameterization $\theta-v\delta_{i}$. The top $m$ pairs are
then used to update $\theta$ using the following update rule:

$\theta\leftarrow\theta+\frac{\alpha}{m\sigma}\sum_{i=1}^{m}[(r_{\theta+v\delta_{i}}-r_{\theta-v\delta_{i}})*\delta_{i}]$

where $\alpha$ represents the learning rate.

As the agent learns and the average reward size increases, the summation over
$r_{\theta+v\delta_{i}}-r_{\theta-v\delta_{i}}$ will increase as well,
resulting in proportionately larger step sizes as the policy improves
overtime. To mitigate this inconsistency in step size, the learning rate
$\alpha$ is normalized by the standard deviation of the rollout rewards
$\sigma$. This normalization is the primary difference between basic random
search (BRS) and augmented random search (ARS). The reward standard deviation
$\sigma$ overtime is plotted in Figure 3.

Figure 3: ARS reward standard deviation $\sigma$ over 1500 episodes.

Lastly, each state observation is normalized before it is fed through the
model to determine an optimal action, eliminating unintentional bias towards
state variables with larger absolute bounds.

Algorithm 2 ARS Algorithm

Initialize hyperparameters: $\alpha$, $\theta$, number of noise matrices $h$,
noise scalar $v$, number of best performing noise matrices $m$

for episode in training duration do

$\delta=[\delta_{1}:\delta_{h}]$ where each $\delta_{i}$ is a randomly sampled

zero mean noise sample matrix the size of $\theta$

for each $\delta_{i}$ in $\delta$ do

$r_{\theta+v\delta_{i}}=$ reward collected from rollout using

$\theta+v\delta_{i}$ to determine action $a$ at each step

observation $s$

$r_{\theta-v\delta_{i}}=$ reward collected from rollout using

$\theta-v\delta_{i}$ to determine action $a$ at each step

observation $s$

$r\\_list\leftarrow$ store
$(r_{\theta+v\delta_{i}},r_{\theta-v\delta_{i}},v\delta_{i})$ in list

end for

$\sigma=$ std of
$\\{r_{\theta+v\delta_{1}},r_{\theta-v\delta_{1}},...,r_{\theta+v\delta_{h}},r_{\theta-v\delta_{h}}\\}$

sort $r\\_list$ using criterion
$\max(r_{\theta+v\delta_{i}},r_{\theta-v\delta_{i}})$

discard $h-m$ worst performing tuples with lowest

rewards

$\theta\leftarrow\theta+\frac{\alpha}{m\sigma}\sum_{i=1}^{m}[(r_{\theta+v\delta_{i}}-r_{\theta-v\delta_{i}})*\delta_{i}]$

end for

## Analysis

The primary metric indicating how well the agent learns is the change in the
average collected reward per episode. The ultimate goal is to achieve high
reward collecting performance in as few training iterations as possible. The
performance of each algorithm is also evaluated compared to three naive
approaches: a random policy model, a manually encoded policy which consists of
a set of actions that inch the walker forward, and a policy that has the
walker take one step forward and stay still for the duration of the episode
iteration.

With both Q-learning and ARS, there are many hyperparameters that can be
adjusted when training the bipedal robot how to walk. For instance, in the
above implementation of Q-learning the following need to be specified:
learning rate $\alpha$, exploration parameter $\epsilon$, discount factor
$\gamma$, action space discretization bin sizes, neural network structure
(number and size of hidden layers), layer activation functions, batch size,
and total number of episode iterations. With ARS, the following
hyperparameters must be specified: learning rate $\alpha$, noise scalar $v$,
the number of random noise matrices $h$, the number of highest ranked matrices
to use $m$, and the total number of episode iterations. It is infeasible to
comprehensively explore all combinations of hyperparameters. Through
experimental trial runs, however, adjusting the learning rate proved to result
in large improvements in learning.

## Results

As benchmarks for the ARS performance, we compare accumulated rewards per
episode with other agents. The random agent selects a random sample along a
uniform distribution from the action space, regardless of the state.
Evaluating the random agent for 1000 episodes gives clusters of accumulated
reward just below 0 and just above -100 in Figure 4. Thus, the agent either
falls over after some effort forward, or gradually loses points by twitching
about, but stays upright. The random agent rarely accumulates a net positive
reward in any individual episode, let alone the success threshold of 300
reward points.

Figure 4: Random policy rolled out for 1000 episodes.

However, the random agent still outperforms the predetermined agent shown in
Figure 5, which periodically repeats an inching forward action regardless of
state. This agent has fewer episodes with zero reward and a higher percentage
of episodes with approximately -100 reward, as the periodic agent either falls
over or expends excess control effort. Because the periodic agent tends to
stand for longer than the random agent, it spends longer expending control
effort without making up for it by covering linear distance, thus resulting in
a lower average reward than the random agent.

Figure 5: Predetermined crawling agent rolled out for 1000 episodes.

A more advanced benchmark is the deep Q-learning agent with decaying
$\epsilon$-greedy exploration. Although this algorithm takes significantly
more episodes (and thus run-time) to obtain results, Figure 6 shows the
majority of episodes accumulate rewards closer to zero than the random or
periodic agents, outperforming both these agents.

Figure 6: Deep Q-Network performance over 10000 episodes with decaying
$\epsilon$-greedy exploration and learning rate $\alpha=0.001$.

Implementing ARS for the same learning rate gives better performance than the
Q-learning approach, in that the average reward is greater than zero, and the
robot rarely accumulates rewards less than zero, as shown in two separate
experiments in Figure 7. This is indicative of the robot taking a single step
to stay upright, then staying still to avoid losing reward points. However,
there is not enough exploration in the policy space for the agent to take
sequential steps. Aggressively increasing the learning rate should remedy this
issue.

Figure 7: ARS learning curve with $\alpha=0.001$, averaged over two separate
experiments.

The ARS agent with a learning rate of $\alpha=0.02$ demonstrates significantly
better performance in accumulated reward than the other agents we have covered
so far. Results over three experiments are shown in Figure 8, illustrating the
agent’s ability to take multiple steps and maintain balance. The agent’s
learning curve varies broadly over these experiments because the learning
process is stochastic, so the episode in which the agent learns to take a
certain stride can happen in a wide range of possible times.

Figure 8: ARS learning curve with $\alpha=0.02$, averaged over three separate
experiments.

Increasing the learning rate to $\alpha=0.06$ demonstrates even higher
performance, making sharp increases in reward accumulation in approximately
100 episodes, and clearing the reward threshold in less than 400 episodes.
Results from two experiments are shown in Figure 9.

Figure 9: ARS learning curve with $\alpha=0.06$, averaged over two separate
experiments.

Algorithms and Run Times
---
Algorithm |  Total Run Time (s) |  Episodes
Random |  181 |  1000
Periodic |  298 |  1000
Stationary |  354 |  1000
Deep Q |  3861 |  10000
ARS $\alpha=0.001$ |  1212 |  1500
ARS $\alpha=0.02$ |  8428 |  1500
ARS $\alpha=0.06$ |  6864 |  1500

## Discussion

When choosing a reinforcement learning method to implement in a system, there
is no all-encompassing algorithm that universally produces the best results.
It is necessary to not only consider the structure of the problem, but also
experiment with various algorithms and various hyperparameters for those
algorithms.

In context of the bipedal walker, the realm of possible algorithms to choose
from is narrowed due to the large continuous state space and action space,
producing the best results with a model-free method. While with some
continuous problems it is possible to discretize the state space and action
space, doing so is computationally intractable for the problem at hand since
producing a sufficiently granular state space and action space would result in
memory overflow as the discretized space scales exponentially with bin size.
This motivated the decision to use deep Q-learning and ARS.

The advantage of deep Q-learning lies in the ability to have a continuous
state space; however, our implementation still requires discretizing the
action-space. The robot did not yield successful learning results likely
because of this discretization. Due to the exponential growth in discretized
action space with each additional bin, the possible actions the agent could
execute were very limited compared to the full continuous action space. Many
of the trial runs produced policies which have the robot take one step forward
and stay still for the duration of the episode. The robot learned that
attempting to move posed a high risk of falling (which is penalized with a
large negative reward), and as a result chooses to remain still and receive a
higher reward, albeit still negative or close to zero.

Perhaps a more aggressive exploration strategy or a finer action space
discretization could nudge the agent out of this local optimum. Furthermore,
future work could involve experimenting with different training batch sizes
and overall length of training.

ARS, however, yielded consistently reliable learning. According to OpenAI, the
`BipedalWalker-v3` is considered “solved” when the agent obtains an average
reward of at least 300 over 100 consecutive episodes. Several runs of ARS in
this study achieved this metric, in particular the runs using learning rates
$\alpha$ in the range of $[0.01,0.06]$. One of the benefits of ARS is that it
allows for a continuous state space and continuous action space, an advantage
likely contributing to the remarkable learning.

Future work could involve applying the algorithm to the
`BipedalWalkerHardcore-v3` environment, which poses a more challenging problem
introducing obstacles and pitfalls. Because the hardcore version of the
problem still requires the robot to learn how to walk on flat ground, it may
be beneficial to initialize future training models with the output model
parameterization $\theta$ of the original non-hardcore version of the problem.
Doing so will allow the walker to spend more time exploring the expanded state
space imposed by the new obstacles rather than relearning how to walk as well.

## Conclusion

This study presented an application of deep Q-learning and ARS to control a
bipedal robot in a simulated environment. While deep Q-learning yielded poor
results, ARS demonstrated quick and efficient learning for this application.
ARS can be extended beyond a two-dimensional simulated environment, and can
ultimately be deployed for real world applications of robot autonomy.

## Group Member Contributions

Jacob:

  * •

Deep Q-Learning code using `pytorch` (?)

  * •

ARS code (?)

  * •

Collected raw training data

  * •

Wrote the following sections of this paper:

Abstract

Approach

Analysis

Discussion

Conclusion

Jack:

  * •

Deep Q-Learning code using `TensorForce` (?)

  * •

ARS code (?)

  * •

Produced data visualizations using `pgfplots`

  * •

Wrote the following sections of this paper:

Problem Statement

Other Work

Results

References

## References

  * [2016] Brockman, G.; Cheung, V.; Pettersson, L.; Schneider, J.; Schulman, J.; Tang, J.; and Zaremba, W.  2016\.  Openai gym.
  * [2017] Kuhnle, A.; Schaarschmidt, M.; and Fricke, K.  2017\.  Tensorforce: a tensorflow library for applied reinforcement learning.  Web page.
  * [2021] Li, Z.; Cheng, X.; Peng, X. B.; Abbeel, P.; Levine, S.; Berseth, G.; and Sreenath, K.  2021\.  Reinforcement learning for robust parameterized locomotion control of bipedal robots.  CoRR abs/2103.14295.
  * [2018] Mania, H.; Guy, A.; and Recht, B.  2018\.  Simple random search provides a competitive approach to reinforcement learning.  CoRR abs/1803.07055.
  * [2004] Morimoto, J.; Cheng, G.; Atkeson, C.; and Zeglin, G.  2004\.  A simple reinforcement learning algorithm for biped walking.  In IEEE International Conference on Robotics and Automation, 2004\. Proceedings. ICRA ’04. 2004, volume 3, 3030–3035 Vol.3.
  * [2019] Paszke, A.; Gross, S.; Massa, F.; Lerer, A.; Bradbury, J.; Chanan, G.; Killeen, T.; Lin, Z.; Gimelshein, N.; Antiga, L.; Desmaison, A.; Kopf, A.; Yang, E.; DeVito, Z.; Raison, M.; Tejani, A.; Chilamkurthy, S.; Steiner, B.; Fang, L.; Bai, J.; and Chintala, S.  2019\.  Pytorch: An imperative style, high-performance deep learning library.  In Wallach, H.; Larochelle, H.; Beygelzimer, A.; d'Alché-Buc, F.; Fox, E.; and Garnett, R., eds., Advances in Neural Information Processing Systems 32. Curran Associates, Inc.  8024–8035.
  * [2018] Skow, C.  2018\.  Coding demos from the school of ai’s move37 course.  Web page.