Emulator-based Bayesian inference on non-proportional scintillation models by compton-edge probing

Domain: earth_science
Published: 2023-11-28
Editorial outcome: minor_revision

Abstract
Scintillator detector response modeling has become an essential tool in various research fields such as particle and nuclear physics, astronomy or geophysics. Yet, due to the system complexity and the requirement for accurate electron response measurements, model inference and calibration remains a challenge. Here, we propose Compton edge probing to perform non-proportional scintillation model (NPSM) inference for inorganic scintillators. We use laboratory-based gamma-ray radiation measurements with a NaI(Tl) scintillator to perform Bayesian inference on a NPSM. Further, we apply machine learning to emulate the detector response obtained by Monte Carlo simulations. We show that the proposed methodology successfully constrains the NPSM and hereby quantifies the intrinsic resolution. Moreover, using the trained emulators, we can predict the spectral Compton edge dynamics as a function of the parameterized scintillation mechanisms. The presented framework offers a simple way to infer NPSMs for any inorganic scintillator without the need for additional electron response measurements.

Introduction
Inorganic scintillation detectors are a prevalent tool to measure ionizing radiation in various research fields such as nuclear and particle physics, astronomy, or planetary science1–7. Other applications include radiation protection, medical diagnostics, and homeland security8,9. In almost all applications, the measured signal needs to be deconvolved to infer the properties of interest, e.g. the flux from a gamma-ray burst or the elemental composition on a comet. This deconvolution requires accurate detector response models and consequently detailed knowledge about the scintillation mechanisms themselves.

Detector response models can either be derived empirically by radiation measurements or numerically using Monte Carlo simulations10. Regarding the numerical derivation, the most common approach to simulate the detector response is to use a proportional energy deposition model. In this model, the scintillation light yield is assumed to be proportional to the deposited energy6,11. Consequently, the detector response characterization is reduced to a comparably simple energy deposition problem, which can be solved by any standard multi-purpose Monte Carlo code.

However, thanks to the development of the Compton coincidence measurement technique12, recent studies could conclusively confirm the conjecture reported in earlier investigations13–15 that not only organic but also inorganic scintillators exhibit a pronounced non-proportional relation between the deposited energy and the scintillation light yield16–18. The origin of this scintillation non-proportionality seems to be linked to the intrinsic scintillation response to electrons and the different mechanisms associated with the creation and transport of excitation carriers in the scintillation crystal19,20. Nonetheless, our understanding of these phenomena is still far from complete and, thanks to the advent of advanced experimental techniques and the development of new scintillator materials, interest in scintillation physics has steadily grown over the past years16–24.

Regarding the detector response modeling, the scintillation non-proportionality has two major implications. First, it leads to an intrinsic spectral broadening and thereby sets a lower limit on the spectral resolution achievable with the corresponding scintillator1,25–28. Second, various studies stated the conjecture that specific spectral features such as the Compton edges are shifted and distorted as a result of the non-proportional scintillation response1,14,15,29,30. Furthermore, additional studies revealed a complex dependence of the scintillation non-proportionality on various scintillator properties including the activator concentration, the temperature, and the crystal size, among others1,21,22,25,28,31–34.

Based on these findings, we conclude that non-proportional scintillation models (NPSM) should be included in the detector response simulations to prevent systematic errors in the predicted spectral response. Non-proportional effects are known to increase with increasing crystal size25,28,31. NPSMs are therefore particularly relevant for scintillators with large crystal volumes, e.g. in dark matter research, total absorption spectroscopy or remote sensing1–7,30. In addition, especially due to the sensitivity of the activator concentration and impurities34, NPSMs need to be calibrated for each individual detector system. In case the scintillator properties change after detector deployment, e.g. due to radiation damage or temperature changes in space, this calibration should be repeated regularly.

Currently, K-dip spectroscopy, the already mentioned Compton coincidence technique as well as electron beam measurements are the only available methods to calibrate NPSM12,35–38. Moreover, only a very limited number of laboratories are able to perform these measurements. Therefore, these methods are not readily available for extensive calibration campaigns of custom detectors, e.g. large satellite probes or scintillators for dark matter research. Additionally, they cannot be applied during detector deployment, which, as discussed above, might be important for certain applications such as deep space missions.

In this study, we propose Compton edge probing together with Bayesian inversion to infer and calibrate NPSMs. This approach is motivated by the already mentioned conjecture, that the Compton edge shifts as a result of the scintillation non-proportionality1,14,15,29,30. We obtain the spectral Compton edge data by gamma-ray spectrometry using a NaI(Tl) scintillator and calibrated radionuclide sources for photon irradiations under laboratory conditions. We apply Bayesian inversion with state-of-the-art Markov-Chain Monte Carlo algorithms39 to perform the NPSM inference with the gamma-ray spectral data. In contrast to traditional frequentist methods or simple data-driven optimization algorithms, a Bayesian approach offers a natural, consistent, and transparent way of combining prior information with empirical data to infer scientific model properties using a solid decision theory framework40–42. We simulate the detector response using a multi-purpose Monte Carlo radiation transport code and solve the required inverse problem by using machine learning trained polynomial chaos expansion (PCE) surrogate models to emulate the simulated detector response43. This approach offers not only a simple way to calibrate NPSMs with minimal effort—especially during the detector deployment—but it also allows additional insights into the non-proportional scintillation physics without the need for additional electron response measurements.

Results
Compton edge probing
To obtain the spectral Compton edge data, we performed gamma-ray spectrometry under controlled laboratory conditions30. The adopted spectrometer consisted of four 10.2 cm × 10.2 cm × 40.6 cm prismatic NaI(Tl) scintillation crystals with individual read-out. We used seven different calibrated radionuclide sources (57Co, 60Co, 88Y, 109Cd, 133Ba, 137Cs, and 152Eu) for the radiation measurements. However, only 60Co, 88Y, and 137Cs could be used for Compton edge probing. For the remaining sources, the Compton edges were obscured by additional full energy peaks (FEPs) and associated Compton continua. We used those remaining sources for energy and resolution calibrations. A schematic depiction of the measurement setup is shown in Fig. 1a.

