Atmospheric Neutrino Charged-Current Interactions at Large Liquid-Scintillator Detectors: I. Physics of Neutrino-Antineutrino Discrimination

Over the past several decades, the discovery of neutrino oscillations [1,2] has established the existence of neutrino masses and lepton flavor mixing. In the standard three-flavor neutrino framework, the neutrino flavor eigenstates $ \nu_\alpha $ ($ \alpha=e,\mu,\tau $) can be written as linear superpositions of the neutrino mass eigenstates $ \nu_i $ ($ i=1,2,3 $) via the lepton flavor mixing matrix [3−5], which is parametrized as

where $ c_{ij} \equiv \cos\theta_{ij} $ and $ s_{ij} \equiv \sin\theta_{ij} $ ($ ij=12, 13, 23 $), and δ is the Dirac CP-violating phase. The three neutrino oscillation frequencies are determined by the three mass-squared differences:

where $ (i,j = 1,2,3,\, i>j) $ and $ m_i $ is the mass of the i-th mass eigenstate $ \nu_i $.

In the standard three-neutrino oscillation framework, the mass-squared difference $ \Delta m_{21}^2 $ and the mixing angle $ \theta_{12} $ have been determined by the KamLAND reactor antineutrino experiment and several solar neutrino experiments [6]. Meanwhile, the parameters $ |\Delta m_{31}^2| $ and the mixing angle $ \theta_{23} $ are primarily measured in atmospheric and accelerator neutrino experiments, including the T2K and Super-Kamiokande experiments [7,8]. The third mixing angle, $ \theta_{13} $, has been precisely measured in a series of reactor antineutrino experiments, particularly the Daya Bay experiment [9], and later by long-baseline accelerator neutrino experiments. Although these parameters have been measured with high precision, two significant unknowns remain in the current standard three-neutrino oscillation model: the neutrino mass ordering (NMO), which addresses the distinct patterns of the neutrino mass spectrum, and leptonic CP violation. Understanding these unknowns is crucial for advancing our knowledge of particle physics and the evolution of the universe.

As a fundamental property of neutrinos, NMO is a critical parameter for many other important measurements. On the one hand, it influences the precision measurements of oscillation parameters and the determination of leptonic CP violation in neutrino oscillations. On the other hand, it provides constraints on the fundamental nature of neutrinos (i.e., whether they are Majorana or Dirac) and on their absolute masses. The experimental determination of NMO can be categorized into two methods. The first uses the Mikheev–Smirnov–Wolfenstein (MSW) matter effect [10,11], which is the foundation for long-baseline atmospheric and accelerator neutrino experiments, such as DUNE [12], Hyper-Kamiokande [13], PINGU [14], and ORCA [15]. Due to the MSW effect, neutrinos in the energy range of 3 to 10 GeV experience resonant flavor conversion in matter for normal mass ordering, while antineutrinos undergo the conversion for inverted mass ordering. This effect is especially pronounced in the electron flavor channel. In accelerator neutrino experiments, the differentiation between neutrinos and antineutrinos can be accomplished by adjusting the configuration of the accelerator beam. Conversely, in atmospheric neutrino experiments, the simultaneous presence of both neutrino and antineutrino types requires their clear separation, which is essential for accurately determining the mass ordering.

The Jiangmen Underground Neutrino Observatory (JUNO) [16,17], with a baseline of 52.5 km and a large liquid scintillator (LS) detector of 20 kton, adopts a unique approach by leveraging the interference of quasi-vacuum oscillations of reactor antineutrinos, setting it apart from other neutrino oscillation experiments. Meanwhile, JUNO is also capable of detecting atmospheric neutrinos, and the combination of reactor and atmospheric neutrino measurements is expected to enhance the sensitivity to the NMO. To achieve this goal, robust capabilities in neutrino directionality and type discrimination are essential. In Refs. [18,19], a preliminary attempt using machine learning techniques was developed. However, the underlying mechanisms and physical origins remain to be fully understood.

In this work, we investigate the characteristics of events from GeV atmospheric charged-current (CC) neutrino interactions in LS detectors and focus on the corresponding neutrino-antineutrino identification capabilities. Our study delves into the fundamental physics involved in identifying atmospheric neutrinos in large LS detectors through comprehensive Monte Carlo simulations. Due to high light yields and low energy thresholds, these detectors are particularly effective at capturing final-state hadronic components, which provide a distinctive statistical approach to differentiate between neutrinos and antineutrinos. Neutron tagging has proven to be particularly effective in LS detectors, in which the neutron multiplicity resulting from atmospheric neutrino interactions serves as a key factor for neutrino-antineutrino discrimination. The inelasticity of neutrino interactions, defined as the ratio of the energy carried by final-state hadrons to the total energy transfer, is another important parameter that contributes to this differentiation. Finally, we introduce a boosted decision tree (BDT) method to demonstrate the effectiveness of these two physical parameters in the neutrino-antineutrino discrimination.

The remainder of this paper is organized as follows. In Sec. 2, we will introduce the interaction tools used to generate final-state particles and perform detector simulations for this analysis. Sec. 3 will focus on the characteristics of both final-state leptons and hadrons in the context of detector response and physics. The effects of the detector size on the characteristics of final-state particles will be discussed in Sec. 4. Finally, the performance of neutrino and antineutrino discrimination using the BDT method will be presented in Sec. 5.

For this study, the simulation of neutrino CC interactions with ¹²C (~88%) and ¹H (~12%) in LS detectors is performed using the GENIE generator (version 3.2.0), in which the dominant contributions to the cross section arise from quasi-elastic scattering (QEL), nuclear resonance production (RES), and deep inelastic scattering (DIS) across different energy ranges. The secondary interactions of final-state particles in LS are also simulated with Geant4 (version 10.4.2). In this calculation, four neutrino types ($ \nu_{e} $, $ \overline{\nu}_{e} $, $ \nu_{\mu} $, $ \overline{\nu}_{\mu} $) are generated, with one million events per type and isotropic directions, and provided to GENIE. The tune G18_10b_02_11b is employed in this study, and a detailed description of the tuned model in GENIE V3 can be found in Ref. [20]. The neutrino fluxes used as input to GENIE are set to be uniform in order to eliminate energy dependence. The final-state particles generated by GENIE are injected into a standalone Geant4 framework to simulate their secondary interactions in LS. For each category of neutrino interaction, event vertices are sampled uniformly within the LS volume. A customized Livermore physics list, which includes the positronium process, is used to simulate electromagnetic interactions [21]. Additionally, a QGSP_BERT_HP-based physics list with a modified neutron-capture process is employed to simulate the secondary hadronic interactions [22]. As an illustrative framework, a simplified JUNO-like detector has been constructed in Geant4 to capture the primary features of the LS medium and basic detector geometry [16,23]. This model includes 20 kton of LS contained within an acrylic sphere with a radius of 17.7 m. A total of 17,612 20-inch photomultiplier tubes (PMTs) are uniformly installed on the spherical surface of a stainless-steel structure with a radius of 20 m. The space between the acrylic sphere and the PMTs is filled with purified water. No other detector components are included in this simplified framework. In the simulation, the geometries of the PMTs are modeled explicitly, and the number of incident optical photons on each PMT is recorded. However, detailed simulations of the front-end electronics, as well as the subsequent waveform and event reconstruction processes, are not considered in this work; accordingly, no energy threshold is applied. A schematic diagram of the simulation workflow for atmospheric neutrino interactions in LS is shown in Figure 1, including the atmospheric neutrino interactions, the production of charged leptons and final-state hadrons, subsequent secondary interactions, and final neutron captures on hydrogen.

Figure 1.  (color online) A schematic diagram of atmospheric neutrino interactions in LS is presented, including charged-current interactions, sequential secondary interactions, and final neutron captures.

In Figure 2, we illustrate two examples of charged-current interactions by atmospheric neutrinos, each with an energy of approximately 10 GeV. The solid lines depict the tracks of the final-state particles and their secondary products, while the dashed lines represent the incoming atmospheric neutrinos. Triangles mark the locations of neutron captures. For electron neutrinos, the final-state electron initiates an electromagnetic shower that produces secondary electrons and gamma rays. A hadronic shower from the $ \nu_e $ interaction is produced concurrently, contributing additional energy deposition and potentially yielding neutrons and other hadrons. For muon neutrinos, the final-state muon loses energy primarily through ionization, leaving a clear track within the detector medium. In this case, the muon track is typically much easier to distinguish from the hadronic shower.

Figure 2.  (color online) Typical examples of charged-current interactions of $ \nu_e $ (left panel) and $ \nu_\mu $ (right panel) with ¹²C. The solid lines show the tracks of final-state particles and their secondary particles, and the dashed line indicates the incoming atmospheric neutrino. The triangles indicate the positions of neutron captures.

In this section, we present detailed characteristics of atmospheric-neutrino interaction events using simulated data samples. For charged leptons, electrons and muons exhibit distinct signatures, namely electromagnetic showers and tracks, respectively. The resulting topological, temporal, and spatial patterns of photoelectrons detected by the PMTs can be used for flavor identification. Meanwhile, the detection of hadronic components primarily enables discrimination between neutrinos and antineutrinos. Inelasticity and neutron multiplicity are two essential parameters for this purpose and are discussed below. These parameters provide valuable insights into the dynamics of neutrino interactions and secondary interactions occurring within the LS, and thus play a crucial role in distinguishing between neutrinos and antineutrinos.

When an electron interacts in liquid scintillator (LS), it initiates an electromagnetic shower—a cascade of secondary particles, including gamma rays and additional electrons. This process produces copious scintillation light over an extended volume, resulting in a spatially diffuse light distribution. By contrast, a traversing muon primarily ionizes the medium, depositing energy in a more localized manner along a straight, narrow track, yielding a more concentrated scintillation signal. Hadronic showers, on the other hand, are characterized by complex, clustered patterns arising from interactions among multiple secondary particles in the LS. Consequently, PMTs record spatially extended, multi-site patterns for both electron-induced and hadronic showers, whereas muons produce a well-defined linear track. The distinct spatiotemporal profiles of these light patterns, as reflected in PMT hit-time distributions, provide a robust basis for particle identification and neutrino-interaction discrimination in the LS medium. In addition to these primary light patterns, delayed signatures of neutrino interactions—such as neutron capture and delayed beta decays of nuclei—offer further information for identifying neutrino interactions.

Figure 3 illustrates the distinct signatures of hit-time distributions for electrons, muons, protons, and pions in different angular regions, where α denotes the angle between the incoming lepton momentum direction and the PMT position. Colored lines correspond to electrons, muons, protons, and pions with momenta of 3 GeV (upper panels) and 1 GeV (lower panels), respectively. For each panel, the plots are arranged with increasing α values from upper left to lower right. Note that all particles originate at the detector centre, and the distributions are averaged and normalized over PMTs within each selected α region. From the figure, we observe that, at multi-GeV energies, the hit-time profiles of muons and electrons differ markedly. Higher momentum shifts the first hits to earlier times—a trend more pronounced for muons than for electrons and stronger at smaller α values. Moreover, 3 GeV muons exhibit more late-time hits than electrons of the same momentum at larger α, an effect that diminishes as the momentum decreases. Regarding hadrons, protons and $ \pi^{+} $ mesons with higher momentum (3 GeV) induce secondary interaction processes, such as inelastic scattering or decays, which generate secondary particles that contribute to hadronic showers. As a result, PMTs near these secondary particles may register signals much earlier, as illustrated in Figure 3. Unlike the charged leptons and $ \pi^{+} $ mesons, protons produce fewer photons due to stronger quenching. Additionally, at lower momenta, protons can be approximated as single, point-like sources of light.

Figure 3.  (color online) PMT hit-time distributions for different angular regions (α), where α is the angle between the PMT position vector and the charged-particle momentum. The lines correspond to electrons, muons, protons, and pions with momenta of 3 GeV (a) and 1 GeV (b). In the GEANT4 simulation, secondary interactions and their resulting products are included.

To quantify the differences between the aforementioned leptonic and hadronic signatures, we introduce three PMT-level observables to capture their distinct hit patterns [19]:

● FHT: the hit time of the first photoelectron detected by the PMT.

● NPE: total number of photoelectrons collected by the PMT.

● Slope: the slope of the rising edge of the PMT hit distribution, defined over the interval from the FHT to the peak time.

Notably, in LS detectors that use PMTs as detection units, custom flash analog-to-digital converters are widely employed in the readout electronics [24]. This approach enables precise acquisition of all three PMT-level observables under realistic detection conditions.

For illustration, these PMT observables are shown quantitatively in Figure 4 and Figure 5 using the Mollweide projection. This equal-area pseudo-cylindrical projection preserves relative areas while mapping the spherical surface onto a plane. The left, middle, and right panels correspond to the Slope, FHT, and NPE variables, respectively. The injected particle momenta in Figure 4 and Figure 5 are 3 GeV and 1 GeV, respectively, with momentum directions perpendicular to the plane and pointing outward.

Figure 4.  (color online) Mollweide projections of PMT-level observables for electrons, muons, protons, and pions with a momentum of 3 GeV. Columns from left to right show Slope, FHT, and NPE. All particles originate perpendicular to the plane and point outward. In the GEANT4 simulation, secondary interactions and their products are enabled.

Figure 5.  (color online) Mollweide projections of PMT-level observables for electrons, muons, protons, and pions with a momentum of 1 GeV. From left to right, the columns show Slope, FHT, and NPE. All particles are emitted normal to the plane and directed outward. In the GEANT4 simulation, secondary interactions and their resultant products are included.

Compared to electrons, a 3 GeV muon can travel roughly 15 m in LS, comparable to the ~20 m radius of JUNO's stainless-steel structure that houses the inward-facing PMTs. When the muon stops near the PMT array, PMTs aligned with its trajectory receive a bright, prompt light signal. This yields a distinct slope pattern (left column of Figure 4) and a rapid fall-off of NPE with increasing α relative to the muon direction. These features enable excellent discrimination between electron and muon neutrinos. At 1 GeV, however, the muon track is shorter, and the detector response behaves more like a point-like source. The hit-time distributions of 1 GeV muons resemble those of electrons, resulting in less distinctive PMT-level signatures (Figure 5). Consequently, as shown in Ref. [19], flavor-identification performance degrades at lower neutrino energies.

The Mollweide maps for hadronic events are characterized by distinct patterns from multiple hadronic showers, which exhibit clustered structures that differ from those produced by electrons and muons. These clusters arise from ionization-induced scintillation in the LS caused by multiple secondary particles. When the momentum decreases to 1 GeV, protons deposit energy almost exclusively through ionization. Owing to their larger mass and greater ionization energy loss, they leave very short tracks in the LS, effectively producing a point-like light source near the detector center. Consequently, Figure 5 shows relatively small variations in the number of photons received by PMTs at different positions. In contrast, for $ \pi^{+} $ particles, in addition to ionization from the pion itself, charged leptons produced in their decay also contribute to ionization-induced scintillation. Therefore, even at lower energies, PMT hit patterns generated by pions still exhibit multiple clustered structures.

From the discussion above, it is evident that different charged leptons and hadrons elicit distinct responses in the LS detector. This variability enables discrimination between charged-current and neutral-current interactions, as well as flavor identification. These PMT-level observables constitute a two-dimensional data array for each physics event, which can be exploited by deep-learning models employing multilayer neural networks to extract features effectively and reduce data dimensionality. In Ref. [19], it has been demonstrated that excellent discrimination performance can be achieved using point-cloud-based machine-learning methods under conditions that closely resemble actual detector operations. This makes the separation of neutrinos and antineutrinos feasible in realistic experimental observations.

Having established the fundamental patterns of electromagnetic showers, tracks, and hadronic showers associated with electrons, muons, and hadronic activity, we are now prepared to discuss the parameters that can be used to discriminate between neutrinos and antineutrinos. In general, neutrino–nucleus CC interactions can be described in a semi-inclusive manner as:

where N denotes the initial nucleus, M the daughter nucleus, $ l^{\pm} $ the outgoing charged lepton, and X the hadronic final state. The inelasticity y for such interactions is defined as

where $ E_{\nu} $ and $ E_{l} $ denote the energies of the incoming neutrino and outgoing lepton, respectively.

The CC interactions of neutrinos and antineutrinos exhibit distinct inelasticity distributions (Figure 6), due to differences in their chiral couplings and in the parton distributions within nucleons. Neutrinos tend to transfer a larger fraction of their energy to the hadronic system than antineutrinos, resulting in greater hadronic energy deposition. At low energies, quasi-elastic (QEL) scattering dominates and is characterized by small y values; in this regime the neutrino interacts with a single nucleon, with most of its energy carried by the outgoing lepton.

Figure 6.  (color online) The inelasticity distributions for neutrino (upper panel) and antineutrino (lower panel) interactions are shown for multiple energy ranges. For comparison, the separate contributions from QEL, RES, and DIS are also shown. Each distribution is normalized to the number of events in the corresponding energy range.

As the neutrino energy increases, deep inelastic scattering (DIS) becomes prevalent, with neutrinos scattering off quarks within nucleons, yielding larger y values and a more complex hadronic final state. The figure also shows the separate contributions from QEL, resonance (RES), and DIS processes, each normalized per energy bin. A clear asymmetry is observed between neutrinos and antineutrinos: antineutrinos systematically transfer less energy to hadrons. This difference provides a key handle for distinguishing the two interaction types.

Recent studies have demonstrated the reconstruction of the inelasticity parameter in water- and ice-Cherenkov detectors (Ref. [25]). This capability enhances the separation between neutrinos and antineutrinos and improves sensitivity to the NMO (Refs. [26,27]). In this work, we focus on LS detectors, which, owing to their lower energy thresholds, offer distinct advantages for detecting the hadronic components of neutrino interactions. Consequently, we anticipate improved reconstruction of the inelasticity parameter in LS-based detection; detailed results will be reported in a separate publication.

The primary interactions of atmospheric neutrinos in LS produce hadronic final states, which mainly consist of neutrons, protons, and $ \pi^{\pm} $ mesons. Figure 7 illustrates the average multiplicity of these hadrons as a function of the neutrino energy. Antineutrino interactions produce more neutrons and $ \pi^{-} $ mesons than neutrino interactions, whereas neutrino interactions yield more protons and $ \pi^{+} $ mesons. Additionally, even within the same hadronic isospin multiplet—such as protons and neutrons—an asymmetry is observed between neutrinos and antineutrinos. At the same incident energy, neutrino interactions produce more protons than antineutrino interactions produce neutrons. As discussed above, neutrino interactions transfer more energy to the final-state hadronic system than antineutrino interactions, leading to generally higher hadron multiplicities for neutrinos.

