A universal energy response model for determining the energy nonlinearity and resolution of e±\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$e^\pm $$\end{document} and γ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma $$\end{document} in liquid scintillator detectors

Energy nonlinearity and resolution in liquid scintillator (LS) detectors are correlated and particle-dependent. A unified energy response model for liquid scintillator detectors has been presented in details. This model has advanced a data-driven approach to calibrate the particle-dependent energy response, using both the monoenergetic γ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma $$\end{document}-ray sources and the continuous β\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta $$\end{document} spectra of 12B\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^\textrm{12}\textrm{B}$$\end{document} and Michel e-\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$e^-$$\end{document} induced by cosmic muons. Monte Carlo studies have demonstrated the effectiveness and robustness of the proposed model, in particular, the positron energy resolution can be extracted in the absence of positron sources. This work will provide a feasible approach of simultaneous calibration of energy nonlinearity and resolution for the running and future LS detectors.


Introduction
Liquid scintillator (LS) detectors, composed of LS and photosensors, such as photomultiplier tubes (PMTs), have been widely used in neutrino experiments since the first detection of reactor antineutrinos [1]. Among the reactor neutrino experiments, KamLAND has measured the θ 12 driven oscillation [2], and θ 13 was measured at Daya Bay [3], RENO [4] and Double Chooz [5]. Besides, Borexino has made remarkable contributions in solar neutrino measurements [6,7]. LS detectors also play important roles in the searches for neutrinoless double-beta (0νββ) decays, such as KamLAND-Zen [8] and SNO+ [9]. In the future, JUNO will build the largest LS detector to determine neutrino mass ordering (NMO) by a e-mail: miaoyu@whu.edu.cn (corresponding author) b e-mail: wenlj@ihep.ac.cn c e-mail: xiangzhou@whu.edu.cn precisely measuring the fine oscillation patterns with reactor neutrinos, which requires the uncertainty of the positron kinetic energy scale better than 1% and an unprecedented energy resolution of 3% at 1 MeV of visible energy [10][11][12][13][14].
For LS detectors used in reactor neutrino experiments, e + from the inverse beta decay (IBD) interaction is the target signal, however, γ s instead of βs are the most common calibration sources. Charged particles e ± lose energy mostly by ionizing and exciting the solvent molecules, while γ s generate e − (and at sufficiently high energies e ± pairs) firstly. [15]. The deexcited molecules release scintillation photons. Charged particles moving faster than light in the medium produce Cherenkov photons. Photons propagate in the detector and undergo series of optical processes until they hit PMTs and get converted into photoelectrons (PEs). The mapping between the number of photoelectrons (N pe ) and the particle deposited energy reflects the energy response of the LS detector, which is crucial to the spectral analysis of neutrino oscillations.
Due to different energy deposition processes, energy response in LS is particle-dependent. The nonlinearities in energy releases of e ± and γ have been thoroughly studied [14,16]. For LS detectors with relatively poor energy resolution, the reported energy resolution is usually calibrated with mono-energetic γ sources without considering particle types, e.g., see Ref. [17]. LS detectors with high energy resolution will have the potential to distinguish the energy resolution curves of different particle types, and the spectral analysis for physics topics like NMO, solar neutrino and 0νββ decays, etc, will benefit from more comprehensive knowledge of energy resolution. Thus, it is necessary to construct a unified energy resolution model.
The N pe of mono-energetic charged particles like e ± is determined by the underlying energy deposition and optical processes in LS. γ s can be described by e ± because they deposit energies by generating primary e ± firstly. However, each γ has a different deposition mode, among which the multiplicities and energies of primary e ± are different. Therefore, the N pe distribution of mono-energetic γ s needs to be derived from the collection of all deposition modes instead of a one-dimensional probability density function (PDF) of the kinetic energy of the primary e ± (e.g., see Fig. 7 in Ref. [16]). All in all, the nonlinearity and resolution are intrinsically correlated.
In this paper, we develop a unified energy model to describe particle-dependent energy response for LS detectors. Following a similar strategy as Ref. [16], we construct the energy response model for e − firstly, then derive the γ and e + models from it and calibrate the model with γ sources and continuous β spectra. By taking into account the dispersions in γ energy deposition modes mentioned above, the model is capable to describe nonlinearity and resolution simultaneously. The structure of this paper is as follows: in Sect. 2 details of the Geant4-based [18] Monte Carlo simulation are presented. In Sect. 3, we elaborate on the connections among particles in LS. The methods of model construction are described in Sect. 4. Then in Sect. 5 the calibration procedures and performances of the model are presented. Furthermore, we discuss the implications of calibration inputs and various energy resolution scenarios in Sect. 6. Finally we give some further remarks and summarize our main conclusions in Sect. 7.

