Calibration strategy of the PROSPECT-II detector with external and intrinsic sources

As described above, P-I calibration data can be used for an external calibration study by analyzing the response of a subset of segments, as shown in figure 3. There were four calibration campaigns throughout the operation of the P-I detector. Each calibration campaign deployed various $\gamma$-ray and neutron sources, see table 1. All $\gamma$-rays and neutron sources provided $\gamma$-rays of known energy used for calibration. In addition, the positron annihilation in the ²²Na source allows study of annihilation $\gamma$-rays transport. Neutron production from spontaneous fission in the ²⁵²Cf source serves as a source of tagged, internally produced,  2.22 MeV $\gamma$-rays from the p(n,$\gamma$)d reaction. The AmBe neutron source provides 4.4 MeV $\gamma$-rays via $\alpha$-particle capture on ⁹Be. The PROSPECT collaboration acquired the AmBe source near the end of data-taking and only used it for data-MC validation purposes in P-I analysis. Additionally, intrinsic cosmogenically produced ¹²B (predominantly $\beta$-decays with an energy up to 13.4 MeV) was used to establish detector response in the higher energy region.

During P-I calibration campaigns, the sources were deployed at various source tube locations inside the detector at seven positions 20 cm apart along the length of the segments. The center source deployment typically ran for 10 minutes mainly for energy scale analysis, while off-center deployments were usually 5 minutes for edge effect studies. The AmBe source was deployed more than 10 hours overnight to accumulate sufficient event statistics. The layout of segments and source locations evaluated here are shown in figure 4. Segments represented with white squares in figure 4(a) were non-operational at the time of calibration due to PMT voltage divider instabilities. Data from many operational segments are excluded from the calibration analysis depicted in figure 4(b) to mimic the proposed P-II external calibration scheme. Note that particle interactions in all non-operational or analysis-excluded segments are tracked in the supporting MC simulations.

Calibration events are identified via selection criteria based on PSD particle identification, the number of segments with energy depositions above a 90 keV energy threshold (event multiplicity), and spatial and temporal correlation of particles. Electron-like events from the $\gamma$-ray sources are selected within 3$\sigma$ of the mean electron-like PSD band as a function of energy [27]. A typical integral $\gamma$-ray rate of $\mathcal{O}$(1000) counts per second for the entire energy range is observed, after all selection cuts have been applied. ²⁵²Cf events are selected by identifying a correlated pair of a prompt $\gamma$-ray from ²⁵²Cf fission followed by a 2.22 MeV $\gamma$-ray from the p(n,$\gamma$)d thermal neutron capture interaction. The time coincidence cut requires the n-H capture $\gamma$-ray to be within 200 $\mu$s of the prompt $\gamma$-ray. ¹²B event tagging, via the ¹²C(n,p)¹²B interaction, requires a correlated prompt single-hit proton-like recoil event within the energy range from 0.7 MeV_ee to 10 MeV_ee, followed by a delayed signal that has electron-like PSD with energy less than 15 MeV and multiplicity $\leq{3}$. The time difference between the prompt and the delayed signals is required to be between 3 and 30 ms, and the distance separation is required to be less than 12 cm.

Birks’ law [28] is an empirical description of the light yield per path length as a function of the energy loss per path length for a particle traversing a scintillator. Birks proposed a simple expression for the light output which took into account the quenching of light output for cases of very large energy loss. Although other P-I analyses employ a second order Birks’ model, it was found that the current data are well reproduced by a first order model.

A Geant4-based Monte Carlo simulation package was developed for the P-I detector. The MC package has been extensively benchmarked and reproduces detector response well [22]. The MC simulates particles traversing the liquid scintillator, converting the deposited energy to scintillation light via a nonlinear relation. Comprehensive calibration data collected with neutron and $\gamma$-ray radioactive sources and intrinsic backgrounds, spanning the entire PROSPECT energy scale, has been utilized to tune MC simulation and accurately reproduce P-I scintillator energy response nonlinearities.

In Geant4, the energy deposition is accumulated with each simulation step (G4Step) along particle tracks. We calculate dE/dx in each G4Step and then apply effective nonlinear corrections to calculate the resulting quenched scintillation emission. The fractional conversion rate of true deposited energy in the scintillator to scintillation light is calculated and summed at each step [27]:

where $E_{scint,i}(k_{B})$ is the scintillation light produced in the scintillator at each simulation step $i$, following Birks’ empirical model of light yield, and $E_{c,i}(k_{C})$ is light yield from the Cerenkov process, and $A$ is an overall normalization factor. Both $k_{B}$ and $k_{C}$ are effective parameters that incorporate scintillator properties as well as other detector effects that produce non-proportional response. Birks’ law is parameterized as follows [28]:

where $k_{B}$ is the effective first-order Birks’s constant and $dE/dx$ is the true deposited energy in that step. Cerenkov light production, absorption, and subsequent scintillation photon re-emission in simulation step $i$ is modelled as

where the integrand is the number of Cerenkov photons emitted per unit wavelength, $\alpha$ is the fine structure constant, $z$ is the particle’s electric charge, $\beta$ is the speed of the particle, $n(\lambda)$ is the index of refraction of the medium, $E_{\lambda}$ is the energy of those Cerenkov photons, and $k_{C}$ is an effective normalization parameter that scales Cerenkov light production with respect to a default estimate based on constant scintillator refractive index assumptions for wavelengths in 200-700 nm range.

We simulated $1\times{10}^{6}$ and $4\times{10}^{5}$ events for the $\gamma$-ray sources (⁶⁰Co, ¹³⁷Cs, ²²Na) and ²⁵²Cf neutron source, respectively, and $1\times{10}^{6}$ events for cosmogenic ¹²B, for each set of nonlinearity parameters. Simulation statistics are over four times larger than the event statistics of the calibration data to enable energy scale model optimization while keeping computation time moderate. Identical event selection rules are applied to simulations, as previously described in section 2.1. MC spectra using nonlinear energy model parameter $k_{B}$ and $k_{C}$ from eq. (2.1) are subsequently rescaled by the overall scaling factor $A$. This allows the parameter $A$ to be varied without re-running the simulation. The $\chi^{2}$-based energy scale model optimization uses this simplification in the following section.

The energy response model of the P-I detector simulation is fit to the calibration data by minimizing the $\chi^{2}$ function [22]:

The $\chi^{2}$ function consists of contributions from both the energy spectra and the event multiplicity distributions from radioactive sources described in section 2.1. The $\chi^{2}$ contribution from the individual $\gamma$-ray spectra is denoted as $\chi^{2}_{\gamma}$. The event multiplicity distributions for each source are included in $\chi^{2}_{\mathrm{multi}}$. The intrinsic cosmogenic ¹²B $\beta$-spectrum contribution is denoted as $\chi^{2}_{{}^{12}\mathrm{B}}$. The $\gamma$-ray event multiplicity distribution captures the detector response to low energy deposits and threshold effects and thus plays an important role in constraining the energy model for highly-segmented inhomogeneous liquid scintillator detectors. The ¹²B $\beta$-dominated spectrum constrains the energy response model in the high energy region up to $\sim$14 MeV.

A point-wise $\chi^{2}_{\mathrm{data-MC}}$ is calculated in a three-dimensional detector response parameter space ($A$,$k_{B}$,$k_{C}$) using weighted histogram comparison as follows [29]:

where $W_{E}$ and $W_{S}$ are the integrated bin contents from experimental and simulation histograms, respectively, $w_{i}$ is the individual bin height, $\sigma_{i}$ is the individual uncertainty from each bin, $n_{\mathrm{min}}$ and $n_{\mathrm{max}}$ is the lowest and highest energy bin of the range of interest. The $\chi^{2}$ function is calculated in the energy range covering all peak areas and regions between the peaks. More specifically, all peaks in the energy spectrum from the data are identified by a peak-finding algorithm, and their lower and upper edges are located at $\pm 3\times$ energy resolution of 4.5 %$/\sqrt{E}$ around each peak energy, which is the measured energy resolution from the P-I detector [22]. The energy range of $\chi^{2}$ calculation for each type of source spectra is defined from the lower energy edge of the lowest energy peak to the upper edge of the highest energy peak. For example, figure 5 shows ²²Na calibration data in comparison to MC generated by an energy scale model corresponding to a set of nonlinearity model parameters (the best fit set is used here). For the energy distribution, the $\chi^{2}$ contribution from ²²Na associated with this point is calculated in the energy range indicated by the gray lines in the top panel, with 0.01 MeV bin size. The event multiplicity range for $\chi^{2}_{\mathrm{multi}}$ calculation is between 1 and 9, beyond which event counts are negligible, as illustrated in the bottom panel of figure 5.

The three-dimensional $\chi^{2}_{\mathrm{data-MC}}$ map in detector response parameter space ($A$,$k_{B}$,$k_{C}$) can be further simplified by optimizing the overall energy scaling factor $A$ once the other two parameters are set. This simplification does not change the fundamentals of the model optimization and significantly decreases the computation time for a point-wise simulation in axis-$A$. With $A$ decoupled from the three-dimensional response parameter space, a two-dimensional $\chi^{2}_{\mathrm{data-MC}}$ map can be easily visualized in the parameter coordinate ($k_{B}$,$k_{C}$) (figure 6). An elliptical contour with a $\Delta\chi^{2}=30$ relative to the minimum value is fitted to determine each parameter’s best fit value and uncertainty. Best fit values of $k_{B}$ and $k_{C}$ parameters are the ellipse’s center value, and their associated uncertainties are assigned as half of the projected width in each parameter axis. The best fit value of $A$ is the average of all values of $A$ within the contour, and the uncertainty is half of the $A$ variation. The choice of $\Delta\chi^{2}=30$ is intentionally conservative, chosen to be significantly larger than fluctuations due to the limited size of the simulation samples (approximately 10 units). The remaining 20 units corresponds to about 4$\sigma$ confidence intervals. Most importantly, the same $\Delta\chi^{2}=30$ is used for both internal and external energy scale model uncertainty analysis.

Following the prescription described in section 2.3, table 2 summarizes the best fit values of each parameter and associated uncertainties. While the external calibration energy model is compatible with the internal calibration for all parameters within the assigned uncertainty, it is less well constrained by approximately a factor of two. We will discuss the implications of this model parameter uncertainty in terms of systematic contribution to spectrum uncertainty in section 3.3 and how to reduce parameter uncertainties with additional high-energy $\gamma$-rays generated from an AmBe neutron source in section 3.2. Given the reduced statistics of calibration events in the external configuration, the results presented here are conservative and improvements, such as higher activity sources and longer calibration acquisition time, are expected from an actual P-II detector setup.

Using the best fit values in table 2 for internal and external calibration, both models demonstrate good agreement between data and MC in the reconstructed energy spectrum and event multiplicity, as shown in figure 7. Noticeable changes in both the spectral shape and multiplicity are apparent between internal and external configurations. The detector active volume only covers half of the solid angle of the external sources and thus is less likely to fully capture multiple $\gamma$-rays, causing the change of energy spectrum and downshift in event multiplicity shown in figure 7. nH capture $\gamma$-ray rates are reduced due to less neutron captures internally, similar to the $\gamma$-ray source. The event multiplicity of nH event will suffer a similar downshift, depending on how deep the nH reaction happens in the fiducial volume. On the other hand, ¹²B, generated intrinsically via the ¹²C(n,p)¹²B reaction in the liquid scintillator, is independent of the calibration setup. In the external calibration setup, both the $\gamma$-ray calibration sources and the 2.2 MeV $\gamma$-ray from neutron captures on hydrogen see a reduction of event statistics that could be compensated by higher activity sources in the actual P-II calibration.

