First calculation of cosmic-ray muon spallation backgrounds for MeV astrophysical neutrino signals in Super-Kamiokande

# First calculation of cosmic-ray muon spallation backgrounds for MeV astrophysical neutrino signals in Super-Kamiokande

## Abstract

When muons travel through matter, their energy losses lead to nuclear breakup (“spallation”) processes. The delayed decays of unstable daughter nuclei produced by cosmic-ray muons are important backgrounds for low-energy astrophysical neutrino experiments, e.g., those seeking to detect solar neutrino or diffuse supernova neutrino background (DSNB) signals. Even though Super-Kamiokande has strong general cuts to reduce these spallation-induced backgrounds, the remaining rate before additional cuts for specific signals is much larger than the signal rates for kinetic energies of about 6 – 18 MeV. Surprisingly, there is no published calculation of the production and properties of these backgrounds in water, though there are such studies for scintillator. Using the simulation code FLUKA and theoretical insights, we detail how muons lose energy in water, produce secondary particles, how and where these secondaries produce isotopes, and the properties of the backgrounds from their decays. We reproduce Super-Kamiokande measurements of the total background to within a factor of 2, which is good given that the isotope yields vary by orders of magnitude and that some details of the experiment are unknown to us at this level. Our results break aggregate data into component isotopes, reveal their separate production mechanisms, and preserve correlations between them. We outline how to implement more effective background rejection techniques using this information. Reducing backgrounds in solar and DSNB studies by even a factor of a few could help lead to important new discoveries.

###### pacs:
25.30.Mr, 25.40.Sc, 24.10.Lx, 26.65.+t

## I Introduction

Neutrinos are powerful probes of the universe and its contents. They are abundantly produced by nuclear fusion processes that convert protons into neutrons, through the decays of unstable particles and nuclei created in high-energy processes, and through pair production in hot, dense environments. They can reach us unattenuated and undeflected from vast distances or from behind enormous column densities of matter, directly revealing the energies and timescales of the processes that made them. Even in a core-collapse supernova, where the neutrinos are thermalized by scattering, they emerge at energies MeV over about 10 s, compared to photons, which emerge at energies eV over months. The detection of astrophysical neutrinos allows us to probe physical conditions and neutrino properties beyond the reach of laboratory experiments.

The first great challenge of neutrino astronomy is the fact that the small interaction cross sections that make the above possible make detection difficult. This can only be solved by brute force — building large enough detectors to ensure adequate event rates. We focus on Super-Kamiokande (Super-K), the world’s largest low-energy neutrino detector, which has a fiducial mass of 22.5 kton of water and a total mass of 50 kton of water Fukuda et al. (2003); Abe et al. (2014). (For comparison, neutrinos were first detected in the Reines-Cowan reactor experiment with a detector using less than 1 ton of scintillator Cowan et al. (1956).) Even with such a large detector, the measured rates of low-energy astrophysical neutrinos are very small: about 15 solar neutrino events (all flavors of neutrinos elastically scattering electrons) detected per day Hosaka et al. (2006); Cravens et al. (2008); Abe et al. (2011) and an upper limit of several events (primarily inverse beta decay) detected per year from the diffuse supernova neutrino background (DSNB) Malek et al. (2003); Horiuchi et al. (2009); Beacom (2010); Bays et al. (2012); Zhang et al. ().

The second and far greater challenge of neutrino astronomy is reducing detector backgrounds to isolate these rare signals. Immense care and sophistication is required, and continual progress with existing detectors is possible. The primary backgrounds for solar and DSNB signals are MeV electrons and positrons from the decays of nuclei and muons. Below about 6 MeV detected electron kinetic energy, intrinsic radioactivities are the dominant background in Super-K Gando et al. (2003); Hosaka et al. (2006); Cravens et al. (2008); Abe et al. (2011), and these are controlled through the selection and purification of materials, choice of water circulation pattern to minimize radon ingress, and software processing (e.g., reconstruction quality and fiducial volume cuts). From about 6 to 18 MeV kinetic energy, induced radioactivities produced by cosmic-ray muons are the dominant background Gando et al. (2003); Hosaka et al. (2006); Cravens et al. (2008); Abe et al. (2011), and there is great potential to reduce these with the help of theoretical work.

To reduce cosmic-ray backgrounds, Super-K was built under 1000 m of rock (2700 m water equivalent) in the Kamioka mine in Japan Fukuda et al. (2003); Abe et al. (2014). As cosmic-ray particles interact with the rock and lose energy, their flux is reduced. The only high-energy particles that reach the Super-K detector are muons and neutrinos. The muon flux is m hr at sea level, and is reduced to 9.6 m hr at Super-K Galbiati and Beacom (2005), which corresponds to a muon rate in the detector of about 2 Hz Hosaka et al. (2006). It is easy to veto the muons themselves, but they frequently produce relatively long-lived radioactive isotopes through the breakup (“spallation”) of stable nuclei directly or, more commonly, through secondary particles produced through muon energy-loss processes. The spallation rate is large, interaction per through-going muon in Super-K, though many of the daughter nuclei are stable or decay in ways that do not produce Cherenkov signals.

Super-K has cuts to reduce backgrounds from the decays of spallation products, but these have to be limited to not overly discard signal events. Many of the unstable isotopes produced have half-lives of order 1 s, comparable to the time between successive muons. It is easy to estimate that a simple cut of all events in a cylinder of radius even a few meters around each muon track for a few seconds leads to a detector deadtime of . The real algorithm used by Super-K is more complex, and is based on a likelihood analysis that takes into account distance and time from the preceding muon as well as a variable related to muon energy loss, but a similar deadtime is achieved Koshio (1998); Hosaka et al. (2006). Even though the Super-K spallation cuts have a rejection efficiency of  Gando et al. (2003), the remaining background rate is still times greater than the solar neutrino signal rate above several MeV (this is then reduced by another factor 10 by the solar direction cut, leaving a background comparable to the signal) Abe et al. (2011). For the DSNB search, a higher energy threshold can be used to dramatically reduce backgrounds, but spallation decays are still overwhelming below about 18 MeV Malek et al. (2003) (16 MeV with new techniques Bays et al. (2012)).

