New methods and simulations for cosmogenic induced spallation removalin Super-Kamiokande-IV

I. INTRODUCTION

Spallation from cosmic-ray muons produces radioactive 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, cosmic ray muons and the showers they produce sometimes interact with 16 O nuclei within the detector volume, producing radioactive isotopes.

Showers induced by the muons are primarily electromagnetic in nature (γ rays and electrons) as a result of deltaray 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 isotopes can be found up to several meters away from the muon track [1][2][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 interval 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 challenging.

In previous SK analyses, spallation reduction algorithms characterized muon signatures solely by considering their reconstructed tracks and prompt light deposition patterns. This method had important limitations due to its reliance on the muon track reconstruction quality. Moreover, while the muon prompt light contains information about produced showers, the shower signatures are partially obscured by muon Cherenkov light. As a consequence, cylindrical cuts around the entire muon track are often necessary.

For the solar neutrino analysis in SK, applying a likelihood cut based on time difference, distance to the muon track, and muon light yield removed 20% of the signal while rejecting 90% of the background in the 6.0-19.5 MeV kinetic energy range. The remaining background is dominated by decays of 16 N, which is not only the most abundantly produced isotope [4], but also is particularly difficult to identify [5].

N is primarily produced through ðn; pÞ interactions on 16 O involving neutrons from muon-induced hadronic showers. Since such neutrons can reach GeV-scale energies, 16 N can be found up to several meters away from the muon track. This large distance, together with the long halflife of this isotope, 7.3 s, makes it particularly difficult to correlate 16 N decays with their parent muon. Moreover, these decays occur either through the βγ (66%) or the β (28%) channel, producing particles with energies ranging between 3.8 and 10. 4 MeV, well within the solar neutrino energy range. Hence, removing 16 N using only muon track information is particularly difficult and results in significant reduction in neutrino signal efficiency.

Showers producing 16 N typically contain many neutrons (Oð100Þ) [6] which capture on hydrogen (H) after thermalizing, as follows:

A single 2.2 MeV γ is difficult to see in SK, but the large neutron multiplicity makes it possible to directly tag the showers and use them for spallation identification purposes.

In this paper, we develop a framework to characterize muon-induced spallation processes. In particular we describe new methods to improve spallation identification at SK by tagging the hadronic shower components, identifying clusters of spallation isotopes, and expanding the previously developed spallation cut. Additionally we present a complete simulation of cosmic muon spallation in SK, inspired by the FLUKA simulations [7,8] developed by [1].

The characteristics of the SK detector and its trigger system are described in Sec. II, and the reconstruction algorithms targeting muons and low-energy events are explained in Sec. III. Then, in Sec. IV we present a FLUKA-based simulation framework modeling muon propagation and shower generation in water. In Sec. V we then show how to take advantage of a recently deployed trigger system at SK to characterize muon-induced hadronic showers and use them to reduce spallation backgrounds. We demonstrate the ability of our simulation to accurately model the shapes and sizes of these hadronic showers in Sec. VI. Finally, in Sec. VII and Sec. VIII, we discuss how to use the insight gained through our approach to design new spallation reduction tools for the solar neutrino analysis. We first describe the solar neutrino analysis strategy and the spallation cut used for previous searches. Then we propose an update of this spallation cut and discuss how this new approach leads to a significant increase in the signal acceptance. Although we focus on the solar neutrino analysis for this paper, the techniques we present here can be readily applied to a wide range of lowenergy neutrino searches, targeting, e.g., astrophysical transients and the diffuse supernova neutrino background. We measure the yields of spallation isotopes in Sec. IX and demonstrate the FLUKA-based simulation's ability to successfully predict these yields within hadronic model uncertainties. Sec. X discusses the paper's conclusions.

II. SUPER-KAMIOKANDE EXPERIMENT

Located in Gifu Prefecture, Japan, the SK experiment is a 50 kton cylindrical water Cherenkov detector, with 2700 m water-equivalent overburden. Major changes in the detector or its electronics define the phases of the experiment [9]. This analysis is based on the SK-IV data-taking phase, the longest running phase of SK (August 2008-May 2018). The experiment is composed of an optically separated inner detector (ID) with 11,129 20-inch photomultiplier tubes (PMTs) and an outer detector (OD) with 1885 8-inch PMTs from Hamamatsu, with the ID having a height of 36.2 m and radius of 16.9 m (32.5 kton). The OD acts as a buffer for backgrounds emanating from the surrounding rock and a cosmic-ray veto. The information from the OD is sometimes used to fit and categorize cosmic rays as well. Inside the ID a fiducial volume (FV) is defined by selecting events with vertices located at least 2 m away from the ID boundary in order to reduce backgrounds arising from radioactivity in the surrounding rock, the PMTs, and their stainless steel support structure as well as the tank.

The detector response has been modeled using various calibration procedures, a detailed description of which can be found in Ref. [10]. The two sources most relevant for lowenergy analyses are a linear accelerator (LINAC) and a deuterium-tritium (DT) neutron generator. The LINAC calibration checks the absolute energy scale as well as vertex and direction reconstruction by injecting monoenergetic electrons into SK at various locations, as described in Ref. [11]. The DT generator produces 14 MeV neutrons from DT fusion; they in turn produce 16 N isotopes from 16 O via ðn; pÞ reactions, also described in Ref. [11]. The 16 N decays are used to check the absolute energy scale. The DT calibration is much easier to perform compared to LINAC calibration, and therefore the DT calibrations are done at regular time intervals while LINAC calibrations are only done once every few years. Also, LINAC electrons are limited to the downward direction while the 16 N decays isotropically.

SK-IV introduced electronics using charge to time converters and ethernet readout (QTC based electronics with Ethernet or QBBEE) and new data acquisition (DAQ) computers allowing for enhanced data processing capability [12]. The new DAQ substantially increases the bandwidth of the information collected from the PMTs: each PMT triggers by itself and then its integrated charge and trigger time are digitized and processed. The QBEEs are equipped with three different amplifiers for the charge digitization. These amplifiers enhance the dynamic charge detection range, allowing a ∼5 times increase in the maximum detected charge for an individual PMT.

Owing to these upgrades the former hardware trigger is replaced by software triggers: Data is acquired in 17 μs segments controlled by a 60 kHz clock. The 17 μs segments, still referred to as hardware triggers, connect seamlessly to each other. The standard software trigger applies a simple coincidence criterion: at least 31 triggered PMTs (or hits) within about 200 ns of each other. The dark noise rate per 200 ns of all PMTs is about 11, so a coincident signal from about 20 PMTs is needed. PMT data are saved in a window starting from 500 ns before and 1000 ns after trigger time. If the coincidence is larger than 47 triggered PMTs, the trigger window is expanded to stretch from 5 μs before and 35 μs after the trigger time. While the standard SK trigger is well-suited for events with kinetic energies down to ∼3.5 MeV, it is not sensitive to neutron captures on hydrogen which are the key to identifying hadronic showers in SK-IV and before.

A neutron capture on hydrogen emits a single 2.2 MeV γ ray. Such events produce little Cherenkov light; in SK-IV 2.2 MeV γs result in only about seven detected photoelectrons on average, so the standard software trigger efficiency is very small. After an event of at least ∼8 MeV electron-equivalent energy, the standard software trigger will automatically issue an after-trigger (AFT) that record all hit PMTs within the next 500 μs. A cosmic muon easily fulfills this condition, and the AFT triggers would catch most 2.2 MeV γ rays from subsequent neutron captures, although reliable identification of the 2.2 MeV γ signal is possible for only ∼20% of them. However, to not unduly strain the standard software trigger, these AFT triggers are disabled after cosmic muons and thus do not allow one to identify neutrons from muon-induced showers.

The wideband intelligent trigger (WIT) receives the full copy of all PMT data before triggering and runs in parallel to the standard software trigger [13,14]. It is designed to trigger on electrons of at least 2.5 MeV kinetic energy and has some capability to trigger 2.2 MeV γ rays without relying on a previously triggered event. WIT is implemented on a computer farm consisting of ten machines connected via 10 GBit Ethernet lines. Seven machines are dual-CPU with eight cores per CPU (16 hyperthreaded cores per CPU) while three newer machines have 28 cores per CPU (56 hyperthreaded cores). WIT receives blocks of 1344 consecutive hardware triggers (∼23 ms). These blocks overlap by 64 hardware triggers (1.1 ms) and are distributed among nine online machines where they are processed separately in one of the hyperthreaded cores. The remaining one of the ten computers sorts all processed files and assembles "subruns": 4,000 consecutive blocks gathered into one file (∼90 s of data). All data processing and sorting is handled without disk write operations; fast RAM memory of the WIT computers replaces the usual hard drives.

The trigger criterion of WIT is more complex than the standard software trigger as coincidence is applied to PMT hit time residuals Δt i with respect to a list of possible vertices v α : Δt i ¼ t i -j v αp i j=c watert 0 with where for the i th hit p i is the PMT hit vector, t i is the PMT hit time, t 0 is the light emission time, and c water is the speed of light in water. The applied trigger condition is

The time uncertainty is chosen to be σ ¼ 5 ns. "Max" refers to the list of possible vertices. In order to create this list, four hit combinations among a set of selected hit PMTs are chosen, and these four hit time residuals are required to be exactly zero in order to define a possible vertex. To improve the speed of the algorithm, a time-varying coincidence condition of more than 11 hits above the current dark noise rate estimate within an absolute time window of 230 ns is applied, and the hits in that window are required to obey relations δx ij > c water δt ij where δt ij (δx ij ) is the time (spatial) difference between hit i and j. The raw PMT times and pulse heights are converted from the digitized counts to calibrated times and photoelectrons in real time.

If a trigger is found, a fast vertex fit to the set of hits used for the construction of four-hit combinations is done. Only events reconstructed further than 1.5 m from any PMT are passed to the online version of the standard SK vertex fitter. For that, a 1.5 μs window is formed around the trigger time (500 ns before and 1000 ns after). If the event reconstructs at least 2 m from any PMT, and if the number of hits with time residuals between -6 ns and þ12 ns is larger than ten, the event is saved. The trigger efficiency of 2.2 MeV γ rays is 13% averaged over the entire detector and 17% averaged over the fiducial volume.

The 1.5 μs events are stored in ROOT format output files [15] containing raw PMT time and pulse height digitized counts, calibrated PMT times and pulse heights, the trigger time and position, the reconstructed vertices by both vertex fits, and the number of PMT hits within 18 ns. Also stored are GPS time stamps, pedestal measurements, and calibration measurements. Storing these low-energy events will allow the resolution of the structure of muoninduced hadronic showers and the design of an efficient spallation reduction strategy at SK-IV. Although WIT records (but not reconstructs) muons as well, we correlate WIT-triggered neutron capture candidates with muons triggered by the standard software trigger. This correlation is done after entire data taking runs are processed completely.

A. Cosmic rays and track fits

Muons pass through SK at a rate of approximately 2 Hz and their tracks are reconstructed from the PMT hits within the ID. The muon reconstruction used for this analysis, outlined in Ref. [16,17], accurately fits tracks as well as categorizes muons. After removing PMT noise hits, the fitter makes an initial guess on the track using the earliest hit PMT with at least three neighboring hits as an entry point and time, and the largest cluster of charge as the exit point. The track parameters 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 characteristics 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. This 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 projected exit point of the muon and nearby OD information.

(3) Multiple (∼7%): Bundles of muons passing through the detector simultaneously, identified by light inconsistent 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 corner clipping, ∼13% of multiple, and ∼30% of stopping muons were found to be something else by eye scan. Almost all of the ∼30% mistagged stopping muons were identified as through going muons by the eye scan. Resolution studies found the entry point resolution to be 100 cm for all types of muons, except multiple muons with more than three tracks, and a directional resolution of 6° [16]. In this analysis, if multiple 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 (vertex), 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 :

where t i is the time of the PMT hit, t TOF is the time of flight from the test vertex to the PMT, and t 0 is the event time. This likelihood is unbinned and therefore it cannot be used to evaluate the goodness of the vertex fit (g t ).

To create a goodness of fit (g t ) from the timing residuals τ i the weighted sum of Gaussian functions Gðτ i jσÞ ¼ exp ½-0.5ðτ i =σÞ 2 is used:

where σ ¼ 5 ns is the width of the Gaussian, obtained by combining the single photoelectron PMT timing resolution of 3 ns to the effective time smearing due to light scattering and reflection. The weights W i are W i ¼ Gðτ i jωÞ= ð P j Gðτ j jωÞÞ with a "weight width" of ω ¼ 60 ns. Event direction reconstruction is a maximum likelihood method comparing data and MC simulation of the PMT hit pattern caused by Cherenkov cones. This likelihood is dependent on reconstructed energy and the angle between the event direction and the direction to individual PMTs. Using the reconstructed direction, the azimuthal symmetry of the PMT hit pattern is probed with the goodness g p , a Komolgorov-Smirnov (KS) test:

;

where ϕ uni is the angle for evenly spaced hits around the event direction and ϕ data is the actual hit angle around event direction. Like g t , g p also tests the quality of the vertex reconstruction: a badly misplaced vertex often presents the direction fitter with a Cherenkov cone pattern 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 arrival times. We then define the effective number of hits N eff as the maximal number of hits in a 50 ns coincidence window. This number is then corrected for water transparency, the angle of incidence to the PMTs, 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 polynomial dependent on N eff for energies in the solar neutrino range. The energy reconstruction assumes an electron interaction. 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.

Muon generation and travel

Most of the muons reaching the Earth's surface are produced at an altitude of around 15 km from the interactions of primary cosmic rays in the atmosphere, principally the decays of charged mesons [21]. The shapes of the meson production energy and angular distributions reflect a convolution of the production spectra, the energy loss, and the decay probability in the atmosphere. In this study we model the muon flux at the surface using a modified Gaisser parametrization optimized for detectors at shallow depth such as SK [18]. 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 observables between these two configurations will be entirely 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 detector we now need to propagate muons through the rock surrounding the detector. We simulate muon propagation using the MUSIC [22][23][24] propagation code. MUSIC integrates models for all the different types of muon interactions with matter leading to energy losses and deflections, such as pair production, bremsstrahlung, ionization, 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 10 7 GeV. Here, we used the rock composition model described in Ref. [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 constitutes a particularly efficient shield, effectively blocking muons with energies lower than 600 GeV. The cosmic muon flux is reduced from 6.5×10 5 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][29][30].

Although MUSIC allows one to evaluate the effect of muon transport through rock on the muon directional distribution, it does not account for the detector's cylindrical geometry. We account for these effects by assuming that cosmic muons arriving from a given direction are uniformly distributed across the detector's cross sectional area and that their tracks inside SK are quasiparallel. 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 distributed 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 straightforwardly accounting for geometrical effects. This procedure allows us to convert the directions generated by MUSIC into a sample of entry points distributed on the surface of the SK inner detector.

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 interactions 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, simulating all the relevant physics processes that lead to energy losses and creation of secondary particles: ionization and bremsstrahlung, gamma-ray pair production, Compton scattering, and muon photonuclear interactions. Hadronic processes such as pion production and interactions, lowenergy neutron interactions with nuclei, and photodisintegration are also modeled.

The 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 Ref. [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 setting that is critical for our study.

Crucial options complement the default setting: EVAPORAT(ion) and COALESCE(nce) give a detailed treatment of nuclear deexcitations 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 (θ,φ) 2R cosθ H sinθ FIG.

3. Spatial distribution of trajectories for muons produced in the same area of the sky. These muons can be considered almost parallel when reaching SK, and the intersection of their trajectories with a plane perpendicular to their direction 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 intersections of the muon trajectories with the plane.

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 structure, 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

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 AE 0.06) measured 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%.

Detector response and event reconstruction

We model the detector response using SKDetSim. This simulation is based on GEANT 3.21 [20] for detector modeling and uses a customized model for light propagation and collection. It covers all aspects of event detection, from the initial interaction to the light collection on the PMTs and event reconstruction.

For this purpose, SKDetSim models in detail the entire geometry of the detector, the particle propagation in water, the emission of Cherenkov photons, reflection and absorption of light on materials, photoelectron production, and electronic response. Simulated data and real data are processed similarly. For this reason, particles simulated in SKDetSim are recorded in trigger windows corresponding to the ones applied to data. For each event, the detector dark rate is also simulated and can be added or not to the outputs depending on the specific needs.

Since FLUKA already models the shower development in water, special care must be taken when interfacing with SKDetSim, which is used for the detector response. It is important to avoid having SKDetSim do a parallel generation of muon-induced showers, which would lead to a double counting of isotopes. In this study, we therefore focus on observables that are particularly robust against possible mismodeling of the shower Cherenkov light pattern: the yields of the spallation-produced isotopes and the characteristic neutron clouds discussed in Sec. V B 2. Isotope decays and neutron captures can indeed be simulated in isolation from their parent muons in SKDetSim, and will hence not lead to unwanted interactions. Muon-induced showers still need to be simulated, as they affect the reconstruction of the muon track and hence neutron identification; however this effect on reconstruction is limited. The impact of shower mismodeling can therefore be mitigated using a few simple steps. The following points describe the interfacing between FLUKA and SKDetSim for a typical event, which consists of a primary muon, a hadronic shower, potential neutron captures, and isotope decays.

An essential change is the deactivation of muon-nucleus interactions (GEANT-MUNU is set to zero) to prevent the muon from inducing extra showers. Thus, the muon will on average behave only as an ionizing particle. Radiative losses through interactions with atomic nuclei are deactivated, but processes such as bremsstrahlung and pair production are still possible, and their importance increases with muon energy.

Together with the muon, particles produced in hadronic and electromagnetic showers in FLUKA are injected at t ¼ d=c vac , where d is the distance from the muon entry point to the shower particle and c vac is the speed of light in a vacuum. This typically corresponds to a few nanoseconds.

Here again, in order to avoid double counting, we deactivate photofission as well as secondary particle generation for all hadronic interactions in SKDetSim, which is used only to model electromagnetic (EM) processes, including the emission of Cherenkov light. We only inject particles from FLUKA that are commonly produced in hadronic showers and can lead to the prompt emission of Cherenkov light: γ rays, pions, and kaons. The only exception occurs if a radioactive isotope is produced in a shower initiated by electromagnetic processes, e.g., by interactions involving electrons, photons, or positrons. In this case, the particles initiating the shower are also injected into SKDetSim, provided their energy is larger than 0.1 GeV since no isotope production has been observed below this threshold in previous simulation studies [2]. This scenario is however particularly rare, and the impact of these extra particles on the muon light pattern will be extremely limited. Table I summarizes the processes treated in FLUKA and SKDetSim.

While the Cherenkov light emission from pions and γ rays coincides with the one from muon ionization, neutron capture typically occurs over much larger timescales. For a given muon, we therefore treat each neutron capture separately in SKDetSim. Since neutrons have already been propagated by FLUKA, we directly simulate a 2.2 MeV γ ray at each capture vertex.

Finally, using the isotope production vertices and decay times given by FLUKA, we simulate isotope decay products in SKDetSim. Each isotope decay can produce either a γ ray or a γ ray and a β ray, and is sometimes followed by a neutron. In what follows we will consider only products leading to prompt Cherenkov light emission, as their reconstructed energy distribution will affect isotope yield measurements.

Each of the steps described above requires modeling not only the signal, but also the noise in the detector. For muons, shower particles, and isotope decays, we use the modeling of the PMT dark noise from SKDetSim, based on regular measurements made over the SK-IV period. The treatment of neutron capture signals is more complex; this signal in water is particularly weak, and the predicted performance of the neutron tagging algorithm is highly sensitive to dark noise modeling. Therefore, we inject noise samples from data into the signal simulation results. These samples were taken at different times over the whole SK-IV using a random trigger, and can hence be used to reflect the time variations of the noise in the detector. Finally, the muon, neutron capture, and radioactivity events undergo the same reconstruction and reduction steps as the ones described in Sec. II and V B 2 for data. WIT triggers were used to reconstruct neutron captures.

Neutron candidate hit time and light pattern

The amount of Cherenkov light deposited by an event in SK, as well as the shape of the associated ring, allow the identification of well-reconstructed neutron capture events. Using the tools outlined in Sec. III, cuts are placed on neutron capture candidates based on these variables to increase sample purity.

a. Reconstructed vertex: Intrinsic radioactivity (except for radon daughters) events typically occur on the walls of the detector or in the surrounding rock. Hence, their reconstructed vertices will typically lie either outside the ID or near its walls. In what follows we will therefore require the reconstructed vertices of the candidate neutrons to lie inside the ID. Note that for events with energies lower than ∼3.5 MeV the WIT trigger already requires online reconstructed vertices to lie in the FV. WIT triggers above that energy only need to be reconstructed inside the detector.

b. Reconstructed energy E rec : When a neutron is captured on hydrogen, a single 2.2 MeV γ ray is released. The SK energy reconstruction for this analysis assumes an interaction resulting in an electron. Thus, the reconstructed energy is the equivalent to that of an electron interaction with the same amount of Cherenkov light production. As the energy is calculated based on the number of effective hit PMTs, it is subject to Poisson fluctuations of that number. The reconstructed energy is hence expected to, on average, underestimate the true photon energy; however its distribution extends far beyond 2.2 MeV due to Poissonian fluctuations in the light yield. To account for these fluctuations we will typically require E rec < 5-6 MeV.

c. Reconstruction goodness g t , g p : These observables have been defined in Sec. III B and measure the goodness

TABLE I. The table summarizes the main processes activated and deactivated in SKDetSim and the particles generated or ignored in FLUKA and injected in SKDetSim. Since neutrons have already been propagated by FLUKA, we directly simulate a 2.2 MeV γ ray at each of the capture vertices found by FLUKA. To avoid double counting we do not record the captures of the neutrons generated by SKDetSim. More description is included in the text. SKDetSim Deactivated Activated Photofission Electromagnetic interactions Secondary particle Cherenkov radiation Generation in inelastic interactions Decays FLUKA Generated Ignored Electromagnetic showers with spallation Electromagnetic showers without spallation (If E > 0.1 GeV) γ, π, K from inelastic interactions γ from n captures Isotope decay production

of the position and direction reconstruction of the neutron candidate. In order to determine suitable spallation cuts, 2.2 MeV γ rays were generated using SKDetSim. The events were simulated without dark noise, overlaid with real SK online raw data, and then processed with WIT. The resulting sample was then split into "good" and "bad" events, whose reconstructed vertices lie less and more than 5 m away from the true vertices respectively. Figure 4 shows the two-dimensional distributions of g t and g p for both event categories, with "good" events being associated with higher goodness values. In the simple case studies presented in this section, we will only require g t > 0.5. In the more refined strategy discussed in Sec. VIII we will impose cuts on both g t and g p in order to improve the neutron cloud localization. The resulting refined cut is overlaid on the goodness distributions shown in Fig. 4.

Correlations with muons

Neutrons produced in muon-induced showers will be found close in time and space to muon tracks. The observables we describe here allow one to identify these space and time correlations. Moving forward, we make the reasonable assumption that muons follow a straight path inside the detector. We then parametrize their tracks using an entry point and a direction, which are determined using the fitter described in Sec. III A. We then define the following observables: the time difference between a neutron candidate and a muon, Δt, and the transverse and longitudinal (with respect to the muon entry point) distances from the candidate to the muon track, l t and l LONG .

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 one to unambiguously link a neutron capture event to its parent muon. Additionally, for high-quality neutron cloud samples, 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, l t and l LONG : l t defines the distance between a neutron candidate and the closest point on the muon track while l LONG is the longitudinal distance, the distance between that closest 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 one to characterize the shapes and sizes of the neutron clouds. In this figure, hl LONG i refers to the average l LONG 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 l 2 t instead of l t to reflect the amount of phase space available.

Incorporating the observables defined above in spallation analyses allow one to identify neutrons and accurately locate neutron clouds and the showers that generated them. In what follows we will showcase studies highlighting the validity of our neutron identification and localization procedure.

Identifying individual neutrons

We build a high-purity neutron sample using WITtriggered events verifying g t > 0.5, l t < 1.5 m, and E rec < 5 MeV. We estimate the number of neutrons in this sample by comparing the distribution of time differences between neutron candidates and their parent muons, Δt, to the one expected from calibration studies 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:

where N is the number of events, τ is the exponential decay time constant, and the constant C absorbs remaining background contributions. For the WIT data the time constant was measured to be τ ¼ 211.8 AE 1.7 μs. In comparison, in AmBe calibration studies, the neutron capture time on hydrogen was measured to be τ ¼ 203.7 AE 2.8 μs. This results in about a 2.5σ difference between the two measurements. This discrepancy is believed to be due to missing neutrons within higher multiplicity showers from pileup. 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 dead time [36] for captures close together in time, biasing for a longer capture time constant. The energy of showering neutrons will have a subdominant contribution to capture time as well due to their effect on thermalization times. This discrepancy also is reduced/ increased when fitting to samples comprised of low/high multiplicity showers.

neutron/ spallation candidate l t (transverse distance) l LONG (longitudinal distance) θ μ track FIG. 6. Diagram showing the l t and l LONG observables associated with neutron identification, as described in Sec. VA 2. . Time distribution of WIT-triggered events with l t < 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.

Identifying neutron clouds

After having successfully identified neutrons using WIT triggers, we will now investigate correlations between 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. Additionally, we require at least one of these WIT events to have Δt > 20 μs, g t > 0.5, and E < 5 MeV. Note that fewer events per unit volume are reconstructed outside than inside 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 one 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 proximity to the center of a neutron cloud. Spallation candidate 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 spallation 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 16 N decays, and the 5 m value reflects the general size of muoninduced showers, as seen in Fig. 7. In addition to this signal sample, we build a background sample using candidates found within 5 m of and up to 60 s before neutron clouds. This background sample allows one to estimate the fraction of spurious pairings between spallation candidates and uncorrelated neutron clouds in the signal sample, and subtract off the corresponding effects.

A strong spatial correlation between spallation candidate 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 spallation candidates relative to the others at low multiplicities is due to accidental pairings between candidates and neutron clouds. Conversely, high-multiplicity neutron clouds can be more easily located and allow one 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 allows one to remove 55% of spallation events with a little more than 4% dead time. Owing to the low detection efficiency for the 2.2 MeV γs resulting from neutron capture on hydrogen, small showers are missed, and therefore this procedure 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 traditional 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 g p to FIG. 8. Δt (black dots) and resulting fit (red 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 AE 1.7 μs. The horizontal axis is distance 3 from the neutron cloud center so that random coincidence events would be more or less uniformly distributed. For each multiplicity these random coincidence events are not subtracted, but the distribution is normalized by the total number of spallation events (i.e., the difference between the total events and the random coincidence events). As multiplicity increases, the steepness in the tail of the distribution increases as expected from a more accurate parametrization of the cloud and from a larger fraction of spallation events for larger hadronic showers. Note that the vertical axis is on a logarithmic scale.

better identify well-reconstructed neutrons and will account for the elongated shape of neutron clouds when investigating spatial correlations with isotope vertices.

Although the neutron cloud cuts presented here cannot be used as a standalone spallation reduction technique in SK-IV, their impact is expected to drastically increase in the new SK-Gd phase of the detector. Estimating this new impact and redesigning spallation cuts accordingly requires extensive simulation studies. In the next section, we compare neutron cloud measurements in pure water to predictions from the simulation described in Sec. IV and show that this simulation can be reliably used to optimize future neutron cloud cuts.

VI. SIMULATED HADRONIC SHOWERS

When comparing SK-IV data to the simulations described in Sec. IV we neglect effects associated with muon bundles and consider only single muons traveling through the whole detector, as previously mentioned. Bundles make up to 5% of the total number of muons, and only one of them is expected to produce a shower; therefore, for moderate-size bundles the neutron multiplicity is not expected to change. In order to accurately model muon showers and reconstruction effects in the simulation we use the following two samples: a. Hadron producing muons: We define hadron-producing muons as muons that undergo an inelastic interaction with a nucleus in water leading to the production of one or more particles of the following types: pions, kaons, or neutrons. These particles often initiate a cascade. These showers can be observed in about 11% of the muons generated using the procedure described in Sec. IV, and 5% of them lead to the production of isotopes contributing to backgrounds in low-energy SK searches. Spallation backgrounds are hence induced by about 0.5% of the muons passing through SK. For this study we generate a sample of 2.7 × 10 5 muons, 97% of which are reconstructed as single through-going muons [37].

b. Electromagnetic only muons: This category includes muons that either undergo only minimum ionizing interactions or induce purely electromagnetic showers, with no pions, kaons, or neutrons. Here we generate a sample of 7.5 × 10 5 muons. Since we do not need to model hadronic showers and isotope production, we do not use FLUKA to propagate particles in water and directly combine MUSIC with the detector simulation in the pipeline described in Sec. IV. The fraction of single through-going muons in the final sample is only 87%, as muons in this sample are typically less energetic than spallation-inducing muons and more likely to stop inside the detector.

We use the analysis procedure detailed in Sec. V to extract observables characterizing the neutron cloud shapes and the neutron multiplicity in muon-induced showers. Comparing the simulation results to the SK-IV data in pure water will allow us to estimate the systematic uncertainties associated with the neutron modeling and motivate the use of a similar simulation to optimize spallation cuts for the current SK-Gd phase [38].

Neutron cloud shapes

We define the shape of a neutron cloud using the l t and l LONG observables defined in Sec. VA and shown in Fig. 6. For the latter, in order to compare multiple clouds, we use the average l LONG of the cloud as a reference and consider

In what follows we define a neutron cloud as a cluster of two or more reconstructed neutrons. In the simulation these reconstructed clouds represent only about 5% of all real neutron clouds, due to the low neutron identification efficiency.

In spite of the low neutron mistag rate, due to the large fraction of minimum ionizing muons, the data sample will be contaminated by non-negligible contributions from fake neutrons. Estimating the contribution from these fake neutrons requires determining the fraction of muons not leading to hadronic showers (referred to as electromagnetic-only muons at the beginning of this section), which is determined by nuclear effects that are difficult to model accurately [39]. This fraction, however, can be readily extracted from data since fake neutrons will populate the tails of the l t and Δl LONG distributions while contributions from real neutrons, namely neutrons produced in muon induced showers in the simulation, will dominate at small distances. At SK, with a muon rate of about 2 Hz, only a few months of data taking are needed for this estimate.

For this study, we evaluate the fraction of EM muons by fitting the predicted l t and Δl LONG distributions to the SK-IV data. We perform a separate χ 2 fit for each observable in order to evaluate the robustness of our model and find that the best-fit fractions of muons without hadronic showers are 96% (χ 2 =NDF ¼ 36.3=33, with NDF being the number of degrees of freedom) and 97% (χ 2 =NDF ¼ 158.4=99) for l t and Δl LONG respectively. These values are compatible with each other, but larger than the FLUKA prediction of 89%. For this analysis we use a fraction of 96.5% and treat the difference between the fit results for the l t and Δl LONG as a systematic uncertainty. Table III shows the resulting muon rates for different tagged neutron multiplicities.

The distributions of l t and Δl LONG for all neutrons are shown in Fig. 12 for the simulation. Here the uncertainties combine both the statistical uncertainties and the systematics associated with the different optimal EM muon fractions described in the previous paragraph. As expected, fake neutrons are associated with larger distances as they can be observed anywhere within 5 m of the muon track. The real neutron component overwhelmingly dominates for Δl LONG < 5 m and l t < 3 m.

Figure 12 shows the l 2 t and Δl LONG distributions for the SK-IV data and the simulation. The observed neutron clouds have an elongated elliptical shape, with average transverse and longitudinal extensions of 3 m and 5 m respectively. For distances of less than about 5 m, where contributions from real neutrons dominate, the predictions differ from the data by at most 15%. This excellent agreement motivates the use of a FLUKA-based simulation to predict neutron cloud shapes and optimize the associated cuts for future spallation analyses, notably at SK-Gd and Hyper-Kamiokande [40]. At SK-Gd in particular, these simulationbased studies will allow one to significantly reshape the spallation reduction procedure, as gadolinium doping will sizably increase the neutron identification efficiency.

Neutron multiplicity

We finally estimate the number of reconstructed neutrons 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 possible mismodeling of this fraction as a systematic uncertainty. The neutron cloud multiplicities for both simulation and data are shown in Fig. 13. There are of course additional uncertainties associated with the details of the simulation of neutron production in FLUKA. These are not discussed here, since the focus of this work is the presentation of the experimental observation and comparing them to the nominal FLUKA prediction. The abundance of low-multiplicity clouds is due to both fake neutron contributions 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

TABLE II. Neutron tagging efficiencies for different selection criteria applied sequentially: WIT trigger, timing cut to remove the primary muon signal and the associated after-pulse, proximity to the muon track, fiducial volume cut, goodness-of-fit cut, and reconstructed energy (E n ) cut. We consider a simulated neutron to be tagged if its true capture time lies within 50 ns of the reconstructed capture time of a neutron candidate. Selection requirements Efficiency values (%) Total efficiency (%) WIT trigger 13.3 AE 0.1 13.3 20 μs < Δt < 500 μs 89.4 AE 0.4 11.9 l t < 500 cm 76.7 AE 0.3 9.1 FV 93.3 AE 0.5 8.5 Fit quality 76.9 AE 0.4 6.5 E n < 5 MeV 99.9 AE 0.6 6.5 TABLE III. Muon rates in SK vs number of tagged neutrons. A neutron cloud is found only if at least two neutrons are tagged. The muon rates are given separately for hadron producing and electromagnetic only muons, where the fraction of the second sample is extracted from a minimum χ 2 fit as described in the main text. Only through-going muons are considered. Muon rate [Hz] if 0n tagged ≥ 1n tagged ≥ 2n tagged Hadronic muons 0.051 0.0088 0.0023 EM muons 1.48 0.066 0.0018

simulate the tails of the data distribution. Note, however, that such large clouds are typically associated with showers producing hundreds, sometimes up to thousands of neutrons. Muons associated with these high-multiplicity showers are not only rare but also deposit a high amount of light in the detector and are hence easier to identify using other observables, such as the residual charge Q res that will be introduced in Sec. VIII C. In any case, neutron multiplicity distributions 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 lowenergy analyses at SK target antineutrinos, whose interactions 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, t (top) and Δl LONG (bottom) distributions for neutrons belonging to a cluster of multiplicity larger than or equal to 2, for simulations (solid black) and data (dashed red). The fraction of muons that do not produce hadronic showers is set to 96.5%; fake neutrons from those are included (see Fig. 11). All distributions are normalized to an area of 1. The dotted lines denote the average transverse and longitudinal cloud extensions of 3 m and 5 m respectively. 11. l t (top) and Δl LONG (bottom) distributions for neutrons belonging to a cluster of multiplicity larger than or equal to 2, for simulation results. The fraction of muons that do not produce hadronic showers is set to 96.5%. The black histograms with a solid line represent contributions from all reconstructed neutrons while the gray filled and black filled distributions show the true and fake neutron components respectively. 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 analysis and, in particular, of the previous spallation reduction strategy.

B and hep solar neutrinos from the Sun scatter elastically off electrons in SK, with the recoiling electron producing a single Cherenkov ring. The produced electrons have energies up to about 20 MeV; however, intrinsic radioactivity 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 cuts on the full SK-IV solar neutrino analysis will be described in [44].

After high-efficiency noise reduction cuts are applied, notably rejecting events outside the FV, energy-dependent cuts are imposed on the quality reconstructed vertices and directions of the events, the goodness observables described in Sec. III B. The major quality cuts select events based on their reconstruction goodness (g R ¼ g 2 tg 2 p ) and their effective distance to the ID wall-obtained by following their reconstructed direction backwards. We then check for consistency with a single electron ring (rather than, e.g., light deposited by βγ decays of 16 N). We parametrize the amount of multiple Coulomb scattering of the electron (the "fuzziness" of the ring) with the "multiple scattering goodness" to separate lower energy β decays. Most of the remaining background 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 considering event directions. Indeed, electrons recoiling from neutrino elastic scattering will be almost collinear with the incoming neutrinos, that is, the angle between the direction 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 quality cuts described above is paired with all 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 transverse distance of the candidate to the muon track l t (Fig. 6), and the residual charge Q res , defined as the excess charge deposited by the muon in the detector compared to the expectation from minimum ionization. For each observable, probability distribution functions (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 us to define a log likelihood function log 10 L, whose functional form is as follows: Neutron multiplicity distributions for simulation (solid black) and data (dashed-plus red) with the EM muon contribution shown for the simulation (filled gray) which is significant only for less than 4 neutrons. 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). The vertical axis is the percentage of neutrons of a given multiplicity compared to all neutron clouds passing the cuts.

where PDF spall;i and PDF ran;i designate the PDFs associated 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 candidates 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 detector. When paired with preceding muons, this random sample estimates the random coincidence contribution, so the corresponding PDFs are extracted. A statistical subtraction of the random sample distributions from the data sample distributions yields the spallation PDFs. Alternatively, we invert the time sequence and pair solar candidates with all 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 observables, the l t PDFs are computed for seven different Q res bins. Since the muon fitter used for these cuts considered 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.

A cut on the likelihood [Eq. ( 6)] used previously in the solar analysis [5] removed 90% of spallation events with a position-averaged 20% dead time. This dead time 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 one to further reduce this dead time 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. the muon track reconstruction was replaced. Previously we used a simple, fast muon track fitter developed at the beginning of SK. It assumes through-going single muons and misreconstructs 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 background by dE=dx reconstruction in the search for diffuse supernova neutrino interactions in SK-I, II, III [45], which inspired the development of FLUKA-based simulation studies [1], highlighted the importance of muon-induced hadronic showers for isotope identification, and allows the characterization of the shapes and sizes of the neutron cloud. Finally, the improved detector electronics associated with the SK-IV phase has a larger dynamic range and therefore allows detection of higher values of the total charge (and therefore Q res ), as well as identify neutron clouds as described in Sec. V.

The new spallation reduction strategy proceeds as follows. First, we apply two sets of preselection cuts in order to remove a sizable fraction of spallation events with minimal harm to the signal efficiency. These cuts aim at removing solar neutrino interaction candidates 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 remaining 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. VA, we define a set of cuts to reliably identify neutron clouds and investigate 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 E rec < 5 MeV. We then assign weights to these candidates depending on their vertex and direction goodness g t and g p . Specifically, we consider three regions of weights, 0, 1, and 2 in the g tg p space shown in Fig. 4.

Once clouds are identified and their barycenter is defined, their positions and detection times can be compared to the location and times of solar event candidates. For this analysis we consider SK-IV low-energy events that passed all the cuts defined for the solar analysis in Sec. VII except spallation cuts. This sample is expected to be largely dominated by spallation isotope decays. We consider 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 specific 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 the 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 transverse and longitudinal distances of low-energy events, l t and Δl LONG ¼ l isotope LONGl n-cloud LONG , defined in Fig. 14, to the origins of their respective coordinate systems. The distribution of l t and Δl LONG 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 origin significantly reduces the spread of this distribution in l t , a spread that is primarily driven by contributions from low multiplicity clouds. We hence choose this definition 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, l t and Δl LONG for each low-energy event-neutron cloud pair, with Δt defined 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 neutrons. Then, we define multiplicity bins of 2, 3, 4-5, 6-9, and ≥ 10 neutron candidates and, for each bin, define a specific ellipsoidal cut on l t and Δl LONG . Since clouds with only two candidates are often not associated with hadronic showers, as shown in Sec. VI, we require Δt < 30 s. For higher multiplicities we consider all muons up to 60 s before the low-energy event. The different cuts are summarized in Table IV. If a low-energy event-neutron cloud pair correlation is within the required Δt, l t , and l LONG values, the corresponding low-energy event is rejected.

Since the sample used for this analysis is largely dominated by spallation isotope decays, the background rejection rate of the neutron cloud cut can be readily estimated. To evaluate the dead time, we use a sample of events with reconstructed energies between 3.5 MeV and 5 MeV whose vertices have been replaced by randomly generated ones. We then pair these events with muons observed up to 60 s after them similar to the steps outlined in Sec. VII B and apply the reduction steps outlined above, with the sign of Δt inverted. The dead time is then given by the fraction of remaining low-energy events and has been found to be 1.3%.

B. Multiple spallation

Since neutrons in a hadronic shower are a good indicator of the production of spallation nuclei, the observation of the

μ track iso tope x y z θ n μ track x y z isotope Δl LONG l t μ track x y z isotope Δl LONG l t FIG. 14. Possible options for the neutron cloud coordinate systems. Left: the origin is the neutron cloud center, and the axes are the neutron cloud axes. θ n is the angle with respect to the muon track. Center: the origin is the neutron cloud center, and the z axis is aligned with the muon track. Right: the origin is the projection of the neutron cloud center onto the muon track, and the z axis is aligned with the muon track. We also show the definition of the transverse and longitudinal distances l t and Δl LONG , between a spallation isotope and the origin. decay of a spallation nucleus is likewise an indicator for a hadronic shower. We therefore apply, in addition to the neutron cloud cut, a preselection removing clusters of isotopes produced by the same muon. Instead of pairing candidate events with possible parent muons, this multiple spallation cut identifies clusters of low-energy events observed within a few tens of seconds and a few meters of each other. Here, we consider a sample composed of all SK-IV events passing the first reduction and quality cuts defined for the solar analysis, as discussed in Sec. VII. Since we need to take all spallation isotope decays into account we apply neither the old spallation cut nor the pattern likelihood cut, as the latter targets 16 N βγ decays. We find that the optimal cut removes candidates found within 4 m and 60 s of any event from this sample. This cut allows one to remove 45% of spallation background events with a dead time of 1.3%. The solar angle distribution of the rejected and the remaining events is shown in Fig. 16. The absence of a peak around cos θ sun ¼ 1 for the rejected events confirms the low dead time for this cut.

Longitudinal distance Δl LONG

Defining PDFs for the difference in longitudinal length Δl LONG ¼ l isotope LONGl dE=dx peak LONG allows one to include information about the muon-induced hadronic shower into the likelihood cut. Here, we developed a new method to reconstruct the energy loss of the muon along its track Residual Charge [p.e.] 0 1000 2000 3000 4000 5000 3 10 × (dE=dx) based on a previously published method [45]. Defects and possible improvements to this method were suggested in Ref. [3]. The method presented here remedies 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 verifying the time correlation equation for each PMT hit:

where t PMT and t entry 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 emitted from, d is the distance from the emission point to the PMT, and c vac and c water are the speeds of light in vacuum 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 estimating 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 Ref. [1], we spread the charge from each hit across multiple bins to account for the PMT timing resolution. More specifically, we take the contribution, g ij , of the ith PMT to the jth bin to be

where

σ i is the timing resolution for the observed charge by the ith PMT, and τðl j Þ is the t PMT that solves Eq. (11). For the jth bin boundary's d and l, Erf is the standard error function, Q j 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 nonlinear, 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Þt i ¼ 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 location of these showers. We hence define Δl LONG ¼ l isotope LONGl dE=dx peak LONG as the longitudinal distance of an isotope to this segment along the muon track. Distributions of this observable for spallation and uncorrelated pairs are shown in Fig. 19. To define the corresponding PDFs, the Δl LONG distribution for low-energy events found within 2 m and 10 s of a muon is fitted by a sum of three Gaussians (parameters are listed in Table XI):

where A i , x, and σ i are the fit parameters for each PDF. For uncorrelated pairs, one of the Gaussian fits was degenerate and dropped from the final form. For minimum ionizing muons, the dE=dx peak is more likely to be at the end of the track, resulting in the background distribution being slightly shifted away from 0. Preliminary studies exploring the reasons for this shift were inconclusive and left for future work. Although improvements can be made, consistent application of this approach allows for trusted results.

Using the PDFs defined above, we define a log likelihood 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-energy event sample that we are considering is dominated by spallation and solar events. We can hence readily estimate 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 random 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 neutron cloud information is available. The availability of neutron cloud information notably allows one to loosen the likelihood cut.

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 reduces 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, which shows the cos θ sun distribution for events passing the new spallation cut and failing the one described in Sec. VII. The clear peak around cos θ sun ¼ 1 shows the new procedure retrieves a sizable number of solar events. In the final sample, this cut increases solar neutrino signal events by 12.6%, with a 0 500 1000 1500 2000 2500 3 10 × 2 in cm 2 radius 1500 -1000 -500 -0 500 1000 1500 z in cm 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 [9]. 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 counts/bin 100 -0 100 200 300 400 500 600 FIG. 22. Figure showing the difference in events between final solar samples using new spallation cut compared to the previous cut for 5.99 to 19.5 MeV. The peak at the right contains 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 interaction 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.

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.

A. Selection of spallation isotopes

To build a spallation-rich sample, we select events passing the noise reduction steps outlined in Ref. [44] and apply the following quality cuts:

(1) Events with reconstructed vertex more than 200 cm from the ID wall (FV cut) (2) Events less than 50 μs of the muon were rejected to remove cosmic μ-e decays as well as PMT afterpulsing. In addition to the noise reduction cuts, we apply cuts designed to increase the fraction of spallation events in the sample. First, we require the event kinetic energies to lie between 6 and 24.5 MeV. This window ranges from the energy at which spallation starts dominating over intrinsic radioactivity to the highest possible energy of the isotope decay products.

To fit the time difference between pairs of spallation isotopes and their parent muons, we pair spallation candidates with muons observed up to 100 s before them. Here, we extended the Δt range compared to the solar analysis from Sec. VII in order to include all long-lived isotopes. Additionally, we consider only muons with a reconstructed track of at least 200 cm. For stopping muons, we compute l t by extending the track through the entire detector; for multiple muons we considered only the primary track. Owing to the 200 cm requirement, corner clipping muons will not be taken into account. Their contribution to the spallation background in the fiducial volume is however negligible; the total number of events found within 500 cm of a corner clipping muon track is about 10 -5 of the total number of spallation events.

After pairing spallation candidates with suitable muons as described above, we build two separate samples with a high fraction of spallation pairs. In the first sample, we apply the same selection cut as in the previous study by Ref. [4], requiring l t < 200 cm. In the second sample, we select muons associated with three or more tagged neutrons. In Sec. V B events were required to have at least two events triggered in WIT within 500 μs and 500 cm, and at least one event passing the quality cuts. For this study, we instead require at least three neutrons triggered by WIT to further increase purity. The cloud location is still only defined by the weighted center of quality neutron cloud events. Then, we pair the selected neutron clouds with spallation candidates verifying the conditions listed in Table IV. Within this table, all cuts requiring Δt < 60 s were used for pairs with Δt up to 100 s.

B. Spallation yield fit

To extract isotope yields from data, we fit the decay time distributions of the ten most abundant isotopes found in [1,4] to time difference Δt distribution using the spallation sample described above. Specifically, we parametrized the time dependence of the total event rate R TOT by a sum of exponentials as follows:

where we keep the isotope decay constants τ i fixed and fit the production rates R i . For this study, we perform a χ 2 fit analytically over the whole Δt ¼ 0-100 s range.

To mitigate possible degeneracies, we grouped isotopes together when their decay constants were within 10% of each other. For this study, 8 Li and 8 B as well as 9 C and 8 He were paired together. The decay constant associated with each of these pairs was taken as a weighted average between the two isotopes, using the yields predicted by [1].

The τ i s used were 1.18 s and 0.181 s for the two pairs respectively. For the case of stopping muons, looking at the individual fit contributions, most spallation products were from those which are normally produced from a neutron interacting with a nucleus. The fits are dominated by 16 N, with the next largest contribution being 12 B with an observed raw rate roughly a factor of 6 smaller. For single through-going muons, 16 N is observed roughly twice as often.

In contrast to the previous paper from Ref. [4], the rate distribution (Δt distribution divided by time bin width) was made with logarithmic bins. Also, as mentioned earlier, the maximum Δt considered for the fit was extended from 30 s to 100 s in order to account for long-lived isotopes and better constrain the distribution of accidentally coincident muon/spallation isotope candidate pairs. The result of this fit for all 10 isotopes is shown in Fig. 23. The χ 2 =NDF was 231.5=243 resulting in a p value of 69%. Figure 23 shows the fit to the entire time range as well as the residual of the fits. Overall, an excellent agreement between the data and the fit can be observed, and the rates of the most abundant isotopes can be determined with percent-level precision.

C. Spallation efficiencies

Obtaining isotope yields from the rates estimated above requires computing the efficiencies of the selection cuts described in Sec. IX A. These efficiencies were computed using both MC simulations and data.

The efficiencies of the noise reduction, quality requirements, and energy cuts for the different spallation isotopes have been determined by simulating isotope decays using the SK detector simulation, based on GEANT3.21 [20]. In this simulation, we accounted for the time dependence of the detector properties using measurements of the PMT and water properties taken almost daily. The decay spectra are simulated using GEANT4.10.7 [46,47] for most isotopes, taking into account both β and γ emission. Since the 8 B decay was mismodeled in GEANT4 we modeled its β spectrum using tabulated values from Ref. [48]. The final efficiencies are shown in Table V.

To determine the efficiency of the l t < 200 cm cut introduced in Sec. IX A, we combine estimates from the data and from the spallation simulation described in Sec. IV. The simulation readily gives us efficiencies for individual isotopes paired with single through-going muons, shown in Table V, but does not currently allow one to study other muon categories.

We complement this simulation study by extracting the l t cut efficiency from data for each muon type; we later take the total efficiency to be a weighted average using the observed relative contributions of each muon category. We obtain the efficiencies for each muon type by fitting the Δt distributions observed in 50 cm l t bins using the method described in Sec. IX B. The 50 cm bin Δt distributions were scanned to determine the largest l t where evidence for spallation could be found (the "integration range"). The integration ranges for each muon category, the l t efficiency, and the spallation proportion are found in Table VI. To combine the final l t efficiencies we then compute the relative contributions from each muon category to the spallation sample by considering TABLE V. Spallation sample efficiencies from simulation, for single through-going muons. The second column shows the isotope decay identification efficiencies for the FV, first reduction, and 6 MeV cuts, and the third column shows the efficiencies for the l t ≤ 200 cm cut. The total l t efficiency has been obtained by averaging the efficiencies associated with the different isotopes, weighted by their predicted yields.

Isotope

Fiducial volumeþ preselection þ E ≤ 6 MeV (%) l t ≤ 200 cm (%) 12 N 69.0 94.9 12 B 54.3 89.9 8 He 23.9 94.6 9 C 67.4 94.0 9 Li 44.4 91.8 8 Li 51.8 91.4 8 B 48.1 92.9 15 C 37.1 87.7 16 N 54.8 88.9 11 Be 44.7 86.5 Total 90.3 AE 2.8 t [s] Δ 4 - 23. Left: full Δt distribution from 50 μs to 100 s, with bins corrected for their logarithmic width. The dashed red line corresponds to the entire fit, and the colored lines correspond to the individual isotope fits and background. The black points (partially obscured by the entire fit) are the observed data. The overall event rate for SK-IV using the constant fit here corresponds to .288 mHz for events meeting the cut criteria. Right: fit residuals shown using unscaled data minus the fit value at the bin center multiplied by the width of the corresponding time bin. all pairs within l t < 2000 cm-the maximal distance between spallation isotopes and their reconstructed parent muon track. The resulting total l t efficiency is calculated to be 73.6%.

A second method identical to the one used in [4] was implemented to compare: After applying a precut on Δt, the efficiency is the ratio of the number of events within the l t cut value over all events passing the precut. Eight different precut 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 l t efficiency of 74.2 AE 0.5%. The total l t efficiency is then taken to be 74.0 AE 0.7% to account for the discrepancies between the two methods. The l t 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 l t distribution was included in the measurement by scaling each isotope 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 systematic error. This error was relatively small for most isotopes, with only 15 C and 11 Be having an effect greater than the full l t systematic error. Here, using the simulation allows one to refine the l t efficiency estimate performed in Ref. [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 l t efficiency obtained from the data is combined with the efficiency associated with the noise, quality, and energy 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 isotopes that were paired together for the Δt fit, the contributions of the

TABLE VI. l t cut efficiency and proportion of total spallation production for the different muon categories. This proportion is measured as the ratio of the integral of a category from 0 to 20m to the sum of those integrals over all categories. Category Range [m] Efficiency [%] Proportion [%] Single 0-5 92.3 54.6 Multiple 0-30 49.3 43.5 Stopping 0-4 93.4 1.9 Corner clipping N/A N/A <10 -3 All N/A 73.6 100 TABLE VII. Observed and calculated spallation rates and yields. Observations are based either on a simple l t cut (2m) of muon tracks to subsequent low-energy events or neutron clouds with at least three neutrons with a required vertex correlation according to Table IV.

The first column indicates the considered isotopes. (For the two sets of isotopes that could not be separated, the systematic errors will cover the range of yields corresponding to changing the relative fraction of an isotope from 0 to 1.) The second column is the decay constant for the given isotope. The third (fourth) column shows the calculated (observed via the l t cut total rate) production yield. The fifth (sixth) column is the total rate based on the l t cut with (without) efficiency correction. The seventh column shows the neutron cloud total rate without efficiency correction. The eighth column is the ratio of neutron cloud rate over l t rate. If the isotope mix in both samples were identical, this ratio would be constant (and would represent the efficiency ratio). The last column is the fraction of isotopes as found in simulation, with blank entries corresponding only when an upper limit could be found in data. If an isotope was not observed, a 90% confidence level (C.L.) upper limit is shown using the positive systematic and statistical error. Compared to the previous results, 15 C has been discovered, and a much stronger constraint on the 11 Be measurement has been placed. The 11 Be rate fit exceeds zero by 1.5σ. The 9 C= 8 He and 9 Li fits for neutron clouds were 98.2% anticorrelated. The total number of events for the two fit contributions are 65% of the corresponding l t contributions. The calculated relative fractions of 8 Li and 8 B in the 8 Li þ 8 B sample are 70.1% and 29.9% respectively. The calculated relative fractions of 9 C and 8 He in the 9 C þ 8 He sample are 78.8% and 21.2% respectively. Isotope τ i [s] Yields [10 -7 cm 2 μ -1 g -1 ] Rates [kton -1 day -1 ] Neutron data Fraction of l t data Calculated Isotope Fractions Calculated Observed Total rate (l t data) Raw rate (l t data) Raw rate (neutron data) 12 N 0.0159 0.92 1.72 3.04 AE 0.06 AE 0.028 1.55 1.08 70% 2.3% 12 B 0.0294 8.6 12.9 22.86 AE 0.11 AE 0.21 9.19 5.95 65% 21.1% 9 C= 8 He ∼0.175 0.8 < 0.61 <1.08 0.11 0.20 176% Á Á Á 9 Li 0.257 1.5 0.67 1.19 AE 0.33 AE 0.010 0.39 0.13 34% 3.7% 8 Li= 8 B ∼1.16 13.4 5.11 9.04 AE 0.17 þ0.60 -1.1 3.69 2.59 70% 32.8% 15 C 3.53 0.55 1.57 2.78 AE 0.45 AE 0.032 0.76 0.37 49% 1.3% 16 N 10.3 14.5 27.3 48.43 AE 0.60 AE 0.49 19.64 12.01 61% 35.3% 11 Be 19.9 0.61 <1.05 <1.9 0.33 0.19 56% Á Á Á

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 simulated spallation sample of 1.362 × 10 8 initial muons is generated in order to accumulate enough statistics for the low yield isotopes compared to 4.8 × 10 8 observed initial muons. 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:

where ρ is the density of water, R i 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.

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 detector. We notably developed algorithms locating muoninduced hadronic showers, both by improving the reconstruction of the energy deposited along the muon track and by identifying neutrons using a recently implemented low-energy trigger. New spallation cuts based on these algorithms allow one to reduce the dead time of the solar neutrino analysis by a factor of 2, allowing a gain of the equivalent of one year of exposure at SK-IV. Moreover, the profiles of the neutron clouds produced in muoninduced hadronic showers are well reproduced by the spallation simulation, motivating its use to develop spallation reduction algorithms for future analyses.

In addition to developing new spallation cuts, we computed 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 2 of the observed values, well within the uncertainties associated with hadron production models. This study also demonstrated that identifying neutron captures associated with muons allows one 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 muoninduced showers. At SK-IV, the performance of our neutron tagging algorithm has been limited by the low live time 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 upcoming spallation reduction algorithms. In this context, the simulation presented in this paper will be instrumental 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. the random coincidence-subtracted spallation sample. For contributing variables with only one spallation and random coincidence function (Δl LONG , Δt, and Q res ) there is no normalization factor while for l 2 t a normalization factor is needed due to the multiple bins for the function. The normalization factors are not listed.

Residual charge

Both signal and background PDFs in the positive Q res region are a sum of exponential functions:

where c i and p i are the fit parameters for the exponential functions. For negative Q res , no analytical form was assumed, and a linear interpolation of the sample bins was used between the four points listed in Table X.

TABLE X. PDF points for negative residual charge. Residual charges below -100; 000 are not considered. PDF Point Residual charge [1000 p.e.] PDF value Spallation 1 -100 10 -6 2 -65 0.000713 3 -22.5 0.02107 4 -7.5 0.34053 Random 1 -100 10 -6 2 -65 0.07034 3 -22.5 0.609 4 -7.5 23.29 TABLE XI. PDF fit parameters for Δl LONG . PDF Gaussian Amplitude Offset [m] Width [m] Spallation 1 989.8 -3.090 0.9217 2 1185 26.32 1.852 3 254.6 -58.87 13.49 Random 1 150.0 -1.160 13.41 2 20.18 -25.16 3.976 TABLE XII. Spallation and random PDF fit parameters for l 2 t . Time [s] resQ [1000 p.e.] Spallation Random c 1 p 1 [m -2 ] c 2 p 2 [m -2 ] c 3 p 3 [m -2 ] l 2 t0 [m 2 ] p 0 p 1 [ð100 mÞ -2 ] p 2 [ð100 mÞ -4 ] 0-0.1 < 0 8.023 5.506 7.678 2.082 5.180 0.4476 0 8.035 35.09 140.6 0-25 8.929 4.170 8.145 1.511 5.329 0.3031 64 11.19 16.27 -75.68 25-50 8.354 4.860 7.278 1.348 4.696 0.2588 0 2.638 30.20 11.59 50-100 8.395 2.820 6.036 0.5158 1.774 0.03065 0 2.005 35.58 114.2 100-500 8.770 2.248 6.129 0.3654 2.712 0.01120 0 1.954 35.02 71.72 500-1000 6.869 1.906 4.353 0.2588 1.713 0.00967 0 0.1952 36.25 0 > 1000 5.772 2.542 4.728 0.3412 1.900 0.00755 0 0.2675 61.74 0 0.1-3 < 0 8.065 4.446 5.021 0.4160 7.063 1.773 0 273.2 26.02 -47.17 0-25 8.758 3.067 6.683 0.6694 -16.92 0.02444 16 365.6 10.89 -107.8 25-50 7.543 4.283 7.361 1.656 4.470 0.2228 9 75.61 19.93 -115.7 50-100 8.087 3.467 6.093 0.6448 2.574 0.06516 4 60.11 23.04 -146.0 100-500 8.377 2.684 6.633 0.5569 2.579 0.01071 0 64.97 28.45 -165.6 500-1000 6.533 1.847 4.544 0.3578 1.464 0.01112 0 7.004 45.45 -1.876 > 1000 5.587 1.997 4.139 0.2551 2.013 0.01305 0 3.265 53.02 -199.6 3-30 < 0 8.742 2.810 6.099 0.4985 -100117 349594 0 2557 26.18 -39.10 0-25 2.611 0.00484 7.220 0.6481 8.999 3.075 16 3377 10.11 -115.0 25-50 8.080 3.452 7.078 0.8928 3.131 0.03730 0 719.2 19.35 -104.9 50-100 8.071 2.845 6.805 0.8500 3.284 0.03119 0 571.7 23.95 -103.2 100-500 8.656 2.048 6.302 0.3115 2.164 0.00498 0 630.7 29.57 -120.8 500-1000 6.782 1.465 1.451 0.00314 4.601 0.1978 0 62.35 35.25 -149.7 > 1000 5.885 1.994 2.999 0.01909 4.521 0.03328 64 26.68 56.73 161.6