Abstract
HTML
Reference
Related
PDF
-
The discovery of neutrino oscillations has established that at least two neutrino mass eigenstates are massive. This unambiguous conclusion has provided the first evidence of new physics beyond the Standard Model which assumes neutrinos are massless. In the framework of three neutrino oscillations, six fundamental parameters are involved to describe the oscillation phenomena. The solar and reactor neutrino experiments have determined the mass-squared difference
$ \Delta m^{2}_{21} $ and the mixing angle$ \sin^2 \theta_{12} $ . Atmospheric and accelerator experiments have measured$ |\Delta m^{2}_{32}| $ (or$ |\Delta m^{2}_{31}| $ ) and$ \sin^2 \theta_{23} $ . Reactor antineutrino experiments have provided the most precise determination of$ \sin^2 \theta_{13} $ , and have also achieved a comparable precision of$ |\Delta m^{2}_{32}| $ with the measurements from the atmospheric and accelerator experiments. The two independent mass-squared differences and the three mixing angles have been measured with a precision of a few percent [1]. However, the sign of$ \Delta m^{2}_{32} $ (referred to as neutrino mass ordering, NMO) and the value of the CP violating phase$ \delta_{\rm{CP}} $ remain partially unknown [2−4], thus, several experiments [5−8] have been proposed to determine these two crucial parameters.The Jiangmen Underground Neutrino Observatory (JUNO) [9−11] in southern China is designed to determine the NMO by detecting reactor antineutrinos at baselines of 52.5 km from the Taishan and Yangjiang nuclear power plants. The detector is equipped with a 20 kton spherical volume of liquid scintillator (LS) located deep underground, providing an overburden of 650 m. JUNO aims to resolve the NMO by probing the interference effect of the two fast oscillations induced by
$ |\Delta m^{2}_{31}| $ and$ |\Delta m^{2}_{32}| $ of reactor antineutrinos. Sufficient energy resolution is crucial for accurately measuring the oscillation pattern in the antineutrino spectrum, and it directly determines the sensitivity to the NMO. In this context, the energy resolution refers to the sigma-to-mean ratio of the Gaussian function used to fit the energy distribution of positrons at a fixed energy. The designed energy resolution of JUNO is 3% at 1 MeV. Taking this assumption as a benchmark, JUNO reported an NMO sensitivity in the design study at a confidence level of 3-4σ after 6 years of data taking [10]. In addition, JUNO is expected to enable sub-percent level precision measurements of the oscillation parameters:$ \sin^2 \theta_{12} $ ,$ \Delta m^{2}_{21} $ , and$ \Delta m^{2}_{31} $ [12]. The experiment also serves as an observatory to detect neutrinos originating from both celestial and terrestrial sources [13−16]. Furthermore, the experiment is sensitive to search for rare decays that go beyond the predictions of the Standard Model [17]. Currently, JUNO is in the construction phase and is expected to start data-taking in the near future.Over the past few years several factors related to the NMO sensitivity in JUNO have undergone changes. As a result, the JUNO collaboration recently updated the NMO sensitivity using more realistic inputs [18]. These updates include the incorporation of updated reactor thermal powers and new constraints on the reactor neutrino spectrum we expect to obtain from the TAO [19] satellite detector. Additionally, the accidental backgrounds have been revised based on radioassay results of the detector's raw materials and the estimation of radioactive impurities during the detector installation [20]. The cosmogenic backgrounds have been re-evaluated due to changes in the overburden, and new backgrounds from global reactors and atmospheric neutrinos have been included. Furthermore, the detector energy resolution has been updated through the full detector simulation and event reconstruction, which take into account a better understanding of the detector structure, as well as the performance of the liquid scintillator (LS) and photomultiplier tubes (PMT) [21, 22]. This paper serves as a complementary report to the updated NMO sensitivity paper of JUNO and provides a detailed description of the latest understanding made to the energy resolution prediction in the JUNO experiment, which is essential in the NMO analysis.
In previous LS-based neutrino detectors as the Borexino [23] and KamLAND [24] experiments achieved the energy resolution levels of 5% and 6.5% at 1 MeV, respectively. In these experiments, statistical fluctuations dominate the energy resolution, with the reported photoelectron (PE) yields of 511 PE/MeV and 250 PE/MeV, respectively. The JUNO detector was designed with a focus on achieving high LS light yield, high transparency, and high photon collection efficiency to minimize the statistical fluctuation of the detected PE. The LS recipe has been optimized to maximize the light yield using one of the antineutrino detectors in the Daya Bay experiment [25]. Additionally, 17,612 high detection efficiency 20-inch PMTs (LPMTs) [21] and 25,600 3-inch PMTs (SPMTs) [26] have been deployed, resulting in a total photocathode coverage of 78%. However, due to the unprecedented energy resolution requirement of JUNO, several other factors, which are usually negligible in previous experiments, significantly contribute to the energy resolution and must be carefully evaluated. These factors include the quenching effect in LS, Cherenkov radiation, electronics and PMT responses, among others. All of these effects have been modeled in the full detector simulation and will be thoroughly discussed in this paper, along with the impacts from the vertex and energy reconstruction algorithms.
The remaining parts of this article are organized as follows. We first provide an introduction to the JUNO detector in Sec. 2, followed by a summary of the main factors that impact the energy resolution throughout the processes of light production, propagation, and detection in Sec. 3. The full detector simulation is then presented in Sec. 4, focusing on the optical models of the LS and PMTs. Updates and improvements made to these models based on the latest knowledge and measurements are highlighted. The approaches used to determine key parameters in these models are also discussed. Additionally, the electronics simulation is introduced in Sec. 5, which provides a comprehensive model of the detector response. Section 6 describes the calibration and event reconstruction procedures and reports the newly predicted energy resolution in JUNO, providing in Sec. 7 a decomposition of the energy resolution to identify the major contributing factors. Finally, a summary is given in Sec. 8.
-
A schematic drawing of the JUNO detector is shown in Fig. 1. The detector consists of three main components: the central detector (CD), the water pool (WP), and the top tracker (TT).
Figure 1. (color online) Schematic drawing of the JUNO detector. TT: top tracker, LS: liquid scintillator, CD: central detector, SS: stainless steel.
The CD is located in the center of the WP and is composed of a 35.4 m inner diameter acrylic sphere with a wall thickness of 12.4 cm. The vessel is supported by 590 connecting bars and filled with 20 ktons of LS. The LS recipe consists of linear alkylbenzene (LAB) as the solvent, 2.5 g/L of 2,5-diphenyloxazole (PPO) as the fluor, and 3 mg/L of p-bis-(o-methylstyryl)-benzene (bis-MSB) as the wavelength shifter. Four purification plants have been constructed to ensure the desired radiopurity and target attenuation length. These plants allow for alumina filtration, distillation, water extraction, and gas stripping. The radiopurity of the LS will be monitored by the OSIRIS detector [27] before and during the filling process.
Scintillation and Cherenkov light produced in the LS is collected by PMTs. PMT modules are the basic units used to assemble PMTs on the stainless steel (SS) supporting structure of the CD, accommodating multiple PMTs in each module. Among the LPMTs, 5,000 are dynode PMTs manufactured by Hamamatsu Photonics K.K. (HPK), while the remaining LPMTs are a new type of Micro-Channel Plate Photomultiplier Tubes (MCP-PMTs) [28] manufactured by Northern Night Vision Technology Co. (NNVT). Furthermore, a fraction of the NNVT MCP-PMTs, manufactured with improved technologies during photocathode fabrication, are labeled as NNVT HQE MCP-PMT. The arrangement of these LPMTs has been extensively studied, and their photocathode coverage has been determined to be 75%. Each LPMT is equipped with a protection cover to prevent chain explosions in case of PMT breakage under the high water pressure. The protection consists of a front acrylic cover with a thickness of approximate 10 mm that follows the shape of the photocathode with a minimum 2 mm gap filled with water, working together with a back SS cover to provide full protection [29]. SPMTs are mounted in the gap between LPMTs, providing an additional
$ \sim $ 3% photocathode coverage. The PMT modules made of SS play a role of optical isolation between the outer and inner WP. In addition, a comprehensive calibration system is also implemented in the CD to control systematic errors in the energy scale and energy response of the detector [30].The output signals from each LPMT are duplicated and amplified by two Trans-Impedance Amplifiers (TIA) with different amplification factors. This configuration allows for a good signal-to-noise ratio for low-energy events, such as the inverse beta decay process (IBD), with the high gain configuration. At the same time, it maintains a large dynamic range with the low gain configuration to detect high-energy events. The amplified analog signals from each TIA are then digitized using a custom flash analog-to-digital converter (FADC). The FADC provides a sampling rate of 1 GHz and a resolution of 14 bits. The digitized waveform data are processed in the Field Programmable Gate Array (FPGA) to extract the charge and time information. The charge and time values of triggered events, along with the waveform data, are then sent to the Data Acquisition (DAQ) cluster. In the DAQ cluster, complex online event classification algorithms are employed to decide whether to save the full waveform or just the charge and time information onto the disk. This decision is made in order to optimize the data bandwidth and storage capacity. For IBD events, the waveform data is typically recorded to allow for more precise processing and analysis offline. More details of the LPMT electronics and trigger design are available in [31−33]. Regarding the SPMT electronics, CatiROC chips are used to readout the signals and only charge and time information are saved. More information is outlined in [34].
The WP serves as a water Cherenkov detector and provides shielding from external radioactivity. The TT, located above the CD and WP, is a plastic scintillator detector used to tag cosmic muons [35]. More detailed information on the JUNO detector can be found in [9, 11].
-
Reactor electron antineutrinos with energy exceeding 1.8 MeV can be detected within the LS target from interactions of the products of the IBD process:
$ \bar{\nu}_e+p \to e^+ + n $ . The positron ($ e^+ $ ), being much lighter than the neutron (n), carries away almost all the kinetic energy of the electron antineutrino.It quickly deposits its kinetic energy in the LS and annihilates into two 0.511 MeV gammas, generating a prompt signal within an energy range of 1 to 12 MeV. After scattering in the LS with an average lifetime of approximately$ 200\; \mu $ s, the neutron is captured by a hydrogen or a carbon nucleus, consequently producing a delayed gamma signal of 2.22 MeV or 4.95 MeV, respectively. The coincidence signature of the prompt and delayed signals can effectively distinguish the IBD signals from the backgrounds. An accurate measurement of the positron energy and understanding its energy resolution in the JUNO detector are crucial for the NMO determination. The impact from the neutron recoils will not be discussed in this article, but is included in the NMO paper [18].The energy deposition (
$ E_\text{dep} $ ) of positrons primarily occurs through excitation and ionization processes. Since the concentrations of PPO and bis-MSB in the LS are relatively low, the deposited energy is mainly transferred to LAB molecules. Subsequently, the excited LAB molecules can transfer their energy to fluor molecules through complex molecular-scale processes [36]. During the de-excitation of the fluor molecules, scintillation photons are emitted. However, the number of emitted scintillation photons is not proportional to the deposited energy due to the quenching phenomenon [37], where some excited molecules release energy without radiation emission. Various models exist to describe this quenching effect, the most common one is the semi-empirical formula proposed by Birks [38], which is used in this study,$ \begin{aligned} \frac{dN_ \text{ScintOP}}{dx}=Y\frac{\dfrac{dE}{dx}}{1+kB\dfrac{dE}{dx}}, \end{aligned} $
(1) where
$ N_ \text{ScintOP} $ represents the number of optical photons produced by scintillation process, Y corresponds to the scintillation light yield in unit of photons per MeV without the presence of quenching,$ dE/dx $ denotes the energy loss per unit path length (also known as the stopping power), and$ kB $ is the Birks' coefficient that depends on the particle type. Given a fixed energy loss ($ dE $ ) over a short path length of$ dx $ ,$ dN_ \text{ScintOP} $ follows a Poisson distribution due to the random nature of the molecular-scale energy transfer process of scintillation, where the variance of the distribution is determined by Y. However, the total number of scintillation photons ($ N_ \text{ScintOP} $ ) produced by a particle no longer follows the Poisson distribution due to fluctuations in the energy loss described by a Landau distribution, the quenching effect modeled by the Birks' Law, and the following random secondary particle generation processes:• Energetic δ-electrons may be produced by ionization. They travel a significant distance away from the primary track and typically exhibit stronger quenching due to their higher
$ dE/dx $ compared to the primary positron.• Gammas can be generated through Bremsstrahlung radiation and positron annihilation, and interact through photoelectric effect, Compton scattering, and pair production. Then the secondary electrons and positrons deposit their energies and produce scintillation photons through the aforementioned processes.
When the positron and its charged secondaries move through the LS at a velocity higher than the group velocity of light in the LS, Cherenkov photons are emitted. The Frank-Tamm formula [39] is commonly used to calculate the number of Cherenkov photons produced per wavelength per unit path length travelled by particles of charge
$ ze $ ,$ \begin{aligned} \frac{d^2N_ \text{CherenOP}}{dx d\lambda}=\frac{2\pi\alpha z^2}{\lambda^2}\left(1-\frac{1}{\beta^2 n(\lambda)^2}\right). \end{aligned} $
(2) In the formula, α represents the fine structure constant, β is the ratio of the particle's velocity to the speed of light in vacuum, and n is the wavelength-dependent refractive index of the LS. According to Eq. (2), the number of Cherenkov photons is inversely proportional to the square of the wavelength (λ), resulting in a higher Cherenkov light yield at shorter wavelengths. However, there is very limited information on the refractive index at short wavelengths, which introduces significant uncertainties in estimating
$ N_ \text{CherenOP} $ . Although the Cherenkov process contributes additional light, it has a detrimental effect on the energy resolution. This is because the Cherenkov light yield depends on the path length of the charged particle in the LS, which can vary significantly due to the generation of secondary particles. This variation in path length can cause the Cherenkov photon distribution to deviate from Poisson statistics, leading to a degradation of the energy resolution.Scintillation and Cherenkov photons produced in the LS may undergo various processes as they propagate through the detector, including absorption, scattering, reflection, and transmission at material boundaries. Absorption within the LS can be modeled as a competition among the LAB, PPO, and bis-MSB molecules. The absorption length of each component determines the probability of photon absorption. When photons are absorbed, they may be re-emitted again, and this re-emission probability is determined by the quantum yield (QY) of the fluors in the LS. Scattering within the LS is dominated by Rayleigh scattering rather than Mie scattering, because LS is expected to have an extremely low concentration of large-sized impurities. The probability of Rayleigh scattering depends on the scattering length of the photons. The propagation of optical photons inside the LS can be well modeled with the optical parameters such as absorption length, re-emission probability, Rayleigh scattering length, and refractive indices of the detector components. Only a fraction of these photons will be detected by the PMTs and converted into PE. The total number of PE (
$ N_ \text{PE} $ ) is influenced by factors such as the scintillation light yield, the quenching effect, the LS refractive index, the light collection efficiency during transportation, and the photon detection efficiency (PDE) of the PMTs. A simple energy reconstruction method is to divide$ N_ \text{PE} $ , including contributions from both scintillation PE ($ N_ \text{ScintPE} $ ) and Cherenkov PE ($ N_ \text{CherenPE} $ ), by the energy scale which is defined as the average PE number per MeV calibrated by 2.2 MeV gammas from neutron capture on hydrogen at the center of CD. The average number of PE produced in both scintillation and Cherenkov processes are not proportional to$ E_ \text{dep} $ , leading to LS non-linearity (LSNL) denoted as:$ \begin{aligned} f_ \text{LSNL}(E_ \text{dep})=\frac{E_ \text{vis}(E_ \text{dep})}{E_ \text{dep}}, \end{aligned} $
(3) where
$ E_ \text{vis} $ is the visible energy defined as the expected reconstructed energy assuming a perfect energy resolution. Due to the aforementioned fluctuations in the energy loss and Cherenkov process,$ N_ \text{PE} $ does not follow a simple Poisson distribution but exhibits a larger standard deviation. To account for these fluctuations, a more general formula can be introduced as follows:$ \begin{aligned}[b] \sigma_ \text{PE}=&\sqrt{\sigma_ \text{ScintPE}^2 + \sigma_ \text{CherenPE}^2 + 2\cdot cov[N_ \text{ScintPE}, N_ \text{CherenPE}]}\\ &> \sqrt{\langle N_ \text{PE}\rangle}, \end{aligned} $
(4) where
$ \sigma_ \text{PE} $ is the standard deviation of the$ N_ \text{PE} $ distribution. Equation (4) includes the correlation between the fluctuations in the scintillation and Cherenkov radiation which has been observed in the JUNO simulation and will be further discussed in the following section. The value of$ N_ \text{PE} $ also depends on the position of the IBD event. This position-dependent effect is referred to as the detector non-uniformity with respect to the amount of the collected photons for the same energy release, which arises from geometrical factors. Consequently, the energy resolution is also position-dependent in the JUNO detector.The PMT and electronics responses can introduce additional fluctuations to the detected
$ \langle N_ \text{PE}\rangle $ . The PMT response includes various factors such as the single photoelectron (SPE) resolution, which arises from the inherent fluctuations in the dynode or MCP charge amplification process. Additionally, effects such as dark count rate (DCR) and afterpulses can contribute to the fluctuations if they occur within the signal readout window and mix with the PE signal. The transit time and transit time spread (TTS) of the PMTs can impact the accuracy of vertex reconstruction and further influence the energy resolution, due to the detector non-uniformity. On the electronics side, the presence of electronics noise can further smear the charge measurements. The digitization process can also contribute to the overall fluctuations. All of these factors and parameters have been well characterized during the PMT mass testing [21] and electronics development phases [40]. Their impacts on the resolution depend also on the specific vertex and energy reconstruction algorithms employed.Based on the waveform recorded in each readout channel, the charge and time information can be extracted by the waveform reconstruction algorithm. Then, the charge can be calibrated and converted to the number of PE using PMT calibration algorithms. Eventually, the event position and energy can be obtained from the vertex and energy reconstruction algorithms using charge and time information of each readout channel as inputs. During reconstruction, additional fluctuations, from PMT calibration, residual energy non-uniformity, among others, will be added to the reconstructed energy, where their amplitudes depend on the performance of the algorithms.
The main factors that may impact the energy resolution in the JUNO detector are summarized below, and also illustrated in Fig. 2 with key factors indicated:
Figure 2. (color online) Summary of the key factors impacting the energy resolution throughout the processes of light production, propagation, detection, calibration and reconstruction.
1. Positron energy loss fluctuation in LS and variation in its secondary particle generation in LS.
2. Scintillation light yield and ionization quenching effect during the scintillation photons production.
3. Cherenkov photons emission.
4. Light propagation and detection processes.
5. PMT and electronics responses.
6. Calibration scheme and energy reconstruction algorithm.
All of the mentioned items are carefully considered in the full JUNO detector simulation chain and algorithms of event reconstruction. Items 1-4, which include the energy deposition, scintillation (Sec. 4.2) and Cherenkov light production (Sec. 4.3), photon propagation (Sec. 4.4) and detection (Sec. 4.6), are taken into account in the simulation. Item 5, which involves the PMT and electronics responses (Sec. 5.1 and Sec. 5.2), is modeled in a dedicated electronics simulation. Item 6 is included in the data processing of the calibration (Sec. 6.1) and reconstruction (Sec. 6.2).
-
The JUNO detector simulation software has been developed based on GEANT4 [41−43] version 10.04.p02. The simulation software, along with other offline data processing modules, is implemented within the SNiPER (Software for Non-collider Physics Experiment) [44, 45] framework. Further details are outlined in [46]. Aspects that potentially influence the energy resolution which are modeled by the detector simulation are discussed in the following sections. Aspects include: the detector geometry, particle interactions in the detector media, light production, LS optical model describing the light propagation, and PMT optical model accounting for the reflection and angular-dependent PDE of the PMTs.
-
In detector simulation, mechanical design drawings, survey data and the best-known values for the optical properties of materials and component dimensions are at the basis to enable reliable simulation of light propagation, collection, and the detector non-uniformity response. The key components modelled include:
• LS sphere with a radius of 17.7 m.
• Acrylic sphere with an inner radius of 17.7 m and a thickness of 12.4 cm.
• Inner water buffer located between the acrylic sphere and PMT module with a thickness of 1.8 m.
• Detailed LPMT geometry [21], including photocathode, reflective aluminum film, inner electrode, and supporting structures, which are important for tracking photons inside the LPMTs. The 17,612 LPMTs are positioned in the inner water buffer facing the acrylic sphere, with a minimum distance of 1.42 m between acrylic and the center of the photocathode. In addition, the acrylic and SS protection covers [29] are also modeled, located at the front and the back of the LPMTs, respectively.
• 25,600 SPMTs [26] located within the gaps between LPMTs with a position distribution arranged to guarantee a uniform SPMT density.
• 590 acrylic support nodes with detailed geometry and material definition, with corresponding SS support structure.
• The PMT modules which optically isolate the inner and outer water pools.
• Calibration system [30] with its anchors mounted on the acrylic sphere.
• Latticed SS shell used to support the LS acrylic sphere and mount the PMT modules and its own support structures.
• Outer water pool with a height and radius of 43.5 m.
• Chimney and calibration house on the top of the detector.
• TT consisting of 3 layers of plastic scintillator.
The optical properties of each key component in the CD have been well defined in the simulation. In Figure 3a, the refractive index of the acrylic, water and PMTs' glass (Pyrex) is summarized as a function of wavelength. The refractive index of acrylic and glass is measured at 5 different wavelengths, indicated as markers. These measured values are then extrapolated to other wavelengths using the dispersion relation:
Figure 3. (color online) (a) Refractive index of acrylic (red), PMT glass (green), and water (blue) as a function of wavelength. (b) Attenuation length of acrylic (red) and water (blue) as a function of wavelength. Markers indicate measurements, while solid lines are interpolations or extrapolations.
$ \begin{aligned} n^2 - 1 = \frac{p_0 \lambda^2}{\lambda^2 - p_1}, \end{aligned} $
(5) where λ is the wavelength.
$ p_0 $ and$ p_1 $ are two parameters obtained by fitting the 5 measured data points using Eq. 5. Refractive index below 300 nm is not set in the simulation, because a negligible number of photons from LS are expected to reach the acrylic or PMTs in this wavelength region, due to the strong absorption of LS. The refractive index of water is obtained from measurements in [47]. The attenuation length of acrylic is obtained by analyzing the transmittance data published in [48, 49], shown as the red curve in Fig. 3b. The attenuation length of water, shown as the blue curve in Fig. 3b, is assumed to be 40 m at 430 nm and its dependence with wavelength is provided by the Daya Bay collaboration. The absorption length is not set for the PMTs' glass, since this effect is implicitly included in the PMTs' PDE. The reflectivity of the SS components is assumed to be 53.5% without angular and spectral dependence, which is calculated using Fresnel’s equations with the refractive index obtained from [50]. This quantity will be measured in situ at JUNO and will be updated in the simulation in the future. The optical properties and modelling of the LS and PMTs will be discussed in more detail in Sec. 4.4 and 4.6, respectively. -
The low energy Livermore model [51] is chosen to simulate the electromagnetic processes of electrons, positrons, gammas, hadrons, and ions in the JUNO simulation. These processes include the photo-electric effect, Compton scattering, Rayleigh scattering, gamma conversion, Bremsstrahlung, and ionisation. The Livermore model is expected to possess improved simulation accuracy of energy deposition and secondary particle generation compared to the standard electromagnetic models, because it directly uses shell cross-section data, instead of parameterizations of these data. The Livermore model can handle particle interactions with energies down to 10 eV and is valid for elements with atomic number between 1 and 99. The physics list construction of G4EmLivermorePhysics provided by GEANT4 is directly used in the JUNO simulation. The production cuts are set to 1 mm for gammas and 0.1 mm for electrons and positrons, which are used to determine the energy thresholds of secondary particle generation during GEANT4 tracking. These cuts can significantly influence the determination of the Birks' coefficient. Beside the low energy electromagnetic processes, a complete physics list has been constructed in simulation, including hadronic processes, ion physics, lepton and gamma-nuclear interactions, absorption, short-lived particles, radioactive decays, and optical processes.
During simulation in GEANT4, the trajectory of a particle is tracked step by step. At each step, GEANT4 calculates the energy deposited by the particle, taking into account the fluctuations associated with the energy deposition process. The deposited energy in each step is then used to calculate the number of scintillation photons produced using Birks formula (Eq. 1) with the addition of Poisson fluctuations. The total number of scintillation photons produced in a physical event is obtained by summing over all the steps of both primary and secondary particles.
-
Constraining the Birks' coefficient (
$ kB $ ) through the energy non-linearity response of an LS-based detector, such as those measured by the Daya Bay and Borexino detectors [52, 53], can be challenging. This is because$ kB $ is strongly correlated with the Cherenkov contribution, as discussed in the Daya Bay publication [52]. However, it is possible to determine$ kB $ by fitting the LS non-linearity data obtained from table-top measurements. These measurements typically involve electrons with energies below the Cherenkov threshold, typically around 0.2 MeV. Nevertheless, the$ kB $ values reported in publications are not directly applicable for use in detector simulations, because when extracting$ kB $ using Eq. 1, only the$ dE/dx $ of the primary particle is considered, usually calculated using tools like ESTAR [54] or SRIM [55]. This calculation does not account for the production of secondary particles. However, in the JUNO detector simulation, a significant number of secondary particles can be generated, carrying a fraction of the primary particle's energy. These secondary particles are independently tracked in GEANT4. To address this issue, a new fitting method has been explored to determine$ kB $ , which can be directly used in simulation. This fitting method takes into account the production and distribution of energy among secondary particles during simulation.The electron quenching effect of JUNO LS has been investigated by two groups, the Institute of High Energy Physics (IHEP) and the Technical University of Munich (TUM). Both groups use electrons generated from the Compton scattering of incident gammas emitted by radioactive sources. The gamma sources are placed outside a cylindrical LS container and can be rotated around it. The data set marked with brown color is measured by the IHEP group [56], in which a
$ ^{22} \text{Na} $ source is used and only energies below 0.2 MeV are considered. While the data set shown in blue is taken by the TUM group, in which gamma source of$ ^{137} \text{Cs} $ is rotated to change its position, and four different data sets are obtained, corresponding to four different rotation angles. In Fig. 4, only one TUM data set is shown, while the remaining data sets are used to estimate the systematic errors.
Figure 4. (color online) Combined fitting results on the electron quenching data below 0.2 MeV with production cuts of 1 mm for gammas and 0.1 mm for electrons and positrons. The data set shown in brown color is from [56]. The data set in blue is taken by the JUNO collaborators of TUM, in which the systematic errors are estimated by measuring the quenching curves at four different incident angles of a
$ ^{137} \text{Cs} $ gamma source with respect to the LS sample. The red curve is the result of the fit using Eq. (7).At each data point in Fig. 4, electron interactions in the LS are simulated using GEANT4, with the initial electron energy
$ E_ \text{true} $ derived from the Compton-scattered gamma. The production cuts in simulation are kept the same with those used in JUNO. It is important to note that the choice of production cuts can affect the energy deposition and, consequently, the value of$ kB $ . Therefore, it is crucial to keep the production cuts consistent throughout the study to ensure the reliability of the results. During the simulation, the energy deposition in each step is recorded. The visible energy,$ E_ \text{vis} $ , for each event is calculated by summing up all the steps in the event as indicated by Eq. 6, where S is a fitting parameter and represents a normalization factor. A combined fit is then performed on the collected experimental data on electron quenching by minimizing the value of the$ \chi^2 $ function Eq. 7:$ \begin{align} E_ \text{vis} &= \sum\limits_ \text{step} \frac{S\cdot dE}{1+kB\left(\dfrac{dE}{dx}\right)} \end{align} $
(6) $ \begin{align} \chi^2 &= \sum_n\sum_i\left(\frac{\bar{E}_ \text{vis}/E_ \text{true}-M_i^n}{\sigma_i^n}\right)^2, \end{align} $
(7) where n denotes the number of data sets, i refers to the number of data points in each data set,
$ {E}_ \text{vis} $ is the average visible energy of the simulation samples at the corresponding electron energy$ E_ \text{true} $ , and M indicates the measured ratio of$ E_ \text{vis} $ to$ E_ \text{true} $ with its corresponding error denoted by σ. Without loss of generality, we scale the$ E_ \text{vis}/E_ \text{true} $ ratio to 1 at 0.1 MeV for each data set. From the fitting, Birks' parameter$ kB $ is determined to be$ (12.1\pm0.3)\times 10^{-3} \text{g}/ \text{cm}^2/ \text{MeV} $ , with the uncertainty including both statistical and systematic errors. -
The refractive index of LS is crucial for calculating the number of Cherenkov photons using Eq. 2. In our simulation, we use the LS refractive index data shown in Fig. 5, which were obtained with a precision better than 0.01% at five different wavelengths (indicated as markers) using the V-prism refractometer [57]. Then, we employ the dispersion relation given by Eq. 5 to extend the refractive index down to 200 nm. For the wavelength range between 120 nm and 200 nm, we adopt the refractive index shape from the KamLAND experiment [58] and scale it to match the refractive index of the JUNO LS at 200 nm. Finally, we perform a linear extrapolation to estimate the refractive index down to 80 nm, which has a value close to unity.
Figure 5. (color online) Refractive index of LS as a function of wavelength. The markers represent the data points measured by experiments using the V-prism refractometer [57]. Then, the dispersion relation (Eq. 5) is used to extend the refractive index down to 200 nm (yellow region). For wavelength range between 120 nm and 200 nm (gray region), the refractive index shape is from the KamLAND experiment. For wavelengths below 120 nm (green region), a linear extrapolation is utilized.
Cherenkov photon production is handled by GEANT4, although with some modifications. This is necessary because GEANT4 assumes a refractive index monotonically increasing with photon energy, where the maximum refractive index corresponds to the maximum photon energy. However, this is not the case for the LS refractive index. To address this, we have enhanced the Cherenkov process in GEANT4 to handle more general forms of refractive index vs photon energy curves. Initially, we use the photon energy range of the refractive index that are above the Cherenkov threshold based on the velocity information of incident particles. Subsequently, we calculate the number of Cherenkov photons for each energy range using Eq. 2. Finally, we sample the energies of the emitted Cherenkov photons according to the LS refractive index curve.
It is essential to acknowledge that the refractive index and re-emission probability of the LS carry significant uncertainties, especially at shorter wavelengths such as the vacuum ultraviolet region (
$ <200 $ nm). These uncertainties introduce a potential bias in predicting the Cherenkov light yield. To address this, we introduce a Cherenkov light yield factor, denoted as$ f_C $ , which is used to adjust the Cherenkov light yield in the simulation. The Cherenkov light yield factor is applied as follows:$ \begin{aligned} N_ \text{CherenOP} = f_C \cdot N_ \text{CherenOP}^ \text{G4}. \end{aligned} $
(8) Here,
$ N_ \text{CherenOP}^ \text{G4} $ represents the calculated Cherenkov photons number obtained from GEANT4 using the LS refractive index as input. The factor$ f_C $ is determined by constraining it with the LS energy non-linearity and energy scale measurements performed by the Daya Bay detectors. Further details regarding this constraint will be discussed in Sec. 4.5. -
Photon propagation in the LS is governed by the LS optical model, which takes into account the processes of emission, scattering, absorption, and re-emission. In this model, the three components of the LS, namely LAB, PPO, and bis-MSB, are treated as a collective entity and share a set of equivalent optical parameters. These parameters include the photon emission spectrum, absorption and scattering lengths, and the quantum yield after photon absorption. The LS optical model is illustrated in Fig. 6. During the propagation of photons, they can either be absorbed or scattered by the LS, depending on the absorption and scattering lengths and their respective energy. If a photon is absorbed without undergoing the re-emission process, its trajectory is terminated. However, if re-emission occurs, a new photon is generated, and its energy is sampled from the LS emission spectrum. This newly generated photon continues its propagation within the LS.
Figure 6. (color online) The schematic diagram of light propagation in the LS optical model.
The reflection and transmittance at the interfaces between the LS and the acrylic, as well as between the acrylic and the inner water buffer, or between the water and the PMT glass bulb, are accounted for by GEANT4 using the Fresnel formula and the predefined refractive index of the LS, acrylic, water, and glass (Fig. 3a). Additionally, photons have a probability of being absorbed by the acrylic or the water, determined by their respective absorption lengths. Scattering within water uses Rayleigh scattering lengths calculated from the refractive index by GEANT4. Scattering within acrylic is currently neglected due to the lack of scattering length information. The boundary processes at the surfaces of the other detector components are also described by GEANT4, utilizing predefined optical properties such as reflectivity. After photon propagation, a fraction of the generated photons can impinge upon the PMT photocathodes, where they are further handled by the PMT optical model discussed in Sec. 4.6.
The LS optical properties employed in the LS optical model are obtained either from bench tests or inherited from the Daya Bay experiment. These properties are summarized as follows:
• Emission spectrum: the LS optical model employs the emission spectrum of bis-MSB, considering the large size of the JUNO detector. After undergoing several cycles of absorption and re-emission processes, the scintillation photons are expected to shift towards longer wavelengths and be primarily dominated by the emissions from bis-MSB, rather than PPO fluor. The bis-MSB emission spectrum is measured using a Fluorolog Tau-3 spectrometer, as shown in Fig. 2 of [59]. This spectrum is used to sample the energies of the scintillation photons during their production induced by ionization excitation and the re-emission process.
• Time profile: the time profiles employed in the simulation are obtained from dedicated measurements by exciting the LS with different charged particles, such as electrons, protons from neutron recoils, and alphas. The measured time profile is fitted using four exponential components, in which time constants of electrons/positrons/gammas are found to be 4.5 ns (70.7%), 15.1 ns (20.5%), 76.1 ns (6.0%), and 397 ns (2.8%) [60]. The time profile allows for the sampling of timing information for each scintillation photon. For the re-emission process, a single exponential component with a time constant of 1.5 ns, as measured in [61], is employed.
• Rayleigh scattering length: the Rayleigh scattering length of the JUNO LS is obtained from measurements reported in [57] and is shown in Fig. 7, yielding a value of 27.0 m at 433 nm. This value is extrapolated to other wavelengths using the Einstein-Smoluchowski-Cabannes formula [62].
Figure 7. (color online) Quantum yield (left) and Rayleigh scattering length (right) as a function of wavelength.
• Absorption length: the attenuation length (
$ L_{att} $ ) of the JUNO LS is assumed to have a target value of 20 m at 430 nm, which comprises both the absorption length ($ L_{abs} $ ) and the scattering length ($ L_{sca} $ ). This can be expressed as$ \dfrac{1}{L_{att}(\lambda)} = \dfrac{1}{L_{sca}(\lambda)}+\dfrac{1}{L_{abs}(\lambda)} $ . Subtracting themeasured scattering length from the attenuation length yields an absorption length of 77 m at 430 nm. The wavelength dependence of the absorption length is assumed to be the same as that of the Daya Bay LS [63].• Quantum yield: the spectrum of the quantum yield is assumed from the Daya Bay experiment as shown in Fig. 7, which has been fine-tuned to achieve agreement between the LS energy non-linearity in data and the simulation. The LS replacement experiment at Daya Bay [25] indicates that the differences in PPO and bis-MSB concentrations between the Daya Bay LS (3 g/L PPO and 15 mg/L bis-MSB) and the JUNO LS have a negligible impact on the quantum yield.
The LS optical model, along with the optical properties, has been implemented in the detector simulation. The model is essential for a reliable simulation of photons in the detector, producing an accurate result on the light collection efficiency.
-
To predict the absolute light yield in LS, the two remaining unknown parameters, Y in Eq. 1 and
$ f_C $ in Eq. 8, can be constrained based on the LS energy non-linearity curve and energy scale measured by Daya Bay. This is done under the assumption that the scintillation light yield and energy non-linearity response are the same for both the Daya Bay and JUNO LS.To ensure consistency between the JUNO and Daya Bay simulations, several modifications are made to the Daya Bay detector simulation, including:
• The same Livermore low-energy electromagnetic model is used, with the same production cuts.
• The same quenching effect model and the
$ kB $ parameter discussed in Sec. 4.2.2 are employed.• The LS refractive index and the improved Cherenkov process are the same.
• The LS optical properties, such as the emission spectrum, quantum yield, scattering length, and time profiles, are assumed to be identical. The shape of the LS absorption spectrum is the same, however, the absolute absorption lengths differ. In Daya Bay, the absorption length is 27 m at 430 nm, as measured in the experiment.
• In Daya Bay, each PMT PDE is considered to be the measured quantum efficiency (QE), assuming a 100% collection efficiency (CE). In JUNO, both QE and CE are taken into account, and their product, PDE at normal incidence, is constrained by the PMT mass testing data [21]. The same PMT optical model (more details in Section 4.6) is used to describe the PMT reflection and angular responses. However, different optical properties of the photocathode are applied for the 8-inch PMTs in Daya Bay and the 20-inch PMTs in JUNO.
• The optical properties of other detector components in Daya Bay remain unchanged.
After making these modifications, the determination of the parameters Y and
$ f_C $ in the Daya Bay simulation is carried out as follows:1. Fix the scintillation light yield Y of the LS and tune
$ f_C $ in the Daya Bay simulation to match the gamma energy non-linearity curve in [52]. By tuning$ f_C $ , the ratio of scintillation photons to Cherenkov photons, which determines the shape of the energy non-linearity curve, can be adjusted. Figure 8a shows that good agreement can be achieved between the simulation and calibration data in Daya Bay.
Figure 8. (color online) Comparisons of the LS energy non-linearity and the energy scale curves between the Daya Bay simulation and experimental data with the newly constrained parameters of Y and
$ f_C $ .2. Simulate a
$ ^{60} $ Co radioactive source at the center of the detector in the Daya Bay simulation to obtain the visible energy. The energy scale in Daya Bay is defined as the average PE number per MeV using$ ^{60} $ Co events at the detector center.3. Determine the value of Y and
$ f_C $ by comparing the scaled visible energy to the expected value of$ ^{60} $ Co, as shown in Fig. 8b. In this procedure, the ratio$ f_C/Y $ is kept constant to ensure that the energy non-linearity response remains unchanged. Finally, Y is determined to be 9846/MeV, while$ f_C $ is 0.52.Table 1 summarizes a few parameters used in this work, allowing for the prediction of the detector energy resolution in JUNO.
Y [/MeV] $ kB [ \text{g}/ \text{cm}^2/ \text{MeV} $ ]
$ f_C $ 
Production cut [mm] gamma $ e^+/e^- $ 
9846 12.1×10−3 0.52 1.0 0.1 Table 1. A summary of a few parameters used in this work to predict the energy resolution in JUNO.
-
The PDE responses of LPMTs in JUNO are determined using the results from the LPMT mass testing setups [21] and the developed PMT optical model [22]. Two mass testing setups, the scanning station and the container system, are used for LPMT acceptance tests and performance evaluation. The container system evaluates the characteristics of each LPMT, including PDE, DCR, gain, and features of SPE. A large area pulsed light source with a central wavelength of 420 nm is used to illuminate the entire photocathode of the LPMTs. The scanning station performs more detailed characterizations for LPMTs using 7 LED sources deployed along the longitude. The light beam from each LED is approximately perpendicular to the PMT surface and has a diameter of about 5 mm with a central wavelength of 420 nm. By rotating the PMTs, the whole photocathode can be scanned. However, even if only about 5% of JUNO LPMTs are tested following this full scan procedure, the results are representative of the total JUNO LPMT production.
The number of PMT's detected photons is assumed to follow a Poisson distribution. The average value measured with the container system is converted to the PDE defined by the scanning station system for all measured PMTs, because the light intensities in this system are calibrated by a reference PMT. The conversion coefficients are determined by comparing the PDEs measured by both the scanning station and the container system for the three different types of LPMTs, including HPK dynode-PMT, NNVT MCP-PMT, and NNVT HQE MCP-PMT. The PDE given by the scanning station is defined as the averaged value across the PMT area, and it is determined from PDEs measured by the 7 LEDs and their surface area weights. The area weights, calculated based on the respective positions on the photocathode of the 7 LEDs, are dependent on the LPMT types (HPK and NNVT), and are summarized in Table 2.
LED1 LED2 LED3 LED4 LED5 LED6 LED7 NNVT 4.8% 9.0% 12.6% 17.2% 20.0% 18.0% 18.4% HPK 4.5% 8.8% 13.5% 17.1% 20.5% 18.6% 17.0% Table 2. Surface area weights of 7 LEDs for NNVT and HPK PMTs.
In detector simulation, the position-dependence of the PDE along the latitude of the photocathode is modeled using the CE curves. The PDE is considered as a product of the QE and the CE. Figure 9a shows the CE as a function of the polar angle in the PMT's local coordinates for the NNVT MCP-PMT (red), NNVT HQE MCP-PMT (blue), and HPK dynode-PMT (violet). The CE curves are obtained by averaging the PDEs measured by the 7 LEDs along the zenith angles in the scanning station. The maximum CE values are set to 100% for the NNVT MCP-PMT and 93% for the HPK dynode-PMT, indicated by the electrostatic simulation results from HPK and NNVT private communications.
Figure 9. (color online) (a) Angular dependence of collection efficiency for NNVT MCP-PMT (red), NNVT HQE MCP-PMT (blue), and HPK dynode-PMT (violet); (b) Wavelength dependence of quantum efficiency for NNVT MCP-PMT (red), NNVT HQE MCP-PMT (blue), and HPK dynode-PMT (violet).
To compute the QE of a given LPMT, the following equation is used:
$ \begin{aligned} \text{QE}= \text{PDE}/\sum_{i=1}^{7}( \text{CE}_i\times w_i), \end{aligned} $
(9) where i represents the index of the LED in the scanning station and
$ w_i $ denotes the surface area weight of the i-th LED as summarized in Table 2. This calculation is performed for each LPMT in the CD to ensure that the product of QE and CE is consistent with the measured PDE, as published in [21]. The QE spectral responses are also implemented in the simulation based on laboratory measurements, as shown in Fig. 9b. The NNVT MCP-PMT (red) and NNVT HQE MCP-PMT (blue) have identical QE spectral responses, while the HPK dynode-PMTs (violet) have a different response at higher wavelength ($ \lambda>500\; \text{nm} $ ). It is assumed that LPMTs of the same type share the same CE curve and QE spectral response in the simulation.The PMT optical model [22] takes into account the angle of incidence (AOI) dependence of the PDE, as well as the reflections on the photocathode and the optical processes inside the PMTs. In this model, the PMT window is treated as a multi-layer optical stack, consisting of water, PMT glass, anti-reflective coating, photocathode, and vacuum layers from outside to inside. The anti-reflective coating and photocathode are considered coherent layers due to their comparable thicknesses with the light wavelength. The model incorporates the light interference effect and the multiple reflections between adjacent boundaries using the transfer matrix method.
The refractive index n, extinction coefficient k, and thickness d of the anti-reflective coating and photocathode are determined in the wavelength range of 390 nm to 500 nm by analyzing the reflectance data of NNVT MCP-PMT, NNVT HQE MCP-PMT, and HPK PMT immersed in the LAB liquid. The optical properties of the other components inside the PMTs, such as dynode (MCP), supporting structure, and aluminum film, are constrained by the QE data. These parameters are available in [22] and are directly used in this work. With these inputs, the PMT optical model calculates both the reflectance and absorbance for a given photon and AOI, assuming uniformity for the anti-reflective coating and the photocathode. The reflected photons are then propagated by GEANT4. The absorbance at a specific AOI is converted to the QE using the QE calculated above and absorbance information at normal incidence. This optical model has been integrated into the detector simulation and used in this work. It has also been applied in the Daya Bay simulation, using the optical parameters obtained from the reflectance and QE data of an 8-inch Daya Bay PMT.
-
In the detector simulation, a simple PMT optical model is used for SPMTs. It assumes that photons hitting the photocathode are 100% absorbed and converted to free PE by applying the PMTs' QE. The QE at 420 nm is implemented for each SPMT, using the value obtained from the characterization of SPMTs and published in [26]. The QE dependence on wavelength is considered according to the vendor's datasheet [64].
-
In the updated detector simulation, the predicted PE yield of the JUNO detector is 1665 PE/MeV, calibrated using neutron capture events on hydrogen at the detector center. The PE yield is larger compared to the previously reported value of 1345 PE/MeV [30]. This enhancement can be attributed to three factors:
• Improved PMT's PDE: In previous studies, the LPMT's PDE was assumed to be 27%. However, the LPMT mass testing has shown that the actual average PDE in CD is 30.1% [21]. This higher efficiency accounts for an
$ \sim $ 11% increase in PE yield.• More realistic PMT optical model: Previous simulations used a simplified PMT optical model that neglected the angular dependence of PDE in water and PMT photocathode reflection. By incorporating a more realistic optical model [22], we observed an additional
$ \sim $ 8% increase in PE yield.• Detector geometry updates: The detector geometry was updated based on the final mechanical design, with reflections on several detector components taken into account, leading to a
$ \sim $ 3% increase in PE yield.Furthermore, the Cherenkov to total PE ratio is found to be 1.1% for 1 MeV positrons, and the detector simulation indicates a PE yield of 1785 PE/MeV for uniformly distributed neutron capture events on hydrogen.
-
The LPMT and SPMT electronics simulations are important for predicting the energy resolution of the detector, as they contribute to the overall electronic noise and affect the signal-to-noise ratio of the recorded signals. These simulations are integrated into the SNiPER framework to provide a comprehensive modeling of the detector response.
-
In the LPMT electronics simulation, LPMT operation parameters are obtained from the PMT mass testing systems and are directly assigned to each LPMT using its serial number. These characteristics include the dark count rate (DCR), gain, and SPE charge resolution.
To simulate the LPMT dark noise, SPE pulses are uniformly distributed in a readout time window of
$ \Delta T=1\mu s $ . The number of dark noise pulses is sampled from a Poisson distribution with mean value$ \text{DCR}\times\Delta T $ .When converting a SPE from the detector simulation into a pulse in the electronics simulation, the amplitude of the pulse is modeled by a combined distribution of a Gaussian and an exponential function. The SPE charge resolution of each PMT, obtained from PMT mass testing system [21], is directly used as the sigma of the Gaussian component in the simulation. The exponential component exhibits distinct behaviors between HPK LPMTs and NNVT LPMTs. For the former, the contribution from the exponential component is found to be 1% with an average amplitude of 1.1 PE, resulting in a SPE relative variance of 0.4. While for the latter, the exponential component describes the large signal observed in the measured charge spectrum of MCP-PMT (blue line in Fig. 10a) with an average amplitude of 2.2 PE, resulting in a SPE relative variance of 0.7. The ratio of large signals in the SPE spectrum of MCP-PMT, as a function of the zenith angle on the photocathode, has been characterized by laboratory measurements, as shown in Fig. 10b. Good agreement between Monte Carlo simulation and experimental data is achieved, illustrated in Fig. 10a.
Figure 10. (color online) (a) The charge spectrum of NNVT MCP-PMT, measured with a pulsed light source, exhibits a long tail, which is described by an exponential distribution in the Monte Carlo simulation. Good agreement between data and Monte Carlo is achieved. (b) The ratio of the large signals depends on the positions on the photocathode, and is extracted from experimental data.
The transit time and TTS are also implemented for both HPK and NNVT LPMTs. Since these quantities are only measured for
$ \sim $ 25% of LPMTs [21], a random sampling based on the distributions is used to obtain transit time and TTS values for other LPMTs in the CD. The transit time and TTS values are position-dependent on the photocathode, and this effect is also considered in the electronics simulation based on dedicated measurements made during the PMT mass testing.The after-pulses of LPMTs are modeled based on the results published in [65]. Since the earliest after-pulse component occurs at about 1 μs after the primary signals, the overlap with the signals is expected to be small and thus has a negligible impact on the energy resolution. The non-linearity effect of LPMTs is also taken into account in the electronics simulation. However, it primarily affects high-energy events rather than IBD ones.
-
In the LPMTs electronics simulation, the SPE average waveform is obtained using the PMT testing data for both NNVT and HPK LPMTs. This waveform serves as a template in the electronics simulation to model the SPE response. The amplitude of the waveform is determined by the gain and SPE charge resolution discussed earlier.
To model the overshoot effect, the same model as that used in the Daya Bay experiment [66] is adopted, but with parameters tuned specifically for the JUNO experiment. The overshoot model is parameterized using an exponential function plus a Gaussian function. To capture the onset of the overshoot, the exponential function is multiplied by a Fermi function. The maximum amplitude of the overshoot is set to be 1% of the primary pulse [67].
White noise is applied to the raw waveform by sampling a Gaussian function with a mean of 10% of the SPE amplitude. This value is based on the requirements of the JUNO experiment. The digitization process is also modeled in the simulation, taking into account the effects of the front-end analog-to-digital converter (FADC) resolution, the high-gain and low-gain amplifications of the SPE signal, as well as the baseline offset. The simulation also considers the non-linearity effects of the electronics. However, it is anticipated that these effects will have a negligible impact on the energy resolution for positrons coming from IBD events.
Overall, the LPMT electronics simulation in the CD takes into account various aspects of the LPMT response, including the SPE waveform, overshoot, white noise, digitization, and non-linearity effects. These simulations aim to accurately model the electrical responses of the LPMTs and the readout electronics, and their impact on the energy resolution of the detector.
-
The SPMT electronics simulation aims to accurately model the response of SPMTs and the readout electronics. For each SPMT, a DCR value is randomly sampled and assigned based on the DCR distribution obtained from the SPMT mass production measurements [26]. The mean values of the TTS (1.6 ns) and the SPE charge resolution (33%) are shared by all SPMTs.
The key features of the CatiROC ASIC are implemented in the simulation based on the measurements in Ref. [34], including the two types of dead times, the 70 ns trigger dead time and the 9.3 μs digitization dead time. In cases where a hit is not triggered due to dead time, a fraction of its charge can still be calculated. This fraction is called the charge acceptance and is parameterized based on the time difference with the previous hit. In addition, the charge in picocoulombs is converted into an ADC unit, with two conversion coefficients used responsible for high gain (less than 10 PE) and low gain (more than 10 PE).
-
As described in Sec. 5, the PMT together with the electronics system will convert the photons to waveforms for each LPMT and to the pairs of charge and time for each SPMT. For each LPMT waveform, the charge and hit time of every pulse are reconstructed using the deconvolution algorithm [68], which exhibits smaller charge non-linearity compared to other algorithms. Meanwhile, the hit time of the first pulse, referred to as the first hit time, is of particular importance. A dedicated algorithm with a linear fit of the rising edge of the first pulse is developed to achieve a more accurate first hit time and reduce its charge dependence. Given that the charge and hit time information of PMTs are the inputs to the vertex and energy reconstruction in Sec. 6.2, the performance of the waveform reconstruction will also contribute to the energy resolution. Moreover, the CD has about O(10
$ ^4 $ ) PMTs and their characteristics may differ. Thus the reconstructed charge and hit time of each PMT must be calibrated to account for the different PMT parameters such as gain, TTS, relative PDE, and DCR. A comprehensive calibration strategy was developed in Ref. [30] to extract these parameters for all the PMTs and continuously monitor their time dependence. -
The aim of the event reconstruction is to derive the energy and vertex of the event from the charge and time information of photon hits captured by PMTs. However, for large volume LS detectors such as JUNO, various complex optical processes occur during the photon propagation. It is usually rather difficult to build a comprehensive optical model to precisely describe the photon hits of PMTs.
Given that a non-uniform detector response is one of the main contributors to the energy resolution for large LS detectors, precise vertex reconstruction is needed to correct for the energy response non-uniformity. Several algorithms have been developed for vertex reconstruction in JUNO [69, 70]. These algorithms utilize the time information of the first photon hit of PMTs, together with the residual time probability distribution function (pdf), which is mainly determined by the LS timing profile and PMTs' TTS. Meanwhile, an optical-model-independent method [71] has been developed to reconstruct the event energy in JUNO and was optimized in Ref. [72] to further improve the energy uniformity. The basic principle is to obtain the expected charge for PMTs using calibration data from the automatic calibration unit (ACU) and the cable loop system (CLS), which is then used to build a likelihood function given the observed charge of all the PMTs.
A calibration data driven, simultaneous vertex and energy reconstruction method has been developed for JUNO [73], based on the vertex and energy reconstruction methods described above. The observables used are the charge and time information of the PMTs. Calibration data with known vertex and energy are used to construct the expected charge and time response for each PMT. Given the observed and expected charge and time information of PMTs, a maximum likelihood method is developed to reconstruct the event vertex
$ {\bf{r}} $ =$ \vec{r}(r,\theta,\phi) $ and visible energy E simultaneously, using the likelihood function in Eq. 10.$ \begin{aligned}[b] &{\cal{L}}(\{q_{i}\}; \{t^{j}\}|{\bf{r}},E_{ \text{}}, t_0)\\ = \;&\prod_{i} \left(\sum_{k=1}^{\infty} P_{Q}(q_{i}|k) \times P(k, \mu_i) \right) \\ \times\; &\prod_{j} \frac{\sum^{K}_{k=1} P_{T}(t^j_{res}|r, d_j, \mu^{l}_j, \mu^{d}_j, k)\times P(k, \mu^{l}_j+\mu^{d}_j)} {\sum^{K}_{k=1} P(k, \mu^{l}_j+\mu^{d}_j)}, \end{aligned} $
(10) $ \begin{array}{*{20}{l}} \begin{aligned} &\mu_i({\bf{r}},E_{ \text{}}) = E_{ \text{}} \times \hat{\mu_i}^L(r,\theta,\theta_{\text{PMT},i}) + \mu_i^D \\ &t^{j}_{res} = t^{j} - t^{j}_{tof}(d_j) - t_0 \\ \end{aligned} \end{array} $
(11) The first product on the right side of Eq. 10 corresponds to the charge-based likelihood function. The index i runs over all the PMTs. The term in the brackets simply describes the probability of observing charge
$ q_i $ on PMT i when the expected number of PE is$ \mu_i $ , which strongly depends on both the vertex and energy of the positron. Given that photons emitted from the same particle have strong temporal correlation and usually arrive on PMTs within a few hundred ns, while PMT dark noise occurs randomly in time, a signal window of 420 ns was set to reduce the PE contamination from dark noise.$ \mu_i^D $ represents the residual PE contribution from dark noise within the signal window. As one of the most crucial ingredients of the reconstruction,$ \hat{\mu_i}^L(r,\theta,\theta_{ \text{PMT},i}) $ represents the expected number of LS PE per unit of visible energy, originating from particles within the signal window. And it is obtained using the simulated calibration data in which the calibration source is deployed at various positions in the CD. Here r and θ are the components of the particle vertex$ {\bf{r}} $ , while$ \theta_{ \text{PMT},i} $ is the angle between$ {\bf{r}} $ and the PMT position vector$ {\bf{r}}_{ \text{PMT},i} $ .$ P(k, \mu_i) $ is the Poisson probability for detecting k PE,$ P_{Q}(q_{i}|k) $ is the probability of observing charge$ q_i $ on PMT i given the charge pdf of k PE$ P_{Q}(q|k) $ , which can be constructed by convolving the SPE charge spectrum with$ P_{Q}(q|k-1) $ . Note that for PMTs which do not pass the firing threshold of$ q_i > 0.1\; \text{PE} $ , this term simplifies to$ P(0, \mu_i) $ +$ P_{Q}(q_i<0.1PE|k=1) $ *$ P(1, \mu_i) $ . Moreover, the index k ends when$ P_Q(q_i|k)<10^{-8} $ is met to simplify the calculation.The second product on the right side of Eq. 10 corresponds to the time-based likelihood function. The index j only runs over the fired PMTs satisfying
$ -100\; \text{ns} < t^j_{res} < 500\; \text{ns} $ and$ 0.1\; \text{PE} < q_{j} < K $ . A cutoff value$ K=20 $ is set for the detected nPE k to simplify the calculation. The residual hit time$ t_{res} $ of PMTs is obtained by subtracting the time of flight$ t_{tof} $ and reference time$ t_0 $ from the first hit time t of PMTs. The distance between the vertex and the PMT is denoted as d.$ \mu^l $ and$ \mu^d $ stand for the expected number of PE originating from particle or dark noise within the full electronic readout window respectively. Another crucial ingredient of the reconstruction is the residual time PDF$ P_{T}(t^j_{res}|r, d_j, \mu^{l}_j, \mu^{d}_j, k) $ , which is obtained using the same calibration data used for$ \hat{\mu_i}^L $ . The fine-grained parameterization of$ P_{T}(t^j_{res}) $ takes into account its dependence on the vertex radius as well as the distance between the vertex and the PMT. Meanwhile, the impact from dark noise is also included via an analytical approach. More details of the construction of$ P_{T}(t^j_{res}) $ can be found in Ref. [73].$ P(k, \mu^{l}_j+\mu^{d}_j) $ acts as a weight for different k values.Comparing to previous reconstruction studies, a few important updates should be mentioned. Firstly, two crucial ingredients have been made more realistic: on the one hand, the residual time pdf are derived from calibration data instead of MC simulation, and all PMT electronics effects are considered for the construction of the expected nPE map
$ \mu_i $ . Secondly, the fine-grained parameterization of the residual time pdf as well as the calibration of the time of flight as a function of the photon propagation distance makes the pdf more accurate. Thirdly, the charge and time information of PMTs are combined together to improve the vertex resolution, especially near the acrylic sphere edge. Finally, the vertex and energy are reconstructed simultaneously, which naturally handles the strong correlation between these two quantities.To evaluate the performance of the energy reconstruction in JUNO, a few sets of positron samples with different kinetic energies
$ E_k $ = (0, 0.5, 1, 2, 3, 4, 5, 6, 8, 11) MeV were produced by MC simulation as summarized in Tab. 3. For each data set, the positrons were uniformly distributed in the CD and the total statistics per set is 500 k. In addition, the$ ^{68} $ Ge calibration data listed in Tab. 3 were also produced to obtain the expected nPE map and time pdfs of PMTs. The positions and type of calibration source have been slightly optimized based on Ref. [72] to improve the energy uniformity. A realistic detector geometry with all the latest knowledge on the properties of LS and PMTs from previous sections were implemented in the simulation. For each set of positrons with fixed kinetic energy$ E^i_k $ , the simultaneous vertex and energy reconstruction using Eq. 10 was applied. The distribution of the reconstructed visible energy$ E_ \text{rec} $ was fitted with a Gaussian function ($ E^i_ \text{vis}, \sigma^i $ ) and the results are summarized in Tab. 4. The corresponding energy resolution was defined as the ratio of$ \sigma^i $ /$ E^i_ \text{vis} $ and the energy non-linearity was calculated by$ E^i_ \text{vis} $ /$ {E^i_ \text{dep}} $ where$ {E^i_ \text{dep}} $ =$ {E^i_{k}} + 1.022\; \text{MeV} $ .Type Energy Statistics Position e+ Ek=(0, 0.5, 1, 2, 3, 4, 5, 8, 11) MeV 500 k/set uniform in CD 68Ge 0.511·2 MeV 20 k/point ACU+CLS (293 points) Table 3. List of the MC simulation samples.
$ E_{k} $ [MeV]
0 0.5 1 2 3 4 5 8 11 $ E_ \text{dep} $ [MeV]
1.022 1.522 2.022 3.022 4.022 5.022 6.022 9.022 12.022 $ {E_ \text{vis}} $ [MeV]
0.9205 1.422 1.947 3.007 4.069 5.133 6.197 9.392 12.59 $ E_ \text{res} $ [%]
3.122 2.414 2.046 1.682 1.484 1.354 1.256 1.077 0.9661 $ {E_ \text{vis}} $ /
$ E_ \text{dep} $ 
0.901 0.934 0.963 0.995 1.012 1.022 1.029 1.041 1.047 Table 4. Summary of the energy reconstruction results. All the energy units are in MeV. For each set of positrons with different
$ E_k $ , the reconstructed visible energy is fitted with a Gaussian function, where$ E_ \text{vis} $ and σ represent the Gaussian mean and sigma, respectively. The energy resolution$ E_ \text{res} $ is defined as σ/$ {E_ \text{vis}} $ . In addition, we also report the ratio of the visible energy to the deposited energy.In Fig. 11, the left plot shows the energy resolution as a function of the average visible energy
$ E_ \text{vis} $ . The data points were fitted with a generic parameterization formula of Eq. 12.
Figure 11. (color online) Energy resolution (left) and energy non-linearity (right) for positrons using samples from Tab. 3. A fit with Eq. 12 was done for the points in the left plot, while interpolation was used instead in the right plot.
$ \begin{aligned} \frac{\sigma}{E_ \text{vis}}=\sqrt{\left(\frac{a}{\sqrt{E_ \text{vis}}}\right)^2+b^2+\left(\frac{c}{E_ \text{vis}}\right)^2}. \end{aligned} $
(12) In this equation, a is the statistical term mainly driven by the Poisson statistics of detected PE. The b term is a constant, independent of energy and mostly contributed by the scintillation quenching effect, Cherenkov radiation and energy non-uniformity. The c term accounts for the PMT dark noise and also the positrons' annihilation γs. The best-fit results of a, b, and c were listed as the Default Case in Tab. 5 and the fitted energy resolution is 2.95% at 1 MeV. The right plot in Fig. 11 shows the energy non-linearity curve. These results were used for the NMO sensitivity in JUNO paper [18].
The relative improvement of the energy resolution at 1 MeV with respect to the previous case is also shown in the last column. Case a [%] b [%] c [%] Eres@1 MeV [%] Sequential improvement Default 2.614 0.640 1.205 2.948 - A (- vertex uncert.) 2.581 0.667 1.206 2.925 0.78% B (- dark noise) 2.571 0.671 0.956 2.824 3.45% C (- waveform reco) 2.542 0.647 0.973 2.798 0.92% D (- SPE charge smear) 2.445 0.600 1.079 2.739 2.1% Table 5. Comparison of the energy resolution among all the cases. Here a,b and c correspond to the three parameters from Eq. 12.
A few additional checks were performed to validate the results. The left plot in Fig. 12 shows that the energy non-uniformity is within 0.4% for positrons with different energies. The right plot shows how the energy resolution changes with respect to
$ r^3 $ , which is caused mainly by the change of the total number of detected PE.
Figure 12. (color online) Energy response non-uniformity checks. Normalized average
$ E_ \text{rec} $ (left) and energy resolution (right) with respect to$ r^3 $ for positron samples with different energies. The normalization is done by dividing the average value within the FV. The Red vertical line corresponds to the FV cut. -
The contribution of major factors to the energy resolution budget is estimated at two different phases: the detector simulation phase and the calibration and reconstruction process. In the detector simulation phase, positron samples are generated at the CD center. However, in the reconstruction process, the positron vertexes are generated with a uniform distribution within the LS volume.
-
The positron energy resolution at the center of the CD is obtained from detector simulation. According to the recorded truth information, the contributions from scintillation light, Cherenkov light, and their covariance to the energy resolution are extracted and major effects in the energy resolution budget are shown in Fig. 13. It has been previously discussed that the quenching effect causes the number of scintillation PE to deviate from a Poisson distribution. Therefore, the contribution from scintillation PE can be further decomposed into two parts: the Poisson fluctuation (blue curve) and the quenching effect (red curve). The contribution from Poisson fluctuation is evaluated using the square root of the number of scintillation PE, which is given by
$ \sigma_ \text{stat} = \sqrt{N_ \text{ScintPE}} $ . This component represents the statistical fluctuations in the scintillation process. The contribution from the quenching effect is then obtained by subtracting the Poisson standard deviation from the standard deviation of the scintillation PE number distribution. Mathematically, this is expressed as$ \sigma_ \text{quench} = \sqrt{\sigma_ \text{ScintPE}^2 - \sigma_ \text{stat}^2} $ , where$ \sigma_ \text{ScintPE} $ is the standard deviation of the scintillation PE number distribution. In this figure, the yellow curve represents the contribution from Cherenkov radiation, and the dark red curve represents the correlation between the scintillation and Cherenkov processes. As expected, the gray curve, which is the square root of the quadratic sum of the four components, accurately reproduces the total energy resolution (black curve) obtained directly from the standard deviation of the total PE spectrum.
Figure 13. (color online) Energy resolution decomposition. The black curve refers to the energy resolution obtained with total PE. The blue, red, dark red, and yellow curves represent the contributions from scintillation Poisson statistics, quenching effect, Cherenkov radiation and the covariance between scintillation, and Cherenkov processes, respectively. The gray curve which superimposes with the black curve is the root squared sum of the four components.
It is worth noting that the standard deviation of the total PE is significantly larger than the Poisson standard deviation, indicating the presence of significant systematic effects in the energy resolution beyond statistical fluctuations. This highlights the importance of considering these systematic effects in the analysis of the energy resolution in the detector simulation.
-
There are a few key factors which could affect the energy resolution in the event reconstruction. Due to energy non-uniformity, the vertex uncertainty will propagate to the energy resolution. A more precise vertex could improve the energy resolution. Another important contributor comes from the PMT dark noise, which contaminates the signal photons and worsens the energy resolution. Thirdly, the inputs to the event reconstruction, namely the charge and time of PMTs which are obtained from the PMT waveform reconstruction, the charge non-linearity and waveform reconstruction uncertainty will propagate to the energy resolution as well. Finally, the PMT charge could only provide a rough estimation of the number of detected PE. The intrinsic PMT charge resolution, which includes contributions from both Gaussian and exponential components (Sec. 5.1.1), will lead to charge smearing and potentially worsen the energy resolution for the charge-based energy reconstruction. The impact of these factors on the energy resolution is decomposed in this section. Different cases listed below were compared, in which the factors above were removed sequentially.
• Default Case: vertex and energy are simultaneously reconstructed.
• Case A: true vertex is used in the energy reconstruction.
• Case B: PMT dark noise is removed in the samples and true vertex is used.
• Case C: on top of Case B, the waveform reconstruction is replaced by a toy simulation to provide the PMT charge and time in the reconstruction.
• Case D: in addition to the changes in Case C, SPE charge resolution and other electronics effects are also removed.
In the default case, the positron vertex is simultaneously reconstructed with its energy using the method described in Sec. 6.2. The vertex resolution is around 10 cm at 1 MeV and decreases at higher energies. Meanwhile, the vertex bias is less than 2 cm within the FV. In case A everything is the same as the default case, except that the MC truth vertex is used in the energy reconstruction, which is equivalent to a vertex resolution of 0 mm. Besides using the true vertex, the PMT dark noise is removed in case B. Given that there is no more dark noise, the 420 ns signal window is discarded in both the construction of the nPE maps as well as in the energy reconstruction. In case C, the waveform reconstruction is replaced by a toy simulation that smears the charge of each PE, providing both charge and time information. This smeared PMT charge from the toy simulation is used in both the construction of the nPE maps as well as the energy reconstruction. In comparison, the charge reconstructed from PMT waveforms are used in previous cases. By using these toy electronic simulation samples, we are able to remove the impact from the waveform reconstruction as well as other electronics effects except the PMT SPE charge smearing. The true vertex is used for case C. For case D, the PMT SPE charge smearing is removed on top of case C by directly using the detector simulation samples for the reconstruction. As a consequence, the new observable is the number of PE, instead of the reconstructed charge. This corresponds to the most ideal scenario and leads to the best energy resolution.
Figure 14 shows the comparison of the energy resolution and energy non-linearity among all the cases. We can see that the energy resolution improves monotonically among the cases, indicating that the energy resolution keeps improving by removing each contributing factor. The fitting results of the energy resolution for all the cases are summarized in Tab. 5. The relative improvement of the energy resolution at 1 MeV for each case with respect to the previous case is also presented. We can see that at 1 MeV, the PMT dark noise and SPE charge smearing are the two dominant contributing factors to the energy resolution with a relative improvement of about 3.45% and 2.1%, respectively. Removing the charge uncertainty from the PMT waveform reconstruction together with charge non-linearity leads to a 0.92% relative improvement. The impact of vertex resolution is relatively small, the relative improvement is about 0.78% using the true vertex. The black curve from Fig. 13 is also shown as the dashed curve here, and its corresponding energy resolution is close to case D.
Figure 14. (color online) Comparison of the energy resolution (left) and energy non-linearity (right) among the different cases. The energy resolution for each case is better with respect to the previous one, indicating that the energy resolution improves by moving from the real to the ideal situation. The black curve from Fig. 13 is represented by the dashed curve here, and its corresponding energy resolution is close to case D.
The comparison above shows the average impact of each factor on the energy resolution. In addition, it also enables us to decompose the energy resolution at each discrete energy point. With the energy resolution curves of all the cases, the default energy resolution can be decomposed into five components, as it is shown in Fig. 15. The black solid curve corresponds to the default case, which is used for the NMO analysis [18]. The teal solid curve corresponds to case D and represents the ideal energy resolution we could obtain. The four dashed curves show the contribution of each factor to the energy resolution, and the values at each energy point are calculated as
$ \sqrt{(E^{X}_ \text{res})^2 - (E^{Y}_ \text{res})^2} $ , where Y goes from case A to D and X corresponds to the preceding case. Two features stand out immediately: first, the contribution from dark noise decreases as the energy increases as expected. Moreover, in the energy range below 1.5 MeV, dark noise is the major contributor. For the energy range above 1.5 MeV, particularly the most sensitive region (1.5, 3) MeV for the NMO analysis, the SPE charge smearing is the dominant factor. Preliminary studies indicate that its impact could be partially mitigated by PE counting, especially for PMTs with 1 or 2 PE.
Figure 15. (color online) Decomposition of the energy resolution. The black solid curve represents the default case. It is decomposed into five parts: the four dotted curves correspond to the contribution of vertex resolution (red), dark noise (green), waveform reconstruction (blue), and SPE charge smear (purple) respectively. The light blue solid curve represents case D (ideal case). Taking 1.022 MeV as an example, the values for all the curves are 3.12%, 0.35%, 0.83%, 0.39%, 0.58% and 2.90%, and the decomposition can be written as
$ 3.12^2 = 0.35^2 + 0.83^2 + 0.39^2 + 0.58^2 + 2.90^2 $ . -
The energy resolution is a crucial parameter in determining the sensitivity of the JUNO experiment to NMO. However, in addition to the statistical fluctuations in the detected number of PE, several other effects impact the energy resolution in the detection of the IBD signals. This paper presents a comprehensive study of the energy resolution in the JUNO experiment, incorporating the latest knowledge and updates in the detector construction stage. This includes a better understanding of the detector structures, more precise measurements of the optical properties of the LS, comprehensive evaluations of the characteristics of both LPMTs and SPMTs, improved modeling of the spectral and angular dependencies of the LPMTs' PDE, and better constraints on the absolute scintillation and Cherenkov light yield based on data from the Daya Bay experiment. All of these updates have been implemented in the JUNO detector simulation, which also includes detailed modeling of the PMT responses and readout electronics in the electronics simulation.
By using data samples generated with the full detector and electronics simulation, a full-chain data processing of calibration and reconstruction is performed to evaluate factors that can smear the energy resolution. These factors include residual energy non-uniformity after reconstruction, accuracy of reconstructed energies from reconstruction algorithms, impacts of dark count rate and SPE charge resolution, among others. After considering these factors, it is found that an overall energy resolution of 2.95% at 1 MeV can be achieved for positrons from IBD signals, and the obtained energy resolution curve has been applied in the recent NMO analysis of JUNO. Furthermore, the contribution of major effects in the energy resolution budget is estimated. This study serves as a reference for interpreting future measurements of energy resolution in JUNO data taking and provides a guideline for understanding the energy resolution of LS-based detectors.
After the data-taking phase of the JUNO experiment begins, further updates and improvements are expected in the understanding of the energy resolution. The analysis of the collected data will provide valuable insights and allow for refinements in the modeling and calibration of the detector. These updates and improvements will contribute to the ongoing efforts to optimize the performance of the JUNO detector and enhance its sensitivity to NMO.
-
We acknowledge the Daya Bay collaboration for providing the Monte Carlo simulation software. We are grateful for the ongoing cooperation from the China General Nuclear Power Group.
Prediction of Energy Resolution in the JUNO Experiment
- Received Date: 2024-06-07
- Available Online: 2025-01-01
Abstract: This paper presents the energy resolution study in the JUNO experiment, incorporating the latest knowledge acquired during the detector construction phase. The determination of neutrino mass ordering in JUNO requires an exceptional energy resolution better than 3% at 1 MeV. To achieve this ambitious goal, significant efforts have been undertaken in the design and production of the key components of the JUNO detector. Various factors affecting the detection of inverse beta decay signals have an impact on the energy resolution, extending beyond the statistical fluctuations of the detected number of photons, such as the properties of liquid scintillator, performance of photomultiplier tubes, and the energy reconstruction algorithm. To account for these effects, a full JUNO simulation and reconstruction approach is employed. This enables the modeling of all relevant effects and the evaluation of associated inputs to accurately estimate the energy resolution. The study reveals an energy resolution of 2.95% at 1 MeV. Furthermore, the study assesses the contribution of major effects to the overall energy resolution budget. This analysis serves as a reference for interpreting future measurements of energy resolution during JUNO data taking. Moreover, it provides a guideline in comprehending the energy resolution characteristics of liquid scintillator-based detectors.
DownLoad: