Pulse shape discrimination technique for diffuse supernova neutrino background search with JUNO

Pulse shape discrimination (PSD) is widely used in particle and nuclear physics. Specifically in liquid scintillator detectors, PSD facilitates the classification of different particle types based on their energy deposition patterns. This technique is particularly valuable for studies of the Diffuse Supernova Neutrino Background (DSNB), nucleon decay, and dark matter searches. This paper presents a detailed investigation of the PSD technique, applied in the DSNB search performed with the Jiangmen Underground Neutrino Observatory (JUNO). Instead of using conventional cut-and-count methods, we employ methods based on Boosted Decision Trees and Neural Networks and compare their capability to distinguish the DSNB signals from the atmospheric neutrino neutral-current background events. The two methods demonstrate comparable performance, resulting in a 50\% to 80\% improvement in signal efficiency compared to a previous study performed for JUNO~\cite{JUNO:2015zny}. Moreover, we study the dependence of the PSD performance on the visible energy and final state composition of the events and find a significant dependence on the presence/absence of $^{11}$C. Finally, we evaluate the impact of the detector effects (photon propagation, PMT dark noise, and waveform reconstruction) on the PSD performance.


Introduction
The diffuse supernova neutrino background (DSNB) is a low-energy steady neutrino flux that permeates the universe and originates from the cumulative emissions of all past core collapse supernovae [2][3][4].The study of DSNB neutrinos holds great promise for advancing our understanding of the evolution of the universe, core-collapse SNe and the fundamental properties of neutrinos.The detection of DSNB neutrinos is challenging due to their low flux and energies and the overwhelming background from other neutrino sources.Nevertheless, several experiments, including Super-Kamiokande (SK) [5][6][7][8], KamLAND [9,10] and Borexino [11], have been improving their search techniques to get a hint of this flux.The next generation, which includes experiments like Jiangmen Underground Neutrino Observatory (JUNO) [12,13], Super-Kamiokande Gadolinium experiment (SK-Gd) [14] and Hyper-Kamiokande (Hyper-K) [15], will be key for its discovery.
JUNO will be the largest ever liquid scintillator (LS) experiment and is planned to be in operation for the physical runs by end of 2024 [12].It is loaded with 20 kton of LS in an acrylic sphere with a diameter of 35.4 m.The DSNB flux can be detected via the inverse beta decay (IBD) channel, especially within the energy window around 10-30 MeV, as the reactor neutrinos will dominate below 10 MeV and the atmospheric neutrino charged current (CC) events will become significant above 30 MeV.The largest background within the energy range of 10-30 MeV comes from the atmospheric neutrino neutral current (NC) interactions on 12 C, whose event rate in the whole detector volume is around 60 per year, in comparison with a rate of 2-3 events per year for the DSNB signal.Therefore, powerful background suppression techniques are required to further increase the signal-to-background ratio, and obtain a significant sensitivity for the DSNB observation.
The pulse shape discrimination (PSD) technique is an excellent method used in LS detectors to distinguish different types of particles based on their energy deposition patterns.This technique relies on the fact that different types of particles produce distinct time profiles of the scintillation light in the LS.Therefore, a precise determination of the scintillation emission time for different particles is essential for the PSD performance.In this work, we present a detailed study on the PSD technique that has been used in the context of DSNB sensitivity studies in Ref. [13].
Initially, we employ methods from Ref. [13,16] to predict the energy spectrum distribution of final-state particles for the DSNB signals and atmospheric neutrino NC backgrounds, which serve as inputs for the detector simulation.Subsequently, we perform a comprehensive simulation using the JUNO offline software framework [17,18], incorporating all detector and electronics effects, as well as waveform and event reconstruction.Various factors, including the LS optical model, photon arrival time resolution, and PMT dark noise, are considered in the simulation.Following this, we conduct a systematic analysis of the characteristic photon emission time spectra of positron and non-positron events, extracting them as inputs for different PSD methods.Instead of more conventional cut-and-count methods such as the tail-to-total method [19] and the Gatti method [20], two PSD methods based on a Boosted Decision Tree (BDT) and a Neural Network (NN), respectively, are used for background suppression.Additionally, we evaluate the energy dependence of the PSD performance and, for the first time, evaluate the corresponding efficiency to the DSNB search in JUNO.Finally, we perform a comprehensive investigation of the limiting factors influencing the PSD performance in JUNO by incrementally introducing these factors.
The paper is organized as follows: the DSNB signals and atmospheric neutrino NC backgrounds are reviewed in Sec.2; the pulse shape discrimination methods are introduced in Sec.3; then, the performances of both methods are presented and discussed in Sec.4; a summary will be presented in Sec. 5.

DSNB signals
The DSNB signal prediction, based on Ref. [13], depends on a variety of crucial ingredients [21][22][23].The first one, which serves as a link to the cosmic history of star formation, is the cosmological SN rate as a function of the progenitor mass and redshift.The second ingredient is the average energy spectrum of SN neutrinos.The last one is the contribution of the failed SNe [21,22], which may alter the DSNB signal and will feature a hotter neutrino energy spectrum compared to the neutron-star-forming SNe (i.e., successful SNe).A reference DSNB model is adopted in this paper to explore the performance of the pulse shape discrimination technique as described in details in Sec. 3. In the reference DSNB model, the failed SNe fraction (f BH = 0.27), the SN Rate at red-shift z = 0 (R SN (0) = 1.0 × 10 −4 yr −1 Mpc −3 ) and the average neutrino energy (⟨E⟩ = 15 MeV) are used to predict the DSNB flux.For details see Ref. [13].
Thanks to the relatively large cross section, the IBD reaction, ν e + p → e + + n, is the ideal channel for the detection of the DSNB.The interaction produces a positron prompt signal and a delayed neutron capture signal in the LS.In addition, the time coincidence of the prompt and delayed signals arising from the final-state e + and n, respectively, perfectly eliminates single-event backgrounds such as signals from the radioactivity of detector materials, the decays of cosmogenic isotopes and the recoil electrons from solar neutrino interactions.We need to take into account the IBD cross section, the target mass, and the detector response in order to calculate the expected DSNB energy spectrum that would be observed in JUNO.We take the free proton number in the JUNO LS as 7.15 × 10 31 kton −1 [1], assuming a 20 kton of LS.The differential IBD cross section is taken from Ref. [24].The visible energy of the DSNB signals in JUNO is expected to be below 100 MeV.

Atmospheric neutrino NC backgrounds
In the search for DSNB, the relevant background source includes reactor νe 's, NC and CC interaction of atmospheric neutrinos, fast neutrons and long-lived cosmogenic βn-emitters 9 Li and 8 He.Among these, the dominant backgrounds within the energy window ranging from 10 MeV to 30 MeV are the atmospheric neutrino NC background events.It is essential to have a comprehensive understanding of the NC background in order to develop analysis techniques for its suppression.
As described in Ref. [16], we have used dedicated simulations to predict spectra and rates of background events resulting from atmospheric neutrino NC interactions with 12 C nuclei in the LS.Specifically, we employed a representative nuclear model from the GENIE generator (version 2.12.0) [25] to simulate neutrino-nucleus interactions within the JUNO LS detector.The GENIE model incorporates a nucleon axial-vector form factor parametrization with an axial mass of M A = 0.99 GeV, determined from deuterium measurements [26].To model the nuclear structure, the GENIE model incorporates the relativistic Fermi gas (RFG) model with "Bodek-Ritchie" modifications [27].The package TALYS (1.8) [28] is used to describe n's, p's, α's, γ-rays and other particles from the deexcitations of the final-state nuclei.Table 1 summarizes the rates of DSNB signals and NC interactions between atmospheric neutrinos on 12 C in inclusive and exclusive channels.For the DSNB search, NC interaction channels with single neutron emission are the major background as they can form prompt-delay coincidences that mimic IBD signals.However, while the visible energy (number of photoelectrons) detected is in the same range as that of DSNB-induced positrons, the prompt signals are usually induced by a collection of neutrons, protons, alphas and de-excited gammas (c.f."Prompt signal" in Table 1).The most relevant NC backgrounds are the quasielastic scattering (QEL) interactions, in which only one neutron is produced, causing a prompt signal with energy below 100 MeV.The expected interaction rates are shown as "Raw" in Table 1.A set of selection cuts is used for filtering out IBD events from DSNB signals where the visible prompt energy window is between (12,30) MeV.After the IBD selection, the typical predicted rate of DSNB signals in the optimized window is around 0.14 kt −1 yr −1 , which is about one order of magnitude smaller than the atmospheric NC background rate 2.56 kt −1 yr −1 .Therefore, powerful background suppression techniques for non-positron prompt events are necessary to achieve an unambiguous discovery of the DSNB signals.
3 Background discrimination with pulse shape discrimination The prompt signal in DSNB IBD events is a e + accompanied by its annihilation and the consequent emission of two γ's (0.511 MeV each), compared to the complicated composition of γ, n, α, and p in NC background events.The intrinsic scintillation timing for γ/e ± is different from that for α/p/n as the more heavily ionizing particles produce a higher fraction of slow scintillation photons [29].In addition, the spatial topologies of energy deposition of γ's and α/p/n are different and introduce a different level of broadening effects on the timing spectrum, which can help particle identification.Therefore, the timing information is used to develop PSD methods to discriminate the DSNB signals and NC backgrounds.Detailed descriptions on their implementations are given below.

MC dataset and software framework
The official JUNO offline software [17] is employed to simulate the DSNB signals and atmospheric neutrino NC backgrounds.The simulation process includes event generation, detector simulation, electronics simulation, waveform reconstruction, and event reconstruction.Below, we provide detailed information on the configurations and procedures involved in each step of the simulation chain.
Event generator: The event generator primarily uses the DSNB spectrum described in Sec.2.1.Additionally, the final states of atmospheric neutrino NC background events, detailed in Sec.2.2, are also computed and incorporated.Both data samples are uniformly generated within the JUNO LS volume.
Detector simulation: We perform a comprehensive simulation of the detector utilizing the most up-to-date design geometry [12].The energy deposition of particles is simulated using Geant4 [30] (version 4.10.p02),while the properties of the LS [31][32][33] and scintillation photon emission time profiles [34] are based on experimental measurements.Optical processes, including photon absorption, re-emission, scattering, refraction, and total reflection, are implemented for photons propagating in the LS [35].Based on measurements of the photon detection efficiency [36], for every 1 MeV of deposited energy at the center of the detector, an average of approximately 1350 photo-electrons (p.e.) is detected by the photomultiplier tubes (PMTs).A "hit" refers to the p.e. detection by a PMT.In this study, all hits collected by the PMTs, not just the first one, are utilized.

Electronics simulation:
The electronics response is modeled in the simulation including the transition time spread (TTS) and the dark noise rate of the PMTs.In JUNO, two types of high quantum efficiency 20-inch PMTs [37]   that of dynode PMTs.25,600 3-inch PMTs are installed in the gaps between the 20-inch PMTs.However, their data is not utilized in this study due to a low photo coverage of only 3% [38].
Reconstruction: The waveforms of the PMTs are generated using realistic measurements of single p.e. (SPE) shapes.These waveforms are then digitized and stored for offline data reconstruction.During this process, the time and charge information of PMT hits is reconstructed.Furthermore, a reconstruction algorithm [39] is applied to estimate the total event energy and vertex, leading to in a reconstructed vertex resolution that is greater than 15 cm.
The time-of-flight (TOF) is calculated using the reconstructed vertex, assuming an optical path that follows a straight line and taking into account the effective refractive indexes as described in Ref [40].The photon emission time (PET) is reconstructed by subtracting the TOF from the hit time obtained through the waveform reconstruction.It should be noted that waveform reconstruction can only determine the arrival time of the first photon when multiple single-photon waveforms overlap in a PMT within 1 ns.To mitigate the impact of this effect on the time profile, the hit time is weighted by the charge for the subsequent PSD methods.In Fig. 1, we present the average normalized reconstructed PET distributions for both the DSNB signal and the NC background events.The peak time for each event is aligned to 0 ns.Importantly, the NC event distributions display elongated tails in contrast to the DSNB signal distributions.Additionally, the peaks in the NC event distributions exhibit a relatively shallower shape compared to the DSNB signal distributions.This difference, that will be exploited by the PSD techniques, can be ascribed to the extended travel distance of neutron recoils within the prompt signal of NC events.The PSD algorithms are specifically designed to differentiate between signal and background events based on their timing spectral shape differences.However, several factors can adversely impact the performance of PSD.Further discussion on these factors will be presented in Sec.4.3.