Our goal for this paper is to detail the production processes for spallation backgrounds in Super-K and the physical characteristics of where, when, and with what associated particles these decays occur. With this information, it will be possible to make better cuts to reject backgrounds while preserving signals. For solar neutrinos, such improvements could help improve the significance of the 2.7- hint of the day-night effect from neutrino mixing in Earth Liu et al. (1997); Maris and Petcov (1997); Blennow et al. (2004); Renshaw et al. (2014). They may also help lead to the first detection of the hep neutrino flux, which is likely only a factor of a few away from detection Marcucci et al. (2000); Bahcall et al. (2001); Fukuda et al. (2001); Bahcall and Pinsonneault (2004); Hosaka et al. (2006). Such measurements would improve our knowledge of the Sun and of neutrino mixing parameters Fogli et al. (2011); Haxton et al. (2013); Bellini et al. (). Reduction of spallation backgrounds would also help lower the energy threshold in the DSNB search Malek et al. (2003); Bays et al. (2012) to where the signal is larger Horiuchi et al. (2009); Beacom (2010), which might help lead to a first detection.

Until now, there has been no detailed published study of spallation backgrounds in water. The Super-K cuts have been developed from empirical studies Hosaka et al. (2006); Cravens et al. (2008); Abe et al. (2011); Bays et al. (2012), and not from theoretical calculations. Further, they treat all isotopes together, without taking into account significant differences in their production, properties, and distributions. With Super-K nearly reaching the sensitivity needed for the above discoveries, a more detailed approach is needed. The interactions of muons with scintillator have been studied extensively with underground Aglietta et al. (1989); Hertenberger et al. (1995); Boehm et al. (2000); Menghetti (2005) and accelerator Hagner et al. (2000); Chazal et al. (2002); Lin (2013) experiments, and measurements like these have been incorporated into simulation packages like FLUKA Ferrari et al. (); Battistoni et al. (2007) and GEANT4 Agostinelli et al. (2003); Allison et al. (2006). This gives an opportunity to check our work and to understand the expected uncertainties of the simulations.

This paper is not meant to be a comprehensive study. It is a first step in understanding spallation backgrounds in water-based detectors, beginning with the yields and the average physical distributions of secondaries and isotopes. In two subsequent papers, we will go further, showing how characteristics of the showers of secondary particles that produce isotopes can be used to tailor better cuts Bays et al. (2012) and how those would be improved if Super-K gained the ability to detect neutrons by adding dissolved gadolinium Beacom and Vagins (2004).

This paper is organized as follows. In Sec. II, we describe the setup for our simulation. In Sec. III, general points about muon energy loss and secondary particle production are discussed. Our main results are in Sec. IV, where we calculate the neutron and isotope yields and study the properties of the induced backgrounds. Finally, we present our conclusions in Sec. V.

## Ii Setup of calculations

The Monte Carlo code FLUKA (version 2011.2b.3) Ferrari et al. (); Battistoni et al. (2007) is used for this work. It is a comprehensive code for particle energy loss and interactions with matter. For our purposes, FLUKA simulates all the physics processes relevant for the interactions of muons and their secondaries with water, including electromagnetic processes such as charged-particle ionization and bremsstrahlung, gamma-ray pair production and Compton scattering, and hadronic processes such as pion production and interactions, photo-disintegration, and low-energy neutron interactions with nuclei. It has been extensively used to simulate muon interactions in underground detectors, e.g., Refs. Wang et al. (2001); Kudryavtsev et al. (2003); Galbiati and Beacom (2005); Mei and Hime (2006); Abe et al. (2010); Bellini et al. (2013). The FLAIR interface Vlachoudis () is used when running FLUKA.

Most of the relevant physics processes and libraries are included in the FLUKA defaults. To make the low-energy neutron treatment more straightforward, the PRECISIOn card was chosen. Some muon processes, such as photo-nuclear and bremsstrahlung, were specifically activated. The new ion transport library was used.

The first main input for our simulation is the detector setup. The Super-K detector is a cylinder of water of diameter 39.3 m and height 41.4 m Hosaka et al. (2006). The outer detector (OD) is separated from the inner detector (ID) by a layer of photomultiplier tubes, most inward-facing, some outward-facing. The ID is about 2.5 m away from the edge of the detector Fukuda et al. (2003); Abe et al. (2014). Our results are calculated only in the fiducial volume (FV) region, which is a virtual cylinder with each side 2 m away from the ID (and about 4.5 m from the outer edge of the OD), containing 22.5 kton of water Hosaka et al. (2006). Water is one of the FLUKA pre-defined materials, including the natural abundances of hydrogen and oxygen isotopes. Muons may also interact with the surrounding rock to produce showers that enter the detector and produce isotopes. In the geometry setup, we include 2 m of rock outside the detector to induce secondary production (see Refs. Zbiri (2010); Empl et al. ()), though it has only a modest effect.

The other main input for our simulation is the muon energy spectrum shown in Fig. 1. The curve is the simulated muon flux at Super-K Tang et al. (2006). Because the muon energy is plotted on a log scale, the flux is plotted as Ed/dE 2.3d/dlogE, so that the integrated number of particles per decade (or other interval of fixed multiplicative width) is proportional to the value of this curve (i.e., plotting just d/dE underweights the importance of high-energy bins). The two vertical lines indicate characteristic energies. The one near 6 GeV is the minimum ionization energy loss for muons passing vertically through Super-K. Muons with less energy stop in the detector (as shown in the figure, these are only 5% of all muons). The line near 1000 GeV is the muon critical energy, at which the radiative energy loss equals the ionization energy loss. Muons with higher energies are more likely to produce showers, and thus more isotopes.

By number, most muons are in the range 30 – 700 GeV, with an average energy of 271 GeV Tang et al. (2006). The spectrum drops at high energies due to the falling spectrum of cosmic rays and at low energies due to muon energy loss in the rock above Super-K. Integration of the spectrum gives a muon rate at Super-K of 1.8 Hz Tang et al. (2006), which is consistent with the published values of 2 – 3 Hz Koshio (1998); Fukuda et al. (1999); Habig (); Gando et al. (2003). Specifying the muon rate more precisely requires knowing unpublished details about the muon multiplicity, path length and angular distributions, and stopping fraction. Other studies have shown that the detailed shape of the spectrum, for the same average energy, does not affect the isotope yield much Kudryavtsev et al. (2003); Abe et al. (2010).