Figure 7.  (color online) The averaged multiplicity of final-state hadrons $ \langle N_{h}\rangle $ is shown as a function of the incoming neutrino energy for neutrinos (left) and antineutrinos (right). The solid lines with dots represent the mean multiplicity in each energy interval.

Figure 8 presents the average kinetic energies of neutrons, protons, and $ \pi^{\pm} $ mesons as a function of neutrino energy. Owing to their lower mass, $ \pi^{\pm} $ mesons can acquire greater kinetic energy from interactions, while neutrons and protons exhibit relatively lower average kinetic energies. Furthermore, in charged-current interactions, neutrinos transfer more energy to the hadronic system than antineutrinos. As a result, the average kinetic energies of all four hadron species are significantly higher in neutrino-induced events than in antineutrino-induced events.

Figure 8.  (color online) The average kinetic energy of the four primary final-state hadrons, $ \langle E_k\rangle $, is shown as a function of neutrino energy for neutrinos (left) and antineutrinos (right). The solid lines with dots represent the mean kinetic energy of the hadrons in each energy interval.

Neutron capture is commonly used to detect low-energy antineutrinos via inverse beta decay on free protons. However, the situation becomes considerably more complex for higher-energy neutrinos interacting with bound nuclei. The primary interaction occurs when a neutrino scatters on a nucleus, initiating final-state interactions within the nucleus that can produce neutrons. In this regime, antineutrinos are generally more effective than neutrinos at generating neutrons. In addition, hadrons produced in the primary interaction can undergo secondary interactions in the LS, such as inelastic scattering and $ \pi^{-} $ capture, which further increase neutron production. Finally, the resulting neutrons are captured by nuclei, predominantly on protons in LS, thereby emitting characteristic gamma-ray signals. Understanding this interplay among neutrino interactions, neutron production in LS, and subsequent neutron captures is essential for characterizing the captured neutron multiplicity distribution.

To elucidate the physics underlying the neutron capture multiplicity distributions, we present a detailed analysis of the relevant production and interaction processes. First, the number of neutrons produced in hadronic secondary interactions is strongly influenced by the kinetic energy of the outgoing hadrons, primarily protons, neutrons, and $ \pi^{\pm} $. Figure 9 shows the mean captured neutron multiplicity as a function of kinetic energy for four final-state hadrons. As the kinetic energy increases, the neutron production capability of these hadrons grows markedly. Notably, both neutrons and $ \pi^{-} $ exhibit particularly high neutron yields. The $ \pi^{-} $ is of special interest because it can not only induce secondary interactions that produce additional neutrons, but also be captured by nuclei, further enhancing neutron generation.

Figure 9.  (color online) The averaged captured-neutron multiplicity $ \langle N_{{\rm{n}}-{\rm{captured}}}\rangle $ is shown as a function of kinetic energy for the four primary final-state hadrons. The lines with dots represent the mean value of the captured-neutron multiplicity in each energy interval.

To further elucidate the production mechanisms of captured neutrons, Figure 10 presents the relative contributions of different neutron production processes for each of the four final-state hadrons as a function of kinetic energy. Neutron inelastic scattering is a dominant source of neutron production for all final-state hadrons. In addition to neutron inelastic scattering, inelastic interactions induced by protons and charged pions also contribute substantially to neutron yields. Subdominant contributions arise from photonuclear interactions and decays of unstable hadrons and nuclei. Neutrons produced directly in primary neutrino interactions account for on the order of tens of percent of neutron captures at kinetic energies below 100 MeV, but their contribution declines rapidly with increasing energy.

Figure 10.  (color online) This figure shows the fractions of different neutron-production processes for four final-state hadron species as functions of their kinetic energies. Each fraction is defined as the number of captured neutrons produced by a specific process divided by the total number of captured neutrons. These neutrons are counted separately for protons, neutrons, $ \pi^{+} $, and $ \pi^{-} $ produced in primary neutrino interactions. The x-axis represents the kinetic energy of these four final-state hadron species.