Forward modeling
We simulated the detector response for the performed radiation measurements using the multi-purpose Monte Carlo code FLUKA44. The performed simulations feature fully coupled photon, electron, and positron radiation transport for our source-detector configuration with a lower kinetic energy threshold of 1 keV. As shown in Fig. 1a, the applied mass model includes all relevant detector and source components in high detail. On the other hand, the laboratory room together with additional instruments and equipment are modeled in less detail. For this simplification, care was taken to preserve the overall opacity as well as the mass density.

We used a mechanistic model recently published by Payne and his co-workers to include the non-proportional scintillation physics in our simulations17,18,22. In general, the sequence of scintillation processes in inorganic scintillators can be qualitatively divided into five steps20,45,46. After interaction of the ionizing radiation with the scintillator, the emitted high-energetic electrons are relaxed by the production of numerous secondary electrons, phonons, and plasmons. The low energetic secondary electrons are then thermalized by a phonon coupling mechanism producing excitation carriers, i.e. electron–hole pairs (e−/h) and excitons. These excitation carriers are then transferred to the luminescent centers within the scintillation crystal, where they recombine and induce radiative relaxation of the excited luminescent centers producing scintillation photons. The first two processes, i.e. the interaction of the ionizing radiation with the scintillator as well as the e−–e− relaxation, are explicitly simulated by the Monte Carlo code. The creation and migration of the excitation carriers on the other hand is accounted for by Payne’s mechanistic model.

In this mechanistic model it is assumed that only excitons are capable of radiatively recombine at the luminescent centers. Consequently, e−/h pairs need to convert to excitons by the classic Onsager mechanism47 in order to contribute to the scintillation emission. In addition, creation and migration of the excitation carriers compete with several quenching phenomena. The quenching mechanisms considered in Payne’s model are the trapping of e−/h pairs at point defects20,22 as well as exciton–exciton annihilation described by the Birks mechanism48.

Using this NPSM, the non-proportional light yield L as a function of the differential energy loss dE per differential path length ds for electrons is given by22
L(dE/ds) = 1/ηe=h exp(dE/ds|Ons) exp(dE/ds|Trap) / (1 + dE/ds/dE/ds|Birks) (1)

where ηe/h, dE/ds|Ons, dE/ds|Trap and dE/ds|Birks are the model parameters characterizing the fraction of excitation carriers, which are created as e−/h pairs at the thermalization phase, as well as the stopping power related to the Onsager, trapping and Birks mechanisms, respectively. As a result, all the parameters of the NPSM reflect physical processes after thermalization of the secondary particles, i.e. generation and transport of excitation carriers. Consequently, these processes and thereby also the corresponding parameters can be regarded as statistically independent with respect to the energy of the secondary particles. From a physics perspective, it is important to note that the Onsager and trapping mechanisms are coupled in a nonlinear way, whereas the Birks mechanism can be regarded as independent of the other mechanisms. As discussed in detail by Vasil’ev and Gektin20, we may therefore interpret the trapping of e−/h pairs as a screening mechanism on the Onsager term in Eq. (1). A scheme highlighting the individual scintillation processes included in the present study is presented in the Fig. 1c–e.

Bayesian inversion
We applied Bayesian inversion using Markov Chain Monte Carlo39 to infer the NPSM parameters as well as to predict spectral and resolution scintillator properties from the measured Compton edge spectra and our forward model. To account for the sensitivity of the NPSMs on the activator concentration and other scintillation crystal-specific properties, we developed two separate inversion pipelines. In the first approach, Bayesian inversion is carried out separately for each of the four crystals, using their individual pulse-height spectra. In the second pipeline, we consider all four scintillation crystals as one integrated scintillator and perform the Bayesian inversion on the combined pulse-height spectra (sum channel). Subsequently, we will refer to these two approaches as the single and sum mode inversion pipelines. For both pipelines, we performed the Bayesian inversion on the 60Co (activity A = 3.08(5) × 105 Bq) spectral dataset30 leaving the remaining measurements for validation. 60Co possesses two main photon emission lines at 1173.228(3) keV and 1332.492(4) keV with corresponding Compton edges according to the Compton scattering theory (Methods) at 963.419(3) keV and 1118.101(4) keV, respectively. However, in this study, we will focus on the lower Compton edge at 963.419(3) keV, because the upper edge is heavily obscured by the FEP at 1173.228(3) keV. Furthermore, as suggested by previous investigators18,22, we fixed the Onsager-related stopping power parameter dE/ds|Ons to 36.4 MeV cm−1 in both pipelines.

Because the high-fidelity radiation transport simulations described in the previous section are computationally intense, we emulated the detector response as a function of the NPSM parameters using a machine learning trained vector-valued PCE surrogate model43. The surrogate model has excellent evaluation speed O(10−4 s) on a local workstation compared to O(103 s) required for a single Monte Carlo simulation with sufficient precision on a computer cluster. Correspondingly, the surrogate model provides a significant acceleration of our Bayesian inversion computations, reducing their processing time by a factor of 107. Considering the minimum number of forward model evaluations needed for a single Markov chain (Methods), the evaluation time can be reduced from O(102 d) to O(1 s).

Following the sum mode inversion pipeline, we present the solution to our inversion problem as a multivariate posterior distribution estimate in Fig. 2. We find a unimodal solution with a maximum a posteriori (MAP) probability estimate given by ηe=h = 5.96 + 0.17 − 0.11 × 10−1, dE/ds|Trap = 1.46 + 0.30 − 0.13 × 101 MeV cm−1 and dE/ds|Birks = 3.22 + 0.44 − 0.46 × 102 MeV cm−1, where we used the central credible intervals with a probability mass of 95% to estimate the associated uncertainties. Combining the individual multivariate posterior distribution estimates from the single mode inversion pipeline, we obtain statistically consistent estimates, i.e. ηe=h = 5.87 + 0.22 − 0.20 × 10−1, dE/ds|Trap = 1.41 + 0.15 − 0.17 × 101 MeV cm−1 and dE/ds|Birks = 3.17 + 1.11 − 0.82 × 102 MeV cm−1.