We adopt several simplifications for the primary muons. All muons in our simulation are vertically down-going. In reality, most muons are down-going, but not perfectly Guillian et al. (2007); Tang et alTang et al. (2006) show that about 75 of muons have down-going zenith angle for KamLAND, which is at the same depth and location as Super-K. A complete 2D map of the simulated angular distribution of muons at Super-K is given in Ref. Tang et al. (2006); Guillian et al. (2007). Muons are sent only along the cylinder center. These two simplifications do not affect our results. Super-K has very good reconstruction for muon tracks, and all our secondary and isotope yields are calculated per muon path length. For muons coming in at an angle or a different spot, it would be easy to rescale our results by the actual muon track length. Besides single through-going muons, there are also muon bundles and muons that only go though a detector corner. We focus on single through-going muons, because they are the most common and because the other cases are easily identifiable. We simulate only ; there are also , but the isotope yields from and differ very little Abe et al. (2010); Bellini et al. (2013), except for nuclear captures of stopping , which we discuss below.

A similar setup was adopted for the spallation study by KamLAND Abe et al. (2010). In their study, spallation yields were measured experimentally and compared to simulation results from FLUKA. The Borexino spallation study Bellini et al. (2013) used both simulation packages FLUKA and GEANT4. Overall, it was found that there are factor of 2 discrepancies between the calculated yields and also between those and the measured values, which is reasonable, given the hadronic uncertainties and that yields for different isotopes vary by orders of magnitude.

## Iii Muon energy loss and secondary production

The average muon energy loss rate is Lipari and Stanev (1991); Dutta et al. (2001); Groom et al. (2001); Bulmahn (2010); Beringer et al. (2012)

 dEdx=α(E)+β(E)E. (1)

The term corresponds to the continuous energy losses due to the ionization (and excitation) of atomic electrons. It has a typical value of 2 MeV cm g and does not change much with muon energy. The ionization can be separated into a restricted ionization energy loss, which is the ionization with soft collisions and small fluctuations, and delta-ray production, which has hard collisions and large fluctuations Beringer et al. (2012). The term corresponds to the energy losses due to radiative processes through interactions with atomic nuclei. For muons at hundreds of GeV, pair production and bremsstrahlung are the most important radiative processes, while photo-nuclear has a small contribution Groom et al. (2001). Pair production is a nearly continuous energy loss, but bremsstrahlung and photo-nuclear energy losses have large fluctuations. Ionization and radiation losses are equal at about 1000 GeV for muons in water, which defines the muon critical energy  Groom et al. (2001).

Figure 2 shows the energy loss distribution for vertical (path length 32.2 m) through-going muons in the Super-K FV. The restricted ionization energy loss is about 6 GeV and the pair production loss is about 1 GeV. These two terms have almost no fluctuations and correspond to the minimum energy loss of 7 GeV shown in Fig. 2. On average, muons lose about 11 GeV, which means 4 GeV for the total of the delta-ray production, bremsstrahlung, and photo-nuclear processes. Bremsstrahlung energy loss is primarily responsible for producing the high energy loss tail Beringer et al. (2012).

Muons lose energy to the production of secondary particles, and there is a lot of energy available to make many of them, as shown in Fig. 2. These interactions do not appreciably affect the parent muon, as the energy loss in the detector is small compared to the muon energy. The muon interaction cross sections then do not change much as muons lose energy traveling through the detector Groom et al. (2001). The muon tracks have only minor deflections, with 90 of muons having less than 30 cm transverse displacement when they exit the FV.

Figure 3 shows the average production of secondaries by muons in Super-K. The plotted path length spectrum is the sum of distances traveled by all secondary particles of the same species at certain energy. It is similar to the particle multiplicity times the mean free path. The difference is that here a particle contributes to the path length at low energies after it travels some distance at high energies, so there is a pileup of path length from high energy to low energy. This path length spectrum is the most useful quantity for calculating interactions by these particles. These results do not depend on density because they are calculated per muon path length (here the vertical distance through the Super-K FV).

As shown in Fig. 3, the dominant secondaries are gammas, followed by electrons (and positrons). This makes sense because the primary ways for muons to lose energy other than ionization are delta-ray production, pair production, and bremsstrahlung, all of which are electromagnetic. In Fig. 2, the average radiative muon energy loss is 5 GeV. The accumulated path length of the secondary electrons and positrons should be 5 GeV / (0.2 GeV/m) 25 m, and the integral of their curve in Fig. 3 is close to this.

A similar figure in Ref. Galbiati and Beacom (2005), which is based on independent calculations, shows secondaries produced by muon interactions in scintillator. Detailed comparison between Fig. 3 and Ref. Galbiati and Beacom (2005) (taking into account the different plotting scales) shows consistent results. As expected, there is not much difference between muon interactions in water or scintillator for muon energies of hundreds of GeV. A minor discrepancy is that there are more than in Fig. 3, whereas it is the opposite in Ref. Galbiati and Beacom (2005). To check this, we ran a separate simulation without hydrogen and found that the slight difference in our Fig. 3 between and is due to scattering of on free (hydrogen) protons. Our best guess is that the and curves in the figure of Ref. Galbiati and Beacom (2005) are mislabeled.

All of the results presented here are averaged over many muon path lengths. In fact, secondaries are made primarily in electromagnetic and hadronic showers, not uniformly along muon tracks. In our simulation runs, we see significant correlated variations in the muon energy loss, secondary production, and isotope production along the muon paths. This is hinted at by the high particle energies in Fig. 3. In our follow-up papers, we will discuss the shower nature of secondary production and how taking it into account can help improve background rejection in Super-K.

Muons interact with oxygen nuclei directly to produce isotopes, but the dominant mechanism to make isotopes is through secondaries breaking up oxygen nuclei. The most important secondaries in this regard are neutrons, pions, and gammas. Of all spallation-induced isotopes that cause backgrounds in Super-K, only 11 are made by muons (7 are N from stopping muons plus 4 other isotopes); the rest are made by secondary particles.

The physical distributions of the secondaries tell us where the isotopes are being made. The differences reflect how the different secondaries lose energy. Figure 4 shows the normalized distribution of secondary particle absorption distances to the muon track. The distribution is dN/dr [cm], i.e., the area factor d is included. Compared to Fig. 3, electrons (and protons) are not shown because they are not major parent particles for spallation products. The gammas have a short mean free path and are mostly forward. Most gammas are destroyed by pair production, and the Moliere radius (9.8 cm in water Beringer et al. (2012)) sets a scale for gamma distances from the muon. The mean free path for pions at these energies is about 1 m Beringer et al. (2012). Assuming pions are destroyed after only one interaction (e.g., absorption on p), the falling distribution corresponds to a typical forward direction of . This is consistent with Fig. 3, where most pions are relativistic. Among muon secondaries, neutrons travel the furthest from the muon track, with 98 of neutrons contained within 3 m. The neutron mean free path is 10 cm above a few MeV, and less at lower energies; neutrons go much farther than this because many scatterings are required to stop them McLane et al. (1988). The result is very similar to the neutron distance distribution in scintillator Bellini et al. (2013). The carbon number density in scintillator and the oxygen number density in water are comparable, but the cross section for neutrons on oxygen is slightly higher than that on carbon Chadwick et al. (2011). As a result, neutrons travel a bit less far in water. Compared to the average distance of 74 cm in water, the average distance in scintillator is 81.5 cm Bellini et al. (2013). Most neutrons are absorbed by capture on hydrogen at non-relativistic energies; we also count the reactions of energetic neutrons on oxygen, e.g., (n,p), though this is a small effect. The Borexino Bellini et al. (2013) measurement counts only gamma-ray producing captures on hydrogen (mostly) and carbon.

## Iv Isotope and neutron production and distributions

Using the muon and secondary data, we calculate the isotope and neutron yields in Super-K using FLUKA. The isotope counts are read from the RESNUCLEi card. Neutron counts and production channels are taken from a modified mgdraw.f subroutine. For neutron counts, processes like (n, 2n) are carefully taken into account.

We began our study by reproducing all of the relevant KamLAND results Abe et al. (2010), and extending the isotope yields to include stable isotopes for comparison to the yields of analogous (stable or unstable) nuclei in Super-K. Consistent results, within a factor of 2, validate our approach. The results show interesting differences in the physics of spallation in water and scintillator, as discussed in detail below.

### iv.1 Predicted Yields

Table 1 shows the neutron and isotope yields per muon along with associated details. Almost all isotopes made by muons and their secondaries are listed (we skip isotopes with small yields or small mass numbers). Since Super-K can only detect relativistic charged particles, only betas and gammas (through pair production or Compton scattering) can be seen, while decay products such as neutrons, protons, and alpha particles are invisible (neutron captures on protons are very hard to detect Zhang et al. ()). The top part of the table contains isotopes that decay and thus are backgrounds in Super-K (referred to as background isotopes); the bottom part of the table contains isotopes that are stable, have long half-lives, or decay invisibly.

The half-lives of the unstable isotopes range greatly, from 0.008 s to 13.8 s. A timescale to compare to is the average separation between muons, about 0.5 s. The beta decay spectra are complicated and have various branches. Here only the dominant decay modes are listed, though our calculations take all modes into account. Unsurprisingly, many of the spallation isotopes are short-lived and high-energy compared to intrinsic radioactivities. The half-lives and decay modes are taken from National Nuclear Data Center (). The isotope decay spectra are taken from Ref. Blaufuss et al. (2001) for N, Ref. Winter et al. (2006) for B, and Ref. Tuli () for all other isotopes.

The fourth column shows the isotope yields calculated with FLUKA. These span five orders of magnitude, which is an important point. As noted, the accuracy of the isotope production rate is only about a factor of 2. Yet, because the yields among different isotopes are so different, we can still get a good understanding of their relative importance. Another point is that the production of beta-decaying isotopes is relatively rare. The sum of unstable isotopes is 58 in the units of the table, corresponding to about 0.02 unstable isotopes per muon (i.e., multiplying by the vertical distance of 3220 cm). The sum of the stable or invisible isotopes is around 2950, or about 0.9 isotopes per muon. Neutrons are produced with a yield comparable to that of all isotopes.

The current Super-K solar neutrino analysis has a kinetic energy threshold of 3.5 MeV Sekiya (), and taking this into account changes the importance of different isotopes. The fifth column shows the production rate of isotopes with decay energy larger than 3.5 MeV. Of unstable isotopes with high yields, N is cut the least. For N decay, 66% of the time there is a 6.1 MeV gamma ray, which leads to an electron-equivalent energy reduced by a factor 1/4 Blaufuss et al. (2001). As a result, the beta spectrum is shifted to higher energies, making it unaffected by the 3.5 MeV cut. The sum of the yields of background isotopes is reduced to 50 in the units of the table.

The last column shows the most important production channel for each isotope. For most isotopes, there are several production channels, with different parent particles, often of comparable importance. Statistically, the assignments of parent particles in the simulation are correct. For low energy neutrons ( 20 MeV), FLUKA uses a multi-group treatment, so the correlations among daughter particles are not accurate. In cases where production by neutrons is important, the results provide a good first understanding, but are not accurate descriptions of the actual interactions.

The final states of the production channels for each isotope indicate particles that could possibly be detected in association with creation of the isotope. (In addition, there will frequently be prompt gamma rays from the de-excitation of daughter nuclei Langanke et al. (1996); Nussinov and Shrock (2001); Kolbe et al. (2002); Ankowski et al. (2012), but the Cherenkov light from their subsequent signals will be buried under that from the muon.) It may be possible to identify pion decays in some cases. Protons and alpha particles will almost always be non-relativistic and hence non-detectable. At present, it is very difficult to detect neutrons in Super-K Zhang et al. (), though that would change with the addition of gadolinium Beacom and Vagins (2004); neutron captures are prompt (about 200 s in pure water and about 10 times shorter if gadolinium is added), so they are efficiently removed by even a short time cut following a muon. An important application could be identifying the production of He and Li, the decays of which can mimic an astrophysical inverse beta signal because there is a beta followed by a neutron capture. We find that there is frequently a neutron produced in association with these isotopes, so there would be a neutron capture preceding the He or Li decay, unlike for a real astrophysical signal event. However, we caution that further study of the contributing channels is needed.

The production of N was independently calculated in Ref. Galbiati and Beacom (2005). This is the most abundant background isotope from muons and it has a long half-life. The dominant way to make N is , which has a yield of g cm, to be compared to the value found by Ref. Galbiati and Beacom (2005), g cm. Sudbury Neutrino Observatory has an upper limit on the N yield of 20 – g cm Aharmim et al. (2005). All of these are consistent.

