Neutral-current background induced by atmospheric neutrinos at large liquid-scintillator detectors: I. model predictions

The experimental searches for diffuse supernova neutrino background and proton decay in next-generation large liquid-scintillator (LS) detectors are competitive with and complementary to those in the water-Cherenkov detectors. In this paper, we carry out a systematic study of the dominant background induced by atmospheric neutrinos via their neutral-current (NC) interactions with the $^{12}{\rm C}$ nuclei in the LS detectors. The atmospheric neutrino fluxes at the location of Jiangmen Underground Neutrino Observatory (JUNO) are used, as the JUNO detector is obviously a suitable representative for future LS detectors. Then, we implement the sophisticated generators GENIE and NuWro to simulate the neutrino interactions with the carbon nuclei, and the package TALYS to deal with the deexcitations of final-state nuclei. Finally, the event rates for the production of additional nucleons, $\gamma$'s, $\alpha$'s, pions and kaons are obtained and categorized, and the systematic uncertainty of the NC background represented by a variety of data-driven nuclear models is estimated. The implications of the NC background from atmospheric neutrinos for the detection of diffuse supernova neutrino background and proton decay are also discussed.


Introduction
part of the background is mainly composed of the IBD interactions of the atmospheric ν e as well as the charged-current (CC) interactions of atmospheric neutrinos with 12 C nuclei in LS, where copious neutrons, protons, γ's and α's are generated and can contaminate the IBD signals. The CC background induced by atmospheric neutrinos can be essentially removed by shrinking the energy window. Therefore, the energy range of interest is usually restricted to the region between two intrinsic ν e backgrounds from reactor and atmospheric neutrinos. It is worthwhile to mention that the exact window depends to a large extent on the chosen detector and the control of various backgrounds.
• Unlike the wCh detectors [15][16][17], the cosmic-ray muon spallation related backgrounds are well under control in the LS detectors by implementing the muon veto. Fast neutron background generated by muon spallation outside the detector can be removed by cutting the outer layer of the detector, implying a slight reduction of the fiducial volume. Thus a balance between the signal loss and the reduction of backgrounds should be made.
• In the chosen energy window, we are finally left with the dominant NC background caused by atmospheric neutrinos. Different from the CC interactions, where the associated finalstate charged leptons carry away most of the initial-state neutrino energies and deposit their energies in the LS detectors, the final-state neutrinos in the NC interactions escape from detection. Meanwhile, the produced neutrons, protons, α's and residual light nuclei, whose deexcitations result in high-energy gamma rays, contribute to the main background.
Different from the DSNB, that is undoubtedly presented but yet to be discovered, proton decay is actually predicted by the ultimate dream of elementary particle physics, namely, the grand unified theories (GUTs) [25,26]. The discovery of proton decay will offer a robust evidence for such fundamental theories. Historically, the experimental searches for proton decay p → e + +π 0 in the wCh detectors [27,28] actually accelerated the development of neutrino physics. The latest limits for proton decay at the SK can be found in Refs. [29][30][31], while the studies of p → K + + ν in the LS detectors have been carried out in Refs. [32,33]. Compared to those of the DSNB, the experimental signals of proton decay and the relevant backgrounds occur at higher energies. This requires us to study the NC background over a broad energy range.
Motivated by the prominent importance of the experimental searches for DSNB and proton decay, in the present work we perform a systematic study of the NC background induced by atmospheric neutrinos. For this purpose, several ingredients are required. First comes the precise calculation of the fluxes of atmospheric neutrinos, which originate from the interactions of highenergy cosmic rays with the nuclei in the Earth's atmosphere [34] and the subsequent decays of the produced mesons and secondary muons. However, depending on the distribution of the Earth magnetic fields and the structure of the Earth's atmosphere, the neutrino fluxes should be calculated specifically for the detector site under consideration. To be concrete, we choose the JUNO site, but the details of the JUNO detector are not required. Secondly, the interactions between atmospheric neutrinos of energies ranging from MeV to 10 GeV and the target 12 C nuclei in the LS detectors are to be modeled. We employ the widely-used generators GENIE [35] and NuWro [36] for the neutrino interactions, and the package TALYS [37] for the deexcitations of the final-state nuclei. By using a variety of data-driven nuclear models, we are able to estimate the systematic uncertainties in the prediction of the NC background. For the DSNB search, the NC background consists mainly of the processes with a single final-state neutron and the rate is found to be (3.1 ± 0.5) kt −1 yr −1 in the energy range of 11 MeV E vis 30 MeV, where E vis is the visible energy in the LS detector. For the proton decay, the primary contribution to the background comes from the production of single π ± or K ± in the neutrino-12 C interactions. By restricting the energy transfer in relevant processes, one can find that the single π ± production rate within [150,650] MeV is (0.61 ± 0.14) kt −1 yr −1 and the single K ± rate is (4 +5 −4 × 10 −4 ) kt −1 yr −1 for the range of [150,1000] MeV. For both DSNB and proton decay, the uncertainty in the background prediction has been estimated by averaging over different models for neutrino-12 C interactions. Finally we would like to mention that NC interactions investigated here are also a significant background for searches of neutrinos from the dark matter annihilation [38,39].
The remaining part of this paper is organized as follows. In Sec. 2, we outline our strategy for the practical calculations of the NC background in the LS detectors, where all the necessary ingredients are offered and explained. Then, the final results of the NC background are given in Sec. 3, and the resultant energy spectra of protons, neutrons, γ's, α's and others are also provided, which may be useful for other general-purpose studies. In addition, the NC background rates are estimated for the searches for DSNB and proton decay. Finally, we summarize our main results and conclude in Sec. 4.