It is worth noting that, considering the uncertainty estimates, we observe only minor differences between the different posterior point estimators for both inversion pipelines (Fig. 2 and Supplementary Figs. S11–S14). However, we find statistically significant differences between the posterior point estimators for the individual scintillation crystals (Supplementary Table S2). Furthermore, our results significantly differ from best-estimate literature values, which we obtained using linear temperature interpolation on a dataset provided by Payne and his co-workers, i.e. ηe=h = 4.53 × 10−1, dE/ds|Trap = 1.2 × 101 MeV cm−1 and dE/ds|Birks = 1.853 × 102 MeV cm−1 for an ambient temperature of 18.8 °C22.

Compton edge predictions
We can use the trained PCE surrogate models to predict the spectral Compton edge as a function of the NPSM parameters and consequently the parameterized scintillation and quenching phenomena. In Fig. 3a–c, we present the spectral response of the PCE surrogate model for the sum channel as a function of the Birks-related stopping power parameter dE/ds|Birks, the free carrier fraction ηe/h and the trapping-related stopping power parameter dE/ds|Trap. We observe a shift of the Compton edge towards smaller spectral energies for an increase in dE/ds|Birks and ηe/h as well as a decrease in dE/ds|Trap.

We leveraged the analytical relation between the polynomial chaos expansion and the Hoeffding-Sobol decomposition49 to perform a global sensitivity analysis of the NPSM. Using the sum mode inversion pipeline, we present total Sobol indices ST for the model parameters dE/ds|Birks, ηe/h and dE/ds|Trap in Fig. 3e. We find that the total Sobol indices can be ordered as ST(ηe/h) > ST(dE/ds|Birks) > ST(dE/ds|Trap) over the entire spectral Compton edge domain indicating a corresponding contribution to the total model response variance. We get consistent results for a Hoeffding-Sobol sensitivity analysis of the individual scintillation crystals (Supplementary Fig. S17).

In addition, we can also predict the spectral Compton edge using the prior and posterior predictive density estimates obtained by the two inversion pipelines. A comparison of these densities for the sum mode inversion pipeline indicates that our methodology successfully constrains the adopted NPSM (Fig. 3d). However, we find also some model discrepancies, especially around the Compton continuum at the very low end of the investigated spectral range (<920 keV). We get consistent results using the single-mode inversion pipeline (Supplementary Fig. S15). Furthermore, by comparing the posterior Compton edge predictions for the sum channel, we find no statistically significant difference between the two inversion pipelines (Supplementary Fig. S16). From a modeling perspective, it is interesting to add that we observe no significant difference for Compton edge predictions using the various point estimators discussed in the previous section for both inversion pipelines.

Intrinsic resolution
As already mentioned in the introduction, the scintillation non-proportionality not only distorts the spectral features in the pulse-height spectra but deteriorates also the spectral resolution of a scintillator detector. This contribution to the overall resolution due to the scintillation non-proportionality will be referred to as intrinsic resolution σintr in accordance with previous studies25–28,50–52. The intrinsic resolution is of great importance for two key reasons.

First, it sets a fundamental lower limit on the achievable spectral resolution for a given scintillator material, making it a crucial factor in the development of new scintillators. As an example, in 1991, the scintillator Lu2SiO5 (LSO) was developed as an alternative to other available options at that time, such as Bi4Ge3O12 (BGO). However, the performance of LSO led to considerable confusion within the research community as LSO exhibits a light yield more than four times greater than that of BGO, yet their energy resolutions are comparable16. Consequently, the energy resolution for LSO was not dominated by counting statistics but some other factor. Thanks to the development of the Compton coincidence measurement technique in 199412, subsequent experimental studies have conclusively shown that the pronounced scintillation non-proportionality of the LSO scintillator was indeed the underlying factor responsible for the observed resolution degradation28,53. This example showcases the need for a better understanding and prediction of the intrinsic resolution in the development of new scintillators54.

Second, from a more technical perspective, the intrinsic resolution is a key component in the postprocessing pipeline for Monte Carlo simulations including NPSMs (Methods). In the forward model discussed before, the transport of scintillation photons, signal amplification by the photomultiplier tube and subsequent signal postprocessing in the multichannel analyzer are not included. As a result, to account for the additional resolution degradation by these processes, we need to perform a spectral broadening operation using a dedicated energy resolution model based on the measured total energy resolution as well as the intrinsic contribution1.

Since our forward model explicitly accounts for the non-proportional scintillation physics by adopting an NPSM, we can use this numerical tool not only to predict pulse-height spectra but also to characterize the intrinsic resolution. We adopted a set of multiple monoenergetic Monte Carlo simulations to quantify the intrinsic resolution for different photon energies (Methods). Using this dataset, we then trained a Gaussian process (GP) regression model to predict the intrinsic resolution characterized by the standard deviation σintr for a given photon energy Eγ. The resulting GP model predictions together with the intrinsic data are highlighted in Fig. 3f. In the same graph, we include also the empirical model to describe the overall energy resolution σtot as well as the corresponding empirical dataset30.

Comparing the intrinsic and overall spectral resolution, we find an almost constant ratio σ2intr/σ2tot ≈ 0.35 for Eγ ≳ 1500 keV. Around Eγ ≈ 440 keV, there is a pronounced peak with σ2intr/σ2tot ≈ 0.42 and for Eγ ≲ 440 keV, we observe a significant decrease in σ2intr/σ2tot with decreasing photon energy Eγ. Moreover, we find a more complex behavior in σintr for Eγ ≲ 110 keV. For 28 keV ≲ Eγ ≲ 60 keV, the K-absorption edge for iodine K[I] at Eγ = 33.1694(4) keV55 alters the resolution significantly. On the other hand, at even smaller photon energies, there is again a pronounced increase in σintr with decreasing energy compared to the mere moderate increase for 60 keV ≲ Eγ ≲ 110 keV.