As previously noted, neutrinos and antineutrinos exhibit distinct inelasticity distributions, reflecting differences in the energy transferred to the hadronic system. Because the captured-neutron multiplicity is closely tied to the hadronic kinetic energy, inelasticity and neutron multiplicity are correlated observables that help distinguish neutrinos from antineutrinos. This relationship is illustrated in Figure 11, which shows two-dimensional distributions of inelasticity versus captured-neutron multiplicity for neutrinos (upper panels) and antineutrinos (lower panels), for energies below 5 GeV (left) and above 5 GeV (right). A clear positive correlation emerges above 5 GeV. Compared with neutrinos, antineutrino interactions generally transfer less energy to the hadronic system but produce more neutrons through primary processes. Consequently, at lower energies antineutrino interactions yield more captured neutrons, whereas at higher energies neutrino interactions dominate the neutron yield. This energy-dependent behavior leads to clearly separable two-dimensional patterns in inelasticity and neutron multiplicity, providing a key handle for neutrino–antineutrino discrimination. The same feature is further illustrated in Figure 12, which shows the average captured-neutron multiplicity as a function of energy for neutrinos (left) and antineutrinos (right), with individual contributions from different final-state hadrons indicated for comparison. At low energies, inelastic processes that produce neutrons dominate, with primary neutrons contributing more strongly in antineutrino interactions. At higher energies, charged-pion-induced processes become the most significant source of neutron production.

Figure 11.  (color online) Two-dimensional distributions of inelasticity versus captured-neutron multiplicity for neutrinos (upper panels) and antineutrinos (lower panels) are shown. The left and right panels correspond to (anti)neutrino energies below and above 5 GeV, respectively. For each energy range, flat (anti)neutrino energy spectra are used.

Figure 12.  (color online) The average captured-neutron multiplicity $ \langle N_{{\rm{n}}-{\rm{captured}}}\rangle $ is shown as a function of neutrino (left panel) and antineutrino (right panel) energy; the individual contributions from different final-state hadrons are also shown for comparison.

Based on the preceding discussions, the characteristics of charged leptons, inelasticity, and neutron-capture multiplicity have thus far been examined under idealized physical conditions. In practice, however, realistic detector configurations can significantly alter these distributions and affect discrimination between neutrino types. Taking the JUNO detector as an example, which uses a 17.7 m-radius acrylic sphere containing 20 kton of LS as its target medium, we concentrate on fully-contained (FC) events, defined as interactions that deposit all their energy within the target. This selection enables the study of the relevant observables in a well-controlled environment. In contrast, partially-contained (PC) events are those that begin within the central detector (CD) but deposit a portion of their energy outside the CD.

For electron neutrino interactions, the spatial extent of the final-state electron is generally small because the electromagnetic shower deposits energy rapidly. As a result, most of these events are FC, and their reconstructed properties are only weakly affected by the detector size. The situation is notably different for muon neutrino interactions. As outlined in Section 3.1, the final-state muon usually produces an extended track inside the detector. Hence, FC muon neutrino events are strongly influenced by geometric constraints, which can modify key observables, including energy reconstruction, track reconstruction, and flavor discrimination.

The finite detector size distorts the measured muon kinetic-energy spectra because of the physical cutoff on observable track lengths, especially at high energies. This geometric effect is illustrated in Figure 13, which compares spectra from neutrino and antineutrino interactions under different containment conditions. The fraction of events that are fully contained within the active volume serves as a direct measure of geometric acceptance and decreases with muon energy. Consequently, although antineutrinos intrinsically produce higher-energy muons because they transfer less energy to hadrons, this difference is substantially reduced in the spectrum of FC events. At high incident energies, the geometric cutoff strongly suppresses the high-energy tail of the antineutrino spectrum, causing the observed distributions for the two flavors to converge.

Figure 13.  (color online) The influence of detector size on the kinetic-energy distributions of muon-neutrino samples. Dashed lines represent the combined distributions of all FC and PC events; solid lines represent FC events only. Blue and green lines denote neutrino and antineutrino samples, respectively. The red solid lines indicate the FC fraction as a function of kinetic energy and are identical for both $ \nu_{\mu} $ and $ \overline{\nu}_{\mu} $ events. The left and right panels correspond to different ranges of incident neutrino energy, simulated using flat neutrino-energy spectra.