Atmospheric Neutrino Fluxes
In the first place, the fluxes of atmospheric neutrinos ν µ , ν µ , ν e and ν e must be calculated as precisely as possible. Since we are concerned about the NC interactions of atmospheric neutrinos with the carbon nuclei in the LS detectors, which are insensitive to neutrino flavors, the flavor conversions of atmospheric neutrinos will have no impact on the final results. In our calculations, these fluxes have been taken from the latest results by the Honda group [40], and the associated uncertainties have further been reduced by using the accurately measured atmospheric muon flux and estimated in Ref. [41]. A brief summary of the flux calculations in Refs. [40,41] is as follows.
The calculation of atmospheric neutrino fluxes is carried out by simulating realistic interactions of primary cosmic rays with the atoms in the Earth's atmosphere, the propagation of primary and secondary particles in the geomagnetic fields, the decays and interactions of produced mesons and subsequent muons. First, the temperature and air density profiles of the Earth's atmosphere are provided by the NRLMSISE-00 model [42], which improves the density profile in U.S.-standard 1976 [43] by taking account of the time variation and the position dependence around the Earth. Second, the model of primary cosmic rays has been constructed by incorporating recent precision measurements by AMS02 [44,45] and other experiments [46,47]. Third, the hadronic interaction models [48,49] are implemented to deal with the collisions between the primary cosmic rays and the atoms in the atmosphere, and then the realistic IGRF geomagnetic model [50] is utilized to perform three-dimensional simulations of the propagation of cosmic rays and their secondaries [51]. Finally, the neutrino fluxes are obtained by following their motion in the atmosphere and selecting of atmospheric neutrinos for ν = ν µ , ν µ , ν e , ν e at the JUNO site, as calculated by the Honda group in Ref. [53].
only those registered in the virtual detectors at the chosen experimental site [52].
In Fig. 1 , namely, the predicted fluxes φ ν (E ν ) of atmospheric neutrinos multiplied by E 3 ν for ν = ν µ , ν µ , ν e , ν e , at the JUNO site, as calculated by the Honda group [53]. Notice that the fluxes in Fig. 1 have been averaged over all directions, and there is no mountain assumed at the JUNO site 1 . In addition, the neutrino energies are ranging from 100 MeV to 10 4 GeV, and the uncertainty in the predictions for atmospheric neutrino fluxes varies from about 20% at low energies to less than 10% in the range of (1 − 10) GeV. Starting from the neutrino energy around 10 GeV, the uncertainty is gradually increasing mainly from the variations of hadron interaction models for the π and K production. As already pointed out in Ref. [41], the uncertainty of neutrino fluxes could be further reduced by more precise measurements of cosmic muon flux at the same experimental site.

Neutrino Interaction Generators
Next, given the fluxes of atmospheric neutrinos, we proceed with the neutrino interactions with 12 C in the LS detectors. The cross section for the NC interactions between neutrinos and nuclei suffers from various uncertainties mainly in the nuclear structure and many-body effects in the nuclei. First, the individual nucleon participating in the NC interactions is actually confined in the nucleus, so the nuclear structure of the latter should be taken into account. Second, as the energy of the incident neutrino increases from 100 MeV to GeV or even higher, the dominant contribution to the cross section comes roughly from quasi-elastic scattering (QEL), coherent and diffractive production (COH), nuclear resonance production (RES) and deep inelastic scattering (DIS) in different energy ranges. In general, all these processes should be considered. Third,

Models
Generator  [35] and NuWro [36], where the adopted models differ in the input values of the axial mass M A , whether the relativistic Fermi gas (RFG) model of or the spectral function (SF) approach to nuclear structure is used, or whether the transverse enhancement model (TEM) in the two-body current contribution is considered.
the nuclear effects, such as the final-state interaction, the meson exchange current (MEC) and multi-nucleon correlation, must be included as well in a more complete study [55][56][57].
In order to take account of different contributions and handle nuclear effects in a proper way in the neutrino interactions with nuclei, one usually relies on some sophisticated Monte Carlo generators, which have been carefully constructed and finely adjusted according to the available experimental data on the interaction cross sections. In addition to GENIE and NuWro, which will be used in our calculations to illustrate the dependence on nuclear models, other generators NUANCE [58], NEUT [59] and GiBUU [60] are available. See, e.g., Ref. [61], for a recent review on the neutrino event generators. Depending on incident neutrino energies, the NC interactions with the 12 C nuclei will be dominated by the QEL for the energy range of several hundred MeV, or by the DIS for energies above 100 GeV. For neutrino energies in between, the NC interactions will be complicated by multiple production processes and various nuclear effects [57]. In the present work, we utilize the software version of GENIE (2.12.0) and NuWro (1.7.10) to calculate the cross sections. Both generators provide some options for the physics parameters and nuclear models, among which one key parameter is the axial mass M A in the parametrization of the nuclear axial-vector form factor. In GENIE, the default setting is M A = 0.99 GeV, which has been determined from the deuterium measurements [62]. In NuWro, this parameter is tunable, and thus we take the following three different values.
• M A = 0.99 GeV, i.e., identical to that in GENIE, in order to study the systematic differences arising from the modelling in GENIE and NuWro.
• M A = 1.35 GeV, which is mainly motivated by the latest analysis of the MiniBooNE data on neutrino-carbon interactions [63]. The MiniBooNE data indicate a surprisingly large cross section, which has triggered an intense discussion on possible explanations [64,65].
• M A = 1.03 GeV, which is actually taken the world average from Ref. [66]. Since this value is very close to M A = 0.99 GeV, the final results are expected to be rather similar.
Regarding the models of nuclear structure, GENIE uses the relativistic Fermi gas (RFG) model as   Table 1 have been shown. a default setting, while NuWro provides several choices for the description of the target nucleus, including both RFG and the spectral function (SF) approach. Furthermore, to illustrate the twobody current effects in QEL, we also take the transverse enhancement model (TEM) of the meson exchange current from NuWro, which has been obtained in Ref. [67] by fitting the electron scattering data. By including TEM, the authors of Ref. [67] have demonstrated that the MiniBooNE results can be reproduced with a smaller value of the axial mass, which turns out to be consistent with the world average value obtained from other measurements. It should be emphasized that only the axial mass M A in the treatment of QEL in NuWro has been changed. For both generators, we shall employ their default setting for all other processes.
In Table 1, we summarize the main features of six typical models, which have been implemented in our calculation to illustrate the model variations of neutrino interactions. Some comments on these models are in order. Just one model from GENIE (i.e., Model-G) has been adopted, while Model-N1 from NuWro has been selected with the same input for the nuclear model and the axial mass. The motivation for this setup is to make a comparison between GENIE and NuWro. Within the same generator NuWro, Model-N1, Model-N2 and Model-N3 share exactly the same setup, except for different values of the axial mass M A for the QEL process, which is intended to illustrate the impact of the axial mass on the NC cross section. Additionally, two other models in NuWro are introduced. First, Model-N4 is the only one to include TEM, so that we can examine the two-body current effect by comparing between Model-N1 and Model-N4. Second, the SF approach is incorporated in Model-N5, whereas the RFG model is used in all others, offering a possibility to study the difference between these two nuclear models within NuWro.
In Fig. 2, we have extracted the inclusive cross sections for the NC interactions of neutrinos and antineutrinos with 12 C from GENIE and NuWro, and shown the final results of σ(E ν )/E ν and σ(E ν )/E ν for six representative models in Table 1. A number of important observations can be made from Fig. 2. First, for both neutrino and antineutrino cross sections, Model-N5 gives the smallest value of the inclusive cross section. This is the only model that uses the SF approach for Daughter Nuclei Shell Hole Configuration Probability Excitation Energy Table 2: The probabilities of the configurations for the nuclei 11 C and 11 B with one nucleon disappearing from 12 C in the nuclear model and those for 10 C, 10 Be and 10 B with two nucleons less, where the corresponding excitation energies E * are given in the last column and E * = 0 MeV actually refers to the ground state.
the target nucleus, instead of the RFG model. As previously shown in Refs. [68,69], the inclusion of many-body nuclear effects in the SF approach generally reduces the total cross sections by 20% or so at the energy around 1 GeV. This effect can be clearly seen by comparing between the cross section at E ν = 1 GeV (or E ν = 1 GeV) in Model-N1 and that in Model-N5, for which the same axial mass M A = 0.99 GeV is input. For even lower energies, the relative difference between these two models becomes more significant. However, the impulse approximation made in the SF approach will be invalidated for low momentum transfers [68], rendering such a comparison to be problematic. Second, the predictions from Model-N1 and Model-G are well consistent with each other for the energies above 300 MeV, but remarkable deviations can be found for lower energies. This discrepancy might be attributed to the different treatments of nuclear effects (e.g., the final-state interaction) in NuWro and GENIE, and will be taken as a systematic uncertainty in our calculations. Third, the axial mass M A = 1.35 GeV is set in Model-N3, leading to the largest cross section for both neutrinos and antineutrinos. In the neutrino sector, though M A = 0.99 GeV is small, Model-N4 predicts even larger cross section below 1 GeV due to TEM. Both of them are motivated by the MiniBooNE data, which does not necessarily mean more accurate description of microscopic physics. Notice that the TEM effect is absent in the antineutrino sector, so the predicted cross sections in Model-N1, Model-N2 and Model-N4 essentially coincide. Finally, it is evident that the cross section for neutrinos turns out to be larger than that for antineutrinos in all models under consideration. The main reason is that the opposite sign of the axial-vector coupling for neutrinos and antineutrinos, which have the opposite helicities, results in the constructive interference in the transverse-and axial-vector amplitudes for neutrinos but the destructive interference for antineutrinos [70,71].

Deexcitation of Final-state Nuclei
The last step is to deal with the deexcitation of final-state nuclei from the NC interactions and extract the total energy spectra of γ's, light mesons, protons, neutrons, and α's, which may contribute to the irreducible backgrounds for the experimental searches for rare events in the LS detectors. However, the Monte Carlo generators of neutrino interactions usually do not provide the exact state of the residual nucleus, which may reside in one of its various excited states, and the nuclear deexcitation with additional γ rays, protons, neutrons or other heavier projectiles is obviously important. Thus it is desirable to combine the neutrino interaction generators with the nuclear structure model and the deexcitation tool to complete our calculations.
To this end, the widely-used TALYS software [37] will be employed to treat the deexcitation of daughter nuclei. In general, TALYS is a very useful tool for both nuclear structures and nuclear reactions in the energy range from 1 keV to 200 MeV. For our purpose, the branching fractions of the deexcitation of the excited nucleus and the energy spectra of all the relevant products from the deexcitation are needed. More explicitly, after the excited state of a nucleus N * with the excitation energy E * is chosen, we implement TALYS to simulate all possible deexcitation channels of the given nucleus and follow the further deexcitation of the residual nuclei until all the final-state particles are in their ground states, either stable or beta decaying. When the residual nuclei are unstable, we figure out their decay processes based on the decay types, endpoints, and lifetimes from the nuclear database. Hence the energy spectra of all deexcitation products are obtainable.
Before doing so, we have to specify the nuclear structure of the target nucleus 12 C. In order to simplify our discussions, we make use of the statistical configuration from the nuclear shell model of 12 C [72][73][74]. Since there are six protons and six neutrons in the 12 C nucleus, the s 1/2 and p 3/2 shells of the lowest energy level for both protons and neutrons will be occupied in the simplest configuration. Moreover, we neglect the potential configuration of the 12 C nucleus due to the nucleon pairing correlation between the p 3/2 and p 1/2 shells because of the relatively small energy gap of around 3-4 MeV [72,73]. This effect will be included in our future work with more sophisticated shell model calculations.
Using this statistical model of the ground state of 12 C, we can figure out all possible excited states for those lighter daughter nuclei resulting from the neutrino interactions with 12 C. The excited states of daughter nuclei are obtained by considering the disappearance of one or more nucleons (either protons or neutrons) from the 12 C ground state. For instance, we show the cases of one or two nucleons disappearing from 12 C in Table 2, where the corresponding probabilities of possible configurations are provided together with the excitation energies of these daughter nuclei. The information of these daughter nuclei will be further input in TALYS to obtain their deexcitation and the energy spectra of the associated final-state particles. In a similar way, one can discuss other possibilities of more nucleons disappearing from 12 C. In our calculations, all the daughter nuclei with a mass number larger than five have been taken into account.