Monte Carlo configurations
A set of Geant4-based (version 4.10.p02) Monte Carlo simulation software is developed. For simplicity, a liquid scintillator sphere is implemented as the target which is surrounded by photosensors. PMT geometry and its response are not accounted in this work in order to exclude the corresponding contribution to the LS energy model. The Geant4 physics packages and custom codes jointly describe the physics processes in the detector. The low energy electromagnetic processes are described by the default Livermore model. A custom scintillation process covers the quenching effect, photon emission, absorption and re-emission. Other optical processes including Cherenkov process, Rayleigh scattering and boundary interactions are depicted with the official Geant4 codes. Besides, the optical parameters of the liquid scintillator are key to the optical propagation. Wavelength-dependent refractive indices are derived from measurement values and the dispersion relation [19], and the attenuation length and Rayleigh scattering length are taken from Refs. [19][20][21]. Scintillation properties are taken from the measurements including emission spectrum [22], fluorescence quantum yield [23][24][25] and time pro-file [26,27]. Flexible studies with varied Monte Carlo inputs are feasible using our simulation software.
By tuning the detection efficiency of sensitive detectors, the photon statistics and hence the energy resolution vary among different simulation samples. Series of different detector configurations have been simulated under this framework by covering the range of energy resolution of several typical large-scale liquid scintillator detectors, for example JUNO (∼ 3.0%) [10], Borexino (∼ 6.0%) [7], Kam-LAND (∼ 6.5%) [2] and Daya Bay (∼ 8.4%) [3]. A baseline light yield is set as 1400 photoelectrons/MeV which corresponds to an energy resolution of electron-positron annihilation with zero kinetic energy around 2.9% in the current simulation framework. The other three configurations with light yields scaled by 1/2, 1/4 and 1/8 have the corresponding energy resolutions 3.8%, 5.5% and 7.8%. It is worth mentioning that the real corresponding relation between the light yield and the energy resolution depends on the specific experimental details, the values mentioned above are obtained under the simulation framework in this study.

Connections among e ± and γ in liquid scintillator
For generality and simplicity, the start vertices of all particles are assumed at the detector center and true N pe detected by PMTs is used. The detector-dependent instrumental effects require specific studies for each experiment and can be added to this model.
To convert N pe into the energy dimension, a general definition of "visible energy" is introduced as, where Y is the photoelectrons yield per MeV and is calibrated by the neutron capture on hydrogen using neutron calibration sources. The energy-dependent nonlinearity f and resolution R for the mono-energetic particles are defined as, where E vis and σ are the expected value and standard deviation of E vis , respectively. Similar to Ref. [16], the e − is regarded as a basic particle, which has energy response described by Eqs. (2) and (3) directly. In our model, we consider the e + deposits its kinetic energy T equivalently as e − approximately, followed by two identical 0.511 MeV annihilation γ s. Therefore, the energy response of e + can be predicted with that of e − directly, where the kinetic energy and annihilation energies are considered as two independent parts. One way to describe the two 0.511 MeV annihilation γ s is employing the 68 Ge source. The 68 Ge source is a positron source while the kinetic energy is mostly deposited in the enclosure, and the two annihilation γ s drift out and deposit energy in LS. The explicit formulae can be found in Sect. 4. There are two cases that might introduce systematic bias, annihilation in flight and formation of positronium. Annihilation in flight will generate two γ s with higher energies, and the fraction of this effect is around 1% per MeV of the kinetic energy from the simulation [16]. If positronium is formed, there are two spin states: singlet (p-Ps) and triplet (o-Ps). Around 2% triplets decay to three γ s with the total energy of 1.022 MeV, while two 0.511 MeV γ s are released in the other cases. To include these effects, one way is to parameterize the energy response of γ . If one has the full knowledge of γ energy distributions from above two effects, the energy response of e + can be corrected furthermore. In the simplified consideration of Eqs. (4)-(5), these two effects are ignored due to their relatively small contributions.
The γ s deposit energy mainly by generating primary e ± through three interactions: photoelectric effect, Compton scattering, and pair conversion. The multiplicities and energies of the primary e ± from mono-energetic γ s are different event by event (named "deposition mode") resulting from the competition of those three interactions. According to the energy-dependent response functions in Eqs. (2) and (3), each collection of primary e ± in a deposition mode corresponds to a specific visible energy distribution: where j specifies the given deposition mode, while l j and n j enumerate all the primary e − and e + in the j-th deposition mode, respectively. The positron term in Eq. (6) can be also described by the e − response according to Eqs. (4)-(5). The expected value of E γ vis j and its variance are calculated assuming that primary e ± are independent: The total visible energy distribution is supposed to be a mixed distribution for all the deposition modes. In previous nonlinearity studies, only the expected value E γ vis is considered thus a one-dimensional PDF averaging over all deposition modes is sufficient for calculation [16]. However, to involve the resolution into consideration, each deposition mode requires to be calculated as an independent distribution so that the fluctuations among different modes are not washed out. Consequently, the nonlinearity and resolution for γ s can be calculated with the knowledge of all the possible deposition modes: where w j is the weight of the j-th deposition mode. The two terms in the numerator of Eq. (10) are redefined as σ 2 nonl and σ 2 ave in Eq. (11), respectively. The second term σ 2 ave is the weighted sum of variance for all primary e ± . The first term σ 2 nonl is the weighted sum of the square of mean value deviation which reflects the dispersion of all deposition modes. Apparently, this term requires an individual calculation of each event rather than the one-dimensional PDF, because different collections of the primary e ± generated by the same E γ will transfer into different visible energy due to nonlinearity. The σ 2 nonl term will disappear if the energy response is linear in LS detectors. Therefore, the coupling of nonlinearity and resolution of e − is embedded in the γ resolution.