Bayesian calibrated NPSM simulations
In addition to the insights into the Compton edge dynamics as well as the intrinsic resolution, the Bayesian inferred NPSM in combination with our forward model offers also the possibility to predict the full spectral detector response for new radiation sources accounting for non-proportional scintillation effects over the entire spectral range of our detector system. We used the 88Y (A = 6.83(14) × 105 Bq) and 137Cs (A = 2.266(34) × 105 Bq) radiation measurements to validate our calibrated NPSM. For the Monte Carlo simulations, we applied the posterior point estimators xMAP obtained by the sum mode inversion pipeline in combination with the intrinsic and total resolution models discussed in the previous sections.

In Fig. 4, we present the measured and simulated spectral detector response for 88Y and 137Cs together with 60Co, whose Compton edge domain was used to perform the Bayesian inversion. For the simulations, we adopted a standard proportional scintillation model as well as the Bayesian inferred NPSM presented in this study. In line with the Compton scattering theory (Supplementary Methods S1.4), we find an enhanced shift of the Compton edge towards smaller spectral energies as the photon energy increases. For all three measurements, we observe a significant improvement in the Compton edge prediction for the NPSM simulations compared to the standard proportional approach. However, there are still some discrepancies at the lower end of the Compton edge domain. Moreover, we find also some deviations between the Compton edge and the FEP for 88Y and 137Cs. It is important to note that these discrepancies are smaller or at least of similar size for the NPSM simulations compared to the proportional approach indicating that the former performs statistically significantly better over the entire spectral domain. Additional validation results for 57Co, 109Cd, 133Ba, and 152Eu together with a detailed uncertainty analysis for each source are attached in the Supplementary Information File for this study (Supplementary Figs. S18–S25).

Methods
Gamma-ray spectrometry
We performed gamma-ray spectrometric measurements in the calibration laboratory at the Paul Scherrer Institute (PSI) (inner room dimensions: 5.3 m × 4.5 m × 3 m). The adopted spectrometer consisted of four individual 10.2 cm × 10.2 cm × 40.6 cm prismatic NaI(Tl) scintillation crystals with the associated photomultiplier tubes and the electronic components embedded in a thermal-insulating and vibration-damping polyethylene foam protected by a rugged aluminum detector box (outer dimensions: 86 cm × 60 cm × 30 cm). The spectrometer features 1024 channels for an energy range between about 30 and 3000 keV together with automatic linearization of the individual scintillation crystal spectra30. We used seven different calibrated radionuclide sources (57Co, 60Co, 88Y, 109Cd, 133Ba, 137Cs, and 152Eu) from the Eckert & Ziegler Nuclitec GmbH. We inserted these sources consisting of a radionuclide carrying ion exchange sphere (diameter 1 mm) embedded in a 25 mm × 3 mm solid plastic disc into a custom low absorption source holder made out of a polylactide polymer (PLA) and placed this holder on a tripod in a fixed distance of 1 m to the detector front on the central detector x-axis. To measure the source-detector distances and to position the sources accurately, distance as well as positioning laser systems were used. A schematic depiction of the measurement setup is shown in Fig. 1a.

Between radiation measurements, background measurements were performed regularly for background correction and gain stability checks. For all measurements, the air temperature as well as the air humidity in the calibration laboratory was controlled by an air conditioning unit and logged by an external sensor. The air temperature was set at 18.8(4) °C and the relative air humidity at 42(3)%. The ambient air pressure, which was also logged by the external sensor, fluctuated around 982(5) hPa.

During measurements, additional instruments and laboratory equipment were located in the calibration laboratory, e.g. shelves, a workbench, a source scanner or a boiler as shown in Fig. 1a. The effect of these features on the detector response was carefully assessed in ref. 30.

After postprocessing the spectral raw data according to the data reduction pipelines described in ref. 30, we extracted the Compton edge spectral data from the net count rate spectra. The spectral domain of the Compton edge DE was defined as
DE := {E: ECE − 3 · σtot(ECE) ≤ E ≤ ECE + 2 · σtot(FEP)}, where E is the spectral energy, σtot the energy dependent total resolution characterized by the standard deviation30 and EFEP the FEP associated with the Compton edge ECE. Neglecting Doppler broadening and atomic shell effects, we compute ECE according to the Compton scattering theory10 as follows:
ECE := EFEP [ 1 − 1 / (1 + 2EFEP / mec2) ] (2)

where mec2 is defined as the energy equivalent electron mass. In this study, we consulted the ENDF/B-VIII.0 nuclear data file library61 for nuclear decay-related data as well as the Particle Data Group library62 for fundamental particle properties.

To investigate the sensitivity of the selected Compton edge domain DE on the Bayesian inversion results, we performed a sensitivity analysis on DE. Within the uncertainty bounds, the inversion results have proven to be insensitive to small alterations in DE (Supplementary Table S3).

It is important to note that, if not otherwise stated, uncertainties are provided as 1 standard deviation (SD) values in this study (coverage factor k = 1). For more information about the radiation measurements and adopted data reduction pipelines, the reader is referred to the dedicated study30 as well as the Supplementary Information File for this work (Supplementary Methods S1.3).

Monte Carlo simulations
We performed all simulations with the multi-purpose Monte Carlo code FLUKA version 4.2.144 together with the graphical interface FLAIR version 3.1-15.163. We used the most accurate physics settings (precisio) featuring a high-fidelity fully coupled photon, electron and positron radiation transport for our source-detector configuration. In addition, this module accounts for secondary electron production and transport, Landau fluctuations as well as X-ray fluorescence, all of which are essential for an accurate description of non-proportional scintillation effects16,18,23,52. Motivated by the range of the transported particles, lower kinetic energy transport thresholds were set to 1 keV for the scintillation crystals as well as the closest objects to the crystals, e.g. reflector, optical window and aluminum casing for the crystals. For the remaining model parts, the transport threshold was set to 10 keV to decrease the computational load while maintaining the high-fidelity transport simulation in the scintillation crystals. All simulations were performed on a local computer cluster (7 nodes with a total number of 520 cores at a nominal clock speed of 2.6 GHz) at the Paul Scherrer Institute utilizing parallel computing.