Results and Discussions
Sec. 2 describes the ingredients necessary for a numerical calculation of the background induced by NC interactions of atmospheric neutrinos: the atmospheric neutrino fluxes at the JUNO site in Fig. 1, the inclusive cross sections of (anti)neutrino-12 C interactions in Fig. 2, the possible excited states of daughter nuclei (e.g., those in Table 2) and the corresponding deexcitation processes. In this Section, we proceed to carry out a detailed analysis of the NC backgrounds for the detection of DSNB and proton decay. We shall first present the general calculation of NC interactions in Sec. 3.1, and then the specific results on the background in searches of DSNB and proton decay in Sec. 3.2 and Sec. 3.3 respectively. However, it is worthwhile to emphasize that the background analysis for any practical search of rare events depends very much on the properties of signals, the performance of the LS detectors and other advanced techniques. Thus we concentrate on the main features of the NC backgrounds and leave the intricate strategy for background reduction for future and better works by the experimental collaborations.

NC interaction rates
We calculate the event rate of atmospheric neutrino NC interactions as a function of the energy of the incident neutrino or antineutrino based on where ν = ν µ , ν µ , ν e , ν e runs over neutrinos and antineutrinos of both electron and muon flavors, the factor 4π comes from the integration over the full solid angle, and N ≈ 4.4 × 10 31 is the total number of target 12 C nuclei per kiloton LS, where the contribution of 13 C at the level of 1% has been neglected for simplicity. Note that the cross section σ ν (E ν ) could be either the inclusive one as shown in Fig. 2 or the exclusive one with specified final states listed in Fig. 4. The distribution of the total event rate with respect to the neutrino energy has been calculated by convolving the atmospheric (anti)neutrino fluxes in Fig. 1 Table 1 are shown in Fig. 3. For comparison, in the lower panel of Fig. 3, the ratio of the result in each model to that in Model-N1 is presented. The event rate essentially follows the main features of the inclusive cross section in each representative model. For instance, there are more events for lower neutrino energies but less for higher neutrino energies in Model-G compared to other models using NuWro. Moreover, Model-N5 using the SF approach predicts the least number of NC events. All these features can be easily understood by recalling the inclusive cross sections in Fig. 2.

and the inclusive cross sections in
In addition, we have checked the individual contributions from different physical processes (i.e., QEL, RES and DIS) to the inclusive event rates in six representative models. Roughly speaking, QEL, RES and DIS respectively contribute around 60%, 20%, 10% of the total event rate. The QEL is the dominant process in the energy range of the DSNB signal, while RES and DIS are the major processes in the range for the proton decay search. Note that the inclusion of the meson exchange current in Model-N4 leads to a remarkable impact, which compensates the relatively smaller event rate of QEL in this model due to a smaller value of M A . On the other hand, the individual contributions from ν µ , ν e , ν µ , ν e turn out to be comparable, namely, they account respectively for about 50%, 20%, 20% and 10% of the total events in all six models. This finding is directly related to the differences in the fluxes and the inclusive cross sections. Firstly, the inclusive cross section for antineutrinos is smaller than that for neutrinos, as shown in Fig. 2. Secondly, the atmospheric neutrino flux of the muon flavor is dominant. Ratio to Model-N1 Figure 3: The distribution of the total event rate induced by atmospheric neutrino NC interactions with respect to the neutrino energy for each of six representative interaction models in Table 1.
The atmospheric neutrino fluxes in Fig. 1 and the inclusive cross sections in Fig. 2 have been input in the calculations, and the ratio of the event rate in each model to that in Model-N1 has been shown in the lower panel.
In order to examine how the NC interactions of atmospheric neutrinos affect the detection of DSNB and proton decay, we have to go one step further and analyze the final-state products. For this purpose, we calculate the event rates of atmospheric neutrino NC interactions on 12 C nuclei in the exclusive channels, where different final-state nuclei and the production of one or more nucleons should be taken into account. Our final results are summarized in Fig. 4 and Fig. 5.
The relevant energy range of atmospheric neutrinos in question is 100 MeV < E ν < 10 GeV. Although the atmospheric neutrino fluxes φ ν (E ν ) multiplied by E 3 ν in Fig. 1 are shown up to E ν = 10 4 GeV, the fluxes φ ν (E ν ) themselves will be highly suppressed at high energies. For the differential event rates of the inclusive process of NC interactions shown in Fig. 3, the maximum rates appear at E ν = (200 − 300) MeV. Within this energy range, the QEL process is most important and one or more nucleons will be knocked out from the carbon nucleus. The event rates for the same exclusive processes but with one or more pions turn out to be negligibly small and can be safely ignored.
In Fig. 4, we summarize the event rates of the relevant exclusive processes, which have been categorized by the residual final-state nuclei and the associated superscript " * " reminds us of possible excited states of these nuclei. The results after taking account of the subsequent deex- citation are presented in Fig. 5. We divide the daughter residual nuclei into two groups. The first group includes all the isotopes with the atomic number Z ≥ 3 and the mass number A ≥ 6 except for 12 C. The event rates for this group have been shown in the first three panels (from left to right) of Fig. 4, where for each isotope the predictions from Model-G, Model-N1 and Model-N5 are represented by solid filled histograms, the striped histograms with slashes, and the striped histograms with backslashes, respectively. As Fig. 3 indicated in the previous section, the event rates in the series of models using NuWro (i.e., Model-Ni for i = 1, 2, 3, 4) are quite similar, therefore, only Model-N1 has been chosen as representative for clarity. For a given model, the individual contributions from the subset of processes with accompanying nucleons can be recognized via the colored bars, where the colors refer to different combinations of the neutron and proton multiplicities. The total event rate is given by the height of the whole histogram, which is the sum of the heights of all colored bars. Notice that the scales of the vertical axes in all the panels are different from each other. The second group contains all the isotopes with Z < 3 or A < 6 or 12 C, for which the results are given in the rightmost panel labelled by "Others". For this group, instead of the final-state nuclei, the subsets of processes have been characterized by the neutron multiplicity, which has been represented by different colored bars. As one can observe from Fig. 4,  Figure 5: The event rates for the NC interactions of atmospheric neutrinos (with energies ranging from 100 MeV to 10 GeV) with 12 C nuclei in the exclusive channels, which have been categorized by the daughter residual nuclei in their ground states and the deexcitation of the excited residual nuclei has been simulated by using TALYS. The notations and the patterns of the histograms in each panel exactly follow those in Fig. 4.
Model-G generally predicts a higher rate for the exclusive processes with extra nucleons, which can be ascribed to more low-energy events from the simulations with GENIE. In addition, the event rate is dominated by the processes with the production of 11 C * , 11 B * , 10 B * , 10 C * , 10 Be * , 9 B * and 9 Be * , whose deexciation will be further processed by using TALYS.
In Fig. 5, we summarize the final event rates in Model-G, Model-N1 and Model-N5, after taking into account the deexcitation of final-state nuclei from TALYS. As we have explained in the previous section, the probabilities and excitation energies of different configurations for a specified nucleus have been extracted from a simple statistical model of 12 C, as partially summarized in Table 2. With this input information, one can simulate all possible channels of deexcitation via TALYS until the daughter nuclei are in their ground states. Therefore, all the exclusive processes in Fig. 5 are now categorized by the final isotopes that are left in the ground states, where the notations and the patterns of the histograms follow exactly those in Fig. 4. Comparing between Fig. 4 and Fig. 5, one can see that the event rates for the processes associated with 9 B, 9 Be, 8 Be, 7 Be, 7 Li and 6 Li increase significantly, so does that for "Others". Meanwhile, the event rates for the processes associated with 11 C, 11 B, 10 C are largely reduced because of the further knockout of one or more nucleon in the deexcitation. Moreover, the combinations of neutron and proton multiplicities related to each of these nuclear isotopes become more complicated due to the deexcitation processes.