Boosted Decision Tree (BDT) method
In our study, we have employed the BDT method [41] for particle identification, recognizing that multivariate machine learning techniques offer a promising alternative to conventional cut-andcount methods [19].Tab. 2 provides a summary of the input variables used in the BDT method, based on the features observed in time profiles of different particles.Fig. 2 displays the input variable distributions along with the correlation matrix among these variables.The definitions of the variables are explained below.
• w r and w f respectively denote the rising and falling times within the interval [-20, 20] ns.These ascending and descending segments within the waveforms contribute to 10% of the total area encompassed by this time range.These two variables effectively capture the distinctive features found in the peak regions of the temporal profiles.The calculations for w r and w f are based on hits from the dynode 20-inch PMTs.The utilization of these specific PMTs preserves intrinsic peak profile characteristics, facilitated by their better time resolution.
• R p and R t are defined as the fractions of the pulse area within [20,1100] ns and [400, 1100] ns, respectively.These two variables are used to further distinguish peak shape and tail shape of different particles, and are also calculated using hits from dynode PMTs.
• The tail distributions of the time spectra follow a multi-exponential probability density function (PDF).To model this distribution, we have developed a two-exponential fit as shown in Eq.(3.1), using the un-binned maximum-likelihood method.The fit is applied to the time spectra within the range of [40,1100] ns.
The exponential terms in Eq.(3.1) represent the fast and slow components of the time spectrum.The number of hits associated with these components is denoted as N , with corresponding decay times τ 1 and τ 2 .The parameter η represents the fraction of the fast component within the time spectrum range of [40,1100] ns.The recorded time spectrum includes dark noise hits from PMTs that can not be distinguished from physics hits.In order to accurately estimate the dark noise rate (n dark ), a fitting parameter is introduced, using the time spectrum within the range of [-200, -100] ns, where dark noise hits are present.
Given the restricted statistical dataset, the fit value of n dark displays correlation with other input variables, particularly the effective slow component τ 2 .Nonetheless, the incorporation of n dark contributes to a more accurate representation of the time spectrum.To improve the statistical precision, all hits from both dynode and MCP PMTs are used in the fitted time profile.
• R 3 is defined as the radius cubed of the reconstructed event position, R being the distance from the reconstructed interaction vertex to the center of the detector.This parameter is also used as an input variable in the BDT method to mitigate the impact of the TOF bias on the PSD performance.As shown in Fig. 2, it is evident that the variables associated with the peak area or those in close proximity to the peak area exhibit correlations with the positions.Hence, the variable R 3 is included to enhance the discrimination capabilities of w f , R p , τ 1 , and η.
The BDT method, integrated into the Toolkit for Multivariate Analysis (TMVA1 ) [42], is a widely used supervised learning algorithm in particle physics for classification tasks.It involves combining multiple decision trees to create a strong predictive model.Each decision tree makes a series of binary decisions based on input features (discriminating variables, as shown in Tab. 2), eventually leading a prediction or classification (e.g., signal or background).The BDT model is trained using 1 million signal events and 1 million background events.According to the importance ranking report from TMVA, the variables η, τ 2 , and τ 1 are identified as the most important ones.Additionally, other variables are introduced to further enhance the PSD performances.The distribution of the test samples is consistent with that of the training samples, indicating no evidence of over-training.The results will be discussed in Sec.4.1.

Neural Network (NN) method
An alternative discriminator has been developed based on NN methods.It is worth noting that a similar NN method has been implemented for α/β pulse shape discrimination in Borexino [43].Instead of extracting features from the event time distributions, the latter are directly used as the input.All the hits in an event are binned in time and an array of the number of hits in each bin is taken as the input.Due to the exponential decay characteristic of LS and the need to gather sufficient statistics within each bin, an uneven binning strategy has been implemented.The binning scheme incorporates a carefully designed range of bin widths, spanning from finer intervals to broader ones.The event time distributions are weighted by the reconstructed charge of each hit to accurately accommodate the influence of multiple p.e. hits, following the same methodology as the BDT method outlined in Section 3.2.To mitigate the broadening effect caused by the relatively broad SPE response of the MCP PMTs, unweighted time distributions are also provided as a supplementary input.Moreover, to capture potential vertex-dependent attributes, the reconstructed vertex information R 3 is appended to the input array.These three components are then concatenated into the one-dimensional input for the network.
A simple NN with one fully connected layer is chosen for event classification.The network  is implemented with the MLPClassifier (multi-layer perceptron classifier) module from the scikitlearn toolkit [44].The fully connected layer consists of 100 nodes.An illustration of the input and network structure is shown in Fig. 3.The training dataset matches the one employed in the BDT method.Approximately seven hundred thousand samples of DSNB signals and atmospheric neutrino NC backgrounds are used for training the model, while the remaining events are reserved for testing purposes.

Results and discussions 4.1 Comparison of BDT and NN methods: PSD performance
The left panel of Fig. 4 presents the PSD distributions obtained from the BDT and NN methods.Both methods exhibit similar trends in the one-dimensional PSD distributions, demonstrating clear separation between signal and background events.In accordance with the requirements of the latest DSNB analysis, we choose a residual fraction of 1% for the NC background.Under this condition, within the prompt energy range of [11,30] MeV, a signal efficiency of around 80% is achieved.The right panel of Fig. 4 shows the area under the Receiver Operating Characteristic (ROC) curve distributions as a function of the prompt energy for the BDT and NN methods.The background inefficiency (left) and signal efficiency (right) as a function of the prompt energy with a PSD cut retaining 1% total background.The legend in the right plot is identical to the one in the left plot.In each prompt energy bin, the cut on NN discriminator is chosen to achieve the same background inefficiency as BDT method.The shaded band in the BDT and the error bar in the NN represent statistical errors.The BDT and NN methods give similar energy dependence and separation.Atm-ν NC w/ 11 C events with prompt energy below 18 MeV are more difficult to reject.
Both methods exhibit a similar energy-dependent trend, with the NN method demonstrating better PSD performance at energies below 20 MeV.
The PSD efficiency is energy-dependent, as illustrated in Fig. 5.The figure includes both BDT and NN.Again, a PSD using either BDT or NN discriminators is applied, choosing to maintain an NC background fraction of 1%.The residual NC background events are then categorized into two types based on their final states: one type corresponds to interactions with 11 C in the final state (atm-ν NC w/ 11 C), the other represents interactions without 11 C (atm-ν NC w/o 11 C).
As expected, the PSD performance improves with rising visible energy due to the large photon statistics available.Moreover, both discriminators are more efficient identifying atm-ν NC w/o 11 C events, since apart from n, the prompt energy of these events also originates from α, p, and other particles (c.f.Tab. 1).The atm-ν NC w/ 11 C events are more difficult to distinguish because the single high-energy neutron in the final state can produce gamma rays in inelastic scattering with carbon nuclei [45].Inelastic scattering of high energy neutrons on 12 C excites the nucleus to excited states, with high energy gamma emission from de-excitation.The electron-like gamma contribution to the prompt energy makes it indistinguishable from the DSNB signal events, which explains the feature of increased residual background below 18 MeV for atm-ν NC w/ 11 C in Fig. 5.More details will be discussed in Sec.4.2.
To facilitate the comparison, varying PSD cuts are implemented on the NN discriminator, ensuring that the PSD inefficiency for the atm-ν NC w/ 11 C background aligns consistently with the BDT method across different energy bins.Under these conditions, the NN method exhibits a similar pattern for atm-ν NC w/o 11 C background, but features an improved signal efficiency for the DSNB signal sample.

Energy dependent efficiency for atm-ν NC with 11 C
A clear energy dependence in the PSD efficiency is observed for the atm-ν NC w/ 11 C channel.To investigate this, the eγ ratio R eγ is defined as the fraction of energy contributed by e ± and γ particles over the total visible energy.The prompt energy of the NC background events is a composite of e ± , γ, α, p signals.When more photons are generated from e ± , γ, the hit time distributions get closer to those of DSNB signals.Fig. 6 illustrates the distributions of R eγ as a function of energy for atmospheric neutrino NC background events with and without 11 C.
The left panels of Fig. 6 display the scenarios involving the primary particles from the initial neutrino interaction, where R pre eγ is defined as the ratio of the kinetic energy of e ± and γ particles to the total kinetic energy of all final state particles, excluding the neutrino.According to the deexcitation model used in predicting NC background events, the atm-ν NC w/ 11 C events have no deexcitation γ's that results in R pre eγ being equal to 0. Atm-ν NC w/o 11 C events could have deexcitation γ's originating from excited residual nuclei leading to nonzero R pre eγ , but the majority of events have R pre eγ close to 0. The secondary particles from subsequent interactions of product particles in the detector, are shown in the right panels of Fig. 6.The R post eγ is defined as the ratio of the visible energy contributed by e ± and γ particles to the total visible energy.This calculation is performed using the MC data obtained in the post-simulation stage.In comparison to the scenarios involving the primary particles, for the atm-ν NC w/ 11 C events, the distribution of R post eγ exhibits two additional bands, which are caused by two types of deexcitation γ rays with energies of 4.439 and 9.641 MeV.This phenomenon arises because the energetic fast neutrons from the atm-ν NC w/ 11 C events can excite the 12 C nucleus to an excited state, and the subsequent deexcitation γ emission contributes a significant portion of the emitted photons.Similar event typologies have been observed for neutron interactions with oxygen in water Cherenkov detectors [46,47].
To gain a deeper understanding of the differences in R eγ before and after the simulation, we  use TALYS to calculate the cross sections of exclusive n− 12 C reactions at various incident neutron energies.These results are depicted in Fig. 7.At kinetic energies of the fast neutron below approximately 10 MeV, the dominant interaction between n and 12 C is elastic scattering (n 12 C).However, as the energy exceeds approximately 10 MeV, the inelastic scattering cross section of n 12 C increases significantly, competing with the elastic scattering process.The inelastic process 12 C(n, nγ) 12 C dominates the cross section up to around 20 MeV.Subsequently, processes such as 12 C(n,npγ) 11 B become important, and the emission of kicked-out p and α particles begins to dominate the photon emission process.This explains the presence of a substantial population at around 15 MeV.
For events with R post eγ < 0.6, they can be effectively identified as background with a PSD score (BDT/NN) close to -1/0 (indicating a background-like event).However, for the remaining events, a positive correlation is observed between R post eγ and the PSD score.

Limiting factors for the performance of discriminators
To explore the theoretical limits of separating the DSNB signals from atmospheric neutrino NC backgrounds, we systematically study the PSD capability using the photon emission time distributions obtained from different simulation stages to perform PSD.Detector and reconstruction effects can cause blurring of the reconstructed hit time, leading to a loss of differentiation between the two signal types.Therefore, we divide the complete simulation chain into five different simulation stages and provide detailed comments on these stages to better understand the simulation results: Stage 1: Events are simulated from the beginning at the event generator to the point of initial photon emission.The initial PET value (t s1 PET ) is obtained using detector simulation truth information.
Stage 2: This stage corresponds to the complete detector simulation.Optical processes of a photon from its primary position to a PMT are simulated, and the true PMT hit time is obtained.The associated PET value (t s2 PET ) is evaluated by subtracting the TOF from the true PMT hit time.The degradation of PSD performance in this stage demonstrates the impact of the calculated vs. the true TOF.deviation from the average.The lower panels display the time profiles around the peak using hits from dynode PMTs, following the same format as the top panels.The two rightmost panels illustrate the impact of dark noise on time profiles.
Fig. 8 shows the photon emission time distributions for the DSNB signals and atmospheric neutrino NC backgrounds at the five simulation stages.The curves in the panels represent the averaged time profile, progressing from left to right as Stage i for i = 1-5.During these stages, the introduction of detector, electronics, and reconstruction impacts causes changes in the average time shape.The shaded band along each curve represents a 1σ deviation from the averaged time profile.Shape fluctuations primarily arise from statistical uncertainty.The time range shown in the upper panels is -50 to 800 ns.The panels offer interesting insights.Firstly, the averaged PET profiles from Stage 2 exhibit a delay between around 50 and 300 ns compared to the true time profiles.Since the actual photon path is unknown, a straight line connecting the event vertex and the fired PMT is assumed, which can introduce a bias towards small values of TOF for scattered or reflected photons.Also, optical refraction on the vessel is neglected.Secondly, the average tail time spectra remain consistent throughout all simulation stages, but their shape fluctuation gradually increases.
The lower panels focus on the peak shape of the corresponding time profiles, specifically using hits from dynode PMTs.In the lower panels, the time range is zoomed in -20 to 20 ns.The time profiles from Stage 1 illustrate that the prompt recoil neutron of NC events travels a longer distance and takes more time to deposit its energy.This characteristic primarily impacts the distribution of starting times, thereby providing another discriminating factor in the peak time profiles between DSNB signal and NC background events.As the simulation progresses, the expected blurring of the distinguishing characteristics of this difference is observed in Stage 2 and Stage 3.However, during the reconstruction stage, a sudden sharpening of the peak shape is noticed.This is attributed to the generation of a higher number of photons at the initial stage within a short time frame.The overlapping of these multiple original time waveforms results in the waveform reconstruction detecting only the time of the first photon and assigning it the combined charge of the closeby p.e.'s.Subsequently, when the PET is obtained from the waveform reconstruction and weighted by charge, the peak becomes sharper.The larger impact on DSNB events actually enhances the differences between peaks to some extent, which is reflected in the performance of the PSD.
The BDT method is employed to assess the theoretical limitations of distinguishing the DSNB signals from atmospheric neutrino NC backgrounds.By using the time profiles obtained from these five stages (depicted in Fig. 8) as inputs for the BDT method, the performance of the PSD is evaluated.Note that for the BDT discriminator, n dark is fixed to 0 in Stages 1-4.Fig. 9 illustrates the PSD efficiency for DSNB signals as a function of the PSD inefficiency for atmospheric neutrino NC background within the energy range of 11-30 MeV at Stages i for i = 1-5.The results indicate that the inclusion of additional detector and electronic effects degrades the PSD performance.However, the PSD performance in Stage 4 is comparable to or slightly better than that in Stage 3. As mentioned earlier, the sharper peak in the time profiles of Stage 4 provides additional discrimination power.If an average NC background PSD inefficiency of 1% is required, the DSNB PSD efficiency is approximately 100% starting from Stage 1.However, it decreases to 95% due to the impact of calculated TOF (Stage 2), further reduces to 90% due to the uncertainty of electronics effects (Stage 3), and finally drops to 84% due to the pollution of dark noise (Stage 5).

Conclusion and outlook
In this paper, we presented the details and understanding of the PSD variables developed for suppressing atmospheric neutrino NC backgrounds, dominant for the DSNB search in JUNO.The particle-type dependent scintillation time profiles of liquid scintillators allow for an efficient particle identification.Using state-of-the-art parameters of the JUNO LS recipe [34], we developed two PSD discriminators based on machine learning methods.The BDT and NN methods give consistent results, with the NN method exhibiting superior PSD performance.The NC background is suppressed by two orders of magnitude, expanding JUNO's potential to detect DSNB signals significantly [13].We also found the atm-ν NC with 11 C events are more difficult to reject since these events feature a large γ energy fraction from neutron inelastic processes on carbon in the LS.Similar event topologies on oxygen have been observed as background in water Cherenkov detectors.We also studied how much the detector timing response and reconstruction effects degrade the PSD performance.
The powerful background rejection for atmospheric neutrino NC events provided by the PSD is crucial for the detection of the DSNB in LS experiments.It is important to validate the PSD variable performance with benchmark data in the future.In particular, the dominant NC background channel emits an energetic neutron with 11 C. Therefore, fast neutron sample tagged by the water pool veto detector could be used to constrain the PSD distributions.An AmBe neutron calibration source [48] emits neutrons of several MeVs and an accompanying 4.4 MeV γ.Events with similar R eγ as atm-ν NC w/ 11 C can be utilized to validate the PSD performance.PSD for atm-ν NC w/ 11 C events can be evaluated based on events with tagged 11 C decay.

Figure 1 :
Figure 1: The average reconstructed photon emission time for the DSNB signals (blue) and the atmospheric neutrino NC backgrounds (red).The peak structures of both the DSNB signal and the NC background events are zoomed in on the sub-panel.The dark noise is subtracted from the distributions for demonstration purposes.

Figure 2 :
Figure 2: The area normalized distributions of the input variables for the BDT method from the DSNB signals (blue curve) and atmospheric neutrino NC backgrounds (red curve) have been presented in the diagonal panels.In the lower panels, the scatter matrix of the correlations between the different input variables is shown, with the DSNB signals and atmospheric neutrino NC backgrounds denoted by blue and red dots, respectively.

Figure 3 :
Figure 3: The structure of the neural network.The input is a one-dimensional array encoded with charge and time information of PMT hits and the event vertex R 3 .The network is a fully connected layer composed of 100 nodes and a 0 represent constant term.The output is a signal-like score ranging from 0 to 1.

Figure 4 :
Figure 4: The left panel illustrates the BDT (filled area) and NN (solid lines) distributions.The right panel displays the areas under the ROC curve as a function of the prompt energy from both methods, with the shaded band representing the statistical error.

Figure 5 :
Figure5: The background inefficiency (left) and signal efficiency (right) as a function of the prompt energy with a PSD cut retaining 1% total background.The legend in the right plot is identical to the one in the left plot.In each prompt energy bin, the cut on NN discriminator is chosen to achieve the same background inefficiency as BDT method.The shaded band in the BDT and the error bar in the NN represent statistical errors.The BDT and NN methods give similar energy dependence and separation.Atm-ν NC w/ 11 C events with prompt energy below 18 MeV are more difficult to reject.

Figure 6 :
Figure 6: The R eγ distributions of the exclusive final state channels for atmospheric neutrino NC background events are depicted as a function of energy in the pre-simulation (left panels) and post-simulation (right panels) scenarios.The total kinetic energy (E tot k ) of all final state particles, excluding the neutrino, represents the energy before the detector simulation, whereas the visible energy of the prompt signal (E vis ) represents the energy after the simulation.The exclusive final state channels are categorized into atm-ν NC w/ 11 C in top panels and atm-ν NC w/o 11 C in bottom panels.

12 Figure 7 :
Figure 7: The cross sections of exclusive n− 12 C reactions, extracted from the results of TALYS, as a function of the incident neutron kinetic energy.

Figure 8 :
Figure 8: The photon emission time distributions for the DSNB signals (blue) and atmospheric neutrino NC backgrounds (red) have been illustrated at five different simulation stages.Each curve represents the averaged time profile, and the colored band along each curve indicates a 1σdeviation from the average.The lower panels display the time profiles around the peak using hits from dynode PMTs, following the same format as the top panels.The two rightmost panels illustrate the impact of dark noise on time profiles.

Figure 9 :
Figure 9: DSNB PSD efficiency as a function of atmospheric neutrino NC background PSD inefficiency using the BDT method in the energy range of 11-30 MeV at five different simulation stages.

Table 1 :
Rate [yr −1 kt −1 ] Fraction [%] Rate [yr −1 kt −1 ] Fraction [%]Summary of DSNB signal rates and NC interactions of atmospheric neutrino on12C in the JUNO sintillator.The NC interactions are listed distinguishing final states, only listing those with neutrons.

Table 2 :
A set of input variables for the BDT method.

Stage 3 :
Electronics simulation is performed to incorporate more realistic electronics effects, including TTS.The hit time obtained from this simulation is used to determine the associated PET value (t s3PET ), which is used to evaluate the impact of TTS on PSD performance.Stage 4 and Stage 5: Waveform and event reconstructions are conducted.Stage 4 doesn't include dark noise, while Stage 5 includes dark noise.These settings are used to study the effect of reconstruction and dark noise on particle discrimination.The corresponding reconstructed PET values, t s4 PET and t s5 PET , are obtained after the full simulation chain.