In our simulation, we consider only primary . Other studies have shown that isotope production by and typically differs by only a few percent Abe et al. (2010). One exception is stopping , which can capture on oxygen and make N by . Stopping make 17% of N, for which Super-K has a separate cut Hosaka et al. (2006). Consequently, if we take primary into account, the N yield would change by about 8%. For most subsequent calculations and comparisons to Super-K measurements, we ignore the correction to isotope production.

### iv.2 Comparison to Super-K Measurements

In the following, we focus our comparisons on data above 6 MeV. At lower energies, detector backgrounds from intrinsic radioactivities are dominant. The largest intrinsic radioactivity background in the water itself is due to the Bi beta decay following Rn ingress; though its endpoint is 3.26 MeV, energy resolution smears the spectrum to higher energies Mitsuda et al. (2003); Hosaka et al. (2006); Abe et al. (2011) (see also Ref. Aharmim et al. (2013)). There are also radioactivities in the photomultiplier and other detector elements, and these are largely reduced through the fiducial volume cut Hosaka et al. (2006); Abe et al. (2011). This dividing line of 6 MeV is in good agreement with the demonstrated effectiveness of the spallation cut above this energy Hosaka et al. (2006); Cravens et al. (2008); Abe et al. (2011), as well as by the results of a dedicated spallation study Gando et al. (2003).

Super-K has given a likelihood function of decay time t after the primary muon for decays in a cylinder around the muon path Koshio (1998); Hosaka et al. (2006). This time is well defined because the muon takes only about 100 ns to cross the detector. The likelihood function is an empirical fit to the sum of all spallation backgrounds, and isotopes with similar half-lives are grouped together. With the simulated yields from FLUKA, we have each component of this separately.

Figure 5 (left panel) shows our combined spallation product decay rate compared to the Super-K fit. The normalization is chosen so that the integrated event numbers are the same between the simulation and the Super-K fit. Overall, the total decay rate and the Super-K fit agree well, up to a factor of 2. The four most abundant isotopes have very different half-lives. This figure shows how each contributes to the total decay rate on different timescales. Below about 0.1 s, B is dominant (with a smaller contribution from N, which has a comparable half-life and decay energy); between 0.1 s to 3 s, Li contributes most; and, after about 3 s, N is dominant. We also show Be, which has the longest half-life, 13.8 s. All of the curves in Fig. 5 (left panel) have a kinetic energy cut of 6 MeV.

Figure 5 (right panel) shows a similar result with a 10 MeV kinetic energy cut to the calculation (the similar Super-K measurement is not available). The main effect of the energy cut is to decrease N compared to other isotopes. A relatively high energy cut works well for N because of its low endpoint energy.

Another comparison we can make with Super-K results is the energy spectrum of spallation backgrounds in the FV. Similar to above, Super-K has the total decay energy spectrum from all background isotopes Gando et al. (2003). With the simulated yields, adding up the component spectra from all isotopes gives a total spectrum that can be compared to data.

Figure 6 (left panel) shows that the simulation and the measurement agree quite well above 6 MeV. For this comparison, the isotope yields were multiplied by the average muon rate at Super-K (1.88 Hz) and the average muon track length in the FV (32.2 m). Both numbers have uncertainties because we do not know the precise definitions used by Super-K. This, together with the limitations of the simulation, introduce the biggest uncertainties. Taking energy resolution into account is important: the high energy events seen in the detector are mainly from imperfectly reconstructed lower energy events. The agreement validates our results, especially because the absolute scale is predicted, not fit.

Figure 6 (right panel) shows the isotope spectra after a t 0.3 s cut, which is about an order of magnitude less than our estimate of the time needed for a simple cylinder cut around each muon (see Sec. I). This is chosen to be short enough to not introduce significant deadtime and long enough to eliminate many short-lived isotopes. The total spectrum decreases by about a factor of 2. It also affects the relative contributions of isotopes at different energies. The dominant component at high energy without a time cut is N; after 0.3 s time cut, it is B. The fewer isotopes that contribute, the more effective isotope-specific cuts will be (see below).

The Super-K DSNB analysis of Ref. Malek et al. (2003) has a lower energy threshold of 18 MeV total energy. The total background rate is 0.2 events per day in the 18 – 20 MeV energy bin. The rate in Fig. 6 is consistent because the measured data in Ref. Malek et al. (2003) include an increasing contribution from the decays of invisible muons.

The Super-K N calibration study reports that the production rate of N by stopping muons is 11 per day in an 11.5 kton volume Blaufuss et al. (2001). The rate from our calculation is g cm. Taking into account the fraction in primary muons and the detector efficiency, we predict 22 events per day. The origin of the discrepancy is unknown, but the Super-K study reported problems with their measurement Blaufuss et al. (2001), so we view this factor of 2 as adequate agreement.

The fact that our same FLUKA predictions match both the energy spectrum and the time profile of the Super-K data is a powerful indication that they are accurate. In the energy spectrum, the components are largely overlapping because of the width of the beta spectra and the effects of energy resolution smearing. In the time profile, the components are better separated because of the wide range of half-lives. In combination, these provide strong tests of both the overall production rate of spallation products and the amplitudes of the many components.

### iv.3 Comparison to Yields in Scintillator

A major difference between spallation in scintillator and water is in the absolute background isotope yield. It is 0.3 of the neutron yield in scintillator, whereas it is only 0.03 of the neutron yield in water. (In scintillator and water, the neutron yields are similar to each other.) The reason is that there is a greater fraction of stable or invisibly decaying isotopes produced by muons in water. The neutron number is comparable to the total yield of all isotopes, both in scintillator and water. It is about 0.7 neutron per muon in the Super-K FV.

The production channels allow us to understand the different spallation processes better. The isotope yields between scintillator and water are similar if the production mechanisms are similar. Some of the most abundant isotopes are made by the (,n) reaction, which corresponds to O in water and C in scintillator. They have yields of 351 and 416 Abe et al. (2010) in the units of g cm. Luckily, O has a low beta-decay energy; in scintillator, C is a serious background. The most abundant background isotope in water is N, which corresponds to B in scintillator, which has a comparable yield for the same muon path length.

### iv.4 Parent Particle Energy Spectrum

