Showering Cosmogenic Muons in A Large Liquid Scintillator

We present the results of FLUKA simulations of the propagation of cosmogenic muons in a 20 kton spherical liquid scintillator detector underneath 700 to 900 meters of rock. A showering muon is one which deposits at least 3 GeV in the detector in addition to ionization energy. We find that 20 percent of muons are showering and a further 10 percent of muon events are muon bundles of which more than one muon enters the detector. In this range the showering and bundle fractions are robust against changes in the depth and topography, thus the total shower and bundle rate for a given experiment can be obtained by combining our results with an estimate for the total muon flux. One consequence is that a straightforward adaptation of the full detector showering muon cuts used by KamLAND to JUNO or RENO 50 would yield a nearly vanishing detector efficiency.


Introduction
In the next decade, the experiments JUNO [1] and RENO 50 [2] and perhaps LENA [3] will employ the largest liquid scintillator detectors ever constructed to detect antineutrinos from a variety of sources. While these detectors will be underground to shield them from cosmogenic muons, the first two will be at the relatively modest depths of 700 and 900 meters respectively. The showers induced by some of these muons will react with the carbon in the liquid scintillator and create 9 Li and 8 He some of whose decays yield the same double coincidence [4] used to identify antineutrinos via inverse beta decay (IBD). In previous experiments, such as KamLAND, the resulting background was largely eliminated by vetoing events occurring soon after the passage of cosmogenic muons [5]. However in the case of these two new, large detectors such a veto is delicate as the time between muon events is comparable to both the lifetimes of 9 Li and of 8 He. In this paper we present the results of a series of simulations of muons in such detectors. The results of these simulations on the one hand indicate that KamLAND's cuts cannot be straightforwardly applied to JUNO and RENO 50 but on the other hand can be used to determine the efficiency of any new veto scheme.
We used the FLUKA simulation package [6] to simulate the propagation of cosmogenic muons and antimuons inside of a 20 kton spherical detector consisting of a LAB-based liquid scintillator. This corresponds to the detector proposed for the experiment JUNO, while the RENO 50 detector will be 18 kton and cylindrical. We considered several different muon spectra, corresponding to various overburdens, topographies and cosmogenic muon distribution models. We will report results obtained using the cosmogenic muon distribution of Ref. [7] which is illustrated in Fig. 1. As µ − and µ + interactions in the detector result in different isotope production rates [8], the ratio of µ + to µ − will affect our results. We will assume an energy independent ratio of 1.37 µ + per every µ − corresponding to the value measured by Kamiokande for 1.2 TeV muons in Ref. [9]. In Sec. 2 we will describe the distribution of the muon's energy deposition in the scintillator, subtracting the deposition due to ionization which does not contribute significantly to the 9 Li and 8 He yield. In our next publication we will describe the 9 Li and 8 He yields themselves and their effect on the physics goals of these experiments, which will allow an optimal choice of veto strategies for these experiments given certain assumptions regarding the yet unknown tracking abilities of the detectors.
The unprecedented size of these detectors leads to a second consequence. As shown in Ref. [7], interactions of cosmic rays with the atmosphere often yield muon bundles consisting of multiple muons which travel in nearly the same direction. At the relevant depths the separations of muons are typically of order 10 meters. This means that in the case of a relatively small detector like KamLAND typically only a single muon will be observed in each event. As a result less than 5% of the muon events at KamLAND resulted in the detection of multiple muons [5]. On the other hand JUNO, RENO 50 and LENA are all much larger than 10 meters and so the vast majority of muon bundles will result in multiple muons impacting the detector. This will present several challenges. First, it will be nearly impossible to determine how much energy was deposited by which muon, and so to determine whether one or both muons is showering. Next, it will be more difficult to track the individual muons, although this tracking is important because full detector cuts in such a large detector weigh heavily upon the detector efficiency. Third, while we will find that the probability of an individual muon showering is about 20%, the probability that at least one muon is showering in a two muon event is about 30%. In Sec. 3 we will estimate the muon bundle frequency in various cases of interest.
Finally in Sec. 4 we will determine the expected spectrum and rates of the spallation isotopes 9 Li and 8 He. We will see that these background rates are comparable to the expected reactor neutrino signals, and so some veto scheme is necessary.

