Neutron Identification and Correlation with Spallation Isotopes in SK Detector Simulation, Tesis de Física Médica

¡Descarga Neutron Identification and Correlation with Spallation Isotopes in SK Detector Simulation y más Tesis en PDF de Física Médica solo en Docsity! New Methods and Simulations for Cosmogenic Induced Spallation Removal in Super-Kamiokande-IV S. Locke,7 A. Coffani,11 K. Abe,1, 43 C. Bronner,1 Y. Hayato,1, 43 M. Ikeda,1 S. Imaizumi,1 H. Ito,1 J. Kameda,1, 43 Y. Kataoka,1 M. Miura,1, 43 S. Moriyama,1, 43 Y. Nagao,1 M. Nakahata,1, 43 Y. Nakajima,1, 43 S. Nakayama,1, 43 T. Okada,1 K. Okamoto,1 A. Orii,1 G. Pronost,1 H. Sekiya,1, 43 M. Shiozawa,1, 43 Y. Sonoda,1 Y. Suzuki,1 A. Takeda,1, 43 Y. Takemoto,1 A. Takenaka,1 H. Tanaka,1 T. Yano,1 K. Hirade,1 Y. Kanemura,1 S. Miki,1 S. Watabe,1 S. Han,2 T. Kajita,2, 43 K. Okumura,2, 43 T. Tashiro,2 J. Xia,2 X. Wang,2 G. D. Megias,3 D. Bravo-Berguño,4 L. Labarga,4 Ll. Marti,4 B. Zaldivar,4 B. W. Pointon,6, 47 F. d. M. Blaszczyk,5 E. Kearns,5, 43 J. L. Raaf,5 J. L. Stone,5, 43 L. Wan,5 T. Wester,5 J. Bian,7 N. J. Griskevich,7 W. R. Kropp,7, ∗ S. Mine,7 A. Yankelevic,7 M. B. Smy,7, 43 H. W. Sobel,7, 43 V. Takhistov,7, 43 J. Hill,8 J. Y. Kim,9 I. T. Lim,9 R. G. Park,9 B. Bodur,10 K. Scholberg,10, 43 C. W. Walter,10, 43 L. Bernard,11 O. Drapier,11 S. El Hedri,11 A. Giampaolo,11 M. Gonin,11 Th. A. Mueller,11 P. Paganini,11 B. Quilain,11 A. D. Santos,11 T. Ishizuka,12 T. Nakamura,13 J. S. Jang,14 J. G. Learned,15 L. H. V. Anthony,16 A. A. Sztuc,16 Y. Uchida,16 D. Martin,16 M. Scott,16 V. Berardi,17 M. G. Catanesi,17 E. Radicioni,17 N. F. Calabria,18 L. N. Machado,18 G. De Rosa,18 G. Collazuol,19 F. Iacob,19 M. Lamoureux,19 N. Ospina,19 M. Mattiazzi,19 L. Ludovici,20 Y. Nishimura,21 Y. Maewaka,21 S. Cao,22 M. Friend,22 T. Hasegawa,22 T. Ishida,22 T. Kobayashi,22 M. Jakkapu,22 T. Matsubara,22 T. Nakadaira,22 K. Nakamura,22, 43 Y. Oyama,22 K. Sakashita,22 T. Sekiguchi,22 T. Tsukamoto,22 Y. Nakano,24 T. Shiozawa,24 A. T. Suzuki,24 Y. Takeuchi,24, 43 S. Yamamoto,24 Y. Kotsar,24 H. Ozaki,24 A. Ali,25 Y. Ashida,25 J. Feng,25 S. Hirota,25 A. K. Ichikawa,25 T. Kikawa,25 M. Mori,25 T. Nakaya,25, 43 R. A. Wendell,25, 43 K. Yasutome,25 P. Fernandez,26 N. McCauley,26 P. Mehta,26 K. M. Tsui,26 Y. Fukuda,27 Y. Itow,28, 29 H. Menjo,28 T. Niwa,28 K. Sato,28 M. Tsukada,28 P. Mijakowski,30 J. Lagoda,30 S. M. Lakshmi,30 J. Zalipska,30 C. K. Jung,31 C. Vilela,31 M. J. Wilking,31 C. Yanagisawa,31, † J. Jiang,31 K. Hagiwara,32 M. Harada,32 T. Horai,32 H. Ishino,32 S. Ito,32 Y. Koshio,32, 43 W. Ma,32 N. Piplani,32 S. Sakai,32 H. Kitagawa,32 G. Barr,33 D. Barrow,33 L. Cook,33, 43 A. Goldsack,33, 43 S. Samani,33 D. Wark,33, 38 F. Nova,34 T. Boschi,23 F. Di Lodovico,23 M. Taani,23 S. Zsoldos,23 J. Gao,23 J. Migenda,23 J. Y. Yang,35 S. J. Jenkins,36 M. Malek,36 J. M. McElwee,36 O. Stone,36 M. D. Thiesse,36 L. F. Thompson,36 H. Okazawa,37 K. Nakamura,53 S. B. Kim,39 I. Yu,39 J. W. Seo,39 K. Nishijima,40 M. Koshiba,41, ∗ K. Iwamoto,42 N. Ogawa,42 M. Yokoyama,42, 43 K. Martens,43 M. R. Vagins,43, 7 K. Nakagiri,43, 41 M. Kuze,44 S. Izumiyama,44 T. Yoshida,44 M. Inomoto,45 M. Ishitsuka,45 R. Matsumoto,45 K. Ohta,45 M. Shinoki,45 T. Suganuma,45 T. Kinoshita,45 J. F. Martin,46 H. A. Tanaka,46 T. Towstego,46 R. Akutsu,47 M. Hartz,47 A. Konaka,47 P. de Perio,47 N. W. Prouse,47 S. Chen,48 B. D. Xu,48 Y. Zhang,48 M. Posiadala-Zezula,49 B. Richards,50 B. Jamieson,51 J. Walker,51 A. Minamino,52 K. Okamoto,52 G. Pintaudi,52 and R. Sasaki52 1Kamioka Observatory, Institute for Cosmic Ray Research, University of Tokyo, Kamioka, Gifu 506-1205, Japan 2Research Center for Cosmic Neutrinos, Institute for Cosmic Ray Research, University of Tokyo, Kashiwa, Chiba 277-8582, Japan 3Institute for Cosmic Ray Research, University of Tokyo, Kashiwa, Chiba 277-8582, Japan 4Department of Theoretical Physics, University Autonoma Madrid, 28049 Madrid, Spain 5Department of Physics, Boston University, Boston, MA 02215, USA 6Department of Physics, British Columbia Institute of Technology, Burnaby, BC, V5G 3H2, Canada 7Department of Physics and Astronomy, University of California, Irvine, Irvine, CA 92697-4575, USA 8Department of Physics, California State University, Dominguez Hills, Carson, CA 90747, USA 9Institute for Universe and Elementary Particles, Chonnam National University, Gwangju 61186, Korea 10Department of Physics, Duke University, Durham NC 27708, USA 11Ecole Polytechnique, IN2P3-CNRS, Laboratoire Leprince-Ringuet, F-91120 Palaiseau, France 12Junior College, Fukuoka Institute of Technology, Fukuoka, Fukuoka 811-0295, Japan 13Department of Physics, Gifu University, Gifu, Gifu 501-1193, Japan 14GIST College, Gwangju Institute of Science and Technology, Gwangju 500-712, Korea 15Department of Physics and Astronomy, University of Hawaii, Honolulu, HI 96822, USA 16Department of Physics, Imperial College London , London, SW7 2AZ, United Kingdom 17Dipartimento Interuniversitario di Fisica, INFN Sezione di Bari and Università e Politecnico di Bari, I-70125, Bari, Italy 18Dipartimento di Fisica, INFN Sezione di Napoli and Università di Napoli, I-80126, Napoli, Italy 19Dipartimento di Fisica, INFN Sezione di Padova and Università di Padova, I-35131, Padova, Italy 20INFN Sezione di Roma and Università di Roma “La Sapienza”, I-00185, Roma, Italy 21Department of Physics, Keio University, Yokohama, Kanagawa, 223-8522, Japan 22High Energy Accelerator Research Organization (KEK), Tsukuba, Ibaraki 305-0801, Japan 23Department of Physics, King’s College London, London, WC2R 2LS, UK 24Department of Physics, Kobe University, Kobe, Hyogo 657-8501, Japan ar X iv :2 11 2. 00 09 2v 1 [ he p- ex ] 3 0 N ov 2 02 1 2 25Department of Physics, Kyoto University, Kyoto, Kyoto 606-8502, Japan 26Department of Physics, University of Liverpool, Liverpool, L69 7ZE, United Kingdom 27Department of Physics, Miyagi University of Education, Sendai, Miyagi 980-0845, Japan 28Institute for Space-Earth Environmental Research, Nagoya University, Nagoya, Aichi 464-8602, Japan 29Kobayashi-Maskawa Institute for the Origin of Particles and the Universe, Nagoya University, Nagoya, Aichi 464-8602, Japan 30National Centre For Nuclear Research, 02-093 Warsaw, Poland 31Department of Physics and Astronomy, State University of New York at Stony Brook, NY 11794-3800, USA 32Department of Physics, Okayama University, Okayama, Okayama 700-8530, Japan 33Department of Physics, Oxford University, Oxford, OX1 3PU, United Kingdom 34Rutherford Appleton Laboratory, Harwell, Oxford, OX11 0QX, UK 35Department of Physics, Seoul National University, Seoul 151-742, Korea 36Department of Physics and Astronomy, University of Sheffield, S3 7RH, Sheffield, United Kingdom 37Department of Informatics in Social Welfare, Shizuoka University of Welfare, Yaizu, Shizuoka, 425-8611, Japan 38STFC, Rutherford Appleton Laboratory, Harwell Oxford, and Daresbury Laboratory, Warrington, OX11 0QX, United Kingdom 39Department of Physics, Sungkyunkwan University, Suwon 440-746, Korea 40Department of Physics, Tokai University, Hiratsuka, Kanagawa 259-1292, Japan 41The University of Tokyo, Bunkyo, Tokyo 113-0033, Japan 42Department of Physics, University of Tokyo, Bunkyo, Tokyo 113-0033, Japan 43Kavli Institute for the Physics and Mathematics of the Universe (WPI), The University of Tokyo Institutes for Advanced Study, University of Tokyo, Kashiwa, Chiba 277-8583, Japan 44Department of Physics,Tokyo Institute of Technology, Meguro, Tokyo 152-8551, Japan 45Department of Physics, Faculty of Science and Technology, Tokyo University of Science, Noda, Chiba 278-8510, Japan 46Department of Physics, University of Toronto, ON, M5S 1A7, Canada 47TRIUMF, 4004 Wesbrook Mall, Vancouver, BC, V6T2A3, Canada 48Department of Engineering Physics, Tsinghua University, Beijing, 100084, China 49Faculty of Physics, University of Warsaw, Warsaw, 02-093, Poland 50Department of Physics, University of Warwick, Coventry, CV4 7AL, UK 51Department of Physics, University of Winnipeg, MB R3J 3L8, Canada 52Department of Physics, Yokohama National University, Yokohama, Kanagawa, 240-8501, Japan 53Department of Physics, Faculty of Science, Tohoku University, Sendai, Miyagi, 980-8578, Japan (Dated: December 2, 2021) Radioactivity induced by cosmic muon spallation is a dominant source of backgrounds for O(10) MeV neutrino interactions in water Cherenkov detectors. In particular, it is crucial to reduce backgrounds to measure the solar neutrino spectrum and find neutrino interactions from distant su- pernovae. In this paper we introduce new techniques to locate muon-induced hadronic showers and efficiently reject spallation backgrounds. Applying these techniques to the solar neutrino analysis with an exposure of 2790 × 22.5 kton.day increases the signal efficiency by 12.6%, approximately corresponding to an additional year of detector running. Furthermore, we present the first spallation simulation at SK, where we model hadronic interactions using FLUKA. The agreement between the isotope yields and shower pattern in this simulation and in the data gives confidence in the accuracy of this simulation, and thus opens the door to use it to optimize muon spallation removal in new data with gadolinium-enhanced neutron capture detection. I. INTRODUCTION Spallation from cosmic-ray muons produces radioac- tive isotopes and induces one of the largest backgrounds for the Super-Kamiokande (SK) neutrino signal between ∼6 and ∼25 MeV. Reducing this background is pivotal for the success of many different analyses in this energy range, and has major implications in solar, reactor, and supernova relic neutrino searches. Specifically in SK, cos- mic ray muons and the showers they produce sometimes interact with 16O nuclei within the detector volume, pro- ducing radioactive isotopes. ∗ Deceased. † also at BMCC/CUNY, Science Department, New York, New York, 1007, USA. Showers induced by the muons are primarily elec- tromagnetic in nature (γ-rays and electrons) as a re- sult of delta-ray production, pair production, and bremsstrahlung. However, there is also the possibility for muons to produce secondary particles in the form of neutrons, pions, and others. Recent simulation studies have shown that most spallation isotopes are produced by these secondary particles, with only 11% of isotopes being made directly from muons, and hence these iso- topes can be found up to several meters away from the muon track [1–3]. These isotopes then undergo mostly β or βγ decays, mimicking the expected signal for neutrino interactions. Their half-lives extend from milliseconds to seconds, and thus can be much larger than the time in- terval between two muons in SK, where the muon rate is about 2 Hz. Identifying spallation isotopes by pairing them with their parent muons is therefore particularly 5 III. EVENT RECONSTRUCTION A. Cosmic Rays and Track Fits Muons pass through SK at a rate of approximately 2 Hz, and their tracks are reconstructed from the PMTs hit within the detector. The muon reconstruction used for this analysis (outlined in [16, 17]) accurately fits tracks as well as categorizing muons. After removing PMT noise hits, the fitter makes an initial guess on the track using the earliest hit PMT with at least three neigh- boring hits as an entry point and time, and the largest cluster of charge as the exit point. The track parame- ters are then varied and a likelihood dependent on the expected Cherenkov light pattern is maximized to get a final track fit. Muons are categorized based on character- istics of the observed light, with the four different muon categories described by: 1. Single Through Going (∼82%): Single muons that pass entirely through the detector and is the default fit category for a muon. 2. Stopping (∼7%): Single muons that enter the ID but do not exit it. Identified by low light observed near the exit point of the muon and nearby OD information. 3. Multiple (∼7%): Bundles of muons passing through the detector simultaneously. Identified by light in- consistent with a single Cherenkov cone. 4. Corner Clipping (∼4%): Single through going muons found to have a track length of less than 7 m inside the ID, while also occurring near the top or bottom of the detector. To check the fitter accuracy, ∼2000 events were fit by this method and by hand. For the categories found by the fitter, ∼0.5% of single through going, ∼1.4% of cor- ner clipping, ∼13% of multiple, and ∼30% of stopping muons were found to be something else by eye-scan. The difference in the mistagged stopping muons were hand fit to be through going muons. Resolution studies found the entry point resolution to be 100 cm for all types of muons, except multiple muons with more than 3 tracks, and a directional resolution of 6◦. In this analysis, if mul- tiple tracks were fit, the principal track is used to identify subsequent events. B. Reconstruction of low energy events Events are reconstructed from their hit timing (ver- tex), hit pattern (direction), and the effective number of hits observed (energy). Events below 19.5 MeV only travel several cm within SK and are treated as a point source. The vertex is reconstructed by maximizing a likelihood dependent on the timing residuals of the hit PMTs, τi: τi = ti − tTOF − t0 (2) where ti is the time of the PMT hit, tTOF is the time of flight from the test vertex to the PMT, and t0 is the event time. This likelihood is unbinned, and therefore it cannot be used to evaluate the goodness of the vertex fit (gt). To create a goodness of fit (gt) from the timing residu- als τi the weighted sum of Gaussian functions G(τi|σ) = exp [ −0.5 (τi/σ) 2 ] is used: gt = Nhit∑ i WiG(τi|σ) (3) where σ = 5 ns is the width of the Gaussian, ob- tained by combining the single photoelectron PMT tim- ing resolution of 3 ns to the effective time smearing due to light scattering and reflection. The weights Wi are Wi = G(τi|ω)/( ∑ j G(τj |ω)) with a “weight width” of ω = 60 ns. Event direction reconstruction is a maximum likeli- hood method comparing data and MC simulation of the PMT hit pattern caused by Cherenkov cones. This like- lihood is dependent on reconstructed energy and the an- gle between the event direction and the direction to in- dividual PMTs. Using the reconstructed direction, the azimuthal symmetry of the PMT hit pattern is probed with the goodness gp, a KS test: gp = max[φunii − φdatai ]−min[φunii − φdatai ] 2π , φunii = 2π · i Nhit , (4) where φuni is the angle for evenly spaced hits around the event direction and φdata is the actual hit angle around event direction. Like gt, gp also tests the quality of the vertex reconstruction: a badly misplaced vertex often presents the direction fitter with a Cherenkov cone pat- tern appearing too small (or too big), which implies an accumulation of hit PMTs on “one side” of the best-fit direction. Finally, for the energy reconstruction, we evaluate the photons’ times of flight from the reconstructed vertex to the hit PMTs and subtract them from the measured ar- rival times. We then define the effective number of hits Neff as the maximal number of hits in a 50 ns coinci- dence window. This number is then corrected for water transparency, the angle of incidence to the PMTs and photocathode coverage, dark noise rate, PMT gain over the course of SK-IV, PMT occupancy around a hit, the PMT quantum efficiency, and the fraction of live PMTs. The energy is then calculated from a fifth order polyno- mial dependent on Neff for energies in the solar neutrino 6 Muon flux at surface Modified Gaisser ⇓ Muon travel through rock MUSIC ⇓ Hadronic shower and isotope generation FLUKA ⇓ Detector simulation SKDetSim–GEANT3 ⇓ Event reconstruction FIG. 1. Simulation steps, from the modeling of the muon flux at the surface to event reconstruction in SK. range. The energy reconstruction assumes an electron in- teraction. This is important to note as neutron captures on hydrogen have a single γ which creates less light than a 2.2 MeV electron would. Within this paper the energy of events will be given in terms of the kinetic energy of an electron with equivalent light yield. IV. SPALLATION SIMULATION A. Cosmic muon simulation In order to understand and optimize spallation event removal techniques we simulate the interactions of cos- mic muons and the subsequent production of neutrons and isotopes in the SK water. This muon simulation is composed of five parts. We first model the muon flux at the surface of the Earth using a modified Gaisser param- eterization described in Ref. [18], and propagate muons through the rock to SK using a dedicated transport sim- ulation code. Second, we simulate the production of hadronic showers and radioactive isotopes inside SK us- ing FLUKA [7, 19]. The FLUKA results are then injected into SKDetSim, the official GEANT-3 [20] based detector simulation for SK, that will model detector effects as well as minimum ionization around the muon track. Finally, we reconstruct muon tracks, neutron captures, and iso- tope decays using standard SK reconstruction software as well as the procedure described in Sec. V. The simulation pipeline is summarized in Fig. 1. 1. Muon generation and travel Most of the muons reaching the Earth’s surface are produced at an altitude of around 15 km, from the inter- actions of primary cosmic rays in the atmosphere, prin- cipally the decays of charged mesons [21]. The shape of the meson production energy and angular distributions reflects a convolution of the production spectra, the en- ergy loss, and the decay probability in the atmosphere. In this study we model the muon flux at the surface us- ing a modified Gaisser parameterization, that has been optimized for detectors at shallow depth such as SK [18]. While muons contributing most to spallation background can cross the detector either alone or as part of muon “bundles” –caused by meson decays within cosmic ray showers in the atmosphere– differences in spallation ob- servables between these two configurations will be en- tirely due to track reconstruction issues. In this paper, our primary goal is to evaluate FLUKA’s ability to model muon-induced showers and isotope production in SK. We will therefore consider only single muons and will leave the subject of bundles to a future study. In order to obtain the muon flux entering the SK detec- tor we now need to propagate muons through the rock surrounding the detector. We simulate muon propaga- tion using the MUSIC [22–24] propagation code. MU- SIC integrates models for all the different types of muon interactions with matter leading to energy losses and de- flections, such as pair production, bremsstrahlung, ion- ization and muon-nucleus inelastic scattering. Angular and lateral displacements due to multiple scattering are also taken into account. Muons are transported with energies up to 107 GeV. Here, we used the rock compo- sition model described in [18], with an average density of ρ = 2.70 g cm−3. We compute the muon travel distance within the rock as a function of the incident angle using a topological map of the SK area from 1997 [25, 26]. Figure 2 shows the energy spectrum of the muons that reach SK. The rock above the detector consti- tutes a particularly efficient shield, effectively blocking muons with energies lower than 600 GeV. The cosmic muon flux is reduced from 6.5 × 105 µm−2 h−1 [27] to 1.54× 10−7 cm−2 s−1 which corresponds to a muon rate of 1.87 Hz, as expected from previous measurements [28– 30]. Although MUSIC allows to evaluate the effect of muon transport through rock on the muon directional distribu- tion, it does not account for the detector’s cylindrical geometry. We account for these effects by assuming that cosmic muons are isotropically produced in the atmo- sphere and that all muons produced in the same area of the sky have quasi-parallel trajectories inside SK. For each muon generated by MUSIC with a given direction (θ, φ), with origin of the coordinates at the center of SK, θ representing the zenith angle and φ the azimuthal one set to zero when the final muon travels from east to west, we generate a set of parallel tracks with uniformly dis- tributed intersection points in the plane perpendicular to the (θ, φ) vector, as shown in Fig. 3. We then reject all the tracks that do not cross the detector, thus straight- forwardly accounting for geometrical effects. This pro- cedure allows us to convert the directions generated by MUSIC into a sample of entry points distributed on the surface of the SK inner detector. 7 1−10 1 10 210 310 410 510 610 Muon energy [GeV] 22−10 21−10 20−10 19−10 18−10 17−10 16−10 15−10 14−10 13−10 12−10 11−10 10−10 9−10 8−10 7−10] -1 G eV -1 s -2 M uo n flu x [c m FIG. 2. Muon energy spectrum at the location of SK detector in the mine inside Mt. Ikenoyama. (θ,φ) 2R cosθ H sinθ FIG. 3. Spatial distribution of trajectories for muons pro- duced in the same area of the sky. These muons can be con- sidered almost parallel when reaching SK and the intersection of their trajectories with a plane perpendicular to their direc- tion will be uniformly distributed. Here, R and H are the radius and the total height of SK’s inner detector while θ and φ define the direction of the muons. Here, the stars indicate the muon entry points and the crosses indicate the intersec- tions of the muon trajectories with the plane. 2. Muon interactions in water Propagation and interactions of muons in water are simulated with FLUKA, taking as input MUSIC energy and angular distributions. FLUKA [7, 19] is a general purpose Monte Carlo code for the description of interac- tions and transport of particles in matter. It simulates hadrons, ions, and electromagnetic particles, from few keV to cosmic ray energies. It is built and frequently upgraded with the aim of maintaining implementations and improvements of modern physical models. FLUKA version 2011.2x.7 is used for this work, together with FLAIR (version 2.3-0), an advanced user interface to facilitate the editing of FLUKA input files, execution of the code and visualization of the output files [31]. FLUKA propagates muons into the SK detector, sim- ulating all the relevant physics processes that lead to energy losses and creation of secondary particles: ion- ization and bremsstrahlung, gamma-ray pair production, Compton scattering and muon photonuclear interactions. Hadronic processes such as pion production and interac- tions, low energy neutron interactions with nuclei and photo-desintegration are also modeled. FLUKA code fully integrates the most relevant physics models and libraries. For this work, the simulation was built with the default setting PRECISIO(n). All the specifics related to this setting can be found in [19]. More detail about the models and settings used in this paper can be found in appendix A. In particular, low-energy neutrons, which are defined to have less than 20 MeV energy, are transported down to thermal energies, a set- ting that is critical for our study. Crucial options complement the default setting: EVAPORAT(ion) and COALESCE(nce) give a detailed treatment of nuclear de-excitations while nucleus-nucleus interactions are enabled for all energies via the option IONTRANS. The SK detector is modeled as a cylindric volume of pure water, as described in Sec. II. The PMT structure is not simulated in FLUKA, given that it is fully incorporated in SKDetSim. Since muons can induce showers outside the water tank and secondary products may reach the active part of the detector, previous studies [1] examined the effect of a 2 m thickness of rock surrounding the OD: it was proven that this has a minor effect on the results. Thus the rock, as well as the tank and the support struc- ture, are not simulated in this work. Both negative and positive muons are generated assuming a muon charge ratio, defined as the number of positive over negative charged muons, of Nµ+/Nµ− = 1.27 [32]. Note that the measured values of the charge ratio at SK depth can vary by about 20%, with the highest value (1.37±0.06) mea- sured at Kamiokande [33]. However, since the isotope yield depends only weakly on the muon charge, these variations have a negligible impact on this analysis, with an effect on the predicted yields of less than 1%. 3. Detector response and event reconstruction We model the detector response using the official de- tector simulation for Super-Kamiokande, referred to as SKDetSim. This simulation is based on GEANT 3.21 [20] for detector modeling and uses a customized model for light propagation. It covers all aspects of event detection, 10 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 t g 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1p g 0 2 4 6 8 10 12 14 16 18 "Bad" Quality Events 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 t g 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1p g 0 5 10 15 20 25 30 35 "Good" Quality Events FIG. 4. MC generated 2.2 MeV γ-rays processed by WIT software. The left distribution shows “poorly” reconstructed events as defined as reconstructed vertex being more then 5 m from MC truth and the right distribution shows the “good” reconstructions. The lines separate the different weight regions used to parameterize the neutron cloud cuts described in Sec. VIII A. tances from the candidate to the muon track, lt and lLONG. a. Time difference ∆t: the characteristic neutron capture time in pure water is τcap∼205 µs. We hence expect to find most neutrons within about 500 µs from their parent muon crossing time. Since the muon rate in SK is about 2 Hz, this time window alone allows to unambiguously link a neutron capture event to its parent muon. Additionally, for high-quality neutron cloud sam- ples, we will require ∆t to be larger than 20 µs to account for PMT afterpulsing [34] The afterpulsing features are shown in Fig. 5 for events found within 5 m of a muon track. b. Transverse and longitudinal distances, lt and lLONG: lt defines the distance between a neutron candidate and the closest point on the muon track while lLONG is the longitudinal distance, the distance between that clos- est point on the muon track and the muon track entry point. The definitions of these two observables are shown in Fig. 6. Their two-dimensional distribution, shown in Fig. 7, allows to characterize the shapes and sizes of the neutron clouds. In this figure, 〈lLONG〉 refers to the av- erage lLONG of all the neutrons in the cloud. Note that neutron clouds have an elongated shape along the muon track and can extend up to about 5 m. Here, we showed l2t instead of lt to reflect the amount of phase space avail- able. Incorporating the observables defined above in spalla- tion analyses thus allow to identify neutrons and accu- rately locate neutron clouds and the showers that gen- erated them. In what follows we will show case studies highlighting the validity of our neutron identification and localization procedure. B. Neutron cloud identification: case studies Here, we present examples of how to identify neu- trons and reconstruct clouds using the observables de- fined above. We first demonstrate our ability to iden- tify individual neutrons by associating the WIT trigger with simple goodness and position cuts. Then we use well-reconstructed neutron candidates to build neutron clouds, and exploit their spatial correlation with spalla- tion isotopes. 1. Identifying individual neutrons We build a high-purity neutron sample using WIT- triggered events verifying gt > 0.5, lt < 1.5 m, and Erec < 5 MeV. We estimate the number of neutrons in this sample by comparing the distribution of time dif- ferences between neutron candidates and their parent muons, ∆t, to the one expected from calibration stud- ies using an americium beryllium (AmBe) source [35]. This ∆t distribution is shown in Fig. 8 and was fitted from 50 µs to 500 µs by the following equation N(∆t) = A · e−∆t/τ + C (5) where N is the number of events, τ is the exponential de- cay time constant, and the constant C absorbs remain- ing background contributions. For the WIT data the 11 s]µt [∆ 0 5 10 15 20 25 30 s µ C o u n ts /0 .1 0 10000 20000 30000 40000 50000 60000 FIG. 5. Time distribution of WIT-triggered events with lt < 5 m. The goodness cuts included in the WIT trigger lead to a deficit of events in the 5–10 µs time range while the bumps are contributions from fake neutrons associated with PMT afterpulsing. In order to remove these fake neutrons, quality events are required to have ∆t > 20 µs. neutron/ spallation candidatelt (transverse distance) lLONG (longitudinal distance) θµ track FIG. 6. Diagram showing the lt and lLONG observables asso- ciated with neutron identification, as described in Sec. V A 2. time constant was measured to be τ = 211.8 ± 1.7 µs. In comparison, in AmBe calibration studies, the neu- tron capture time on hydrogen was measured to be τ = 203.7 ± 2.8 µs. This results in about a 2.5σ dif- ference between the two measurements. This discrep- ancy is believed to be due to missing neutrons within higher multiplicity showers from pile up. WIT defines a 1.5 µs window around a triggered event (500 ns before and 1,000 ns after the trigger time). If another neutron capture happens within that window, no new trigger will be issued. This introduces a form of deadtime [36] for 20− 15− 10− 5− 0 5 10 15 20 > (m) LONG ­<l LONG l 0 2 4 6 8 10 12 14 16 18 20) 2 ( m 2 l t 200 400 600800 1000 1500 2000 LONG vs l2lt FIG. 7. Data Vertex correlation of neutron capture events. The vertical axis shows the squared distance to the muon track, the vertical axis the distance along the muon track with respect to the average. Solid contours indicate levels of multi- ple 1,000 events/bin, the thick contours are for 5,000, 10,000, 15,000, and 20,000. The levels of the dotted line contours are marked in the Figure. captures close together in time, biasing for a longer cap- ture time constant. 2. Identifying neutron clouds After having successfully identified neutrons using WIT triggers, we will now investigate correlations be- tween neutron clouds and spallation isotopes. Here, we define neutron clouds as groups of two or more WIT events observed within 500 µs after a muon. Addition- ally, we require at least one of these WIT events to have ∆t > 20 µs, gt > 0.5, and E < 5 MeV. Note that fewer events reconstruct outside the ID (and even the fiducial volume) due to the WIT trigger conditions on the online event reconstruction. As detailed at the beginning of this section, these cuts allow to discriminate neutrons against noise fluctuations, radioactivity, and afterpulsing events. If a cloud was found using the conditions listed above, spallation candidates were searched for in close proxim- ity of the center of a neutron cloud. Spallation candi- date events were preselected using the noise reduction and quality cuts described in Sec. VII and in Ref. [5] for the SK-IV solar neutrino analysis —without the spalla- tion and pattern likelihood cuts. Spallation candidates found within less than 5 m and 60 s after an observed neutron cloud are selected. 60 s was chosen to contain more than 99% of the 16N decays and the 5 m value re- flects the general size of muon-induced showers, as seen in Fig. 7, In addition to this signal sample, we build a back- 12 s]µt [∆ 50 100 150 200 250 300 350 400 450 500 s µ C o u n ts /0 .1 2000 3000 4000 5000 6000 FIG. 8. ∆t (black dots) and resulting fit (grey solid line) using the function from Eq. (5), for neutron events detected after muons by WIT. The time constant obtained from the fit is τ = 211.8± 1.7 µs. ground sample using candidates found within 5 m of and up to 60 s before neutron clouds. This background sam- ple allows to estimate the fraction of spurious pairings between spallation candidates and uncorrelated neutron clouds in the signal sample, and subtract off the corre- sponding effects. A strong spatial correlation between spallation candi- date events and the centers of neutron clouds was found. As shown in Fig. 9, this correlation increases with the multiplicity (the number of WIT triggered events) of the neutron cloud. The slow decrease of the number of spalla- tion candidates relative to the others at low multiplicities is due to accidental pairings between candidates and neu- tron clouds. Conversely, high-multiplicity neutron clouds can be more easily located and allow to reliably identify spallation products. These high-multiplicity clouds are also likely to be associated with multiple isotopes in large hadronic showers, as can be observed from Fig. 10. Identifying neutron clouds and correlating them with spallation candidates using the criteria outlined above al- lows to remove 55% of spallation events with a little more than 4% deadtime. Due to the low detection efficiency for the 2.2 MeV γs resulting from neutron capture on hy- drogen, small showers are missed, and therefore this pro- cedure alone does not suffice to remove spallation in an SK analysis. In Sec. VIII we will show how to optimize neutron cloud reconstruction and associate it to tradi- tional spallation cuts for the SK solar analysis. There, in addition to the observables used for this case study, we will notably use the directional goodness gp to better identify well-reconstructed neutrons, and will account for ] 3 [m 3 Distance 0 20 40 60 80 100 120 C o u n ts /b in /s p a lla ti o n 3− 10 2−10 1−10 Neutron Multiplicities 2 neutrons 3 neutrons 4­5 neutrons 6­9 neutrons 10+ neutrons FIG. 9. Signal distribution for vertex correlation by neutron cloud multiplicity. The figure is normalized by the number of spallation events for each multiplicity (Total events in sig- nal − background distribution). As multiplicity increases, the steepness in the tail of the distribution increases as there is more accuracy in parameterization of the cloud and more spallation is expected in larger hadronic showers. Note that the vertical axis is on a logarithmic scale. Neutron Multiplicity 0 20 40 60 80 100 120 140 S p a lla ti o n M u lt ip lic it y 0 5 10 15 20 25 30 35 40 1 10 210 3 10 410 5 10 6 10 FIG. 10. Two-dimensional distribution showing the corre- lation between neutron capture candidate multiplicity and spallation candidate multiplicity. Note the color scale is on a logarithmic scale. 15 nificantly reshape the spallation reduction procedure, as gadolinium doping will sizably increase the neutron iden- tification efficiency. 0 20 40 60 80 100 120 140 160 180 200 220 240 310× ]2 [cm2lt 0 5 10 15 20 25 30 35 6−10× N or m al iz ed d is tr ib ut io ns MC Data 1500− 1000− 500− 0 500 1000 1500 [cm]ln∆ 0 0.5 1 1.5 2 2.5 3−10× N or m al iz ed d is tr ib ut io ns MC Data FIG. 12. l2t (top) and ∆lLONG (bottom) distributions for neu- trons belonging to a cluster of multiplicity larger than or equal to two, for simulations (black) and data (red). The fraction of muons that do not produce hadronic showers is set to 96.5%. The dotted lines the average transverse and longitudinal cloud extensions of 3 m and 5 m respectively. 3. Neutron multiplicity We finally estimate the number of reconstructed neu- trons associated with the muons in both the data and simulation samples. Here, as in the previous section we use an EM muon fraction of 96.5% and treat pos- sible mismodeling of this fraction as a systematic uncer- tainty. The neutron cloud multiplicities for both simula- tion and data are shown in Fig. 13. The abundance of low-multiplicity clouds is due to both fake neutron con- tributions and the low efficiency of the neutron tagging algorithm. For neutron clouds with multiplicities lower than 10, simulation and data show reasonable agreement. Conversely, for multiplicities larger than 10, FLUKA fails to accurately simulate the tails of the data distribution. Note, however, that such large clouds are typically as- sociated with shower-producing hundreds, sometimes up to thousands of neutrons. Muons associated with these high-multiplicity showers are not only rare but also de- posit a high amount of light in the detector and are hence easier to identify using other observables, such as for ex- ample the residual charge Qres that will be introduced in Sec. VIII C. In any case, neutron multiplicity distribu- tions from data will be used to improve the simulation in future. Our results hence demonstrate the ability of FLUKA-based simulations to accurately model hadronic showers for the types of muons that need most to be studied in future SK analyses. VII. SOLAR NEUTRINO ANALYSIS We will now use the insight gained through the studies described in this paper to design new spallation cuts for a specific SK-IV analysis. For this paper, we will focus on the solar neutrino measurements. Indeed, while many low energy analyses at SK target antineutrinos, whose in- teractions produce neutrons, solar neutrino interactions do not have a signature, so radioactive β decays will look similar to them (except for the direction of the electron). Hence, even with SK-Gd, the solar neutrino analysis will continue to heavily rely on dedicated spallation tagging techniques, in addition to control of intrinsic radioactive backgrounds such as from the radon decay chain [41]. In what follows we present an overview of the current analy- sis, and in particular of the previous spallation reduction strategy. A. Overview Here, we present a brief overview of the cuts used for the SK-IV solar neutrino analysis. A more detailed de- scription of these cuts can be found in Ref. [5]. 8B and hep solar neutrinos from the Sun scatter elasti- cally on electrons in SK, with the recoiling electron pro- ducing a single Cherenkov ring. The produced electrons have energies up to about 20 MeV; however, intrinsic ra- dioactivity is a limiting background that dominates below about 5.5 MeV kinetic energy. In this paper, we therefore consider the 5.49−19.5 MeV energy range. The impact of these new spallation cut son the full SK-IV solar neutrino analysis will be described in [42]. After high-efficiency noise reduction cuts are ap- plied, notably rejecting events outside the FV, energy- dependent quality cuts are imposed on the reconstructed 16 2 3 4 5 6 7 8 9 10 Neutron multiplicity 1 10 210 P er ce nt ag e MC mip+showers MC mip Data C C EM only ata 10 15 20 25 30 35 40 45 Neutron multiplicity 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 P er ce nt ag e MC mip+showers MC mip Data C C EM only ata FIG. 13. Neutron multiplicity distributions for simula- tion (black) and data (red) with the EM muon contribution shown for the simulation (filled gray), this contributes only below 3 neutrons multiplicity. We show these distributions separately for neutron multiplicities from 2 up to 10, in log scale (top) and from 11 to 45 in linear scale (bottom). vertices and directions of the events, the goodness observ- ables described in Sec. III B. The major quality cuts se- lect events based on their reconstruction goodness (gR = g2 t − g2 p) and their effective distance to the ID wall — obtained by following their reconstructed direction back- wards. We then check for consistency with a single elec- tron ring (rather than, e.g. light deposited by βγ de- cays of 16N). We parametrize the amount of multiple Coulomb scattering of the electron (the “fuzziness” of the ring) with the “multiple scattering goodness” to sep- arate lower energy β decays. Most of the remaining back- ground events are removed using a dedicated spallation cut. The number of solar neutrinos can be readily extracted from the sample of events remaining after cuts by con- sidering event directions. Indeed, electrons recoiling from neutrino elastic scattering will be almost collinear with the incoming neutrinos, that is the angle between the di- rection of a given event and the direction from the Sun at its detection time, θSun is small (less than 15◦) and the cos θSun distribution is strongly peaked around 1. The numbers of solar neutrinos and background events are extracted from a fit to the cos θSun distribution. B. Spallation cuts for the previous solar analyses Spallation backgrounds can be significantly reduced by identifying space and time correlations between each isotope decay and its production in a hadronic shower initiated by a muon. To this end, each solar neutrino candidate selected using the noise reduction and qual- ity cuts described above is paired with muons detected up to 100 s before it. For each pair, three observables are then considered: the time difference ∆t between the solar neutrino event candidate and the muon, the trans- verse distance of the candidate to the muon track lt — defined in Fig. 6— and the residual charge Qres, defined as the excess charge deposited by the muon in the detec- tor compared to the expectation from minimum ioniza- tion. For each observable, probability distribution func- tions (PDFs) are then defined for spallation pairs, formed by isotope decay events and their parent muons, and for uncorrelated “random” pairs. These PDFs then allow to define a log likelihood function log10 L, whose functional form is as follows: log10 L = log10 [∏ i PDFspall,i(xi)/PDFran,i(xi) ] (6) where PDFspall,i and PDFran,i designate the PDFs asso- ciated with a given observable i for spallation and random pairs respectively. In the absence of a spallation simulation, the PDFs for spallation and random pairs need to be extracted from data. One sample is built by pairing solar event candi- dates with preceding muons (as described above). It will contain a mixture of spallation and random coincidence pairs. Using the times of events with energies much below 6 MeV, we construct a“random sample” by generating a vertex from a uniform distribution filling the entire de- tector. When paired with preceding muons, this random sample estimates the random coincidence contribution, so the corresponding PDFs are extracted. Also, after subtraction of the random sample, the spallation PDFs are extracted from the data sample. Alternatively, we invert the time sequence and pair solar candidates with muons up to 100 s after them. This inverted sample is used the same way as the random sample. Finally, in order to account for possible correlations between observ- ables, the lt PDFs are computed for seven different Qres 17 bins. Since the muon fitter used for these cuts consid- ered only single through-going muons, a goodness-of-fit cut was also considered for this PDF in order not to be misled by poorly-fitted muon tracks. The final likelihood cut used for the solar analysis re- moved 90% of spallation events with a position-averaged 20% deadtime. This deadtime was measured with the random sample as a function of position. In the next section, we will show how the SK-IV new electronics, the new techniques described in this paper, and better muon reconstruction algorithms allow to further reduce this deadtime for the upcoming analysis. VIII. SPALLATION CUTS FOR THE SOLAR ANALYSIS Here, we present a new spallation cut that improves on the reduction strategy described in Sec. VII B. We take advantage of several improvements and studies that took place within the last decade. First, the muon track reconstruction was replaced. Previously we used a sim- ple, fast muon track fitter developed at the beginning of SK. It assumes through-going single muons and misre- constructs or fails on other muons. The more complex muon track reconstruction of this analysis categorizes as described in Sec. III A, and reconstructs all categories (up to ten tracks). It was used to reject spallation back- ground by dE/dx reconstruction in the search for dif- fuse supernova neutrino interactions in SK-I, II, III [43], which inspired the development of FLUKA-based simula- tion studies [1] and highlighted the importance of muon- induced hadronic showers for isotope identification, and allowed to characterize their shapes and sizes. Finally, the improvements in the detector electronics associated with the SK-IV phase allowed to raise the PMT satu- ration rate and detect higher values of the total charge (and therefore Qres), as well as identify neutron clouds as described in Sec. V. The new spallation reduction strategy proceeds as fol- lows. First, we apply two sets of preselection cuts in or- der to remove a sizable fraction of spallation events with minimal harm to the signal efficiency. These cuts aim at removing events close in space and time to neutron clouds, as well as clusters of low energy events, typically associated with the decays of multiple isotopes produced by the same muons. Then, we remove most of the re- maining spallation events using an updated version of the likelihood cut described in Sec. VII B. A. Neutron Cloud Spallation Cut for Solar Analysis Using the observables defined in Sec. V A, we define a set of cuts to reliably identify neutron clouds and inves- tigate their space and time correlation with solar event candidates (most of which are spallation events before cuts). First, we define neutron candidates as WIT events found less than 500 µs after a muon and within 5 m of its track. The number of these candidates gives the neutron cloud multiplicity. Then, in order to compute the cloud barycenter, we consider a high-purity subsample of these neutron candidates, requiring them to verify ∆t > 20 µs and Erec < 5 MeV. We then assign weights to these can- didates depending on their vertex and direction good- ness gt and gp. Specifically, we consider three regions of weights 0, 1, and 2 in the gt − gp space, that are shown in Fig. 4. Once clouds are identified and their barycenter is de- fined, their positions and detection times can be com- pared to the location and times of solar event candi- dates. For this analysis we consider SK-IV low energy events that passed all the cuts defined for the solar anal- ysis in Sec. VII except spallation cuts. This sample is expected to be largely dominated by spallation isotope decays. We then invert the time sequence similar to the procedure described in Sec. VII B, this time considering neutron clouds observed up to 60 s before each low energy event. In order to take advantage of the expected shower shape of neutron clouds, we then need to define a spe- cific coordinate system for each cloud. Here, we consider three possible options, shown in Fig. 14. First, the axes of the new coordinate system could align with the axes of the best-fit ellipsoid of neutron cloud. This option is not practical, however, due to the large shape uncertainties for the low multiplicity clouds. A second possibility is to use the muon track as the z axis of our coordinate system and the center of the neutron cloud as its origin. Finally, the third option also uses the muon track as a z axis but sets the projection of the cloud center on the muon track as the origin. In order to assess the discriminating power of these last two options, we compute the trans- verse and longitudinal distances of low energy events, lt and ∆lLONG = lisotope LONG − ln-cloud LONG , defined in Fig. 14, to the origins of their respective coordinate systems. The distri- bution of lt and ∆lLONG is shown in Fig. 15 for all neutron cloud multiplicities. We notice that using the projection of the neutron cloud center on the muon track as an ori- gin significantly reduces the spread of this distribution in lt, a spread that is primarily driven by contributions from low multiplicity clouds. We hence choose this def- inition of the origin and set the z axis to be along the muon track for our analysis. Using the coordinate system defined on the right panel of Fig. 14 we then define cuts on ∆t, lt and ∆lLONG for each low energy event–neutron cloud pair, with ∆t de- fined as the time difference between the low energy event and the muon associated with the cloud. We first define spherical cuts, removing events within either 0.2 s and 7.5 m or 2 s and 5 m of clouds with 2 or more neutrons. Then, we define multiplicity bins of 2, 3, 4–5, 6–9, and ≥10 neutrons candidates, and, for each bin, define a spe- cific ellipsoidal cut on lt and ln. Since clouds with only 2 candidates are often not associated with hadronic show- ers, as shown in Sec. VI, we require ∆t < 30 s. For higher 20 0 20 40 60 80 100 ]2 [m2l t 0 100 200 300 400 500 600 c o u n ts /b in 0 100 200 300 400 500 600 700 800 900 ]2 [m2l t 0 10 20 30 40 50 60 70 80 c o u n ts /b in FIG. 17. Distributions with PDF for the spallation (left, normal time sequence) and random coincidence (right, inverted time sequence) samples with corresponding fits (solid lines). The distributions shown here are for the 3–30 s ∆t and 0.5–1.0 Mpe Qres bin. The analytical forms for the fits are shown in Equations 8 and 9. cle (MIP) traveling the same distance inside the detector. The MIP muon charge per unit of track length is defined as the peak value of the distribution of the amount of charge deposited by unit track length for single through- going muons. It is evaluated for each run time period and typically lies around 26.78 photo-electrons (p.e.) per cm. We obtain the PDFs for this observable by using a sam- ple of low energy events within 2 m and 10 s of a muon. The distributions for spallation and random uncorrelated pairs are shown in Fig. 18. Both signal and background PDFs in the positive Qres region are a sum of exponential functions: PDF(Qres) = 5∑ i=1 eci−pi·Qres , Qres > 0 (10) where ci and pi are the fit parameters for exponential functions. For negative Qres, no analytical form was as- sumed and a linear interpolation of the sample bins was used. 4. Longitudinal distance ∆lLONG Defining PDFs for the difference in longitudinal length ∆lLONG = lisotope LONG − ldE/dx peak LONG allows to include informa- tion about the muon-induced hadronic shower into the likelihood cut. Here, we developed a new method to re- construct the energy loss of the muon along its track (dE/dx) based on a previously published method [43]. Defects and possible improvements to this method were Residual Charge [p.e.] 0 1000 2000 3000 4000 5000 3 10× N o rm a liz e d C o u n ts /p .e . 6−10 5−10 4−10 3−10 2−10 1−10 FIG. 18. Signal (red triangles) and Background (blue circles) distributions for the sample used for the Qres PDF construc- tion. suggested in Ref.[3]. The method presented here reme- dies those defects, although the improvements differ from those suggested by Ref. [3]. Using the entry time of the muon in the detector and the PMT hit pattern, dE/dx is estimated by identifying the points of the muon track 21 verifying the time correlation equation for each PMT hit: tPMT − tentry = l · cvac + d · cwater (11) where tPMT and tentry are the PMT time and muon entry time respectively, l is the distance from the muon entry point to the point along the track where the light is emit- ted from, d is the distance from the emission point to the PMT, and cvac and cwater are the speeds of light in vac- uum and water respectively. The dE/dx is computed for 50 cm segments of the muon track, corresponding roughly to the vertex resolution for events in the energy region of 3.49–19.5 MeV in SK. The simplest approach to esti- mating dE/dx is to add the charge of each PMT to each bin containing a solution of Eq. (11). Here, using the method proposed in [1], we spread the charge from each hit across multiple bins to account for the PMT timing resolution. More specifically, we take the contribution, gij , of the ith PMT to the jth bin to be: gij = Qi · ed(lj+25)/λ S(θij , φij) · fij∑ k fik (12) where fij = ∣∣∣∣Erf ( τ(lj)− ti√ 2σi ) − Erf ( τ(lj + 50)− ti√ 2σi )∣∣∣∣ , (13) σi is the timing resolution for the observed charge by the ith PMT, τ(lj) is the tPMT that solves Eq. (11). For the jth bin boundary’s d and l, Erf is the standard error function, Qj is the charge observed by the PMT, the exponential function is water attenuation correction, S is the photocathode coverage correction, and the integral of the sum is normalized to one. This procedure ensures that the charge of a given PMT hit is only counted once. Although error functions are used for the integral, since τj(l) is non-linear, the integral is not easily normalized. Care also has to be taken for the shape of τj(l) as it is not monotonic, therefore special cases are implemented to handle scenarios where τ(l) − ti = 0 and when dτ/dl = 0. For muons that induce particle showers in addition to minimum ionization, the segment of the muon track associated with the largest dE/dx can indicate the lo- cation of these showers. We hence define ∆lLONG = lisotope LONG − l dE/dx peak LONG as the longitudinal distance of an iso- tope to this segment along the muon track. Distribu- tions of this observable for spallation and uncorrelated pairs are shown in Fig. 19. To define the correspond- ing PDFs, the ∆lLONG distribution for low energy events found within 2 m and 10 s of a muon is fitted by a sum of three Gaussians: PDF(∆llong) = 3∑ i=1 Aie −(llong−xi) 2 2σ2 i (14) where Ai, x, and σi are the fit parameters for each PDF. For uncorrelated pairs, one of the Gaussian fits was de- generate and dropped from the final form. For minimum [m]LONG∆l 40− 30− 20− 10− 0 10 20 30 40 F ra c ti o n 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 FIG. 19. ∆lLONG distribution used for PDF fit, for the spal- lation signal (blue triangles) and for spallation accidentals (black circles). The peak of the background distribution is shifted to negative ∆lLONG as a result of non-showering muons being more likely to have a dE/dx peak later in the track. ionizing muons, the dE/dx peak is more likely to be at the end of the track, resulting in the background distri- bution being slightly shifted away from 0. Using the PDFs defined above, we define a log likeli- hood function as shown in Eq. (6). The distributions of this log likelihood for signal and background are shown in Fig. 20. To estimate the impact of the cut on this function, we take advantage of the fact that the low en- ergy event sample that we are considering is dominated by spallation and solar events. We can hence readily es- timate the background rejection rate of our algorithm by computing the fraction of events with cos θsun < 0 removed by the likelihood cuts. Conversely, the signal efficiency can be computed by applying cuts on a ran- dom sample, where low energy events are paired with muons observed after them. We use these techniques to tune the cut point for the log likelihood, maximizing the signal efficiency for a background rejection rate of 90%, similar to the one obtained with the previous spallation cut described in Sec. VII B. Since WIT was running only during a small fraction of the SK-IV period, we apply different likelihood cuts depending on whether the neu- tron cloud information is available. The availability of neutron cloud information notably allows to loosen the likelihood cut. 22 (Likelihood) 10 log 0 1 2 3 4 5 6 7 8 9 10 C o u n ts /b in 0 2000 4000 6000 8000 10000 12000 FIG. 20. Comparison of the logL for the signal and back- ground of the non-neutron data period. This figure shows the distributions for the spallation signal (dashed) and spallation accidentals (solid). The vertical dashed line shows the tuned cut value. Since the multiple spallation cut is already applied, only 82% of remaining spallation is removed to achieve 90% overall spallation removal effectiveness. D. Total Spallation Cut Results After tuning the neutron cloud, multiple spallation, and likelihood cuts described throughout this section, we apply spallation reduction to a sample of SK-IV events passing the noise reduction, quality cuts, and pattern likelihood cut described in Sec. VII. The total dead times for the periods with and without neutron cloud information are 8.9% and 10.8%. Figure 21 displays the position dependence of the dead time with neutron cloud information. The new spallation cut hence allows to reduce the dead time by up to 55% compared to the previous analysis. The effect of this new reduction on the SK-IV solar analysis is shown in Fig. 22, that shows the cos θsun dis- tribution for events passing the new spallation cut and failing the one described in Sec. VII. The clear peak around cos θsun = 1 shows that the new procedure allows to retrieve a sizable number of solar events. In the final sample, this cut allows for an increase of 12.6% solar neu- trino signal events, with a reduction in the relative error on the number of solar events of 6.6%. Compared to the total SK-IV exposure, retrieving this signal corresponds to an increase of roughly a year of detector running. 0 500 1000 1500 2000 2500 3 10× 2 in cm 2 radius 1500− 1000− 500− 0 500 1000 1500 z i n c m 11% 10% 9% 8% 5% FIG. 21. Spallation cut dead time as a function of radius and height. The origin is taken at the center of the detector [44]. This function combines neutron cloud, multiple spallation and likelihood dead time. The fiducial volume is indicated by the shaded area. sun θcos 1− 0.8− 0.6− 0.4− 0.2− 0 0.2 0.4 0.6 0.8 1 c o u n ts /b in 100− 0 100 200 300 400 500 600 FIG. 22. Figure showing the difference in events between fi- nal solar samples using new spallation cut compared to the previous cut for 5.99 to 19.5 MeV. The peak at the right con- tains the additional solar neutrino interactions gained with the change in cut. The small excess in the flat region reflects changes to the other solar neutrino cuts and differing interac- tion with the previous and new spallation cut. The intrinsic radioactivity in the energy region below 6 MeV was rejected by the previous spallation cut due to accidental coincidence. 25 TABLE VI. This table summarizes the proportion of total spallation production for the different muon categories and their respective efficiencies. The proportion of spallation is measured as the ratio of the integral of the category fits at 2000 cm to the sum of all categories. Category Endpoint [cm]Efficiency [%]Proportion [%] Single 500 92.3 54.6 Multiple 3000 49.3 43.5 Stopping 400 93.4 1.9 Corner Clipping N/A N/A < 10−3 All N/A 73.6 100 implemented to validate: After applying a pre-cut on ∆t the efficiency is the ratio of the number of events within the lt cut value over all events passing the pre- cut. Eight different pre-cut values ranging from 10 ms to 30 s were chosen. To estimate these efficiencies, a background sample was made using muons found after spallation candidates, as described in Sec. VIII C, and the ∆t distribution from this sample was subtracted from the spallation sample distribution. This procedure yields an lt efficiency of 74.2± 0.5%. The total lt efficiency is then taken to be 74.0± 0.7% to account for the discrepancies between the two methods. The lt efficiency for single through-going muons obtained from data was compared to the average efficiency from the MC simulation and was found to be about 1σ away. Finally, the isotope-dependence of the lt distribution was included to the measurement by scaling each iso- tope efficiency from MC by the weighted average of all isotopes. This factor was then used to scale the single through going measurement from data, and half of its effect was used to scale the multiple and stopping muon case. Since the multiple muon case was not performed in MC, this allowed for the difference to full or no isotope dependence to be covered in an isotope dependent sys- tematic error. This error was relatively small for most isotopes, with only 15C and 11Be having an effect greater than the full lt systematic error. Here, using the simula- tion allows to refine the lt efficiency estimate performed in [4], where the isotope dependence was covered by a ∼ 4% systematic uncertainties. D. Rate and Yields Calculation To calculate the total rates of the individual isotopes, the lt efficiency obtained from the data is combined with the efficiency associated with the noise, quality, and en- ergy cuts. Raw rates for each isotope obtained from the fits in Sec IX B are then corrected by these efficiencies to obtain the total production rates at SK-IV. For the iso- topes that were paired together for the ∆t fit, the con- tributions of the isotopes were varied from 0 to 1 and used as a systematic error. These results are shown in Table VII. The rates extracted from the data are compared to the FLUKA-based simulation described in Sec. IV. A simu- lated spallation sample of 1.362 × 108 initial muons is generated in order to accumulate enough statistics for the low yield isotopes. The predicted rates are shown alongside the observed results in Table VII. Finally, the isotope yields are obtained by rescaling the total rates computed above for each isotope as follows: Yi = Ri · FV Rµ · ρ · Lµ (16) where ρ is the density of water, Ri is the total rate of the ith isotope at SK-IV, Rµ is the muon rate (2.00 Hz), Lµ is the average length of reconstructed muon tracks, and FV is the fiducial volume of the detector. Table VII shows the final isotope yields for the data. E. Isotope study with neutron clouds In addition to updating the study performed in [4], we investigate the impact of neutron cloud cuts on iso- tope rates. Here, as stated in Sec. IX A, we consider a sample of spallation candidates paired with muons asso- ciated with at least three tagged neutrons. The fitting procedure described in Sec. IX B is then performed for all isotopes, giving a χ2/dof of 252.6/243, which corresponds to a p-value of 0.323. Then, the final rates are scaled to account for the lower live time of the WIT trigger and allow a comparison with the lt < 200 cm sample. The scaled rates without efficiency corrections (also called raw rates) are shown in Table VII for both the neutron cloud and the lt sample. As shown in Table VII, the fitted rates for the neu- tron cloud sample range between 33% and 193% of the rates found for the lt sample. The largest discrepan- cies between the two rates are seen for 8He/9C and 9Li. For these subdominant isotopes, however, the precision of the fit is limited by statistics. Moreover, these isotopes have similar half-lives, which leads to degeneracies; the rates for the 8He/9C and for 9Li are in fact found to be 98.3% anti-correlated when considering the covariance of the fit parameters. For the most abundant isotopes, on the other hand, the ratio between the rates in the neu- tron cloud and lt sample remains around 60–70%. The stability of this ratio highlights the correlation between the neutron cloud and lt cuts, as neutron cloud cuts also make use of the isotope distance to the muon track. X. CONCLUSION In this paper, we presented new techniques to reduce spallation backgrounds for low energy analyses at SK, as well as the first realistic spallation simulation in the de- tector. We notably developed algorithms locating muon- induced hadronic showers, both by improving the recon- struction of the energy deposited along the muon track, 26 TABLE VII. Observed and calculated spallation rates and yields using the lt cut. For the two sets of isotopes that could not be separated, the systematic error covers the range of yields corresponding to changing the relative fraction of an isotope from 0 to 1. The upper limit uses a 90% confidence level (C.L.) using the positive systematic error. Compared to the previous results, discovery for 15C has been made, and much stronger constraint on the 11Be measurement has been made. 11Be was at a 1.5σ excess. The 9C/8He and 9Li fits for neutron clouds were 98.2% anti-correlated. The total number of events for the two fit contributions is 65% the corresponding l2t contributions. The calculated relative fractions of 8Li and 8B in the 8Li + 8B sample are 70.1% and 29.9% respectively. The calculated relative fractions of 9C and 8He in the 9C + 8He sample are 78.8% and 21.2% respectively. Yields [10−7cm2µ−1g−1] Rates [kton−1day−1] Neutron data Fraction of lt data Calculated Isotope FractionsIsotope Calculated Observed Total Rate (lt data) Raw Rate (lt data) Raw Rate (neutron data) 12N 0.92 1.72 3.04± 0.06± 0.028 1.55 1.08 70% 2.3% 12B 8.6 12.9 22.86± 0.11± 0.21 9.19 5.95 65% 21.1% 9C/8He 0.8 <0.61 <1.08 0.11 0.20 176% – 9Li 1.5 0.67 1.19± 0.33± 0.010 0.39 0.13 34% 3.7% 8Li/8B 13.4 5.11 9.04± 0.17+0.60 −1.1 3.69 2.59 70% 32.8% 15C 0.55 1.57 2.78± 0.45± 0.032 0.76 0.37 49% 1.3% 16N 14.5 27.3 48.43± 0.60± 0.49 19.64 12.01 61% 35.3% 11Be 0.61 <1.05 <1.9 0.33 0.19 56% – and by identifying neutrons using a recently-implemented low energy trigger. New spallation cuts based on these algorithms allow to reduce the deadtime of the solar neu- trino analysis by a factor of two, allowing a gain of the equivalent of one year of exposure at SK-IV. Moreover, the profiles of the neutron clouds produced in muon- induced hadronic showers are well reproduced by the spallation simulation, motivating its use to develop spal- lation reduction algorithms for future analyses. In addition to developing new spallation cuts, we com- puted the yields of the most abundant spallation isotopes at SK-IV, updating the study presented in [4] with a 50% increase in exposure. For the isotopes with the highest production rates, the yields can be determined with a precision of a few percent. Overall, the yields predicted by our spallation simulation lie within a factor of two of the observed values, well within the uncertainties asso- ciated with hadron production models. This study also demonstrated that identifying neutron captures associ- ated with muons allows to build spallation-rich samples while keeping the relative fractions of the most abundant isotopes stable. A central piece of the spallation studies described in this paper is the identification of neutrons produced in muon-induced showers. At SK-IV, the performance of our neutron tagging algorithm has been limited by the low livetime of the associated trigger, and the weakness of the neutron capture signal. At SK-Gd, however, the neutron capture visibility will be significantly increased due to gadolinium doping. Hence, algorithms based on neutron clouds will become a key component of the up- coming spallation reduction algorithms. In this context, the simulation presented in this paper will be instrumen- tal in designing future analysis strategies. This paper thus demonstrates that, beyond neutrino-antineutrino discrimination, neutron tagging will impact significantly all low energy neutrino searches at SK. ACKNOWLEDGMENTS We gratefully acknowledge the cooperation of the Kamioka Mining and Smelting Company. The Super- Kamiokande experiment has been built and operated from funding by the Japanese Ministry of Education, Culture, Sports, Science and Technology, the U.S. De- partment of Energy, and the U.S. National Science Foun- dation. Some of us have been supported by funds from the National Research Foundation of Korea NRF-2009- 0083526 (KNRC) funded by the Ministry of Science, ICT, and Future Planning and the Ministry of Educa- tion (2018R1D1A3B07050696, 2018R1D1A1B07049158), the Japan Society for the Promotion of Science, the Na- tional Natural Science Foundation of China under Grants No. 11620101004, the Spanish Ministry of Science, Uni- versities and Innovation (grant PGC2018-099388-B-I00), the Natural Sciences and Engineering Research Council (NSERC) of Canada, the Scinet and Westgrid consor- tia of Compute Canada, the National Science Centre, Poland (2015/18/E/ST2/00758), the Science and Tech- nology Facilities Council (STFC) and GridPPP, UK, the European Union’s Horizon 2020 Research and Innova- tion Programmeunder the Marie Sklodowska-Curie grant agreement no. 754496, H2020-MSCA-RISE-2018 JEN- NIFER2 grant agreement no.822070, and H2020-MSCA- RISE-2019 SK2HK grant agreement no. 872549. Appendix A: Simulation settings In this appendix the main settings chosen to build the FLUKA simulation are described in detail. FLUKA code 27 fully integrates the most relevant physics models and li- braries; it is not possible for the user to modify or adjust them according to their needs. Several default settings are available and must be chosen at the beginning of the simulation depending on the general physics problem the user is dealing with. In addition to this, FLUKA offers several options to customize the default settings enabling or disabling a certain type of processes or changing the treatment of specific type of interactions. For this work, FLUKA simulation was built with the default setting PRECISIO(n) [19]. All the specifics that are particularly important for the scope of this paper are summarized be- low. Low-energy neutron, which are defined to have less than 20 MeV energy, are transported down to thermal energies. The absorption is fully analogue for low energy neu- trons: in a fully analogue run, each interaction is simu- lated by sampling each exclusive reaction channel with its actual physical probability, this allows for event-by- event analysis. In general, this is not always the case, especially concerning low-energy neutron interactions. Muon photonuclear interactions are activated with ex- plicit generation of secondaries. Several options are used to complement the default setting. PHOTONUC option: photon and electron interactions with nuclei are activated at all energies. MUPHOTON option: controls the full simulation of muon nuclear interactions at all energies and the pro- duction of secondary hadrons. EVAPORAT(ion) and COALESCE(nce) options: these two are set to give a more detailed treatment of nuclear de-excitations. Despite the related large CPU penalty, it is fundamental to activate these options when isotope production want to be studied. EVAPORAT en- ables the production of heavy nuclear fragments (A > 1) while COALESCE sets the emission of energetic light- fragments. IONSPLIT option: used for activating ion splitting into nucleons. IONTRANS option: full transport of all light and heavy ions and activation of nuclear interactions. RADDECAY option: activate radioactive decay calcu- lations. Settings are specified in the main input file: FLUKA, unlike other Monte Carlo particle transport codes, is built to get the basic running conditions from a sin- gle standard code. However, due to the complexity of spallation mechanism, standard options do not satisfy to retrieve the problem-specific informations we need to score: customized input and output routines (SOURCE and MGDRAW) are required to be written in order to incorporate non-standard primary particle distributions, the ones calculated with MUSIC simulation, and to ex- tract event-by-event informations for the shower recon- struction. In particular, only primary muons inducing the production of at least one hadron or of an isotope are selected and recorded; the rest are not interesting for this study and are discarded to save computational time. Appendix B: Spallation Likelihood Fit Parameters In the following section of the appendix, all of the fit parameters for the different components of the log10 L will be included. The fit to the random coincidence sam- ple is referred to as “random”, while “spallation” means the fit to the random coincidence-subtracted spallation sample. For contributing variables with only one spalla- tion and random coincidence function (∆lLONG, ∆t, and Qres) there is no normalization factor while for lt2 a nor- malization factor is needed to due to the multiple bins for the function. The normalization factors are not listed. 1. ∆t There is only a spallation function as follows as the constant fit for random coincidences is dropped: PDFsig(∆t) = 7∑ i Aie −∆t/τi (B1) where τi is the decay constant for the isotope and Ai is the fitted amplitude. The fit parameters are listed in Table VIII. TABLE VIII. ∆t fit parameters Exponential i Ai τi [s] Isotopes 1 23500 0.0159 12N, 13O, 11Li 2 83250 0.02943 12B, 13B, 14B 3 234.7 0.2568 9Li, 9C 4 869.2 1.212 8Li, 8B, 16C 5 93.37 3.533 15C 6 468.6 10.29 16N 7 5.400 19.91 11Be 2. Transverse Distance The fit for l2t is carried out over 7 residual charge bins and 3 time ranges, corresponding to 21 different fits. The equations for the spallation and random coincidence fits are as follows: PDFspa,lt(l 2 t ) = 3∑ i=1 eci−pi·l 2 t (B2) PDFrnd,lt(l 2 t ) = { p0 l2t ≤ l2t0 p0e −p1(l2t−l 2 t0)+p2(l2t−l 2 t0)2 l2t > l2t0 (B3)