DSNB
First, let us consider the signals of DSNB and possible backgrounds induced by the atmospheric neutrino NC interactions in an LS detector. As has been mentioned in Sec. 1, the IBD process ν e + p → e + + n is the golden channel for the detection of DSNB ν e , for which the time coincidence between the prompt signal of the positron annihilation and the delayed signal of neutron capture helps reduce greatly the background. However, the irreducible backgrounds come from reactor ν e in the low-energy range and atmospheric ν e in the high-energy range, leaving only a narrow window of visible energies 11 MeV E vis 30 MeV for the DSNB observation [19].
Even in this optimal energy window, the NC interactions of atmospheric neutrinos can mimic the DSNB signals. As the energies of our interest are relatively low, the QEL process is of crucial importance. In Fig. 6, the energy spectra of p, n, γ and α particles have been extracted from the simulations and are shown with respect to their kinetic energies. In the upper panel, the spectra have been divided into two different categories. First, we show the energy spectra of nucleons, either protons or neutrons, which are directly produced from the NC interactions. In accordance with previous observations, the production rate from the Model-G (red dashed curve) is much higher than those from Model-Ni for i = 1, 2, 3, 4 (black dashed curve) and Model-N5 (blue dashed curve) below 100 MeV. Since the predictions from Model-Ni for i = 1, 2, 3, 4 are quite similar, only the average value of them is shown and the gray band along the black dashed curve represents the 1σ deviation. Second, the energy spectra of p, n and α particles from the deexcitation of final-state nuclei that are produced in the NC interactions have been plotted as blue, magenta and orange solid curves, respectively. The colored bands along the solid curves stand for 1σ deviation from the average of all six models. Although the deexcitation processes are the same for these models, the production rates for a given nuclear isotope actually differ. As one can observe from the upper panel, there is significant model dependence of the neutrino interaction generators in the low energy range, which indicates the necessity of using a complete set of neutrino interaction models. Meanwhile, the nucleon knockout from the deexcitation of finalstate nuclei will be comparable to that directly from NC interactions around or below 10 MeV, but decreases rapidly toward high energies. On the other hand, the energy spectrum of γ's has been shown as gray histograms in the lower panel, where discrete lines are superimposed on a spectral continuum reaching to 30 − 35 MeV. The highest-rate line at E γ = 23 MeV can be attributed to the excitation energy of several excited nuclei in Table 2. The p, n, α, and γ-rays will deposit their kinetic energies in LS immediately after their production. If followed by a neutron capture, the prompt scintillation signals can mimic the prompt event of an IBD coincidence. The neutrons are tagged with high efficiency via their captures on hydrogen in the LS detectors. Therefore, the event rates of the exclusive channels of the QEL process have been categorized in Fig. 7 according to the neutron multiplicities. Some comments on Fig. 7 are helpful.
• The event rates for all channels in the series of models (i.e., Model-Ni for i = 1, 2, 3, 4) are quite similar, thus the average value and the standard deviation of the predictions from these four models are shown as squares with error bars, and labelled by "Model-N(1-4)". However, the error bars are too small to be visible. Generally speaking, the event rate decreases as Production rate [kt : deexcitation γ Figure 6: The energy spectra of p or n from the NC interactions of atmospheric neutrinos with 12 C from Model-G (red dashed curve), Model-Ni for i = 1, 2, 3, 4 (black dashed curve) and Model-N5 (blue dashed curve) have been presented in the upper panel, where those of p, n and α from the deexciation of final-state nuclei in the NC interactions have been denoted by blue, magenta and orange solid curves, respectively. The gray band along the black dashed curve represents 1σ deviation from the average of the predictions from Model-Ni for i = 1, 2, 3, 4. The colored band along the solid curve stands for 1σ deviation from the average of all six models, for which the deexcitation processes are the same. In the lower panel, the energy spectrum of γ from the deexciation of final-state nuclei has been shown as gray histograms. The production rates are shown with respect to the kinetic energies of relevant particles, and the QEL with EMC effects dominates the production. the neutron multiplicity increases. To visualize the model dependence, we have also shown the 1σ variation of all the six representative models as the gray shaded band.
• Due to the high neutron tagging efficiency, NC interactions associated without neutrons or  Figure 7: The event rates for the NC interactions of atmospheric neutrinos (with energies ranging from 100 MeV to 10 GeV) with 12 C nuclei in the exclusive channels, which are categorized by the associated neutron multiplicities. The predictions from Model-G, Model-Ni for i = 1, 2, 3, 4 and Model-N5 are denoted as dots, squares and triangles, respectively. Note that the squares with error bars stand for the mean value of the predictions from Model-Ni (for i = 1, 2, 3, 4) and 1σ deviation. In all the panels, the raw rates are extracted from those in Fig. 5, and the gray bands represent the 1σ uncertainties from all six models. In the rightmost panel of the second row for a single neutron, the total rate (16.5 ± 2.8) kt −1 yr −1 (black) and that (3.1 ± 0.5) kt −1 yr −1 for the DSNB-like events (blue) should be read off from the vertical axis with red ticks.
with more than one neutrons can be rejected as the background events in DSNB searches. Only single-neutron events have to be considered, for which the rates are given in the second row of Fig. 7. The rightmost panel of this row provides the total rate of all single-neutron interactions as (16.5 ± 2.8) kt −1 yr −1 . When regarding only events in the prompt visible energy range from 11 MeV to 30 MeV, most relevant for DSNB searches, the rate reduces to (3.1 ± 0.5) kt −1 yr −1 . Notice that the associated uncertainty is about 20%, representing the model variations of neutrino interactions. If the extra uncertainty of 15% from the calculations of atmospheric neutrino fluxes is simply added in quadrature, one may obtain the overall uncertainty of 25% for the NC backgrounds.
• As indicated in the fourth row of Fig. 7, the GENIE generator produces significantly higher rate of the channels with more than two neutrons. In future large LS detectors, the neutron multiplicity distribution can be measured and will be very useful to scrutinize the nuclear models. In addition, the decays of unstable final-state nuclei may provide unique signatures that allow for in situ measurements of the NC backgrounds. In a companion paper, we shall carry out a systematic study of such in situ measurements in the LS detectors [75].
It is worthwhile to emphasize that the conversion of the kinetic energies of final-state particles in the NC interactions to the visible energies of the final events in the LS detectors is nontrivial. Such a conversion is process-dependent and relies on the energy deposition of different types of final-state particles. In our calculations of possible NC backgrounds for DSNB signals, we have employed the method of Monte Carlo simulations based on the software GEANT4 (4.9.4) [76], without including a specific detector geometry and the optical processes. The quenching effect of LS is considered according to the description in the Appendix of Ref. [19]. After converting the event rates in terms of kinetic energies of final-state particles into those of the quenched visible energies, we obtain the total rate of the IBD-like signals of the atmospheric neutrino NC interactions in the energy range from 11 to 30 MeV. Finally, one should notice that the typical rate of DSNB signals in LS detectors is around 0.3 kt −1 yr −1 [19], which is about one order of magnitude smaller than the atmospheric NC background rate (3.1 ± 0.5) kt −1 yr −1 . However, a pulse shape discrimination offers a technique of efficient background suppression for non-positron prompt events, which will be necessary to ultimately achieve an unambiguous discovery of the DSNB signals [19].

Proton Decay
Then, we turn to the detection of proton decay in the LS detectors, for which the signals and relevant backgrounds are quite different from those for DSNB. For the wCh detectors, the decay channel p → e + + π 0 offers the clearest signature for proton decay [30,31]. However, as shown in Refs. [77][78][79], the supersymmetric minimal SU(5) GUTs may predict a highly suppressed decay rate for p → e + + π 0 , but an appreciably large one for p → K + + ν, where the flavors of ν depend on specific models [80] and are irrelevant for our following discussions. With a null signal in this channel, the SK experiment has placed a lower limit on the partial proton lifetime [81,82] and the latest result is τ p /B(p → K + ν) > 5.9 × 10 33 yr at the 90% confidence level [29], where τ p is the proton lifetime and B(p → K + ν) is the branching ratio.
As for proton decay in the channel p → e + + π 0 in the LS detectors, the immediate scintillation light from e + 's and the instantaneous decay of π 0 → 2γ lead to a single prompt signal, which can easily be contaminated by numerous backgrounds. In contrast, the decay channel p → K + + ν serves as the most sensitive probe for proton decay in the LS detectors. The main features of this decay channel and possible backgrounds are summarized as follows.
1. The K + meson decays quickly (with a lifetime τ K + = 12.4 ns) into six channels, namely, µ + ν µ (63.56%), π + π 0 (20.67%), π + π + π − (5.58%), π 0 e + ν e (5.07%), π 0 µ + ν µ (3.35%) and π + π 0 π 0 (1.76%), where the corresponding branching ratios are given in the parentheses [83]. In order to identify the K + signal, one has to analyze these decay products and their signals in the LS detectors. The most important decay modes in LS are K + → µ + ν µ and K + → π + π 0 , because they produce a signal of three-fold coincidence. See, e.g., Refs. [19], [32] and [33], for earlier discussions. As K + is a heavy and highly ionizing charged particle in LS, it will loose its kinetic energy rapidly, producing a first prompt scintillation signal. In either decay mode there is a shortly delayed signal (∼ τ K + ) from the daughter particle(s). In the first decay mode, the kinetic energy of µ + constitutes the shortly delayed signal. Then the final-state µ + decays into e + ν e ν µ with a proper lifetime of τ µ + = 2.2 µs, which is long enough to be separated from the preceding two signals, and the Michel electron from the µ + decay can be reconstructed as the third delayed signal. In the second decay mode, the neutral pion π 0 instantaneously decays into two gamma rays, while the charged pion π + decays primarily into µ + ν µ (with a proper lifetime of τ π + = 26 ns). The deposition of the total energy of the π 0 decay and the kinetic energy of π + are indistinguishable from each other, and will constitute the shortly delayed signal of this mode. Moreover, the daughter µ + has low kinetic energy (∼4 MeV) and its signal will be submerged in the tail of the shortly delayed signal. As a consequence, similar to the first decay mode, the signature of the second decay mode also represents a three-fold coincidence of the prompt signal, a shortly delayed signal and a single Michel electron. In the following discussions, we shall take the proton decay p → K + + ν and the subsequent decays K + → µ + ν µ or π + π 0 as the target signal.
2. Based on the above discussions about the signal of p → K + + ν, we now discuss the possible background in the NC interactions of atmospheric neutrinos with 12 C in the LS detectors. First, for the K ± production, there will be a three-fold coincidence signal from the K ± decay, but the associated production of other particles also contributes to the first prompt signal, which tends to have a relatively higher prompt energy compared to that of the three-fold signature of proton decay. Second, if a single neutral pion π 0 is produced, there will be no three-fold coincidence signals because of the prompt decay of π 0 into two γ's. Third, for the π ± production, the first prompt signal comes from the scintillation light of the π ± kinetic energy and other associated particles, and the third delayed signal is the possible Michel electron from the muon decay. Therefore, to mimic the proton decay signal, the key is to have a suitable shortly delayed signal, which can be produced from the hadronic interactions of high energy π ± and nucleons in LS if they are separable from the prompt scintillation signal. The possibility to achieve the separation heavily depends on the detailed simulation of particle interactions in LS, which will not be further explored in the current paper. Instead, we shall calculate the total rate of single π ± production as a benchmark number. Finally, if one or more neutrons are produced in the NC interactions, high efficiency  rejection power can be achieved using the typical neutron capture signature. Therefore we shall consider the possible production of single π ± or K ± but without neutrons as our focus in the following part. Note that a single negatively-charged pion π − will decay mainly into µ − ν µ . Unlike µ + which only decays, the stopped µ − in LS will be captured into an atomic orbital and then it can either decay or undergo nuclear capture. The muon capture on 12 C in LS may emit one or more neutrons. A recent measurement of muon capture on light nuclear isotopes in LS can be found in Ref. [84].
3. In Fig. 8, neutron multiplicity distributions of the NC interactions of atmospheric neutrinos (with energies ranging from 100 MeV to 10 GeV) with 12 C nuclei have been categorized by the multiplicities of the charged pions in the COH, RES, and DIS processes. Note that these multiplicity distributions can be measured in future large LS detectors and will be very useful to scrutinize the nuclear models. For proton decay, the visible energies of final-state particles will be as high as several hundred MeV, so the corresponding neutrino energies will be much higher, for which the COH, RES and DIS processes will be relevant. However, the COH process produces only π 0 that will not contribute to the background. Moreover, the three-fold coincidence signature consists of only one Michel electron from muon decay, so the processes associated with multiple π ± 's or neutrons are distinguishable from proton decay. For this reason, those with no neutron and only one charged pion are able to contaminate the proton decay signal. Note that (1 − 2)% of the 1π ± 0n processes are accompanied by a charged kaon, but they can be rejected by the criterion of only one Michel electron. As a consequence, the rates for the potential background 0K ± 1π ± 0n 2 from Model-G, Model-N1 and Model-N5 are found to be 2.1 kt −1 yr −1 , 2.4 kt −1 yr −1 and 1.8 kt −1 yr −1 , respectively, where the average and 1σ deviation of these models is (2.10 ± 0.27) kt −1 yr −1 .
4. In Fig. 9, the multiplicity distribution of the produced charged kaons from the COH, RES, and DIS processes is shown. To illustrate the model dependence, we have performed the calculations by using Model-G, Model-N1 and Model-N5, the results for which are represented by red, blue and black histograms, respectively. The average fractions of NC events from the COH, RES, and DIS processes without charged kaons and with a single charged kaon are about 98.9% and 1.0%, respectively. The low rate of kaon production is due to the fact that a relatively high energy for the incident neutrino is required and the event rate above a neutrino energy of 1 GeV decreases rapidly, as shown in Fig. 3. The event rate for a single K ± production (i.e., 1K ± 0π ± 0n) can be estimated as 0.057 kt −1 yr −1 , 3.8 × 10 −3 kt −1 yr −1 and 3.1 × 10 −3 kt −1 yr −1 from Model-G, Model-N1 and Model-N5, respectively. The average and 1σ deviation of these three models is then (0.021 +0.025 −0.021 ) kt −1 yr −1 , where the uncertainty comes from the large variation of model predictions between GENIE and NuWro (see the right panel of Fig. 10).  Figure 10: The differential event rate for the production of π ± and K ± with respect to the energy transfer, where the results of any number of π ± and K ± are given in the left panel and those of the exclusive one single π ± or K ± are shown in the right panel. Note that only Model-G and Model-N1 have been implemented to illustrate the difference between GENIE and NuWro generators. 5. The differential rates of the interactions associated with the production of π ± and K ± are shown in Fig. 10 with respect to the energy transfer, which is defined as the energy difference between incoming and outgoing neutrinos. Some important observations can be made. First, the π ± production is predominantly stemming from the RES process with a relatively low energy transfer in both GENIE and NuWro models, whereas the K ± production needs much higher energy transfers. Second, the energy thresholds for π ± production in the RES and DIS processes turn out to be the same in the GENIE model. But this is not the case for NuWro. In addition, in NuWro, there is no contribution to the K ± production from the RES process, but the differential rates of total K ± production are compatible between GENIE and NuWro. Third, comparing between the curves in the left and right panels, we can see that the single π ± process dominates over the multiple π ± in both GENIE and NuWro. Nevertheless, while the single K ± production without neutrons and π ± 's is most important in the total K ± production of GENIE, it becomes unimportant in NuWro. Such differences between GENIE and NuWro might arise from different treatments of the internal parameters in nuclear models.
As indicated in Fig. 10, the energy transfer for kaon production is at least 500 MeV, so one can apply an energy cut to efficiently reduce the background associated with kaons. In a realistic LS detector, the conversion of the energy transfer to the visible energy is nontrivial and heavily depends on the types of final-state particles. In this paper we directly apply the energy cut on the energy transfer to obtain a rough estimation, and leave the detailed detector response modeling for future works by the experiments. We take the energy transfer range from 150 MeV to 650 MeV for the single π ± production (0K ± 1π ± 0n), and the rate is reduced to (0.61 ± 0.14) kt −1 yr −1 . On the other hand, since a large proportion of the rest mass of K ± will be carried away by invisible neutrinos from kaon decay, we require the energy transfer to be ranging from 150 MeV to 1 GeV for the single K ± production (1K ± 0π ± 0n). Then the corresponding rate will be (4 +5 −4 ×10 −4 ) kt −1 yr −1 .
Here we would like to remind that not all the single π ± events, but only those being recognized as the three-fold coincidence signal will constitute the background of proton decay. In this aspect, we need careful treatments on the particle interaction in LS and detector simulation. The pulse shape discrimination of the single π ± events and even the advanced deep learning techniques would be critical to obtain a reasonable background level for proton decay. Such further exploration will be reported in the future works. Moreover, it should be noted that the energetic neutrons and protons induced by the neutrino interactions with 12 C may produce secondary π ± 's. Such interactions can also contribute to the backgrounds for the searches of proton decay. Additionally, we focus on the NC interactions of atmospheric neutrinos and leave out the CC interactions, which should also be taken into account in more realistic background analysis. It is clear that our calculation strategy of the NC background will play an important role in such an analysis as well.

Summary
In this paper, we have performed a systematic calculation of the NC background induced by atmospheric neutrino interactions with the 12 C nuclei in the LS detectors, which are expected to be crucially important for the experimental searches for DSNB and proton decay. In our calculations, the up-to-date fluxes of atmospheric neutrinos at the JUNO site provided by the Honda group are used. As for the neutrino-nucleus interactions, we have chosen six representative models from the Monte Carlo neutrino event generators GENIE and NuWro and take the model variation as the systematic uncertainty in the background prediction. Then, a statistical configuration model of 12 C is implemented to determine the probability distribution of the excited states of final-state nuclei produced in the NC interactions. The deexcitation processes of these nuclei are handled by using TALYS. Taking account of neutrino interactions and the deexcitation processes, we are able to compute the event rates of the exclusive processes with final-state protons, neutrons, γ's and α's and the corresponding energy spectra. These processes are relevant for the experimental searches for the IBD signals of DSNB ν e 's. In addition, the production of multiple neutrons, kaons and pions in the high-energy region is considered, and its implications for the detection of proton decay are investigated.
For the detection of DSNB in the LS detectors, the golden channel is the IBD process. After rejecting the irreducible backgrounds from reactor ν e 's in the low-energy region and atmospheric ν e 's in the high-energy region, one has to focus on the narrow window of the visible energies 11 MeV E vis 30 MeV. In this energy window, the NC backgrounds induced by atmospheric neutrinos will be dominant. As the visible energy is relatively low, the QEL process of neutrino-12 C interactions turns out to be most important. Though the neutron tagging efficiency is intrinsically high for LS detectors, the NC backgrounds with one neutron production lead to the IBD-like signals and thus are irreducible. With numerical simulations with the chosen neutrino event generators (i.e., GENIE and NuWro) and TALYS for deexcitation, we have found the event rate for the exclusive processes with one neutron production to be (16.5 ± 2.8) kt −1 yr −1 in the whole range of visible energies, where the uncertainty originates from different models for neutrino-12 C interactions. When restricted into the energy window 11 MeV E vis 30 MeV of interest, the event rate is (3.1 ± 0.5) kt −1 yr −1 . Further reduction of the NC backgrounds from atmospheric neutrinos can be achieved by using pulse shape discrimination.
For the proton decay in LS detectors, the decay mode p → K + + ν is most promising. The three-fold coincidence among the energy deposition of K + , the shortly delayed signal of µ + from K + → µ + ν µ (or π + π 0 with π + → µ + ν µ ), and the single Michel electron from µ + → e + ν e ν µ . We identify the relevant backgrounds from NC interactions of atmospheric neutrinos as the production of single charged pion π ± or kaon K ± . When requiring the energy transfer to be ranging from 150 MeV to 650 MeV for the single pion production (i.e., 0K ± 1π ± 0n) and from 150 MeV to 1 GeV for the single kaon production (i.e., 1K ± 0π ± 0n), we find that the background rates are (0.61 ± 0.14) kt −1 yr −1 and (4 +5 −4 × 10 −4 ) kt −1 yr −1 respectively, where the uncertainty arises from the model variation of neutrino interactions. Therefore, further examination of all possible backgrounds for proton decay in the LS detectors must pay a particular attention to the NC events from atmospheric neutrinos, as we have demonstrated.
Apart from the future ordinary LS detectors [19,20], we will have other large-scale detectors with advanced techniques based on water [85], water-based LS [86], liquid-Argon [87,88]. The experimental searches for DSNB and proton decay are also primary physics goals for these detectors, and the NC backgrounds induced by atmospheric neutrinos will be relevant. We believe that the calculations performed in the present work will be not only useful for the LS detectors, but also instructive for the parallel studies for other types of detectors.