Showering Muons
We have used the initial muon energy distribution at several depths of interest as the input of a FLUKA simulation which determined the energy deposited by each muon in a spherical detector containing 20 kton of the liquid scinitillator. The energy deposited via ionization does not contribute to the creation of 9 Li and 8 He and also can be reliably subtracted if the track of the muon is known. Therefore we will not consider the energy deposited by ionization. It may be subtracted in two different ways, which yield essentially the same results. First of all, the ionization energy is proportional to the path length of the muon and so, once the path is known, it can be calculated and subtracted using a constant of proportionally appropriate for the liquid scintillator. This will be the only method which may be employed in the data analysis of the actual experiment. Second, FLUKA automatically separated the energy deposition in categories, of which one is ionization, and so one may simply not include the ionization deposition in the energy loss extrapolated from FLUKA. In the case of multimuon events, the detector cannot separate the energies deposited by the various muons, thus the multimuon event case is quite different from the single muon events. For simplicity, in this section we will consider only single muon events.
The diameter of the detector considered is about 35 meters, which is much larger than the 10 meter typical separation of muons in a muon bundle [7]. This means that if one muon from a bundle strikes the detector, it is very likely that all muons in the bundle strike the detector. As a result, muon bundles lead to multimuon events. Therefore we will assume in this note that the number of muons in a bundle which are detected is equal to the number of muons in a bundle or in other words that either all or none of the muons in a given bundle strike the detector. This is quite different from the case of KamLAND, whose detector diameter is of order a bundle size and so in general at most one muon from each bundle strikes the detector. As a result the single muon rate at KamLAND includes most bundle events whereas in the case at hand it includes essentially no bundle events. Multimuon events will be the subject of Sec. 3.
Following [7], the flux of incident muons arriving from a zenith angle θ is where h is the weighted difference in altitude measured in kilometers between the detector and the surface. The material will always be standard rock, and so we define a weighted altitude h as 2.64 times the true altitude h r . While this note was written for a flat surface in which case h is a constant, we will also consider detectors under mountains in which case h may be a function of θ. This is the correct prescription for generalizing the results of [7] to nontrivial topographies because muons arising from different angles θ are independent of one another, and so h can be chosen independently for each value of θ.
As a result of the factor of cos(θ), Eq. (2.1) yields the flux of muons crossing a horizontal surface of unit area. Our detector is not a horizontal surface, it is a sphere. Therefore the correct angle in our case is not θ but rather the angle between the normal to the detector and the direction from which the muon arrives. However we may simplify our expression for the flux by noting that, for a spherical detector, whatever the angle θ of the incoming muon the total cross sectional area of the detector is πr 2 where r is the detector's radius. We set r to 18 meters although the inner detector radius is actually 17.7 meters, as a result of spill in of the muon showers the fiducial volume for muon events is likely to be greater than 18 meters. Thus the incident flux upon the detector from an angle θ is πr 2 K(h, θ)/cos(θ) where the factor of cos(θ) has been removed as πr 2 is the area of the detector perpendicular to the velocity of the incident muon. Integrating this quantity over the angles from which the muons arrive yields the single muon event rate where we recall that, if the terrain is not flat, h will be a nontrivial function of θ.
The muon flux per unit energy at a given angle θ is up to a normalization term which can be fixed by demanding that the integral of N over E yield the incident flux πr 2 K(h, θ)/cos(θ). E is the energy in units of TeV. Fixing this normalization and integrating over solid angles in the upper half of the detector's surface one finds that the muon rate per unit energy is where again we recall that for a nonflat topography h is a function of θ. The rate per energy is reported Fig. 1 for detectors 700 meters and 900 meters underneath a flat standard rock surface, at the preferred JUNO and RENO 50 sites and also 200 meters below the preferred JUNO site.
Given the initial energy and the trajectory of a muon with respect to the detector, FLUKA yields the energy deposited via various channels. We have simulated muons incident upon the detector with various impact parameters, weighted according to their likelihood such that all muons pass through the inner detector. Each datapoint in Fig. 2 corresponds to the initial energy and the energy deposited not counting ionization of one of our simulated muons. Integrating this distribution over the initial energy one arrives at the upper panel of Fig. 3, which displays the distribution ρ(E) of energies deposited by muons at the preferred site for the JUNO experiment and 200 meters deeper.
Integrating ρ(E) from E to ∞ and normalizing the result to unity we obtain the lower panel of Fig. 3. This displays the fraction of muons for which the deposited energy, not counting that deposited by ionization, is greater than each fixed level E. The definition of a showering muon used by the KamLAND collaboration is one which deposits 3 GeV of energy in addition to that deposited by ionization. As can be seen in Fig 7.18 of Ref. [10], this definition is useful as, in the case of a KamLAND sized detector, the vast majority of 9 Li and 8 He is produced by showering muons. In fact, at JUNO and RENO 50 the contribution of nonshowering muons to 9 Li and 8 He production will have little effect on the science goals of the experiment even if they are not vetoed at all. The lower panel of Fig. 3 shows that the fractional energy deposition distribution of the muons is quite insensitive to the depth in the range considered and in fact the showering fraction is extremely robust. For a 20 kton detector using KamLAND's definition of a showering muon, as had been anticipated in Ref. [11] about 20% of muons are showering at these depths, 200 meter variations in the depth and indeed even 200 meter changes in the topography have little effect on this robustly determined fraction. This is our main result.