Building on this geometric distortion, the inelasticity and captured-neutron multiplicity in FC $ \nu_\mu $/$ \overline{\nu}_\mu $ events are also reshaped by the detector boundary, as shown in Figure 14. At high neutrino energies, the detectable muon energy is limited by the boundary, forcing a larger share of the energy to be assigned to the hadronic system in the reconstruction. This artificially inflates the measured inelasticity and neutron multiplicity while simultaneously suppressing genuine flavor differences that depend on the detailed energy partition, thereby making $ \nu_\mu $/$ \overline{\nu}_\mu $ separation difficult for high-energy FC events. In Figure 14, in addition to the mean values depicted as solid and dashed lines with dots, the bands show the standard deviations of the distributions within each energy interval, representing an effective resolution arising from both the primary neutrino interaction and secondary interactions. This effect is inherently tied to the detector’s physical scale. Given JUNO’s liquid-scintillator diameter of 35.4 m, muons with kinetic energies above ~7 GeV typically escape the active volume. Hence, the majority of events beyond this threshold are necessarily partially contained. Conversely, in the few-GeV energy range most relevant for the NMO determination, the detector-size effect turns from a limitation into an advantage. In this regime, muons are fully contained, and the hadronic system receives sufficient energy to amplify flavor-sensitive signatures in inelasticity and neutron multiplicity. This provides a clear window for effective $ \nu_\mu $/$ \overline{\nu}_\mu $ discrimination, precisely where it matters most for NMO sensitivity.

Figure 14.  (color online) The inelasticity (left panel) and captured-neutron multiplicity (right panel) distributions are shown as functions of the muon-neutrino energy. The dashed lines with dots represent the mean values for all FC and PC events, while the solid lines with dots indicate the mean values for FC events only. The shaded bands indicate the standard deviations of the distributions within each energy interval.

In this section, we quantitatively assess the performance of neutrino-antineutrino discrimination using the event-level observables discussed previously. As an illustration, we employ a simple BDT that takes inelasticity and the multiplicity of captured neutrons as discriminants. More advanced machine-learning methods that exploit PMT-level information have been explored elsewhere for flavor identification (see Ref. [19]) and are not considered here.

We trained BDT classifiers separately for different flavors to evaluate their power to discriminate between neutrinos and antineutrinos in the energy range 0.1–20 GeV, using the area under the receiver operating characteristic (ROC) curve (AUC) as the performance metric. In Figure 15, we show the BDT-score distributions for $ \nu_{e} $/$ \overline{\nu}_{e} $ (left) and $ \nu_{\mu} $/$ \overline{\nu}_{\mu} $ (right) discrimination for the total sample (top) and the FC sample (bottom). The higher the BDT score, the more an event is classified as neutrino rather than antineutrino. From the figure, one finds that the antineutrino samples (blue curves) have smaller BDT scores than the neutrino samples. In addition, the neutrino distributions (red curves) exhibit two peaks at ~0.6 and ~0.9. To understand these features, Figure 16 shows the two-dimensional distributions of inelasticity and BDT score for events with and without a true neutron capture. As illustrated, whether the number of captured neutrons is zero is a crucial discriminator between neutrino and antineutrino events. Samples without a neutron capture show the highest BDT scores, peaking at 0.9. For samples with a captured neutron, higher inelasticity correlates with a higher BDT score. However, this relationship saturates when the inelasticity exceeds 0.8, resulting in a plateau where the BDT score remains unchanged, with a peak at 0.6. This behavior is consistent with the inelasticity curves in Figure 6, where the decreasing inelasticity dependence of antineutrino DIS channels becomes relatively flat.

Figure 15.  (color online) BDT score distributions for $ \nu_{e} $/$ \overline{\nu}_{e} $ (left) and $ \nu_{\mu} $/$ \overline{\nu}_{\mu} $ (right) discrimination, shown for the total sample (top) and the FC sample (bottom).

Figure 16.  (color online) The distributions of inelasticity and BDT score, with and without neutron capture, are shown. The left and right panels correspond to $ \nu_{e} $/$ \overline{\nu}_{e} $ and $ \nu_{\mu} $/$ \overline{\nu}_{\mu} $ discrimination, respectively, while the upper and lower panels correspond to the total and FC samples.

In Figure 17, we show the ROC curve (left panel) and AUC values (right panel) for neutrino–antineutrino discrimination. Figure 17(a) presents the ROC curves for $ \nu_{e} $/$ \overline{\nu}_{e} $ and $ \nu_{\mu} $/$ \overline{\nu}_{\mu} $ discrimination for the total and FC samples. Figure 17(b) shows the AUC as a function of neutrino energy, comparing results for both electron and muon flavors under different geometric containment conditions. The combined use of inelasticity and neutron multiplicity provides effective discrimination, which can exceed 80% at lower energies when neutron multiplicity contributes significant additional separating power. As energy increases, the discrimination power decreases; however, the AUC remains around 75% in the energy range of approximately 3 to 10 GeV (the NMO-sensitive region). For muon neutrinos in the FC sample, performance degrades at higher energies. This decline is a direct geometric effect: as the neutrino energy increases, the muon track length may exceed the detector’s active volume. This containment limit reduces the measurable energy of the final-state lepton and diminishes neutrino-type-sensitive information from the hadronic system. Therefore, a detector with a larger active volume not only improves event statistics but also better preserves the intrinsic separation power, particularly for high-energy muon neutrinos, by allowing longer tracks to remain fully contained.

