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

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 science [1][2][3][4][5][6][7]. Other applications include radiation protection, medical diagnostics and homeland security [8,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 simulations [10]. 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 energy [6,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 technique [12], recent studies could conclusively confirm the conjecture reported in earlier investigations [13][14][15] that not only organic but also inorganic scintillators exhibit a pronounced non-proportional relation between the deposited energy and the scintillation light yield [16][17][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 crystal [19,20]. Nonetheless, our understanding about these phenomena is still far from complete and, thanks to the advent of novel experimental techniques and the development of new scintillator materials, interest in scintillation physics has steadily grown over the past years [16][17][18][19][20][21][22][23][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 scintillator [1,[25][26][27][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 response [1,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 others [1,21,22,25,28,[31][32][33][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 size [25,28,31]. NPSMs are therefore particularly relevant for scintillators with large crystal volumes, e.g. in dark matter research, total absorption spectroscopy or remote sensing [1][2][3][4][5][6][7]30]. In addition, especially due to the sensitivity on the activator concentration and impurities [34], 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 NPSM [12,[35][36][37][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-proportionality [1,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 algorithms [39] to perform the NPSM inference with the gamma-ray spectral data. In contrast to traditional frequentist methods or simple datadriven 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 framework [40][41][42]. We simulate the detector response using a multi-purpose Monte Carlo radiation transport code in combination with parallel computing. To meet the required evaluation speed for the Bayesian inversion solver, we use machine learning trained polynomial chaos expansion (PCE) surrogate models to emulate the simulated detector response [43]. This approach offers not only a novel way to calibrate NPSMs with minimal effort-especially during the detector deployment-but it also allows new insights into the non-proportional scintillation physics without the need for additional electron response measurements.

Compton edge probing
To obtain the spectral Compton edge data, we performed gamma-ray spectrometry under controlled laboratory conditions [30]. 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 ( 57 Co, 60 Co, 88 Y, 109 Cd, 133 Ba, 137 Cs and 152 Eu) for the radiation measurements. However, only 60 Co, 88 Y and 137 Cs 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 FLUKA [45]. 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 modelled in less detail. For this simplifications, 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 nonproportional scintillation physics in our simulations [17,18,22]. In general, the sequence of scintillation processes in inorganic scintillators can be qualitatively divided into five steps [20,47,48]. 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 to radiatively recombine at the luminescent centers. Consequently, e − /h pairs need to convert to excitons by the classic Onsager Fig. 1 Compton edge probing to perform Bayesian inference on non-proportional scintillation models. a Monte Carlo mass model of the experimental setup to perform Compton edge probing with an inorganic gamma-ray scintillation spectrometer under laboratory conditions. The spectrometer consists of four individual 10.2 cm × 10.2 cm × 40.6 cm prismatic NaI(Tl) scintillation crystals with the associated photomultiplier tubes (PMT), the electronic components, e.g. the multi-channel analyzers (MCA), embedded in a thermal-insulating and vibration-damping polyethylene (PE) foam protected by a rugged aluminum detector box. We inserted radiation 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. The mass model figures were created using the graphical interface FLAIR [44]. For better visibility and interpretability, we applied false colors. b Overview of the Bayesian inference framework highlighting the gamma-ray spectrometry based Compton edge probing measurements, the Monte Carlo simulations using the multi-purpose code FLUKA [45] combined with the machine learning trained polynomial chaos expansion (PCE) emulator models supported by principal component analysis (PCA) as well as the Bayesian inference by Markov Chain Monte Carlo (MCMC) itself using UQLab [46]. c Radiation transport mechanisms inside the inorganic scintillation crystal, which is surrounded by a thin reflector layer and a rugged aluminum crystal casing. d Schematic representation of an inorganic scintillation crystal lattice including the activator atoms and point defects. e Mechanistic depictions of the various scintillation and quenching pathways for electron-hole pairs (e − /h) as well as excitons within the inorganic scintillation crystal lattice.
mechanism [49] 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 defects [20,22] as well as exciton-exciton annihilation described by the Birks mechanism [50].
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 by [22]: 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 Gektin [20], 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 Figs. 1c-e.

Bayesian inversion
We applied Bayesian inversion using Markov Chain Monte Carlo [39] 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 60 Co (activity A = 3.08(5) × 10 5 Bq) spectral dataset [30]  −0.82 × 10 2 MeVcm −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×10 1 MeVcm −1 and dE/ds | Birks = 1.853×10 2 MeVcm −1 for an ambient temperature of 18.8 • C [22].

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 the Figs. 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 toward 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 decomposition [51] to perform a global sensitivity analysis of the NPSM. Using the sum mode inversion pipeline, we present total Sobol (3)] keV, respectively. The measured net count rate cexp as well as the simulated net count rate adopting a proportional scintillation model c sim were presented already elsewhere [30]. We obtained the simulated net count rate c corr sim the same way as c sim but accounted for the non-proportional scintillation effects by the Bayesian calibrated NPSM obtained by the sum mode inversion pipeline. For the calibration, we used the 60 Co dataset. For all graphs presented in this figure, uncertainties are provided as 1 standard deviation (SD) shaded areas (coverage factor k = 1). These uncertainties are only visible for cexp.
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 studies [25][26][27][28][53][54][55]. 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 Lu 2 SiO 5 (LSO) was developed as an alternative to other available options at that time, such as Bi 4 Ge 3 O 12 (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 comparable [16]. 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 1994 [12], 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 degradation [28,56]. This example showcases the need for a better understanding and prediction of the intrinsic resolution in the development of new scintillators [57].
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 contribution [1]. Since our forward model explicitly accounts for the non-proportional scintillation physics by adopting a 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 σ 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 dataset [30].
Comparing the intrinsic and overall spectral resolution, we find an almost constant ratio σ 2 intr /σ 2 tot ≈ 0.35 for E γ ≳ 1500 keV. Around E γ ≈ 440 keV, there is a pronounced peak with σ 2 intr /σ 2 tot ≈ 0.42 and for E γ ≲ 440 keV, we observe a significant decrease in σ 2 intr /σ 2 tot with decreasing photon energy E γ . Moreover, we find a more complex behaviour 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) keV [52] 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 88 Y (A = 6.83(14) × 10 5 Bq) and 137 Cs (A = 2.266(34) × 10 5 Bq) radiation measurements to validate our calibrated NPSM. For the Monte Carlo simulations, we applied the posterior point estimators x MAP 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 88 Y and 137 Cs together with 60 Co, 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 scatter theory (Supplementary Methods S1.4), we find an enhanced shift of the Compton edge toward 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 88 Y and 137 Cs. 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 57 Co, 109 Cd, 133 Ba and 152 Eu together with a detailed uncertainty analysis for each source are attached in the Supplementary Information File for this study (Supplementary Figs. S18-S25).

Discussion
Here, we demonstrated that Compton edge probing combined with Monte Carlo simulations and Bayesian inversion can successfully infer NPSMs for NaI(Tl) inorganic scintillators. A detailed Bayesian data analysis revealed no significant differences between standard posterior point estimators and the related spectral detector response predictions for both inversion pipelines. Consequently, the Bayesian inversion results indicate that our methodology successfully constrained the NPSM parameters to a unique solution. However, we found statistically significant differences between our results and best-estimate literature values as well as between the individual scintillation crystals themselves. These results corroborate the experimental findings of Hull and his co-workers [34] and underscore the criticality of the NPSM calibration for every individual detector system.
Various studies reported a distortion of the Compton edge in gamma-ray spectrometry with inorganic scintillators [1,14,15,29,30]. In this study, we presented conclusive evidence that this shift is, at least partly, the result of the scintillation non-proportionality. Moreover, using our numerical models, we can predict the Compton edge shift as a function of the NPSM parameters. We observed a Compton edge shift toward smaller spectral energies for an increase in dE/ds | Birks and η e/h as well as a decrease in dE/ds | Trap . These results imply that an enhanced scintillation non-proportionality promotes a Compton edge shift toward smaller spectral energies. In line with these observations, the non-proportionality is enhanced by a large e − /h fraction, an increased Birks mechanism as well as a reduction in the e − /h trapping rate [20,24,48].
Further, we quantified the sensitivity of the NPSM on the individual NPSM parameters using a PCE-based Sobol decomposition approach. The sensitivity results indicate that η e/h has the highest sensitivity on the Compton edge, followed by dE/ds | Birks and dE/ds | Trap . However, previous studies showed a pronounced dependence of dE/ds | Trap on the ambient temperature [22,33]. In addition, we expect also a substantial change of the crystal structure by radiation damage, i.e. the creation of new point defects in harsh radiation environments [10,58]. Therefore, the obtained sensitivity results should be interpreted with care. dE/ds | Trap might be of significant importance to model the dynamics in the detector response with changing temperature or increase in radiation damage to the crystals, e.g. in deep space missions.
Using the Bayesian calibrated NPSM, we are also able to numerically characterize the contribution of the scintillation non-proportionality to the overall energy resolution. This intrinsic resolution sets a fundamental lower limit on the achievable spectral resolution for a given scintillator material, making it a key factor in the development of new scintillators. At higher photon energies (E γ ≳ 400 keV), we observed a significant contribution to the total spectral resolution (≥ 35%) with a maximum of ≈ 42% around 440 keV. At lower energies (10 keV ≲ E γ ≲ 400 keV), the intrinsic contribution is reduced and shows substantial distortions around the K-absorption edge for iodine at about 33 keV. We conclude that the non-proportional scintillation is a significant contributor to the total energy resolution of NaI(Tl). These observations are in good agreement with previous results [28,55,[59][60][61][62] and thereby substantiate not only the predictive power of our numerical model but showcase also its potential as a novel tool in the development of new scintillators.
Most of the theoretical studies focused on the prediction of NPSMs themselves. In contrast, available numerical models to predict the full detector response are scarce, computational intense and complex due to the adopted multi-step approaches with offline convolution computations [55,56,60]. In this study, we present an alternative way to implement NPSMs and simulate the full spectral detector response to gamma-ray fields by directly evaluating the NPSM online during the Monte Carlo simulations. This approach saves considerable computation time and has the additional advantage of not having to store and analyze large files with secondary particle data. We have used this implementation to predict the full spectral detector response for additional radiation fields accounting for non-proportional scintillation effects. Validation measurements revealed a significant improvement in the simulated detector response compared to proportional scintillation models. However, there are still some model discrepancies, especially at the lower and higher end of the Compton edge domain. These discrepancies might be attributed to systematic uncertainties in the Monte Carlo mass model or deficiencies in the adopted NPSM. Sensitivity analysis performed in a previous study in conjunction with the prior prediction density results might indicate the latter [30].
While we focused our work on NaI(Tl) in electron and gamma-ray fields, the presented methodology can easily be extended to a much broader range of applications. First, it is general consensus that the light yield as a function of the stopping power is, at least to a first approximation, independent of the ionizing particle type [16,31]. Second, the adopted NPSM was validated with an extensive database of measured scintillation light yields for inorganic scintillators, i.e. BGO, CaF 2 (Eu), CeBr 3 , Cs(Tl), Cs(Na), LaBr 3 (Ce), LSO(Ce), NaI(Tl), SrI 2 , SrI 2 (Eu), YAP(Ce) and YAG(Ce), among others [17,18,22]. From this, it follows that, given a gamma-ray field with resolvable Compton edges can be provided, our methodology may in principle be applied to any combination of inorganic scintillator and ionizing radiation field, including protons, α-particles and heavy ions. However, it is important to note that our methodology relies on the observation of Compton edge shifts with a sufficient signalto-noise ratio (SNR). We have shown that these shifts are influenced by the strength of scintillation non-proportionality of a given scintillator. As a result, scintillator materials that exhibit only a mild non-proportional scintillation response, e.g. LaBr 3 (Ce) or YAP(Ce), may present challenges for the calibration of a NPSM due to the reduced SNR in the Compton edge shift. Further investigations are required to assess the applicability of our methodology to such scintillators. That said, the presented methodology can be readily adapted using Bayes's theorem to address low SNR cases more effectively by combining multiple Compton edge domains or by probing additional spectral features distorted by the non-proportional scintillation response.
In summary, we conclude that NPSMs are essential for accurate detector response simulations, especially for scintillators with large crystal volumes [25,28,31], e.g. in dark matter research, total absorption spectroscopy or remote sensing [1][2][3][4][5][6][7]30]. The novel methodology presented in this study offers a reliable and cost-effective alternative to existing experimental methods to investigate nonproportional scintillation physics phenomena and perform accurate full detector response predictions with Bayesian calibrated NPSM. Moreover, this new technique does not require any additional measurement equipment and can therefore be applied for any inorganic scintillator spectrometer, also during detector deployment. This is especially attractive for applications, where the scintillator properties change in operation, e.g. due to radiation damage or temperature changes, but also for detector design and the development of novel scintillator materials. Last but not least, we can use the derived numerical models not only for NPSM inference but also to investigate and predict various scintillator properties, e.g. intrinsic resolution or Compton edge dynamics, and thereby contribute to a better understanding of the complex scintillation physics in inorganic scintillators.

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 thermalinsulating 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 spectra [30]. We used seven different calibrated radionuclide sources ( 57 Co, 60 Co, 88 Y, 109 Cd, 133 Ba, 137 Cs and 152 Eu) 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 [30].
After postprocessing the spectral raw data according to the data reduction pipelines described in [30], we extracted the Compton edge spectral data from the net count rate spectra. The spectral domain of the Compton edge D E was defined as where E is the spectral energy, σ tot the energy dependent total resolution characterized by the standard deviation [30] and E FEP the FEP associated with the Compton edge E CE . Neglecting Doppler broadening and atomic shell effects, we compute E CE according to the Compton scattering theory [10] as follows: where m e c 2 is defined as the energy equivalent electron mass. In this study, we consulted the ENDF/B-VIII.0 nuclear data file library [63] for nuclear decay related data as well as the Particle Data Group library [64] for fundamental particle properties.
To investigate the sensitivity of the selected Compton edge domain D E on the Bayesian inversion results, we performed a sensitivity analysis on D E . Within the uncertainty bounds, the inversion results have proven to be insensitive to small alterations in D E (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 study [30] as well as the Supporting 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.1 [45] together with the graphical interface FLAIR version 3.1-15.1 [44]. 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 effects [16,18,23,55]. 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 eventby-event basis using the custom user routine usreou together with the detect card. The number of primaries was set to 10 7 for all simulations, which guarantees a maximum relative statistical standard deviation σ stat,sim,k /c sim,k < 1% and a maximum relative variance of the sample variance VOV k < 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 [30].
To implement the NPSM described by Eq. 1, we developed an additional user routine comscw. Similar to [1,65], 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 [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 D E 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 models [43]. As shown by [66][67][68], any function Y = M (X) with the random input vector X ∈ R M ×1 and random response vector Y ∈ R N ×1 can be expanded as a so-called polynomial chaos expansion provided that E[∥Y ∥ 2 ] < ∞: where a α := (a 1,α , . . . , a N,α ) ⊺ ∈ R N ×1 are the deterministic expansion coefficients, α := (α 1 , . . . , α M ) ⊺ ∈ N M ×1 the multi-indices storing the degrees of the univariate polynomials ψ α and Ψ α (X) := M i=1 ψ i α i (X i ) the multivariate polynomial basis functions, which are orthonormal with respect to the joint probability density function To reduce the computational burden, 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 N ′ of output variables compared to the original number N of spectral variables, i.e. N ′ ≪ N [43]. Similar to [69], we computed the emulated computational model responsê M PCE (X) in matrix form as: with µ Y and σ Y being the mean and standard deviation of the random vector Y and Φ ′ the matrix containing the retained eigenvectors ϕ from the PCA, i.e. Φ ′ := (ϕ 1 , . . . , ϕ N ′ ) ∈ R N ×N ′ . On the other hand, the vector Ψ (X) ∈ R card(A ⋆ )×1 and matrix A ∈ R N ′ ×card(A ⋆ ) store the multivariate orthonormal polynomials and corresponding PCE coefficients, respectively. The union set A ⋆ := N ′ j=1 A j includes the finite sets of multi indices A j for the N ′ output variables following a specific truncation scheme.
We used a Latin hypercube experimental design X ∈ R M ×K [70,71] 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 ∈ R N ×K for this design was then evaluated using the forward model described in the previous subsection. We adopted a hyperbolic truncation scheme A j := {α ∈ N M : ( M i=1 α q i ) 1/q ≤ p} with p and q being hyperparameters defining the maximum degree for the associated polynomial and the q-norm, respectively. To compute the PCE coefficient matrix A, we applied adaptive least angle regression [72] 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 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 code [46] 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 framework [40], we approximate the measured spectral detector response y ∈ R N ×1 with a probabilistic model combining the forward model M(x M ) and model parameters x M ∈ R M M ×1 with an additive discrepancy term ε, i.e. y := M(x M ) + ε. For the discrepancy term ε, which characterizes the measurement noise and prediction error, we assume a Gaussian model π(ε | σ 2 ε ) = N (ε | 0, σ 2 ε I N ) with unknown discrepancy variance σ 2 ε . On the other hand, as discussed in the previous subsection, we emulate the forward model M(x M ) with a PCE surrogate model M PCE (x M ). Consequently, we can compute the likelihood function as follows: In combination with the prior density π (x), we can then compute the posterior distribution using Bayes' theorem [42]: where we assume independent marginal priors, i.e.
We defined the marginal priors based on the principle of maximum entropy [73] as well as empirical data from previous studies [17,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' theorem [42], we therefore incorporate the results obtained by the sum mode inversion pipeline in 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: All Bayesian computations were performed with the UQLab code [46]. We applied an affine invariant ensemble algorithm [39] to perform Markov Chain Monte Carlo (MCMC) and thereby estimate the posterior distribution π (x | y). We used 10 parallel chains with 2 × 10 4 MCMC iterations per chain together with a 50% burn-in. The convergence and precision of the MCMC simulations were carefully assessed using standard diagnostics tools [42,74]. We report a potential scale reduction factorR < 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 study [30]. 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 a first step, we extracted specific spectral domains containing a singlet or multiplet of FEPs from a set of measured count rate spectra covering a spectral range between 122 and 1836 keV [75]. In a second step, we fitted a spectral peak model based on a sum of independent Gaussian peaks together with a numerical baseline [76] to the selected singlets or multiplets using weighted non-linear least-squares (WNLLS) regression combined with the interior-reflective Newton method [77]. 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 a first step additional Monte Carlo simulations for an isotropic 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 a 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 energy [30]. In a third step, we adopted the extracted σ intr from the individual energy calibrated FEPs to train a Gaussian Process (GP) regression model with [78]: where we applied a polynomial trend function of the second order, i.e. f (E γ ) := 1, E γ , E 2 γ ⊺ and β := (β 0 , β 1 , β 2 ) ⊺ , a homoscedastic noise model with the noise variance It is worth noting that, due to the known asymmetry in the FEPs [1,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 [78]: 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 σ 2 GP 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 [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 [30] with an adapted dispersion σ 2 tot − σ 2 intr . For completeness, we included this adapted dispersion model in Fig. 3f. Data availability. The radiation measurement raw data presented herein are deposited on the ETH Research Collection repository: https://doi.org/10.3929/ethz-b-000528920 [75]. Additional datasets related to this study are available from the corresponding author upon reasonable request.
Note that, in contrast to previous studies [1][2][3], we standardize our model response Y with Y * := diag (σ Y ) −1 (Y − µ Y ) to account for the differences in the variance of the individual response variables. We can then perform an eigenvalue decomposition of the correlation matrix Σ Y with the eigenvalues λ j and eigenvectors ϕ j := (ϕ 1 , . . . , ϕ N ) ⊺ satisfying Σ Y ϕ j = λ j ϕ j for j = 1, . . . , N . Since Σ Y is symmetric and positive definite, the eigenvectors define an orthonormal basis R N = span({ϕ j } N j=1 ) and we can perform an orthogonal transformation of our random vectors Y * as follows: We call the transformed vectors Z := (Z 1 , . . . , Z N ) ⊺ the principal components of Y * . Once we get the principal components, we can transform them back to the original response variable space with To reduce the dimensions of our problem, we retain only N ′ principal components with the highest variance and thereby approximate our random vector Y as where we choose N ′ := min{S ∈ {1, . . . , N } : For the PCE model part, we start again with the polynomial chaos expansion of the model response M (X) with the random input vector X ∈ R M ×1 as described in Eq. 3 in the main study: where a α := (a 1,α , . . . , a N,α ) ⊺ ∈ R N ×1 are the deterministic expansion coefficients, α := (α 1 , . . . , α M ) ⊺ ∈ N M ×1 the multi-indices storing the degrees of the univariate polynomials ψ α and Ψ α (X) := M i=1 ψ i α i (X i ) the multivariate polynomial basis functions, which are orthonormal with respect to the joint probability density function f X of X, i.e. ⟨Ψ α , Ψ β ⟩ f X = δ α,β . For computational purposes, we truncate the PCE series by adopting a truncation set A j for the multi-index α of each individual response variable j = 1, . . . , N resulting in: For the truncation, we can use a hyperbolic truncation scheme defining the multi-index set as A j := {α ∈ N M : ( M j=1 α q j ) 1/q ≤ p} with p and q defining the maximum degree for the associated polynomial and the q-norm, respectively.
To reduce the computational burden, we can now combine these results and perform the PCE not in the original response variable space but in the truncated principal component space. For that, we insert Eq. S5 in Eq. S3: which we can rearrange by introducing the union set A ⋆ := N ′ j=1 A j to: or expressed in a more compact matrix form: as well as the two matrices Φ ′ ∈ R N ×N ′ and A ∈ R N ′ ×card(A ⋆ ) storing the multivariate orthonormal polynomials Ψ α , the retained eigenvectors ϕ j and the PCE coefficients a j,α , respectively.
For model training, we introduce an experimental design with the input matrix X ∈ R M ×K and response matrix Y ∈ R N ×K for K instances, M input variables and N response variables. For the PCA model, we can use the response matrix Y to estimate µ Y , σ Y as well as Σ Y : with Y * denoting the standardized response matrix storing the standardized response variables y * := diag (σ Y ) (y −μ Y ), i.e. Y * := y * (k) , . . . , y * (k) , . . . , y * (K) ⊺ ∈ R N ×K . On the other hand, a rich variety of non-intrusive and sparse methods exist to estimate the PCE coefficient matrix A using both, the input matrix X ∈ R M ×K and response matrix Y [4]. In the main study, we chose the least angle regression algorithm [5] due to its high evaluation speed and its high accuracy even for very small experimental designs.

PCA-PCE based Hoeffding-Sobol decomposition & Sobol indices
One of the major advantages to use PCE emulators for computational intense simulations is the relation between PCE and the Hoeffding-Sobol decomposition and thereby Sobol indices [6]. For completeness, we repeat here some of the theory already discussed elsewhere [3,[6][7][8] and derive the PCA-PCE based Sobol indices accounting for the standardization in the PCA discussed in the previous subsection. We start with the global variance decomposition theory derived by Sobol [8]. It can be shown that for any univariate integrable function M (X) with M mutually independent random input variables X i in D X and i = {1, 2, . . . , M }, there exists a unique functional decomposition, which is often referred to as Hoeffding-Sobol decomposition [8]: 2. All the terms in the functional decomposition are orthogonal: Further assuming that the function M (X) is square-integrable, the functional decomposition in Eq. S10 may be squared and integrated to provide the variance decomposition: with the total variance V and the partial variances V u defined as: Based on these results, Sobol indices S u can be defined as a natural global sensitivity measure of M (X) on the input variables X u : Consequently, S u represents the relative contribution of the set of variables u to the total variance V . First order indices S i indicate the influence of X i alone, whereas the higher order indices quantify possible interactions or mixed influences between multiple variables. In addition, we can also define the total Sobol index S T i to evaluate the total effect of an input parameter X i on M (X): As a result, S T i includes not only the effect of X i alone but in addition the effect induced by all interactions between X i and the other variables. This is also the reason, why the sum of the total Sobol indices i S T i can in fact exceed 1. As an example, if we have an interaction between the variables X 1 and X 2 , their interaction effect on M is counted twice, once in S T 1 and another time in S T 2 . This example in mind, it is easy to see that the peaks in Fig. 3(e) in the main study highlight regions, where the interaction terms between the individual variables significantly contribute to the total effect on M. We have to add that for our study, the absolute values of S T i are of less importance. We are more interested in the relative size of S T i , because the comparison of these values allows us to draw conclusions about the relative importance of the corresponding variables X i for the model response M (X). As shown by [6], S T i can also be computed as: where we use ∼i to denote a set of indices, which do not include i, i.e. S ∼i = S v with v = {1, . . . , i − 1, i + 1, . . . , M }. Suppose now that we have a PCA-PCE surrogate model to emulate the vector-valued model response Y = M (X) with the random input vector X ∈ R M ×1 and random response vector Y ∈ R N ×1 . To derive the S T i,k for each response variable k ∈ {1, 2, . . . , N }, we start with Var X ∼i [E X i [Y k ]] from Eq. S17b by replacing Y k with the k th component of Eq. S8: where we used ϕ row k := (ϕ k1 , . . . , ϕ kN ′ ). We can simplify this expression by expanding the first term and considering that the expectation vanishes for all principal components, i.e. E [AΨ (X)] = 0: As shown by [3], due to the orthonormality of the polynomial basis {Ψ α } α∈A ⋆ , we can further simplify Eq. S19b resulting in: with the subset A ⋆ i=0 := {α ∈ A ⋆ | α i = 0}. Using these results, we can compute the total variance with: In the end, we get the total PCE-PCA based Sobol index S T i,k for the input variable i and the response variable k by inserting Eq. S20 and Eq. S21 into Eq. S17b:

Uncertainty analysis
For completeness, we repeat here the uncertainty analysis pipeline adopted for the measured and simulated pulse-height spectra and highlight some changes to [9]. For the radiation measurements, the statistical uncertainty of the net count rate spectra c exp,k characterized by the standard deviation was computed adopting a probabilistic Poisson model [10]: where C gr,k and C bg,k are the gross and background counts in channel k together with the gross and background measurement live times t gr and t bg , respectively. The small statistical uncertainty in the live time measurement is neglected. To compute the source activity A as a function of the measurement date t, we use the fundamental exponential law of decay, i.e. A = A 0 · 2 −∆t/t 1/2 [10]. The uncertainty induced by the source activity A normalization is quantified using the standard error propagation methodology for independent variables [11]: with the reference activity A 0 and associated uncertainty σ A 0 provided by the vendor, the source half life t 1/2 [12] as well as the time difference ∆t = t − t 0 between the reference date t 0 and the measurement date t. Contributions of the uncertainties in t 1/2 and ∆t to σ A are found to be less than 1% for all performed measurements and are therefore neglected. We then summarize the total experimental uncertainty as follows [11]: For the simulations, we computed the statistical uncertainty of the net count rate spectrum c sim,k characterized by the standard deviation as follows [10]: where c sim,kl are the individual broadened energy deposition events in the detector channel k, N dep the number of recorded events and N pr the number of simulated primaries. It is good practice in Monte Carlo studies to report not only the estimated uncertainty in the sample mean c sim,k using the sample standard deviation σ stat,sim,k but also the so called variance of the sample variance VOV k for the detector channel k to quantify the statistical uncertainty in σ 2 stat,sim,k itself [13]: The propagation of the systematic uncertainties for the simulated detector response was performed by the Monte Carlo sampling technique. We considered the same model parameters for the uncertainty propagation as in [9]. These parameters are the energy calibration factor D 1 keV −1 as well as the empirical resolution parameters B 1 [−] and B 2 [−]. However, we adapted the marginal distributions by introducing truncated normal distributions as summarized in Table S4. In addition, we accounted for the statistical dependence of the model parameters B 1 and B 2 by correlated sampling using the Gaussian copula C N [14]: with the log-transformed variable B * 1 := log (B 1 ), the linear correlation matrix R obtained by the regression analysis, the marginal distribution functions F provided in Table S4, the bivariate Gaussian distribution function Φ 2 associated with the Gaussian copula C N and the inverse cumulative distribution function of the standard normal distribution Φ −1 , respectively. The energy calibration factor D 1 is sampled independently according to the corresponding marginal as in [9]. For more details and relevant literature on the copula theory, the reader is referred to [15,16].
The N MC ∈ N >1 independently drawn input samples X MC = (x (1) , ..., x (m) , ...x (N MC ) ) ⊺ from the probabilistic input model with X := (D 1 , B 1 , B 2 ) ⊺ are then propagated through the postprocessing pipeline described in [9] to obtain the corresponding spectral count rate samples Y MC = (c sim,k ) ⊺ with k ∈ {1, ..., 1024}. These samples can then be used to compute the sample standard deviation σ sys,sim,k similar to Eq. S9b and thereby quantify the systematic uncertainty with respect to the empirical model parameters D 1 , B 1 and B 2 . The total uncertainty characterized by the sample standard deviation can be summarized in the same way as for the experimental uncertainty [11]:

Compton edge shift analysis
To better understand the nature of the Compton edge shift utilized in our study, we quantify here the spectral shift between the measured Compton edge energy and the theoretical value according to the Compton scattering theory (cf. Eq. 2 in the main study). Because the measured Compton edges are obscured by the finite spectral resolution, quantification of the exact position in the measured pulse-height spectrum would require additional coincidence measurements [17,18]. Therefore, similar to a previous study [19], we use an alternative approach by quantifying the spectral shift between the already available Monte Carlo simulations with a proportional scintillation response and the measured pulse-height spectra. It is important to add that, based on the findings reported in the main study as well as due to the improved signal-to-noise ratio, we focus our investigation here on the sum channel.
In a first step, we determine the inflection points as a characteristic measure of the corresponding Compton edges using spline regression [20]. We then compute the Compton edge shift as the spectral difference between the determined inflection points for the measured and simulated spectra. We apply this method to different Compton edges, i.e. 477.334 (3)  In the Fig. S26, we present the results of our Compton edge shift analysis for the sum channel. In general, we can identify a consistent trend in the shift, i.e. an enhanced Compton edge shift toward smaller spectral energies with increasing Compton edge energy. As discussed in the main study, our NPSM predicts the Compton edge shift for all analyzed Compton edges with high accuracy. However, because our approach is based on complex Monte Carlo simulations, interpretability of the results and their connection to the various underlying physical processes is challenging. Therefore, to improve our understanding of the relationship between the NPSM and the resulting Compton edge shift, we develop here a simplified semi-analytical model.
We start by deriving the light yield function L as a function of the initial electron kinetic energy E k by integrating Eq. 1 from the main study from E k down to the mean excitation energy I: To compute the integral in Eq. S30, we first need a model to describe the differential energy loss dE per differential path length ds as a function of the kinetic electron energy E k . Similar to Payne and his co-workers [21,22], we apply a modified Bethe-Bloch model derived by Joy and Luo [23] to described dE/ds as a function of E k in units of eVÅ −1 : with the scintillator related mass density ρ, the atomic number Z and the molecular weight A. In accordance with the results obtained by Payne and his co-workers, we fix the stopping power correction factor c in Eq. S31 to c = 2.8. To account for radiative losses as well as relativistic effects at higher energies, we combine Joy's model with the ESTAR database [24]. We also consulted the ESTAR database for all material related properties of NaI(Tl) (cf. Table S5). The resulting total stopping power model together with the individual model components are shown in Fig. S27(a) for NaI(Tl). We then combine the derived stopping power model with the light yield L (dE/ds) to perform the integration in Eq. S30. For the model parameters in L (dE/ds), i.e. η e/h , dE/ds | Ons , dE/ds | Trap and dE/ds | Birks , we applied the maximum a posteriori (MAP) probability point estimates to compute a mean light yield function as well as the full set of posterior samples to derive the corresponding credible intervals for both, the sum channel and the individual crystals associated with the sum and S7 single mode inversion pipelines, respectively. The resulting relative light yield functions L (E k ) /E k are shown in Fig. S27(b). The characteristic shape of these relative light yield curves in Fig. S27(b) has been extensively documented by numerous previous empirical studies [21,22,[25][26][27][28], illustrating an increase in light yield with increasing energy for E k ≪ 10 keV, a prominent peak around 10 keV, followed by a subsequent decrease in yield for higher energies.
After successfully deriving the light yield function L as a function of the kinetic electron energy E k , we now continue by investigating the connection between the observed negative Compton edge shift and the non-proportional nature of this light yield function. For an ideal detector with a proportional scintillation response, we would have a constant relative light yield function: So, we can easily see that for a detector with a non-proportional scintillation light yield function, we get a spectral shift ∆E, if we convert the produced scintillation light of an electron with energy E k,1 at a different energy E k,2 : It is important to add that in our analysis, we implicitly assume a continuous deceleration of the involved electrons, starting from their initial kinetic energy E k and progressing down to the scintillator specific excitation energy I. In particular. we neglect electrons escaping from the scintillator. From Eq. S33b we can conclude that this shift ∆E will be positive for L (E k,1 ) /E k,1 > L (E k,2 ) /E k,2 and negative for L (E k,1 ) /E k,1 < L (E k,2 ) /E k,2 . Moreover, with an increase in the relative difference |1 − [L (E k,1 ) /E k,1 ] / [L (E k,2 ) /E k,2 ]|, we expect a proportional increase in the magnitude of the spectral shift.
We can now use these insights for our Compton edge shift analysis. In gamma-ray spectrometry with inorganic scintillators, the energy calibration is typically performed using full energy peaks (FEPs), sometimes also called photopeaks [9,10]. In other words, to relate the photon energy with the generated scintillation photons, we use the relative light yield function for j = {1, . . . , N e − } electrons with kinetic energy E k,j generated in the scintillator leading subsequently to a FEP. Consequently, in order to compute the Compton edge shift for inorganic scintillators using our model in Eq. S33b, we need to investigate the light yield for both Compton edge (CE) and FEP events, more specifically the integrated light yield for all electrons generated during these events. We start with the CE events: As already discussed in the main study, in a CE event, a photon enters the scintillator, undergoes a single Compton scattering (COM) event with a deflection angle of 180 • , i.e. full back-scattering, and subsequently escapes the scintillator. During this interaction, some of the photon's energy gets transferred to a single atomic electron. Neglecting Doppler broadening and atomic shell effects [29,30], the transferred energy E CE k is equivalent to the CE energy discussed in Eq. 2 in the main study, i.e.: with E 0 γ being the initial photon energy and m e c 2 the energy equivalent electron mass. Because only one COM event takes place with a deterministic energy transfer, the light yield for a CE event can easily be calculated as follows: FEP events on the other hand are more complex because they involve a variable number of COM events with a subsequent photoelectric absorption (PE) of the photon in the scintillator. For simplicity, we neglect again Doppler broadening as well as atomic shell effects and consider only secondary electrons generated during COM and PE events. In particular, we neglect fluorescence photons and Auger electrons. Using these simplifications, we can calculate the light yield for a FEP event involving j = {1, . . . , N e − } electrons with kinetic energy E k,j as a sequence of N COM COM events followed by a single PE event: where we denote with E PE k,j and E COM k,j the kinetic energies of the electrons generated during PE and COM events, respectively, and with E i γ the photon energy after i subsequent COM events. Both, the number of COM events (N COM ) as well as the transferred energy in these COM events (E i−1 γ − E i γ ) are linked in a complex stochastic process and depend on the photon energy as well as the properties of the scintillator [29][30][31]. To estimate these variables for our specific detector system, we apply once again Monte Carlo methods. More specifically, we adopted the multi-purpose Monte Carlo code FLUKA with the same physics settings as described in the main study [32]. For the semi-analytical model described here, the mass model only included the scintillation crystal embedded in a vacuum environment. To estimate the scintillator response, we irradiated the mass model with an isotropic and uniform monoenergetic photon flux of energy E 0 γ using the FLOOD mode with the BEAMPOSit card. We repeated these simulations for 31 different photon energies E 0 γ in the spectral range [500 , 2000] keV with a spacing of 50 keV. To score N COM as well as E i γ , we applied the user routine mgdraw. In Fig. S28(a), we present the probability density for the scored number of COM events before absorption (N COM ) as a function of the photon energy E 0 γ in a prismatic NaI(Tl) scintillator with dimensions 10.2 cm × 10.2 cm × 40.6 cm, i.e. the same crystal size as for our detector system used in the main study. In line with previous results [33], we find a moderate increase of the mean number of COM events with increasing photon energy ranging from 1.4 at 500 keV up to 2.3 at 2000 keV. Combining these results together with Eq. S35 and Eq. S36b, we can now compute the spectral Compton edge shift according to Eq. S33b as follows: where we used the fact that for a FEP event in our simplified framework, the following special relationship holds: Fig. S28(b), we show the resulting negative Compton edge shift −∆E as a function of the photon energy E 0 γ for the same prismatic NaI(Tl) scintillator. In general, we find a good agreement between the predictions of our simplified semi-analytical model and the experimental results. Deviations can be attributed to the various simplifications and assumptions made during model derivation, e.g.

S9
Landau fluctuations, Doppler broadening, atomic shell effects, detector cross-talk, neglected secondary particles such as fluorescence photons and Auger electrons or escaping electrons.
By discriminating the predicted spectral shift for different number of COM events N COM , we can gain also some further insights in the underlying physics. First, we find a pronounced increase in the Compton edge shift |∆E|, both for an increase in the photon energy E 0 γ as well as for an increase in the number of COM events N COM . This is in line with our predictions discussed above, i.e. we expect a proportional increase in the magnitude of the spectral shift for an increase in the relative difference ]. It is now easy to see that, due to the decreasing trend in L for higher energies (cf. Fig. S27(b)), the relative difference will increase for an increase in E k and subsequently E 0 γ . Furthermore, we find that the relative light yield for CE events is on average smaller than for FEP events with E 0 γ ∈ [500 , 2000] keV, i.e. L(βE 0 γ )/β < ⟨ N e − j=1 L(E FEP k,j )⟩. This explains the negative sign for the spectral shift observed by our NaI(Tl) detector system.
With this newly derived semi-analytical model, we have now also a tool to investigate the relation between the negative Compton edge shift and the size of a scintillation crystal. In Fig. S29(b), we present the predicted mean Compton edge shift as a function of the photon energy E 0 γ alongside the relation between N COM and E 0 γ . We find a pronounced and consistent increase in the negative Compton edge shift for an increase in crystal size over the entire spectral domain [500 , 2000] keV. From our semi-analytical model and the results in Fig. S29(a), it is evident that this trend can be explained by the increase in N COM for bigger scintillation crystals.
In summary, the semi-analytical model derived in this section cannot only successfully predict the trends and sign of the Compton edge shift with increasing photon energy E 0 γ , but it can also be used to investigate the relation between crystal size and Compton edge shift and thereby supports the interpretation of the results and findings in the main study obtained by high-fidelity Monte Carlo simulations. Chain index [ -] Step index [ -]  Chain index [ -] Step index [ -]  Chain index [ -] Step index [ -]  Chain index [ -] Step index [ -]  Chain index [ -] Step index [ -]  x MAP x Mean x Median Step index [-]  x MAP x Mean x Median Step index [-]  x MAP x Mean x Median Step index [-]  x MAP x Mean x Median Step index [-]  x MAP x Mean x Median Step index [-]  x MAP x Mean  x MAP x Mean  x MAP x Mean  x MAP x Mean Methods in the main study). The measured net count rate c exp as well as the simulated net count rate adopting a proportional scintillation model c sim were presented already elsewhere [9]. We obtained the simulated net count rate c corr sim the same way as c sim but accounted for the non-proportional scintillation effects by the sum mode inversion pipeline presented in this study. For the calibration, we used the 60 Co dataset [9]. For all graphs presented in this figure, uncertainties are provided as 1 standard deviation (SD) shaded areas (coverage factor k = 1). These uncertainties are only visible for c exp . Fig. S19: Uncertainty quantification for the 60 Co spectral detector response. The measured and simulated mean net count rates c exp and c sim are shown for the sum channel using a 60 Co calibrated radionuclide source (A = 3.08(5) × 10 5 Bq) together with the corresponding uncertainty estimates, i.e. the combined statistical and systematic measured uncertainty σ tot,exp , the simulated statistical uncertainty σ stat,sim as well as the simulated systematic uncertainty σ sys,sim , using 1 standard deviation values. The measurement results were presented already elsewhere [9]. Two different scintillation models have been used for the simulations: a Proportional scintillation model published in [9]. b Bayesian calibrated non-proportional scintillation model obtained by the sum mode inversion pipeline presented in this study. Distinct spectral regions, i.e. the backscatter peak (BSP), the Compton edge (CE) as well as the full energy peaks (FEP) are highlighted for both graphs. Note that the highlighted CE region refers to the lower Compton edge at 963.419(3) keV associated with the photon emission line at 1173.228(3) keV. The normalized residual level | c exp − c sim | /σ tot with σ tot := σ 2 tot,exp + σ 2 tot,sim for a coverage factor of 2 is marked with the horizontal dash-dotted black line in the lower subfigures. More information on the numerical computation of the uncertainty estimates can be found in Section and [9]. The measured and simulated mean net count rates c exp and c sim are shown for the sum channel using a 88 Y calibrated radionuclide source (A = 6.83(14) × 10 5 Bq) together with the corresponding uncertainty estimates, i.e. the combined statistical and systematic measured uncertainty σ tot,exp , the simulated statistical uncertainty σ stat,sim as well as the simulated systematic uncertainty σ sys,sim , using 1 standard deviation values. The measurement results were presented already elsewhere [9]. Two different scintillation models have been used for the simulations: a Proportional scintillation model published in [9]. b Bayesian calibrated non-proportional scintillation model obtained by the sum mode inversion pipeline presented in this study. Distinct spectral regions, i.e. the backscatter peak (BSP), the Compton edges (CE) as well as the full energy peaks (FEP) are highlighted for both graphs. The normalized residual level | c exp − c sim | /σ tot with σ tot := σ 2 tot,exp + σ 2 tot,sim for a coverage factor of 2 is marked with the horizontal dash-dotted black line in the lower subfigures. More information on the numerical computation of the uncertainty estimates can be found in Section and [9]. i.e. the combined statistical and systematic measured uncertainty σ tot,exp , the simulated statistical uncertainty σ stat,sim as well as the simulated systematic uncertainty σ sys,sim , using 1 standard deviation values. The measurement results were presented already elsewhere [9]. Two different scintillation models have been used for the simulations: a Proportional scintillation model published in [9]. b Bayesian calibrated non-proportional scintillation model obtained by the sum mode inversion pipeline presented in this study. Distinct spectral regions, i.e. the backscatter peak (BSP), the Compton edge (CE) as well as the full energy peak (FEP) are highlighted for both graphs. The normalized residual level | c exp − c sim | /σ tot with σ tot := σ 2 tot,exp + σ 2 tot,sim for a coverage factor of 2 is marked with the horizontal dash-dotted black line in the lower subfigures. More information on the numerical computation of the uncertainty estimates can be found in Section and [9]. i.e. the combined statistical and systematic measured uncertainty σ tot,exp , the simulated statistical uncertainty σ stat,sim as well as the simulated systematic uncertainty σ sys,sim , using 1 standard deviation values. The measurement results were presented already elsewhere [9]. Two different scintillation models have been used for the simulations: a Proportional scintillation model published in [9]. b Bayesian calibrated non-proportional scintillation model obtained by the sum mode inversion pipeline presented in this study. The normalized residual level | c exp − c sim | /σ tot with σ tot := σ 2 tot,exp + σ 2 tot,sim for a coverage factor of 2 is marked with the horizontal dash-dotted black line in the lower subfigures. More information on the numerical computation of the uncertainty estimates can be found in Section and [9]. Bq) together with the corresponding uncertainty estimates, i.e. the combined statistical and systematic measured uncertainty σ tot,exp , the simulated statistical uncertainty σ stat,sim as well as the simulated systematic uncertainty σ sys,sim , using 1 standard deviation values. The measurement results were presented already elsewhere [9]. Two different scintillation models have been used for the simulations: a Proportional scintillation model published in [9]. b Bayesian calibrated non-proportional scintillation model obtained by the sum mode inversion pipeline presented in this study. The normalized residual level | c exp − c sim | /σ tot with σ tot := σ 2 tot,exp + σ 2 tot,sim for a coverage factor of 2 is marked with the horizontal dash-dotted black line in the lower subfigures. More information on the numerical computation of the uncertainty estimates can be found in Section and [9]. Bq) together with the corresponding uncertainty estimates, i.e. the combined statistical and systematic measured uncertainty σ tot,exp , the simulated statistical uncertainty σ stat,sim as well as the simulated systematic uncertainty σ sys,sim , using 1 standard deviation values. The measurement results were presented already elsewhere [9]. Two different scintillation models have been used for the simulations: a Proportional scintillation model published in [9]. b Bayesian calibrated non-proportional scintillation model obtained by the sum mode inversion pipeline presented in this study. The normalized residual level | c exp − c sim | /σ tot with σ tot := σ 2 tot,exp + σ 2 tot,sim for a coverage factor of 2 is marked with the horizontal dash-dotted black line in the lower subfigures. More information on the numerical computation of the uncertainty estimates can be found in Section and [9]. Bq) together with the corresponding uncertainty estimates, i.e. the combined statistical and systematic measured uncertainty σ tot,exp , the simulated statistical uncertainty σ stat,sim as well as the simulated systematic uncertainty σ sys,sim , using 1 standard deviation values. The measurement results were presented already elsewhere [9]. Two different scintillation models have been used for the simulations: a Proportional scintillation model published in [9]. b Bayesian calibrated non-proportional scintillation model obtained by the sum mode inversion pipeline presented in this study. The normalized residual level | c exp − c sim | /σ tot with σ tot := σ 2 tot,exp + σ 2 tot,sim for a coverage factor of 2 is marked with the horizontal dash-dotted black line in the lower subfigures. More information on the numerical computation of the uncertainty estimates can be found in Section and [9].  for collisional losses at low energies derived by Joy and Luo [23] as well as radiative and collisional losses at higher energies predicted by the ESTAR database [24]. b Relative light yield L(E k )/E k as a function of E k for both, the sum channel and the individual scintillation crystals associated with the sum and single mode inversion pipelines, respectively. We applied the maximum a posteriori (MAP) probability point estimates for the individual model parameters, i.e. η e/h , dE/ds | Ons , dE/ds | Trap and dE/ds | Birks , derived in the main study to compute the mean relative light yield function normalized at 1 MeV according to Eq. S30. In addition, we present 99% central credible intervals for each individual relative light yield function using the full set of posterior samples obtained by the sum and single mode inversion pipelines. A list of all material properties for NaI(Tl) used to compute the stopping power predictions as well as the resulting relative light yield functions can be found in Table S5.   Table S2: Posterior statistics summary. This table includes posterior point and dispersion estimators for the Bayesian inverted non-proportional scintillation models obtained by the sum and single mode inversion pipelines. The listed estimators are the maximum a posteriori (MAP) probability estimate x MAP , the posterior mean x Mean and the posterior median x Median together with the 95% credible interval and the posterior standard deviation σ x for the parameters x := dE/ds | Birks , dE/ds | Trap , η e/h , σ 2 ε ⊺ , i.e. the Birks related stopping power parameter dE/ds | Birks , the trapping related stopping power parameter dE/ds | Trap , the free carrier fraction η e/h as well as the discrepancy model variance σ 2 ε .

Supplementary Algorithms
Algorithm S1 COMSCW(dE/ds | Birks , dE/ds | Ons , dE/ds | Trap , η e/h ) COMSCW is a custom user-routine for the multi-purpose Monte Carlo code FLUKA [32] called at each energy deposition event in the scintillation crystal. We adapt this routine by weighting each electron or positron energy deposition event by the adopted non-proportional scintillation model (NPSM) [35]. The algorithm accounts for both continuous as well as local energy deposition events. As described in the Methods section in the main study, we set a kinetic energy threshold of 1 keV below which the electrons and positrons as well as particles generated by atomic deexcitation are no longer transported and their energy is deposited on the spot. We refer to these events as "local" energy deposition events. On the other hand, above this threshold, ionization losses are evenly distributed along the particle step [36,37]. Hence, we call these events "continuous". The pseudo-code added below is a simplified version of the one implemented in our forward model. For more details, we kindly refer to the actual routine deposited on the ETH Research Collection repository: https://doi.org/10.3929/ethz-b-000595727 [38]. load E k ▷ Kinetic particle energy 9: if ParticleID = "electron" then