To understand isotope production mechanisms in more detail, we look at the energy spectra for secondaries making isotopes. Figure 7 (left panel) shows the spectra of parent particles of spallation background isotopes in Super-K. Here the y axis is a histogram of event number per MeV with arbitrary absolute normalization. The relative height reflects how important each parent particle is.

For making spallation backgrounds in Super-K, the most important parent particle is the neutron, as it contributes almost 10 times more than any others. The shape of the spectrum is a convolution of the neutron path length shown in Fig. 3 and the neutron-nucleus cross section. The peak below 20 MeV comes from the cross section Galbiati and Beacom (2005). Due to the nuclear capture of at rest, there is also a huge peak for low energy . Gamma, , and high energy contribute roughly equally, each only about half as much as the first bin. The parent particles of fast neutrons are similar to those for isotopes. Wang et alWang et al. (2001) showed that at GeV, most neutrons are produced by , followed by gamma and neutron.

One interesting feature is that, even though the dominant secondaries produced directly by muons are gammas and electrons, the ones that make background isotopes in water are mainly hadrons. This is consistent with the primary processes shown in Table 1. The fact that the gamma and pion curves initially rise with energy is consistent with the path length spectra in Fig. 3. The fact that these curves continue to high energies indicates the importance of showers for isotope production.

As discussed above, the result is somewhat different from muon spallation in scintillator. A rough count from the KamLAND result tells us that the main parent particle to produce isotopes is gamma, as it is responsible for C and Be production. This is consistent with the result shown in Fig. 7 (right panel). Here we show the parent particle spectra for all isotopes produced in water, including the stable ones and those that decay invisibly. The gamma contribution is significant, comparable to that of the neutron. Also, the relative height between the two panels shows the fraction of isotopes that are dangerous in Super-K relative to all isotopes. The reason for the big difference between the left and right panels is simply that in water, some of the most abundant isotopes made by gammas, e.g., O, N, and C, are invisible in Super-K.

### iv.5 Spatial Distribution of Isotopes

Because spallation products are produced by muons and their secondary particles, there are spatial and temporal correlations between spallation events and the parent muons. The muon itself emits Cherenkov light along its entire path, which makes it easy to detect. Thus, the correlations between muons and isotopes provide an opportunity for physics-motivated cuts.

There are two distances to describe the position of the isotope to the parent muon. One is the perpendicular distance to muon track, which is one of the variables for the Super-K likelihood function for the spallation cut. The other is the isotope position along the muon track.

Once isotopes are produced, they do not move far before they decay. Ions stop in a short distance, and there is no significant bulk motion of the water Abe et al. (2014). This can be seen from the fact that the Super-K likelihood function of isotope distance to the muon track shows a peak at very small distance Hosaka et al. (2006). Figure 8 shows our calculated distribution of isotope distance to the muon track. This shows one of the likelihood functions used for the Super-K spallation cuts. Our results are consistent with those shown in Ref. Koshio (1998) (the Super-K results depend on a variable associated with muon energy loss; we summed over those distributions with appropriate weights). We did not take the Super-K position resolution into account in Fig. 8; it is about 1 m at 5 MeV and about 0.5 m at 10 MeV Abe et al. (2011). We find that 99 of isotopes decay within 3 m.

Each isotope has a different distribution, and we show two examples. The most abundant background isotope is N, and it dominates the low-energy end of the spectrum. On the other hand, B contributes the most at the high-energy end, as shown in Fig. 6. These two isotopes have quite different distributions. The 90 containment distance for B is 1.7 times smaller than that for N, which corresponds to a factor of 3 in cylinder volume. Taking this into account could improve cuts and reduce deadtime. For example, at high decay energies, B but not N can contribute, so a more specific cut could be used.

Figure 8 shows useful features for improving the Super-K likelihood function for the spallation cut. In the Super-K current likelihood function, the isotope distance to the muon track is one variable. However, the distance distribution for each isotope can be appreciably different. As a result, instead of using a combined likelihood function for all background isotopes, a likelihood function for each isotope separately should give a much more accurate description of the physics.

If we consider the isotope distance along the muon track, to first order we would expect a flat distribution when we average over muons (for individual muons, this would have bumps due to showers). The reason is that, on average, muons have hundreds of GeV energy and lose only about 11 GeV during propagation through Super-K. More precisely, the isotope yields decrease smoothly from the top of FV to the bottom of FV by several percent. This is partly due to the decrease of muon energy, and partly due to stopping muons. There is negligible spillover from the rock above Super-K.

## V Conclusions and Future Work

Guided by theoretical understanding and analysis, we use the simulation package FLUKA to study muon interactions with water, the production and properties of secondary particles, and the production and decay of unstable isotopes. Where possible, we compare our results to published measurements from Super-K, finding good agreement on an absolute scale, i.e., a factor of 2, which is reasonable considering the orders of magnitude differences in production rates. The residual discrepancies primarily arise from uncertainties in hadronic interactions and unpublished details of the muon backgrounds, and some of the differences could be reduced by calibration to measured data.

As a check, we also performed similar calculations for scintillator-based detectors, for which there are more extensive theoretical studies and experimental measurements. We focus on comparison to isotope and neutron production in KamLAND Abe et al. (2010) and Borexino Bellini et al. (2013), finding good agreement, within the factors of 2 that have been noted by others between the measurements and calculations and also between calculations with FLUKA versus GEANT4 Wang et al. (2001); Kudryavtsev et al. (2003).

One interesting point for context is how different the spallation backgrounds are in water Cherenkov detectors compared to scintillator detectors. First, for water there is the fortunate point that, although the production rate of all isotopes is comparable to that in scintillator, that of unstable isotopes is about ten times less. Second, many of the unstable isotopes decay without producing Cherenkov light. Scintillator detectors have the ability to detect neutrons through their radiative captures, and this is a significant advantage in identifying spallation products. However, if Super-K adds dissolved gadolinium to enable the detection of neutron captures, it will have a similar capability Beacom and Vagins (2004).

Our calculations for Super-K lead to important new high-level results beyond the details presented here. First, a demonstration that a theoretical calculation of the spallation backgrounds in water is now possible, even though it was not when Super-K began Koshio (1998). Compared to an empirical approach, production mechanisms are revealed, aggregates are separated into components, and correlations are preserved. Second, we show details that were heretofore unavailable. Important examples are differences between the distributions and correlations of each isotope, including temporal distribution after the muon, distance distribution away from the muon, decay energy spectrum, and associated particles.