We scored the energy deposition events in the scintillation crystals individually on an event-by-event basis using the custom user routine usrbin together with the detect card. The number of primaries was set to 107 for all simulations, which guarantees a maximum relative statistical standard deviation σstat,sim,k/csim,k < 1% and a maximum relative variance of the sample variance VOV < 0.01% for all detector channels k. More details on the simulation settings as well as on the postprocessing of the energy deposition data can be found in ref. 30.

To implement the NPSM described by Eq. (1), we developed an additional user routine comscw. Similar to refs. 1,64, we weight each individual energy deposition event in the scintillator, point-like or along the charged particle track, by the scintillation light yield given in Eq. (1) (Supplementary Algorithm S1). The resulting simulated response is then rescaled to match the energy calibration models derived in ref. 30. Using our methodology, we get simulated pulse-height spectra that incorporate non-proportional effects across the entire spectral range of our detector system.

Surrogate modeling
We applied a custom machine learning trained vector-valued polynomial chaos expansion (PCE) surrogate model to emulate the spectral Compton edge detector response over DE for both, the individual scintillation crystals as well as the sum channel. PCE models are ideal candidates to emulate expensive-to-evaluate vector-valued computational models43. As shown in refs. 65–67, any function Y = M(X) with the random input vector X ∈ RM×1 and random response vector Y ∈ RN×1 can be expanded as a so-called polynomial chaos expansion provided that E[|Y|2] < ∞:
Y = M(X) = ∑α∈NM aα Ψα(X) (3)

where aα := [a1,α, ..., aN,α]> ∈ RN×1 are the deterministic expansion coefficients, α := [α1, ..., αM]> ∈ NM are the multi-indices and Ψα(X) are the multivariate orthonormal polynomials corresponding to the probability measure of X, i.e. E[Ψα, Ψβ] = δα,β. To compute the PCE coefficient matrix A, we combined the PCE model with principal component analysis (PCA) allowing us to characterize the main spectral Compton edge features of the response by means of a small number N0 of output variables compared to the original number N of spectral variables, i.e. N0 ≪ N43. Similar to ref. 68, we computed the emulated computational model response M^PCE(X) in matrix form as
Y ≈ M^PCE(X) = μY + diag(σY) Φ0 A Ψ(X) (4)

with μY and σY being the mean and standard deviation of the random vector Y and Φ0 the matrix containing the retained eigenvectors from the PCA, i.e. Φ0 := [ϕ1, ..., ϕN0] ∈ RN×N0. On the other hand, the multivariate orthonormal polynomials and corresponding PCE coefficients, respectively. The union set A* := ∪N0 A includes the finite sets of multi-indices A for the N0 output variables following a specific truncation scheme.

We used a Latin hypercube experimental design X ∈ RM×K69,70 with K = 200 instances sampled from a probabilistic model, which itself is defined by the model parameter priors described in the next subsection. The model response Y ∈ RN×K for this design was then evaluated using the forward model described in the previous subsection. We adopted the q-norm truncation scheme with q = 1, i.e. {α ∈ NM : ∑Mi=1 αi ≤ p} for the associated polynomial degree and the q-norm, respectively. To compute the PCE coefficient matrix A, we applied adaptive least angle regression71 and optimized the hyperparameters p:={1, 2, …, 7} and q:={0.5, 0.6, …, 1} using machine learning with a holdout partition of 80% and 20% for the training and test set, respectively. For the PCA truncation, we adopted a relative PCA-induced error εPCA of 0.1%, i.e. N0 := min{S ∈ {1, ..., N}: ∑Sj=1 λj / ∑Nj=1 λj ≥ 1 − εPCA} with λ being the eigenvalues from the PCA. The resulting generalization error of the surrogate models, characterized by the relative mean squared error over the test sets, are <1% and <2% for the sum channel and the individual scintillation crystals, respectively. All PCE computations were performed with the UQLab code72 in combination with custom scripts to perform the PCA. More information about the PCE-PCA models as well as the PCE-PCA-based Sobol indices including detailed derivations are included in the Supplementary Information File for this study (Supplementary Methods S1.1–S1.2).

Bayesian inference
Following the Bayesian framework40, we approximate the measured spectral detector response y ∈ RN×1 with a probabilistic model combining the forward model M(xM) and model parameters xM ∈ RMM×1 with an additive discrepancy term ε, i.e. y := M(xM) + ε. For the discrepancy term ε, which characterizes the measurement noise and prediction error, we assume a Gaussian model π(ε|σ2ε) = N(ε|0, σ2ε IN) with unknown discrepancy variance σ2ε. On the other hand, as discussed in the previous subsection, we emulate the forward model M(xM) with a PCE surrogate model M^PCE(xM). Consequently, we can compute the likelihood function as follows:
π(y|x) = N(y|M^PCE(xM), σ2ε IN) (5)

with x := [xM, σ2ε]> and xM := [dE/ds|Birks, ηe/h, dE/ds|Trap]> . In combination with the prior density π(x), we can then compute the posterior distribution using Bayes’ theorem42:
π(x|y) = π(y|x) π(x) / ∫D π(y|x) π(x) dx (6)

where we assume independent marginal priors, i.e. π(x) = ∏Mi=1 π(xi). We defined the marginal priors based on the previous studies17,18,22. It should be emphasized that we applied the sum mode inversion pipeline first followed by the single mode inversion pipeline. In accordance with Bayes’ theorem42, we therefore incorporate the results obtained by the sum mode inversion pipeline into the marginal priors used for the single mode inversion pipeline. A full list of all adopted marginal priors for both pipelines is attached in the Supplementary Information File for this study (Supplementary Table S1). Using the prior and posterior distributions, we can then also make predictions on future model response measurements y* leveraging the prior and posterior predictive densities:
π(y*) = ∫D π(y*|x) π(x) dx (7a)
π(y*|y) = ∫D π(y*|x) π(x|y) dx (7b)

All Bayesian computations were performed with the UQLab code72. We applied an affine invariant ensemble algorithm39 to perform Markov Chain Monte Carlo (MCMC) and thereby estimate the posterior distribution π(x|y). We used 10 parallel chains with 2 × 104 MCMC iterations per chain together with a 50% burn-in. The convergence of the MCMC chains was monitored using the scale reduction factor R^ < 1.02 and an effective sample size ESS ≫ 400 for all performed MCMC simulations. Additional trace and convergence plots for the individual parameters x and point estimators (Supplementary Figs. S1–S10), a full list of the Bayesian inversion results (Supplementary Table S2) as well as a sensitivity analysis on the adopted Compton edge domain (Supplementary Table S3) can be found in the attached Supplementary Information File for this study.

Resolution modeling
In this last section, we will discuss the derivation of the energy resolution models adopted in this study. We start with the model to characterize the overall or total energy resolution σtot for our detector system and describe in a second step the derivation of the intrinsic resolution model σintr. It is important to note that in contrast to σintr, we provide here only a short summary of the key aspects involved in σtot. The entire postprocessing pipeline to derive σtot was already thoroughly discussed in a previous study30. For more details, we kindly refer the reader to the dedicated study.

For each scintillation crystal, we quantified σtot by characterizing the spectral dispersion of measured FEPs associated with known photon emission lines from specific radionuclides. The corresponding pipeline can be divided into three steps. In the first step, we extracted specific spectral domains containing a singlet or multiplets of FEPs from a set of measured count rate spectra covering a spectral range between 122 and 1836 keV75. In a second step, we fitted a spectral peak model based on a sum of independent Gaussian peaks together with a numerical baseline76 to the selected singlets or multiplets using weighted non-linear least-squares (WNLLS) regression combined with the interior-reflective Newton method77. In the third step, we extracted the Gaussian standard deviation parameters from the fitted FEPs as a characteristic measure for the spectral resolution. By combining these empirical resolution values with the known emission line energies, we derived an exponential model to describe σtot as a function of the photon energy Eγ adopting again WNLLS. The resulting relative generalization error, characterized by leave-one-out cross-validation, is <0.2% for all scintillation crystals.

To derive a model for σintr, we performed in the first step additional Monte Carlo simulations for anisotropic and uniform monoenergetic photon flux of energy 10 keV ≤ Eγ ≤ 3200 keV. To account for the different spectral scales, we applied a non-uniform experimental design for the photon energy Eγ with a 2 keV spacing below 110 keV and 100 keV spacing above. Moreover, to account for the non-proportional scintillation physics, we ran all simulations with the Bayesian calibrated NPSM, i.e. the derived MAP point estimators. The mass model for those simulations features a 10.2 cm × 10.2 cm × 40.6 cm prismatic NaI(Tl) scintillation crystal embedded in a vacuum environment. In the second step, we extracted the mean light yield values from the simulated FEPs (Supplementary Algorithm S1). Similar to the measured spectra, we can then derive a simple polynomial energy calibration model using WNLLS to convert the simulated light yield to energy30. In a third step, we adopted the extracted σintr from the individual energy-calibrated FEPs to train a Gaussian Process (GP) regression model with78:
σintr(Eγ) ∼ GP(f(Eγ)> β, κ(Eγ, E0γ) + σ2GP δ(Eγ, E0γ)) (8)

where we applied a polynomial trend function of the second order, i.e. f(Eγ) := [1, Eγ, E2γ]> and β := [β0, β1, β2]>, a homoscedastic noise model with variance σ2GP and kernel κ(Eγ, E0γ) := θ exp(− 3|Eγ − E0γ| / θ) with the kernel scale θ. It is worth noting that, due to the known asymmetry in the FEPs1,25,28, we adopted numerical estimates both for the mean and standard deviation parameters associated with the individual FEPs. With the N-dimensional intrinsic dataset {Eγ, σintr}, we can then predict the intrinsic resolution σ*intr for a new set of N* photon energies E*γ using the GP posterior predictive density as follows:
π(σ*intr | E*γ, Eγ, σintr) = N(σ*intr | μGP, ΣGP) (9a)
μGP = F*> β^ + K*> K−1 σintr − F> β^ (9b)
ΣGP = K** − K*> K−1 K* + U> (F K−1 F> )−1 U (9c)
β^ = (F K−1 F> )−1 F K−1 σintr (9d)
U = F* − F K−1 K* (9e)

with the matrices F=f(Eγ) ∈ R3×N, F* = f(E*γ) ∈ R3×N*, K = κ(Eγ, Eγ) + σ2GP IN ∈ RN×N, K* = κ(Eγ, E*γ) ∈ RN×N* and K** = κ(E*γ, E*γ) ∈ RN*×N*. To account for the different spectral scales, we trained two GP models, one for 10 keV ≤ Eγ ≤ 90 keV and the other one for 90 keV ≤ Eγ ≤ 3200 keV, using the MATLAB® code. For both models, we applied 5-fold cross-validation in combination with Bayesian optimization to determine the GP hyperparameters σ2GP and θ. It is important to add that in case of the experimental design X adopted to train the surrogate model, we ran the pipeline for σintr with the corresponding set of NPSM parameters defined by X.