The PROSPECT collaboration acquired an AmBe source near the end of the P-I data-taking period. An AmBe source is of particular relevance to IBD neutrino detectors. Neutrons and a correlated 4.4 MeV $\gamma$-ray are produced in the AmBe source via the ⁹Be($\alpha$,n)¹²C* reaction. The correlated pair of neutron and $\gamma$-ray is similar to IBD events and can be used as a tool to quantify neutron capture efficiency and neutron mobility. They also provide a 4.4 MeV energy calibration feature, further constraining $k_{B}$ and $k_{C}$ over the high-energy region.

The AmBe source was not included in the initial calibrations and was only used as a cross-check for energy scale model data-MC comparison in previously published results [22]. However, we include these data here to test the impact of higher energy calibration features on future calibration campaigns. Inclusion of the AmBe source in the calibration data analysis is accomplished via the equation,

where $\chi^{2}_{\mathrm{data-MC}}$ is defined in eq. (2.4) and the $\chi^{2}$ contributions from the AmBe 4.4 MeV $\gamma$-ray reconstructed energy spectrum and event multiplicity are added to $\chi^{2}_{\mathrm{tot}}$. The best-fit energy scale model with the AmBe source shows good agreement between the AmBe data and MC as shown in figure 8, with $\chi^{2}_{\mathrm{AmBe}}$/ndf = 250.6/245. The fit for other sources shows good agreement as well. For example, ²²Na energy spectrum contributes $\chi^{2}_{\mathrm{Na}}$/ndf = 254.7/106 with the inclusion of the AmBe source, compared to $\chi^{2}_{\mathrm{Na}}$/ndf = 254.3/106 without the AmBe source. The benefit of including an energetic calibration $\gamma$-ray is an improved constraint on the higher energy calibration. In particular, the $k_{C}$ uncertainty is reduced, as shown in table 3, where we see a 31 % improvement. The amount of the improvement is unlikely to increase in P-II calibration run, given that the AmBe data we used contains sufficient events statistics for the purpose of the study.

A key performance metric of the calibration approach is the contribution to the overall systematic uncertainty from the energy scale model. We desire that the systematic uncertainty contribution from the external calibration scheme will be less than or equal to the anticipated statistical uncertainties of P-II. To evaluate the impact of energy scale systematic uncertainty on the oscillation and neutrino spectrum physics goals of P-II, we begin by calculating a covariance matrix.

The covariance matrix is generated by comparing a set of toy model spectra to a reference IBD spectrum. Reference IBD spectra are generated by passing the Huber-Mueller model [30, 31] predicted neutrino spectrum through the detector response with the best-fit value of the parameters, and toy model spectra use Gaussian distributed parameters with sigma values assigned by the corresponding error bar. Toy models are generated with uncorrelated uncertainties. Reference IBD spectra consist of $4\times{10}^{7}$ simulated IBD interactions. Toy model spectra consist of $1\times{10}^{3}$ toy simulations, each of $4\times{10}^{6}$ simulated IBDs. Approximately 10 % of simulated IBDs pass the selection cuts in the fiducial volume.

The diagonal elements of the reduced uncertainty covariance matrix provides a measure of the relative uncertainty in each energy bin. In figure 9, the relative uncertainty is plotted between 0.8 MeV and 7.2 MeV visible energy where most reactor antineutrino events from the IBD process fall. The external calibration uncertainty with the AmBe source (red line) shows a similar level of performance to the P-I internal calibration without the AmBe source (black line). P-II is expected to observe more than $3.5\times{10}^{5}$ IBD events during two years’ deployment, with an estimated signal-to-background ratio of 4.3. Estimated P-II statistical uncertainty with background subtraction is plotted in blue, slightly above the uncertainty of the external calibration. Further reduction of the external calibration systematic uncertainties can be achieved with higher activity sources and longer calibration acquisition time, resulting in smaller systematic uncertainties than the estimated two years’ P-II neutrino statistics with background subtraction. This study suggests that the energy scale model systematic uncertainties will not be dominant for the improved P-II detector [23].

We have demonstrated an external calibration scheme viable for the upgraded P-II detector. Here we continue to briefly address some specific questions related to the design of P-II: the detector-source stand-off distance, increased backgrounds in the location radioactive sources will be placed, and the impact of inter-segment cross-talk. The detector’s size and geometry, and the dimensions of each segment can also impact the performance of the external calibration technique. However, we do not expect significant impact on calibration due to the small change in the segment dimensions.

Having radioactive sources deployed externally introduces a stand-off distance between the source and the boundary of the active liquid scintillator volume, which can be treated as an additional fitting parameter if the transverse source location is not precisely determined. The preliminary P-II design introduces a stand-off distance of approximately 3 cm. We evaluate the precision requirement of stand-off distance by comparing MC spectra at different stand-off distances to the nominal position spectrum. In figure 10, MC simulations of additional source stand-off distances of 1.4 mm further away from the nominal position show negligible differences ($\chi^{2}/$ndf= 250.6/245) in the spectrum. We further increase the offset to 2, 5 and 10 times the unit of 1.4 mm, which corresponds to 1% of the segment width. At those offset positions, apparent spectral deviations are visible with increasing $\chi^{2}/$ndf values, 1.25, 1.74 and 4.19 respectively. A positioning requirement of $\pm$1.4 mm or $\pm$1% of the segment width is achievable through careful engineering and will have minimal impact on the desired energy scale uncertainty.

By definition, an external calibration scheme will be primarily using data from segments on edge, and these are likely to experience higher $\gamma$-ray background rates than those internal to the detector. In the P-II detector, an externally deployed source sits outside the active liquid scintillator but inside P-II exterior shielding packages. The P-I detector observed increased $\gamma$-ray background rates by up to a factor of four in the outer layers of the detector segments compared to the detector interior. It is these outer layers that are primarily calibrated by external sources. By artificially increasing the background rates by a factor of four, we find these increased $\gamma$-ray backgrounds translate to increased uncertainties in the reconstructed energy spectrum and event multiplicity by $\sim$3-4 %. These increased uncertainties are negligible in the $\chi^{2}$-optimization process and could be mitigated by using a higher activity source or better P-II exterior shielding package.

The P-II detector separates the PMTs from the active liquid scintillator volume by an acrylic window. This design change introduces $\sim$0.5 % cross-talk between adjacent segments and $\sim$0.01 % cross-talk between diagonal segments, estimated with optical MC simulations and benchtop testing. To investigate the impact of such cross-talk on calibration performance, various levels of artificial cross-talk are added to the P-I calibration data and simulation, and the energy scale analysis is repeated. The artificial cross-talk takes into account the non-uniform longitudinal dependence along the segment based on the optical simulations. The magnitude of the cross-talk increases non-linearly by about 80 % as the light source moves from the center to the end of a segment. Using an internal calibration setup (similar results for external calibration setup), the effective energy scale model with 2 % cross-talk shows no significant difference in the energy response parameters, as shown in table 4.

A specific concern is that the level of cross-talk will not be exactly known from simulations and in-situ measurements. To investigate the precision required, different levels of cross-talk are applied to the data and MC. We obtain a similar effective energy scale model using an artificial cross-talk of 2 % applied to data and 2.5 % applied to simulation. In this case, a small increase in the $A$ parameter is observed, to compensate for more light lost to neighboring segments (where it falls below the per-segment detection threshold) in the simulation than the data. This indicates a tolerance of up to 0.5 % for cross-talk measurement precision. The energy scale models with various levels of cross-talk are summarized in table 4; they all show good agreement in energy spectrum and event multiplicity.