We demonstrate that there is more information to be gained by having likelihood functions of time and distance for each isotope. Instead of a global likelihood for all spallation decays, our results could be used to construct per-isotope likelihoods that would lead to more precise cuts. Also, a new variable of decay energy can be used in addition to its original three variables of decay distance to the muon track, decay time, and muon energy loss. Even modest improvements, say a factor of a few, could lead to significant gains in the ability to measure signals. This could help lead to first discoveries of the day-night effect and the hep flux in solar neutrinos, as well as the DSNB.

Our results are calculated for Super-K, but they could have wider applicability. The isotope yields per muon vary only moderately with depth, once that depth is appreciable, because they have a modest dependence on the muon average energy, scaling roughly as  Abe et al. (2010). As first estimates, our results would provide useful comparisons for the Sudbury Neutrino Observatory Aharmim et al. (2013), Hyper-Kamiokande Abe et al. (), and the water shields of a variety of neutrino and dark matter detectors.

It would be valuable for Super-K to produce a dedicated study on spallation backgrounds informed by the predictions of this paper. The yields of different isotopes could be identified by a global fit that takes into account the full energy and time information on spallation decays, e.g., energy spectra in different time ranges, as has been done for scintillator detectors Abe et al. (2010); Bellini et al. (2013). Another key observable is the radial distributions of isotopes produced by different types of secondaries. An improved FLUKA simulation could be developed using a more complete description of the detector details, especially the muon distributions. It seems likely that the uncertainties could be reduced to well below a factor of 2 by calibrating the simulation to measured data.

It would also be valuable to have a similar study for the Sudbury Neutrino Observatory Aharmim et al. (2013). The very low muon rate and intrinsic radioactivities would make it easier to identify spallation decays and to avoid confusion over which muon was the parent. In addition, the ability to detect neutrons would help identify isotope production channels. With corrections for the different muon spectrum, detector properties, and the production of neutrons by deuterium photo-disintegration, it would be straightforward to relate these measurements to Super-K results.

In two follow-up papers, we will develop further ways to reduce backgrounds in Super-K and other water Cherenkov detectors. In the first paper, we will study the variations in muon energy loss along the path due to showers, and how this can be used to identify where isotopes are produced. This effect was discovered empirically in Ref. Bays et al. (2012), and our results will provide the first detailed explanation of how it works and how it could be improved. In the second paper, we will show how the ability to detect neutrons using gadolinium in water, as first suggested in Ref. Beacom and Vagins (2004), can be used to improve cuts to reduce spallation backgrounds. These papers will include some surprises that will allow significant gains in sensitivity beyond those enabled by results given here.

## Vi Acknowledgments

S.W.L. and J.F.B. were supported by NSF Grant No. PHY-1101216 to JFB. We thank Sheldon Campbell, Mark Chen, Anton Empl, Cristiano Galbiati, Yusuke Koshio, Pablo Mosteiro, Masayuki Nakahata, Kenny Ng, Itaru Shimizu, Michael Smy, Mark Vagins, and Lindley Winslow for helpful discussions.

### References