Muon Bundles
In this section, we will use the parametric formulae of Ref. [7] to determine the multimuon event rates at 20 kton spherical liquid scintillator detectors in various settings. The m-muon flux is The m-muon rate can be calculated as in Eq. (2.2) In Table 1 we have summarized several relevant quantities. The muon rate is defined to be ∞ m=1 mR(m). The event rate is ∞ m=1 R(m). The bundle rate is ∞ m=2 R(m). The m-muon rate is just R(m) while the mean muon energy is, using the approximations of Ref. [7], independent of m. We have considered a detector 700 and 900 meters underneath a flat surface as well as various topographies corresponding to potential sites for detectors. In particular we have considered the preferred Dong Keng [12] site for JUNO and the preferred Mt GuemSeong [12] site for RENO 50, which are illustrated in the upper and lower panels of Fig. 4 respectively. These are respectively h r = 700 meters and 900 meters underneath the peaks of their corresponding hills, where h r is the rock overburden which is equal to h r = h/2.64. We have also considered a location 200 meters beneath the preferred site for JUNO. Note that the reported mean muon energies are higher than that found in other studies. This increase is caused by a small number of muons in the very high energy tail of the cosmogenic muon distribution, at several TeV.
To determine the expected muon flux ideally one requires a geological survey of the rocks around these sites and a full 3d simulation such as MUSIC. We have instead simply assumed that the rock is standard and have employed a cylindrically-symmetric approximation to the topographies illustrated in Fig. 4. For JUNO we have assumed an overburden of h r = 700 meters of rock in a 100 meter radius circle, followed by an annulus of inner radius 100 meters and outer radius 500 meters with a surface 50 meters lower, than another annulus of outer radius 700 meters with a surface 50 meters yet lower and finally we have assumed that the surface is another 50 meters lower at radii beyond 700 meters. Similarly we have  approximated Mt Guemseong using a series of annuli, beginning with a 200 meter radius circle with h r = 900 meters of rock overburden, followed by a 600 meter annulus with h r = 800 meters of overburden, a 900 meter radius annulus with h r = 700 meters of overburden and then we have assumed a h r = 600 meter overburden at all radii beyond 900 meters.

Isotope Production
We have also estimated the total yield of 9 Li and 8 He at each experimental site. To do this, we ran 10 7 FLUKA simulations of single monoenergetic muons of each energy between 1 GeV and 10 GeV and 10 6 simulations at energy between 11 GeV and 500 with a step size of 1 GeV, with an additional 10 6 simulations at each energy up to 3 TeV with a step size of 50 GeV, for a total of 6.4 × 10 8 simulations. Again the impact parameter was randomized with a distribution reflecting a homogeous distribution of muons. The 9 Li yield and best power law fit, with an exponent of 0.842 ± 0.002 are presented in red in the top panel of Fig. 5. The results of another 6.4 × 10 8 simulations of antimuons are presented in blue, together with a power law fit whose exponent is 0.847 ± 0.002. The same simulations also determined the 8 He rate, as is shown on the bottom panel. The best fit exponents are 0.869 ± 0.006 and 0.861 ± 0.006 for µ − and µ + respectively.
The power law fits reproduce the simulated isotope yields to within the statistical errors from 10 GeV up to 3 TeV, where the vast majority of the spallation isotope production occurs. However below 10 GeV, the isotope production from µ − events increases. Indeed  Figure 5: The mean number of 9 Li (top) and 8 He (bottom) created by a single µ − (red) and µ + (blue) with a given energy. The lines are the best power law fits. The gray band represents the relative statistical uncertainty on the simulated isotope yield. in this range the 9 Li yield is 3.3 times greater for µ − events than for µ + . This is because µ − at these low energies stop and bind with the hydrogen and carbon in the scintillator. The deeper potential well about the carbon nucleii means that any µ − which initially bind to hydrogen are soon transfered to carbon and these quickly decay to 1s or 1p orbitals. Here most of the µ − , like the µ + , simply decay. However 7% of the µ − are then captured by the carbon nucleii, undergoing a charged current interaction which can create spallation isotopes.
To arrive at the rate for each experimental site we have folded the power law fits to the simulated data as a function of muon energy with the expected muon spectra reported in Fig. 1. In contrast with Fig. 1, here we are interested in the total µ rate, not just those arising from single µ events, so we have rescaled this result by the ratio of the total µ rate to the single µ rate reported in Table 1. We have assumed that the ratio of µ + to µ − is 1.37, as was found for 1.2 TeV muons by Kamiokande in Ref. [9]. The resulting 9 Li and 8 He rates per muon energy bin are displayed in Fig. 6. Next, to arrive at the total 9 Li and 8 He rates we have integrated this figure over the µ energy. The resulting total rates at each experimental site are summarized in Table 2. The 9 Li rates are compatible with the rough estimate of Ref. [11] obtained by simply rescaling KamLAND's rate.
As can be seen in Table 2, the total spallation isotope rates are appreciably higher than the expected reactor neutrino IBD signal rates of 3 × 10 −4 Hz for JUNO and 2 × 10 −4 Hz for RENO 50 at Guemseong. However only those decays which produce a neutron yield a false double coincidence signal. These are 51% of 9 Li decays and only 16% of 8 He decays. In the bottom two rows of Table 2 we report the false double coincidence rate expected. Note that this background double coincidence rate is still twice the IBD signal rate at JUNO and RENO 50, although it would be slightly lower than the RENO 50 signal rate at the Munmyeong site of Ref. [14] even with the same overhead burden as the Guemseong site.
In practice the background rate will be reduced, at the expense of dead time, by a veto program. Similarly, the expected background can be subtracted from a signal by a shape analysis. However the statistical fluctuations created by the background will survive this subtraction and obscure the low energy oscillations in the reactor neutrino spectrum whose observation is necessary for a determination of the hierarchy. Thus, both the sensitivity to the hierarchy and the precision with which θ 12 may be measured depend critically on the optimal choice of muon vetoes. As the muon rate is comparable to the 9 Li half life, full detector vetoes are not an option for all but the highest energy showers, thus excellent muon tracking, including the tracking of multimuon events and even muons that arrive horizontally will be necessary to achieve these science goals.    Table 2: Rates in 10 −5 Hz of 9 Li and 8 He production via cosmogenic muon spallation on 12 C under 700/900 meters of rock with a flat surface, at JUNO, 200 meters beneath JUNO and at RENO 50. The resulting 9 Li and 8 He decays only provide false double coincidence signals when their decay produces a neutron, which happens in 51% and 16% of their decays respectively. The last two rows report the rates of these decays. The errors reported reflect only the uncertainty in the muon rate, not the uncertainties in the isotope production per muon.

Concluding Remarks
As can be seen in the lower panel of Fig. 3, for any depth in the range relevant to future large liquid scintillator experiments, about 20% of single muons will be showering in a 20 kton spherical detector. Thus to determine the showering single muon rate, one only needs to find the single muon rate for a given site and multiply it by 20%. Similarly from Table 1 we can see that 17-20% of muon events will occur in muon bundles with multiple muons entering the detector in the vast majority of cases. As, for the large detectors considered here, in most cases muon bundles yield multimuon events in the detector, the multimuon event rate is about 15%. Again, this bundle fraction is fairly robust against changes in depth and topography and so may be applied to a wide variety of candidate experimental sites.
One immediate consequence of our study is that both the bundle and the showering muon rates are about 0.8 Hz (0.5 Hz) at JUNO (RENO 50). This means that KamLAND's veto strategy, employing 2 second full detector cuts for showering and poorly constructed tracks, cannot be applied to these experiments as the events are separated by less than 2 seconds. If JUNO is moved 200 meters lower, KamLAND's cuts would still be problematic although 1 second full detector cuts would be possible with considerable dead time. 1 second cuts are sufficient to reduce the 9 Li and 8 He backgrounds well below the level of the signal. KamLAND's veto strategy has been assumed in all studies of these backgrounds that have so far appeared in the literature, and so our results indicate that these studies need to be redone.
As a result, full detector cuts will be infeasible at these experiments and it is critical that muon tracking be successful for the vast majority of single and multimuon events, so that selective cuts may remove events from cylinders surrounding cosmogenic muon tracks. In our next paper we will investigate the creation of 9 Li and 8 He by these muons and so determine the effectiveness of various veto schemes as well as the cost to the scientific goals both as a result of the increased dead time and as a result of the background which will mask the 1-3 oscillations in the 3-4 MeV range whose detection is essential to determine the neutrino mass hierarchy [13]. However for any proposed veto scheme, the lower panel of Fig. 3 together with the total muon flux may be used to calculate the resulting veto efficiency and so the resulting dead time for the experiment.