Model construction
According to discussions in Sect. 3, the energy response for both e + and γ can be deduced from that of e − . Therefore, it is crucial to construct nonlinearity and resolution models for e − firstly. Details are described in the following.

Nonlinearity
The two luminescent processes in LS are scintillation and Cherenkov processes. It is necessary to discuss the contributions to energy nonlinearity from the two effects separately.
1. Quenching effect. Some fraction of the deposited energy is transferred into heat and only the remaining quenched Fig. 1 The quenching induced nonlinearity as a function of the deposited energy of electrons. Solid lines are derived from the numerical integral with the ESTAR dE/dx and dashed lines are Geant4based simulation results with the production cut of 0.1mm. The red, blue and gray lines are nonlinearity curves with k B = 5.5 × 10 −3 g/cm 2 /MeV, k B = 6.5 × 10 −3 g/cm 2 /MeV (inherent value in simulation) and k B = 7.5 × 10 −3 g/cm 2 /MeV, respectively. The kinks at 80 keV in the simulation curves are caused by the production cut setting in Geant4 energy E q is visible by generating photons [28]. The quenching induced nonlinearity f q is defined as f q = E q /E, and the expected value of visible energy from scintillation photons (E s ) can be expressed by introducing the scintillation light yield Y s , which is the number of scintillation PEs per unit visible energy. As a comparison, Y is the total PEs per unit visible energy: where k B is the Birks' coefficient and larger k B yields larger quenching. Implementations of the quenching effect in Geant4 are based on the Birks' model in discrete steps. Besides, we used a numerical integral of Birks' model with stopping power data from ESTAR [29] for cross-check. The different calculation results of f q (E) are shown in Fig. 1. The curves have apparent discrepancies in the low energy region, while similar best-fit curves of the total nonlinearity and resolution have been obtained in the follow-up analysis, which validates the feasibility and robustness of our model. Hereafter, the results are based on the numerical integral curves. 2. Cherenkov radiation. Cherenkov effect heavily relies on LS refractive index and the absorption re-emission probability, especially in the ultraviolet region. The Frank-Tamm formula [30] is implemented in Geant4 for Cherenkov intensity calculation. The shape of the Cherenkov photoelectron yield as a function of the electron deposited energy f C (E) can be obtained from Geant4 simulation and is shown in Fig. 2. To avoid the heavy dependency on the LS optical properties in simulation, an empirical parametrization of the Fig. 2 The Cherenkov PEs yield f C for e − from simulation (blue markers). The parametrization is shown as the red line Cherenkov light is constructed as, where E = E − E 0 , and E 0 is the energy threshold of the Cherenkov effect. The fitting curve is also shown in Fig. 2, and there is a good agreement between the simulation and our model.
In summary, the expected visible energy is the sum of the scintillation part and the Cherenkov part, as shown in Eq. (14), where the parameter p are the four parameters involved in Eq. (13). A six-parameter energy nonlinearity model has been established for e − based on the physical origins.

Resolution
Similar to the energy nonlinearity model, the contributions to the energy resolution from scintillation photons and Cherenkov photons have been studied separately. Considering the correlation between E s and E C , the total fluctuation of the visible energy can be decomposed as: where σ s and σ C represent the standard deviation of E s and E C , and Cov [E s , Fig. 3a the three parts in Eq. (15) are shown as the three shadow regions. The scintillation process is close to Poissonian as expected (blue inverted triangles). The Cherenkov photons have a rather large smearing due to particle track length fluctuation (red circles). A longer track length leads to more Cherenkov photons and relatively smaller dE/dx, which induces a smaller quenching effect. Thus, a positive correlation between E s and E C degrades the energy resolution furthermore (gray region). The energy resolution for LS detectors is sensitive to the photon composition, and a higher Cherenkov ratio yields worse resolution with the same photon statistics.
Due to the strong correlation, decoupling of the three terms in Eq. (15) is difficult without separate benchmark measurements. Instead, the energy resolution for LS detectors is often empirically modeled as where the three terms in Eq. (16) [14,16,17], while it is released as a free parameter in this analysis for more flexible description extending to the higher energy region. Suppose that PMT dark noises in the readout window obey Poisson statistics, the noise term is particle independent and can be estimated as with the total number of PMTs N PMT , average dark noise rate R DN and readout time window length T . There will be a trade-off between less dark noises in shorter windows and more signal photons in longer windows during the optimization of energy reconstruction, which requires a separate study. In the following, the R noise term is neglected as no PMT noise is considered. But it is straightforward to include dark noise in the future by introducing the c-related term into the model. Figure 3b shows the resolution of e − from simulation, together with the parametrization based on Eq. (16).

Energy response model of γ
The connections between energy resolution of γ and e − have been constructed in Sect. 3. In Eqs. (9) and (10), the calculation of energy response for mono-energetic γ s requires applying e − nonlinearity (Eq. (14)) and resolution (Eq. (16)) on the primary e ± collections of all deposition modes. To avoid complex modeling, the method adopted here is to sample deposition modes from simulation based on a certain physical model. With sufficient sampling statistics, the calculation results are supposed to converge to the true values.
To test the dependency on physical models of energy deposition of γ s, two common low-energy electromagnetic models in Geant4, namely Livermore and Penelope, have been compared. It is claimed by the Geant4 group that these two models can provide reliable results covering electrons and photons physics from 250 eV to 1 GeV. Taking 8 MeV γ s as an example, consistent results of the energy and multiplicity distributions of primary e ± have been obtained between the where L represents the Livermore model and P represents the Penelope model, and displayed in the upper panels two models in Fig. 4. Besides, the interaction cross sections for γ energy deposition are well theoretically calculated and measured in Ref. [31]. An algorithmic calculation of primary e ± distributions has been developed as in Ref. [32] and has excellent agreements with Geant4 simulation. Based on the above cross-validations on the γ s interactions in the reactor antineutrino energy range, we choose the Livermore model in the following simulation.
Decomposition of γ energy resolution as Eq. (11) is shown in Fig. 5. The energy resolution of e − is also displayed as the dark green dotted line for comparison. The excess of γ resolution compared with e − mainly comes from the nonlinearity induced smearing term σ 2 nonl , which contributes more than 20% in the energy range below 10 MeV. The remaining term σ 2 ave (red cross markers) is slightly smaller than e − resolution with the same visible energy, because lower energy e ± in γ energy deposition has fewer Cherenkov photons.

Energy response model of e +
Precision measurement of reactor neutrino oscillations requires a deep understanding of the energy resolution of e + . Meanwhile, knowledge of the individual contributors to the resolution will provide guidelines for future improvement. As discussed in Sect. 3, the energy response of e + in LS can be conveniently expressed by the energy response of e − and calibration source 68 Ge (Eqs. (4) and (5)). Therefore, by in red cross markers. The e − energy resolution curve is shown as the dark green dotted line for comparison. The worse resolution for γ than e − comes from σ 2 nonl = σ 2 − σ 2 ave (gray region) which is caused by the dispersion among deposition modes. The remaining σ 2 ave is slightly smaller than σ e − 2 due to less Cherenkov photons. Upper panel: the ratio of σ 2 nonl /σ 2 , which is larger than 20% for all energies below 10 MeV adding calibration data of 68 Ge upon Eqs. (14) and (16), the detailed expressions of e + energy resolution can be written as, where E 68 Ge vis and σ 68 Ge represent the expected value and standard deviation of the calibrated visible energy spectrum of 68 Ge.

Model inputs
The radioactive gamma sources are the most important ones for LS detector calibration. In our Monte Carlo simulation, the gamma sources from Ref. [14] are used, as shown in Table 1. For simplicity, all γ sources are considered as bare sources without enclosures and simulated at the detector center. To reach 0.01% statistical uncertainty of the energy peaks, 5 × 10 4 events for each source are generated. Besides the regular gamma calibration sources, some isotopes with continuous β spectra induced by energetic cosmic muons can also serve as the model inputs. One important source is 12 B, which decays with a Q value about 13.4MeV and lifetime around 29ms. The visible energy spectrum of 12 B spans from 0MeV to around 14MeV. Another useful sample is Michel e − generated by the decay of stopped cosmic muons. Michel e − has a wide energy spectrum with the cut-off energy at 52.8MeV. To achieve relatively low statistical uncertainties, the 12 B samples and Michel e − events are assumed as 100k and 400k respectively in this study. Taking the JUNO detector as an instance, where the muon rate is estimated as 3.6Hz [33], the charge ratio is approximated as the value at the sea level [34] and the muon stopping rate is around 4%, the above statistics is equivalent to three months running period approximately. All the isotope spectra are simulated at the detector center without considering the detector geometric effects.

Statistical approach via χ 2 minimization
The inputs for model fitting are the energy spectra of the simulated sources mentioned in Sect. 5.1. The whole energy spectra of γ sources are utilized. For 12 B spectrum, an energy window from 3 MeV to 12 MeV is chosen, both for suppression of the natural radioactivity and for covering the high energy range of reactor antineutrinos.
To determine the parameters in our energy model, a χ 2 statistics is defined as Eq. (21) by comparing predictions (P) and measurements (M) of energy spectra for all input sources, , k B , p, a, b, n The first term in Eq. (21) is for γ s. The summation of i enumerates all calibration γ sources, and the summation of j i enumerates all energy bins for the i-th source. M j i and σ j i are the count of j-th energy bin and its statistical uncertainty. Neglecting the energy leakage in large LS detectors via a fiducial volume, the ideal energy spectra of γ sources are approximately expected as normal distributions N E γ vis , σ γ . Thus, the count of j-th energy bin (P j i ) can be predicted according  (Eq. (14)) on the theoretical energy spectrum, and then smearing the spectrum according to the e − energy resolution function (Eq. (16)). The Michel e − spectrum is not used as a fitting input considering uncertainties, more details are discussed in Sect. 6.1.
The description of e − nonlinearity introduces scintillation PEs yield Y s , Birks' coefficient k B in the quenching effect Fig. 6 Comparison of γ energy response and 12 B energy spectrum between the simulated results (red markers) and the best-fit model (black lines) (a gamma energy nonlinearity, b energy resolution, c 12 B visible energy spectrum). The relative residual bias, defined as (F −T )/T (T is the simulation value, F is the best-fit model prediction), is within 0.05% and 1% for nonlinearity and resolution, respectively. For nonlinearity of each multi-γ source, like 60 Co, 68 Ge and n 12 C, the average energy of all its γ s is displayed. The resolution of multi-γ sources deviates from that of the single-γ source, mainly due to different Cherenkov PEs ratios and 4 parameters in the Cherenkov effect f C . Besides, the e − resolution model requires 3 more parameters a, b and n. Therefore, total 9 physical parameters are involved in the model fitting. The χ 2 in Eq. (21) is minimized with TMinuit.

Fitting results
The fitting procedures have been proceeded in the four simulation configurations mentioned in Sect. 2, and the detector with the energy resolution around 2.9% are demonstrated in this section, where the value of Y is determined as 1400 photoelectrons/MeV. The fitted nonlinearity and resolution values for all these calibration γ s are compared with those from simulation in Fig. 6a, b. A residual bias for nonlinearity and resolution of γ s of less than 0.1% and 1%, respectively, can be achieved. The nonlinearity-only calibration in Ref. [14] based on a 4-parameter empirical formula is quoted as a comparison. The 0.1% residual bias level of the nonlinearity in Fig. 6a are comparable with the 0.1% residual bias level there. Among all γ calibration sources, there are three multiγ sources, 60 Co, 68 Ge and n 12 C, which have slightly different behavior of energy response compared to the other single γ sources. The displayed energy E for multi-γ sources in Fig. 6a is calculated as the average value of all the γ s, while the sum of all the γ s are displayed for multi-γ sources in Fig. 6b. The 68 Ge and 60 Co sources have better energy resolution compared with single γ sources with same deposited energy, because the Cherenkov PEs are less for two softer γ s, resulting in smaller non-Poisson fluctuation. While for n-12 C, the worse resolution originates from the E vis discrepancy for the two branches (see Table 1). Furthermore, the best-fit 12 B spectrum has great consistency with the simulation data (see Fig. 6c).
The best-fit parameters are listed in Table 2, and the correlation matrix is shown in Fig. 7. The correlation matrix manifests features of the block matrix, where parameters of the scintillation nonlinearity, Cherenkov nonlinearity and energy resolution are strongly correlated internally. The nonlinearity part and the resolution part are only weakly correlated. However, around 60% correlation is found between the amplitudes of scintillation and Cherenkov components (Y s and p 0 ). Amounts of scintillation photons and Cherenkov photons are negatively correlated due to the total PEs constraints, and may have bias compared to the nominal settings. For example, the best-fit value of the Birks' coefficient k B = (5.76 ± 0.03) × 10 −3 g/cm 2 /MeV deviates from the nominal value 6.5 × 10 −3 g/cm 2 /MeV, while the total nonlinearity fitting performs well as in Fig. 6a. To disentangle the strong correlation, standalone measurements on different components are important.

e + nonlinearity and resolution from data-driven calibration
To validate the model fitting, mono-energetic e + samples are simulated at the detector center covering the reactor antineutrino energy range. According to Eqs. (19) and (20), the e + energy response can be derived from the best-fit e − energy model and the measured data of 68 Ge calibration source. The predictions of the e + nonlinearity and resolution, as well as the uncertainties, are shown in Fig. 8. To calculate the where the relative residual bias for e + nonlinearity and resolution is less than 0.1% and 2%, respectively. Again, the residual bias of the nonlinearity fitting of e + is similar with the 0.2% residual bias level in Ref. [14] and the uncertainty is relatively smaller.

Impacts of fitting data inputs
The fitting inputs for our model are γ calibration sources and cosmogenic 12 B. During the model tuning, we have found Fig. 8 Comparison of e + energy response between the prediction and the simulated data (a nonlinearity, b resolution). The relative residual bias, with the same definition in Fig. 6, is less than 0.1% and 2% for nonlinearity and resolution, respectively. The uncertainty bands are scaled for better visualization Fig. 9 Comparison of Michel e − visible energy spectrum between the model prediction and the simulated data. The relative residual bias, with the same definition in Fig. 6, is shown in the upper panel that the spectra of γ sources are critical to constrain the resolution model due to the mono-energetic nature. However, the bias and uncertainty increase in the high energy region due to lack of calibration data. Unlike the case in Ref. [16], our Monte Carlo data does not include electronics nonlinearity, thus the 12 B spectrum provides limited improvements on the nonlinearity fitting. Moreover, the 12 B spectrum is insensitive to the energy resolution either because of the continuous spectrum property. The fitting performance in the higher energy region is checked by the Michel e − samples. The caveat is that the theoretical spectrum of Michel e − could be significantly affected if muons decay in atomic-bound states [35]. Moreover, there are both μ − and μ + components in the cosmic muons, so the Michel e − spectrum might be contaminated by e + , which introduces additional uncertainties into the model prediction. Therefore, the Michel e − spectrum is not used as a fitting input. As a test, we use the ideal energy spectrum of Michel e − without considering theoretical uncertainties or the e + contamination, and compare it with the model prediction using best-fit values in Table 2. The best-fit spectrum has good agreement with the simulated data, as shown in Fig. 9.
In the energy region of reactor antineutrinos, our model enables a high-precision calibration of the energy nonlinearity and resolution of e + . While extending to higher energy, the energy resolution is less constrained due to the lack of mono-energetic sources. This also motivates the development of new calibration strategies in the specific energy region of interest.

Energy response separation among different particles
In Fig. 10, we compare the best-fit nonlinearity and resolution for γ and e − , as well as the predicted e + responses. The 68 Ge calibration source anchors the starting points of the e + nonlinearity and resolution. Some observations can be made below.
-At the same E vis , the required energy of e − is lower than that of e + , whereas the energy resolution of e − is slightly worse than that of e + , due to more Cherenkov photons and larger non-Poisson fluctuations for e − , as discussed in Sect. 4.1.2. -The nonlinearity curve of γ is below that of e − because its energy deposition is via multiple secondary e ± with smaller energies. At low energies, the deposited energy of e + is dominated by the annihilation part, so e + generates fewer photons than the single γ s with the same total energy. As energy increases, the nonlinearity curve of e + will exceed that of γ s. -Theγ resolution is notably worse than those of e − and e + .
It mainly results from the diversity of energy deposition modes as discussed in Sect. 3. -Similar to previous work [14,16], the current model is able to separate the nonlinearity curves for γ and e ± . Moreover, the new feature is the ability to obtain particledependent energy resolution curves, see Fig. 10b. Highresolution detectors like JUNO are capable of distin- Fig. 10 The best-fit nonlinearity (a) and resolution (b) for γ and e − , as well as the predicted e + responses. Shadow regions on the lines represent the error bands, which have been scaled for better visualization guishing the energy resolution between γ and e + at the reactor antineutrino energy range.

Model results on detectors with varied resolutions
By tuning the detection efficiency of the photosensors in the Monte Carlo software, it is convenient to study the model performances versus different resolutions. Four cases are simulated with the resolution values of zero kinetic energy positron being 2.9%, 3.8%, 5.5% and 7.8%, respectively.

Fig. 11
Fitting results of energy resolution for e + are displayed, and the shadow error regions on the lines have been scaled for better visualization. Calibration data of γ sources are displayed together for comparison. Four different detector configurations have been studied. As the energy resolution of detectors becomes worse, the separation ability between e + and γ decreases The same fitting procedure is applied and the results are shown in Fig. 11. The best-fit curves of e + resolution are shown as colored lines with error bars scaled for better visualization. For comparison, the resolution values of single γ sources are also superimposed. The statistics of N pe dominates the energy resolution variation among the four cases. The discrimination ability between e + and γ decreases as the energy resolution of detectors becomes worse. This can be simply understood that as the LS-induced nonlinearities among all detectors are similar, the fraction of σ 2 nonl in [σ γ ] 2 becomes smaller for detectors with worse energy resolution according to Eq. (11). Results in Fig. 11 indicate that the gamma resolution curve is incapable to approximate e + for high-resolution detectors. The proposed data-driven calibration strategy for e + resolution will strongly advance the precision measurements of reactor antineutrinos.

Summary and prospect
A comprehensive study has been carried out to construct a unified energy model for γ and e ± in LS detectors. Based on Monte Carlo studies, we have demonstrated a promising data-driven calibration approach to predict the e + nonlinearity and resolution simultaneously in the reactor antineutrino energy region. Without considering uncertainties from enclosures of γ sources, backgrounds in the 12 B spectrum and other instrumental induced effects, the relative residual bias in nonlinearity and resolution can achieve within 0.1% and 2%, respectively. For higher energy region, the good agreement with simplified Michel e − spectrum has preliminarily validated the model performance, however, further exploration of dedicated calibration strategies will be needed. The study of energy resolution decomposition will also motivate our future developments on scintillation/Cherenkov photon discrimination algorithms.
In practice, the observed N pe is often influenced by the PMT and electronics response, such as charge smearing and dark noises, and the geometric effects in LS detectors will influence the energy response such as the non-uniformity. Those instrumentally induced effects can be effectively eliminated by the reconstruction algorithms, e.g., see discussions in Refs. [14,36,37] for JUNO. For the application to a realistic LS detector, these effects that impact energy response should be considered carefully and added to our model.
The JUNO-TAO experiment [38], which is a ton-scale satellite experiment of JUNO, aims to precisely measure the reactor antineutrino spectrum. Its targeted effective energy resolution is <2% and a similar calibration strategy has been developed [39] to control the nonlinear energy scale to be <1%. Our energy model would be also applicable to JUNO-TAO, and its high-resolution data could be of great value to disentangle the Cherenkov and quenching effects in energy nonlinearity and resolution.