1. Y. Fukuda et al. (Super-Kamiokande Collaboration), Nucl. Instrum. Meth. A501, 418 (2003).
2. K. Abe et al.Nucl. Instrum. Meth. A 737, 253 (2014).
3. C. L. Cowan, F. Reines, F. B. Harrison, H. W. Kruse,  and A. D. McGuire, Science 124, 103 (1956).
4. J. Hosaka et al. (Super-Kamiokande Collaboration), Phys. Rev. D 73, 112001 (2006).
5. J. P. Cravens et al. (Super-Kamiokande Collaboration), Phys. Rev. D 78, 032002 (2008).
6. K. Abe et al. (Super-Kamiokande Collaboration), Phys. Rev. D 83, 052010 (2011).
7. M. Malek et al. (Super-Kamiokande Collaboration), Phys. Rev. Lett. 90, 061101 (2003).
8. S. Horiuchi, J. F. Beacom,  and E. Dwek, Phys. Rev. D 79, 083013 (2009).
9. J. F. Beacom, Ann. Rev. Nucl. Part. Sci. 60, 439 (2010).
10. K. Bays et al. (Super-Kamiokande Collaboration), Phys. Rev. D85, 052007 (2012).
11. H. Zhang et al. (Super-Kamiokande Collaboration),  arXiv:1311.3738 [hep-ex] .
12. Y. Gando et al. (Super-Kamiokande Collaboration), Phys. Rev. Lett. 90, 171302 (2003).
13. C. Galbiati and J. F. Beacom, Phys. Rev. C 72, 025807 (2005).
14. Y. Koshio, Study of Solar Neutrinos at Super KamiokandePh.D. thesis, University of Tokyo (1998).
15. Q. Y. Liu, M. Maris,  and S. T. Petcov, Phys. Rev. D 56, 5991 (1997).
16. M. Maris and S. T. Petcov, Phys. Rev. D 56, 7444 (1997).
17. M. Blennow, T. Ohlsson,  and H. Snellman, Phys. Rev. D 69, 073006 (2004).
18. A. Renshaw et al. (Super-Kamiokande Collaboration), Phys. Rev. Lett. 112, 091805 (2014).
19. L. E. Marcucci, R. Schiavilla, M. Viviani, A. Kievsky, S. Rosati,  and J. F. Beacom, Phys. Rev. C 63, 015801 (2000).
20. J. N. Bahcall, M. Pinsonneault,  and S. Basu, Astrophys. J. 555, 990 (2001).
21. S. Fukuda et al. (Super-Kamiokande Collaboration), Phys. Rev. Lett. 86, 5651 (2001).
22. J. N. Bahcall and M. H. Pinsonneault, Phys. Rev. Lett. 92, 121301 (2004).
23. G. L. Fogli, E. Lisi, A. Marrone, A. Palazzo,  and A. M. Rotunno, Phys. Rev. D 84, 053007 (2011).
24. W. C. Haxton, R. G. Hamish Robertson,  and A. M. Serenelli, Ann. Rev. Astron. Astrophys. 51, 21 (2013).
25. G. Bellini et al. (Borexino Collaboration),  arXiv:1308.0443 [hep-ex] .
26. M. Aglietta et al.Nuovo Cimento Soc. Ital. Fis., C 12, 467 (1989).
27. R. Hertenberger, M. Chen,  and B. L. Dougherty, Phys. Rev. C 52, 3449 (1995).
28. F. Boehm, J. Busenitz, B. Cook, G. Gratta, H. Henriksen, J. Kornis, D. Lawrence, K. B. Lee, K. McKinny, L. Miller, V. Novikov, A. Piepke, B. Ritchie, D. Tracy, P. Vogel, Y.-F. Wang,  and J. Wolf, Phys. Rev. D 62, 092005 (2000).
29. H. Menghetti (LVD Collaboration), AIP Conf. Proc. 785, 259 (2005).
30. T. Hagner, R. von Hentig, B. Heisinger, L. Oberauer, S. Schonert, F. von Feilitzsch,  and E. Nolte, Astropart. Phys. 14, 33 (2000).
31. V. Chazal, F. Boehm, B. Cook, H. Henrikson, G. Jonkmans, A. Paic, N. Mascarenhas, P. Vogel,  and J.-L. Vuilleumier, Nucl. Instrum. Meth. A 490, 334 (2002).
32. C.-J. Lin, “Proposal to Measure Muon-induced Neutron Background at CERN,” https://zzz.physics.umn.edu/lowrad/_media/lin_cern-neutronprod.pdf (2013), [AARM workshop].
33. A. Ferrari, P. R. Sala, A. Fasso,  and J. Ranft, “FLUKA: A multi-particle transport code,” CERN-2005-10 (2005), INFN-TC-05-11, SLAC-R-773.
34. G. Battistoni et al.AIP Conf. Proc. 896, 31 (2007).
35. S. Agostinelli et al.Nucl. Instrum. Meth. A 506, 250 (2003).
36. J. Allison et al.IEEE Trans. Nucl. Sci. 53, 270 (2006).
37. J. F. Beacom and M. R. Vagins, Phys. Rev. Lett. 93, 171101 (2004).
38. Y.-F. Wang, V. Balic, G. Gratta, A. Fassò, S. Roesler,  and A. Ferrari, Phys. Rev. D 64, 013012 (2001).
39. V. Kudryavtsev, N. Spooner,  and J. McMillan, Nucl. Instrum. Meth. A 505, 688 (2003).
40. D.-M. Mei and A. Hime, Phys. Rev. D 73, 053004 (2006).
41. S. Abe et al. (KamLAND Collaboration), Phys. Rev. C 81, 025807 (2010).
42. G. Bellini et al. (Borexino), J. Cosmol. Astropart. Phys. 1308, 049 (2013).
43. V. Vlachoudis, “FLAIR: A Powerful But User Friendly Graphical Interface For FLUKA,” Proc. Int. Conf. on Mathematics, Computational Methods and Reactor Physics, Saratoga Springs, New York, 2009.
44. K. Zbiri, Nucl. Instrum. Meth. A 615, 220 (2010).
45. A. Empl, R. Jasim, E. Hungerford,  and P. Mosteiro,  arXiv:1210.2708 [astro-ph.IM] .
46. A. Tang, G. Horton-Smith, V. A. Kudryavtsev,  and A. Tonazzo, Phys. Rev. D 74, 053007 (2006).
47. Y. Fukuda et al. (Super-Kamiokande Collaboration), Phys. Rev. Lett. 82, 2644 (1999).
48. A. Habig (Super-Kamiokande Collaboration),  arXiv:hep-ex/0106024 [hep-ex] .
49. G. Guillian et al. (Super-Kamiokande Collaboration), Phys. Rev. D 75, 062003 (2007).
50. P. Lipari and T. Stanev, Phys. Rev. D 44, 3543 (1991).
51. S. I. Dutta, M. H. Reno, I. Sarcevic,  and D. Seckel, Phys. Rev. D 63, 094020 (2001).
52. D. E. Groom, N. V. Mokhov,  and S. I. Striganov, Atom. Data Nucl. Data Tabl. 78, 183 (2001).
53. A. Bulmahn, Lepton pair production in high energy A scattering: cross sections, energy loss, and applications to underground lepton productionPh.D. thesis, University of Iowa (2010).
54. J. Beringer et al. (Particle Data Group), Phys. Rev. D 86, 010001 (2012).
55. V. McLane, C. L. Dunford,  and P. F. Rose, Neutron Cross Sections, Vol. 2 (Academic Press Inc, 1988).
56. M. B. Chadwick et al.Nuclear Data Sheets 112, 2887 (2011), Special Issue on ENDF/B-VII.1 Library.
57. E. Blaufuss et al. (Super-Kamiokande Collaboration), Nucl. Instrum. Meth. A 458, 638 (2001).
58. National Nuclear Data Center, “NuDat 2 database,” http://www.nndc.bnl.gov/nudat2/, [Online version 2.6].
59. W. T. Winter, S. J. Freedman, K. E. Rehm,  and J. P. Schiffer, Phys. Rev. C 73, 025503 (2006).
60. J. Tuli, “Evaluated Nuclear Structure Data File (ENSDF),” http://www.nndc.bnl.gov/ensdf/, [Online; accessed 2013-09-08].
61. H. Sekiya (Super-Kamiokande Collaboration),  arXiv:1307.3686 [astro-ph.SR] .
62. K. Langanke, P. Vogel,  and E. Kolbe, Phys. Rev. Lett. 76, 2629 (1996).
63. S. Nussinov and R. Shrock, Phys. Rev. Lett. 86, 2223 (2001).
64. E. Kolbe, K. Langanke,  and P. Vogel, Phys. Rev. D 66, 013007 (2002).
65. A. M. Ankowski, O. Benhar, T. Mori, R. Yamaguchi,  and M. Sakuda, Phys. Rev. Lett. 108, 052505 (2012).
66. B. Aharmim et al. (SNO Collaboration), Phys. Rev. C 72, 055502 (2005).
67. C. Mitsuda, T. Kajita, K. Miyano, S. Moriyama, M. Nakahata, Y. Takeuchi,  and S. Tasaka, Nucl. Instrum. Meth. A 497, 414 (2003).
68. B. Aharmim et al. (SNO Collaboration), Phys. Rev. C 88, 025501 (2013).
69. K. Abe et al.,  arXiv:1109.3262 [hep-ex] .
180902