As discussed already in the Results section, the intrinsic resolution is also a key component in the postprocessing pipeline for Monte Carlo simulations including NPSMs. Because the Monte Carlo simulations performed for the forward model only inherently include the intrinsic resolution, we need to perform a spectral broadening operation to account for the additional energy resolution degradation due to the transport of scintillation photons, signal amplification by the photomultiplier tube and subsequent signal postprocessing in the multichannel analyzer. Similar to ref. 1, we assume statistical independence between the resolution degradation due to the scintillation non-proportionality and the aforementioned neglected processes in the Monte Carlo simulations. We can then perform the broadening operation as described in ref. 30 with an adapted dispersion σtot = √σ2tot − σ2intr. For completeness, we included this adapted dispersion model in Fig. 3f.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "As discussed in my inifial review, the work reported in this manuscript would advance the state of the art\nfor scinfillator characterizafion. Because inorganic scinfillators are used in many fields, the work would\nbe of broad interest to the scienfific community and is, therefore, suitable for publicafion in Nature\nCommunicafions. The methods are sound and ample detail is provided. The revised version compares\nresults with previous work and discusses possible applicafion of the methods beyond the scope of the\ncurrent study. In my opinion, no further revisions are needed.", "reviewer_2": "On comment #1\nThe derivafion and applicafion of the semi-analyfical model is impressive. The subsequent analysis\ndefinifively allows one to befter understand the physics behind the observed shifts in energy. I thank the\nauthors for adding the secfion S1.4 to the supplementary materials.\nOn comment #2\nI thank the authors for the clarificafion.\nOn comment #3\nI thank the authors for revealing that the sum mode inversion pipeline was the one presented in the\nsubmifted version of the manuscript and supplementary materials, and for adding informafion about the\nsingle mode inversion pipeline in the revised version of the manuscript and supplementary materials. As\nthe authors menfioned, it is interesfing to see that there are no stafisfically significant differences\nbetween the two pipelines, but sfill, understandably, observable differences between the crystals.\nOn comment #4\nI thank the authors for sharing their thoughts on this subject, and I agree with them.\nOn comment #5\nI thank the authors for the explanafion.\nOn comment #6\nI thank the authors for providing insights on the computafional resources required for solving this type of\nproblems. The fime saved by using the surrogate model is indeed substanfial, which certainly mofivates\nits development and usage.\nOn comment #7\nThat’s a good point. I thank the authors for this well thought and demonstrated jusfificafion.\nOn comment #8\nI thank the authors for the explanafion.\nOn comment #9\nThis is doubtlessly a great addifion to the manuscript. The added details on the resolufions assuredly\nallow (1) a smoother reading, (2) a befter understanding of the intrinsic resolufion and of its relevance\nfor this manuscript, and (3) to highlight the power and importance of the tool that the authors\ndeveloped.\nOn comment #10\nI thank the authors for the explanafion and for elaborafing on the energy range of applicability.\nOn comment #11\nI thank the authors for adding the informafion.\nOn comment #12\nInteresfing. It is great to see that the development of the semi-analyfical model led to a befter\nunderstanding of the Compton edge shifts.\nOn comment #13\nI thank the authors for including the pseudo-code in the supplementary materials, which answers the\nquesfions raised in this comment, and for reminding the reviewer about the ETH research data\nrepository.\nTo conclude, I would like to thank the authors one more fime for the fime and efforts that they put\ntowards the revision of the manuscript and supplementary materials. Everything was addressed far\nbeyond expectafions, which aftests the excellence of this work.\nBest of luck with this manuscript and with your future endeavors!\nBest regards.", "reviewer_3": "Please see aftached document, thanks.\nDear Researchers: Thank you for your thoughul responses. I’d like to connue the conversaon on one\nminor topic, with excerpts on my comment and your response below, with a follow-on from me.\n(Original) Referee #3, Comment #3:\nThe comment “Nonproporonality effects are known to increase with crystal size” is not correct.\n(Original) Response: We respecully disagree with the referee that the statement “Nonproporonality\neffects are known to increase with crystal size” is incorrect. As correctly stated by the referee, the\nresoluon or more specifically the intrinsic resoluon, which is a direct consequence of the scinllaon\nnonproporonality (cf. comment #9 for referee #2), increases with crystal size…\nSecond Response: I agree that a Compton scatering event followed by an addional Compton scater\nwill shi the edge to higher energy. But while nonproporonality could have some contribuon to this\neffect, it is NOT required to generate the shi. All that is needed is two Compton events adding up to\nmore than the single (maximum) Compton edge energy, and the release of a gamma photon from the\nscinllator. Consider the simplified gamma spectra below where nonproporonality is NOT included. On\nthe le is the spectrum that includes Compton-followed-by-Compton events, while on the right the\nspectrum only includes Compton singles. The energy shi for the edge from mulple Compton scatering\nevents is apparent – nonproporonality has a minor impact on it and is not required. I hope this makes\nsense… maybe think about it. I leave the final decision of how to handle this up to you.\nFigure: (le) Spectrum with double Compton events; (right) Spectrum with only single Comptons."}, "author_response": "5 Response to Referee #1\nWe would like to thank the referee once gain for the thoughtful feedback and positive evaluation.\nWe believe that the modifications suggested by you have improved the manuscript significantly\nand made it more complete. We are really grateful and appreciate your time and effort in\nreviewing our manuscript. All the best for your future.\nWe would like to thank the referee once again for the thorough review and the detailed feedback.\nWe believe that the modifications suggested by you have improved the manuscript significantly\nand made it more complete. We are really grateful and appreciate your time and effort in\nreviewing our manuscript. All the best for your future.\nWe thank the referee for pointing this out. We recognize that some of our statements in com-\nment #3 in the last peer review report lacked clarity and might have been to a certain extent\nmisleading. Furthermore, we agree that multiple (sequential or independent) Compton scatter\nevents detected in coincidence in an inorganic scintillator such as NaI(Tl) might result in a\nspectral shift of the Compton edge toward higher spectral energies as highlighted in the figure\nprovided by referee #3. That said, we believe that our statements and conclusions remain sub-\nstantiated and are also in agreement with the referee’s observations. In the following paragraphs,\nwe provide a concise discussion of our findings and thereby offering a more in-depth explanation\nof our perspective on comment #3 from the previous peer review report.\nFirst, it is important to note that we observe a shift of the Compton edge not toward higher\nbut toward smaller spectral energies, i.e. a negative Compton edge shift (cf. Fig. 4 in the main\nmanuscript as well as Fig. S26 in the Supplementary Information File). Based on these findings,\nwe conclude that the process highlighted by the referee cannot be the sole factor contributing\nto the observed Compton edge shift.\nSecond, we have to underline that all our Monte Carlo simulations feature a high-fidelity fully\ncoupled photon, electron and positron radiation transport as discussed in the Method section\nin the main manuscript (lines 374–376). Consequently, multiple Compton scatter events are\nproperly simulated for the entire mass model in all performed simulations and the potential shift\ndue to multiple Compton scatter events is included intrinsically in all our results. This is also true\nfor the additional simulations described in the supplementary Section S1.4. Consequently, the\nprocess highlighted by the referee cannot be responsible for the observed discrepancy between\nmeasurements and simulations.\nThird,asdemonstratedinthemainstudyinSection”BayesiancalibratedNPSMsimulations”\nour NPSM predicts this negative Compton edge shift for all analyzed Compton edges with high\naccuracy. Moreover, with the semi-analytical model derived in the supplementary Section S1.4,\nwe are able to show that the relation between photon energy, crystal size and negative Compton\nedge shift can be explained by the distribution of generated secondary electrons in the scintillator\nin combination with the decreasing trend in the relative light yield for higher energies (cf.\nsupplementary Figs. S27–S29). These findings are in line with our statements in comment #3 in\nthe previous peer review report. For more information on this topic, we kindly refer the reader to\nthe dedicated supplementary Section S1.4 in the Supplementary Information File, which provides\na detailed derivation and discussion of the results summarized here.\nMotivated by the comment, we rephrased several statements in the main manuscript and\nthe Supplementary Information File to enhance clarity and prevent potential misinterpretations\nby underlining the negative sign of the measured Compton edge shift. A full list of all minor\nchanges made can be found in Section 8.\nHere, we list all minor changes associated with the second peer review round.\n1. In the main manuscript, we rephrased the sentence ”In line with the Compton scatter theory\n(Supplementary Methods S1.4), we find a moderate increase in the Compton edge shift\nfor higher photon energies.” to ”In line with the Compton scatter theory (Supplementary\nMethods S1.4), we find an enhanced shift of the Compton edge toward smaller spectral\nenergies as the photon energy increases.”.\n2. In the Supplementary Information File (p. S8), we rephrased the sentence ”In general, we can\nidentify a consistent trend in the shift, i.e. an increasing Compton edge shift with increasing\nCompton edge energy.” to ”In general, we can identify a consistent trend in the shift, i.e.\nan enhanced Compton edge shift toward smaller spectral energies with increasing Compton\nedge energy.”.\n3. In the Supplementary Information File, we changed the phrase ”Compton edge shift” to\n”negative Compton edge shift” at multiple locations in the supplementary Section S1.4, i.e.\nonce in paragraph 6 on p. S9 and twice in the second last paragraph on p. S11.\n4. In line with the guide for authors, we revised all vector-graphics, both in the main manuscript\nand the Supplementary Information File, to ensure that all display items are clipped at the\ncorresponding axes outlines."}, {"round": 2, "reviewer_comments": {"reviewer_3": "It appears that the issue was not a technical one, but a mafter of befter wording. I'm pleased that this issue has been cleared up. Thank you."}, "author_response": "9 Response to Referee #3\nWe would like to thank the referee once again for the thorough review. We believe that the\nmodificationsmotivatedbyyourcommentshaveimprovedthemanuscriptsignificantlyandmade\nit more complete. We are really grateful and appreciate your time and effort in reviewing our\nmanuscript. All the best for your future.\nHere, we list all minor changes associated with the third peer review round.\n1. To avoid exaggerated language, we changed the following phrase in the last paragraph of the\nconclusion section in the revised manuscript from ”The novel methodology presented in this\nstudy offers a reliable and cost-effective alternative ...” to ”The methodology presented in\nthis study offers a reliable and cost-effective alternative ...”.\n2. In accordance with the provided guide to prepare final artwork, we changed the label font\nsize in all figures, both in the revised manuscript and the revised Supplementary Information\nFile, from 8 pt to 7 pt. As a result, we also changed the tick font size in all figures from 7 pt\nto 6 pt.\n3. In accordance with the provided guide to format the final article, we changed the font style\nof all mathematical vector terms from bold italic to bold non-italic in both, the revised\nmanuscript and the Supplementary Information File (matrices were already in the correct\nformat). In addition, we ensured that the vector font style in the revised figures is consistent\nwith those changes in the manuscript and the Supplementary Information File.\n4. In accordance with the provided guide to format the final article, we changed the symbol on\nline 20 on page S1 in the revised Supplementary Information File from ε to ε .\nPCA PCA\n5. In accordance with the provided guide to format the final article, we rephrased the Data and\nCodeavailabilitystatementsintherevisedmanuscripttobeinlinewiththeprovidedtemplate\nformat.\n6. In accordance with the provided guide to format the final article, we included a statement\nin the Acknowledgments of the revised manuscript highlighting the persons to which the\ndeclared funding was awarded.\n7. In accordance with the provided guide to format the final article, we changed the title\nof the Supplementary Information File from ”Supplementary Materials” to ”Supplementary\nInformation” as well as the titles of sections S1–S4 from ”Materials & Methods”, ”Supple-\nmentary Figures”, ”Supplementary Tables” and ”Supplementary Algorithms” to ”Materials”,\n”Figures”, ”Tables” and ”Algorithms”, respectively.\n8. For consistency reasons, we changed the term ”Supporting Information File” to ”Supple-\nmentary Information File” in both, the revised manuscript and Supplementary Information\nFile."}]

Decision letter