Figure 17.  (color online) The ROC curve (left panel) and AUC values (right panel) for neutrino-antineutrino discrimination are presented. Discrimination between $ \nu_{e}/\overline{\nu}_{e} $ and $ \nu_{\mu}/\overline{\nu}_{\mu} $ for the total and FC samples is shown separately, using the true values of inelasticity and neutron multiplicity obtained from simulations.. The entire energy range has been integrated for the left panel while displayed for the right panel.

It is important to note that these results are obtained under idealized assumptions, without accounting for reconstruction resolution in inelasticity or the tagging efficiency for captured neutrons. In LS detectors, the low-energy threshold significantly enhances the detection capability for hadronic components, including signals from neutron capture. This capability not only improves neutrino–antineutrino discrimination but also enables deeper insights into the underlying physics of atmospheric neutrino interactions in the LS medium.

Before concluding this work, we highlight the neutrino–antineutrino discrimination potential across different detector technologies. The Super-Kamiokande experiment has been upgraded with gadolinium loading since 2020 to improve neutron-tagging efficiency, enhancing its sensitivity to the diffuse supernova neutrino background (DSNB) [28] and benefiting atmospheric neutrino oscillation studies [8]. Recent studies have also explored the use of final-state hadron detection and machine learning to reconstruct inelasticity in water/ice Cherenkov detectors such as IceCube-Upgrade and ORCA [25,27]. The upcoming DUNE experiment will employ liquid-argon time-projection chamber technology [29], offering excellent tracking and calorimetric precision for both charged leptons and hadrons from neutrino–nucleus charged-current interactions [30]. However, neutron tagging remains challenging in such detectors, which may affect energy resolution and limit the potential for distinguishing neutrinos from antineutrinos.

In this work, we have systematically evaluated the characteristics of charged-current interactions of atmospheric neutrinos with ¹²C and ¹H nuclei in liquid-scintillator detectors, covering both the charged-lepton and hadronic components of the final state. The principal features of these event signatures are summarized as follows.

● Charged leptons and hadrons exhibit distinct light topologies in LS: muons produce clear, track-like ionization patterns; electrons generate diffuse, multi-point electromagnetic showers; and hadrons yield multi-component hadronic showers. These differences provide a strong basis for flavor and neutrino-antineutrino identification at the PMT level.

● Inelasticity, defined as the fraction of energy transferred to the hadronic system, differs markedly between neutrinos and antineutrinos. Neutrino interactions yield an approximately flat inelasticity distribution, whereas antineutrinos preferentially populate events with smaller hadronic energy fractions, resulting in a distribution that falls rapidly with increasing inelasticity.

● Captured neutron multiplicity is dominated by primary neutron production at low neutrino energies; at higher energies, however, secondary interactions of final-state hadrons in the LS medium predominate. Below approximately 5 GeV, antineutrino interactions yield higher neutron multiplicities; above this energy, the trend reverses, with neutrino interactions producing more neutrons.

● Detector size significantly affects muon-flavor events, owing to their long tracks, whereas electron-flavor events remain largely unaffected. For contained muon events, geometric containment enhances both inelasticity and neutron multiplicity, while simultaneously suppressing the intrinsic differences between neutrinos and antineutrinos.

Using the two-dimensional distribution of inelasticity and neutron multiplicity in a BDT classifier, we demonstrate discrimination between neutrinos and antineutrinos with AUC values exceeding 70% under an idealized detection scenario. Performance is better at lower energies, where primary neutron production provides additional discriminating information. For FC muon events, the AUC drops toward 50% around 20 GeV, a direct consequence of geometric containment, which limits the measurable lepton energy and suppresses flavor-sensitive hadronic information. Further studies employing more advanced machine-learning methods that incorporate additional PMT-level features are anticipated to yield improved performance.

The present analysis also enables improved energy and direction reconstruction for atmospheric neutrino charged-current events in LS detectors [18]. Collectively, these results establish a foundation for advanced high-energy neutrino reconstruction and flavor discrimination, enabling more precise studies of neutrino oscillations, particularly measurements of the neutrino mass ordering and CP violation with atmospheric neutrinos. This work establishes the physical basis and capability for neutrino/antineutrino discrimination in large LS detectors using truth-level quantities. However, the effects of detector response on inelasticity and neutron tagging cannot be neglected. For instance, after-pulse effects [31] that follow a GeV-scale event can affect the tagging of captured neutrons produced by atmospheric neutrinos. Additionally, finite inelasticity resolution may further confuse neutrinos and antineutrinos. These effects necessitate further investigation using more accurate and comprehensive simulations, together with realistic data.