Modeling Scintillation Kinetics and Coincidence Time Resolution in Heterostructured Scintillators

In the search for new materials and technologies to push the timing performances of time-of-flight positron emission tomography (TOF-PET) detectors, it is important to have a model capable of predicting the coincidence time resolution (CTR) of the system to be implemented. While for bulk standard scintillators, a model that takes into account the intrinsic properties of the material (and the characteristics of the photodetector) is already well established, it has never been experimentally validated for composite structures. As heterostructured scintillators–i.e., the combination of two or more materials with complementary properties–are emerging as a possible solution to the conflict between fast timing and high detection efficiency for TOF-PET detectors, such validation becomes necessary. In this work, by using a time-correlated single photon counting (TCSPC) setup capable of simultaneously recording the TCSPC signal and the scintillation pulse on an event-by-event basis, we experimentally demonstrate that the scintillation kinetics of heterostructures can be modeled as a linear combination of the scintillation kinetics of the materials that constitute the heterostructure itself. Based on these results, we develop an extension of well-established CTR analytical model which can be applied to heterostructured scintillators.


I. INTRODUCTION
H ETEROSTRUCTURED scintillators are being increas- ingly explored as a potential solution to enhance the performance of time-of-flight positron emission tomography (TOF-PET) detectors [1], [2], [3], [4], [5], [6].In the context of TOF-PET detectors, there is often a trade-off between high sensitivity and fast timing [7].Traditional scintillator materials may offer high sensitivity but lack fast timing capabilities or vice versa.Heterostructures consist of a combination of two different materials with these two complementary properties and rely on the mechanism of energy sharing [2] (see Fig. 1).While the heavy material has higher probability of fully stopping the incoming γ -ray via photoelectric effect, in a fraction of events-the shared events-the photoelectron can escape from it, depositing some of its energy in the fast material, thus generating fast scintillation and improving the overall timing of the detector.Various studies have already assessed the gain in timing performances achievable when combining bismuth germanate (BGO) [or lutetium-yttrium oxyorthosilicate (LYSO)] with an ultrafast scintillator, like BC422 and EJ232 plastic scintillators, or barium fluoride (BaF 2 ) [2], [4], [5].
However, several challenges need to be addressed for the heterostructure to become competitive with L(Y)SO crystals, the current state-of-the-art TOF-PET detectors [8], [9].Plastic scintillators, despite their fast scintillation kinetics (decay time ≈ 1 ns and rise time < 50 ps [5], [10]) have a low density and effective atomic number.Considering the loss of stopping power, the achievable gain in time is not sufficient to lead to an improvement in the detector performances [7].It is therefore necessary to use denser and/or faster materials.One possibility is BaF 2 [4], [11], which has five times the density of plastic and a subnanosecond decay component due to cross-luminescence.The disadvantages in this case are the low-light yield (LY) and ultraviolet (UV) emission of the cross-luminescence process peaked at around 200 nm.While new photodetection technologies are being developed to increase the detection efficiency even in the vacuum UV (VUV) region [11], [12], the interest has also shifted to nanomaterial-based scintillators [13].Benefiting from quantum confinement effects, nanomaterials can feature ultrafast decay time and tunable emission, while also having a high atomic number.For all these reasons, the focus on the applicability of nanomaterials has shifted in recent years from optoelectronics to radiation detectors as well [14], [15], [16], [17], [18], [19].However, their size, instability, self-absorption, and low stopping power (as they are often embedded in a host matrix, usually polymer) still prevent them from exploiting their favorable properties in a full detector [20].
Given the many factors contributing to define the timing performances of a TOF-PET detector, the assessment of an analytic model to describe the coincidence time resolution (CTR) of heterostructured scintillators can simplify and speed up the search for the optimal fast material.For bulk scintillators, a well-established model already exists [10], [21], and in this work, we extend it to heterostructured scintillators providing an experimental validation for a simplified proof of concept.For this purpose, we first demonstrated that the scintillation kinetics of such a composite structure can be modeled as a linear combination of the scintillation kinetics of its constituent materials, weighted by the energy deposited in each.The experimental validation was carried out with a simplified proof of concept (3 × 3 × 3 mm 3 BGO + EJ232 heterostructure) using a time-correlated single-photon counting (TCSPC) [22] setup designed to simultaneously record the TCSPC signal and the scintillation pulse event-by-event [23], [24].Afterwords, knowing the scintillation kinetics of the heterostructure, its effective decay time (τ d,eff ) [10] was evaluated, and by measuring also its CTR, we demonstrated that the relation with LO being the light output (number of detected photons), is also valid for heterostructures.

A. Sample and Experimental Setups
For this study, a 3 × 3 × 3 mm 3 heterostructure made of alternated plates of BGO and EJ232 was used.The plates were optically decoupled, and to keep the structure compact, it was placed in a cubic plastic holder with one face open.We chose a thickness of 100 and 200 µm for BGO and EJ232 plates, respectively.Despite this configuration is not optimal in terms of stopping power [4], [5], [25], a large volume fraction of plastic allows us to maximize the energy sharing probability (i.e., large statistics of fast events) and to simplify this proof of concept study.
The heterostructure was measured in TCSPC mode under 511 keV irradiation ( 22 Na source) using the experimental setup schematized in Fig. 2 (left).The start-signal was given by a reference crystal (3 × 3 × 15 mm 3 LYSO:Ce) coupled to a Hamamatsu S13360-3050PE SiPM and readout by the NINO chip [27].The stop-signal was given by an ID-Quantique (IDQ) ID100-50 sensor [26], detecting in TCSPC mode the light produced by the heterostructure.A hole of 2 mm was drilled in the plastic holder ensuring a sampling over both materials from the IDQ.While the open face of the heterostructure was placed in dry contact with a Hamamatsu S13360-6050CS SiPM (see Fig. 2, top left), to also record the energy information and perform pulse shape discrimination.
For both the start-and stop-detector, the energy information was processed by a linear operational amplifier (AD8000).For a more detailed description of the setup and the assessment of the TCSPC conditions, please refer to [23].
The impulse response function (IRF) of the system was measured performing an analogous measurement by replacing the heterostructure with a PbF 2 crystal, as also described in [24].As PbF 2 is a sole-Cherenkov emitter, and as Cherenkov emission is prompt, this allowed us to accurately measure the IRF accounting also for photon time spread in the crystal [28], resulting in 75 ps sigma.
Afterward, we measured the CTR of the heterostructure with the setup described in [29].For this measurement, we used the Hamamatsu S13360 6050CS SiPM for extracting both the timing and energy information.The energy signal was processed, as before, by a linear operational amplifier, while for the time signal a high-frequency amplifier (two cascade radio frequency BGA616 amplifiers) was used (an updated version of [30]).

B. Methods
With the TCSPC setup, we were able to measure the scintillation time profile of the heterostructure (see Fig. 2 bottom right) and correlate it to the deposited energy.
The energy information was obtained by evaluating both the amplitude and integral of the analog energy signal.We integrated the pulse in its positive part in a time window up to 1.5 µs.
The reason why we recorded both the amplitude and the integral of the energy signal is to be able to classify the events according to the material where the energy is deposited-pure BGO, pure EJ232, or shared events.BGO and EJ232 have similar (LY, number of produced photons per energy unit)-8-10 ph/keV-but the effective decay time of BGO is almost a factor 100 slower than EJ232 [20].As described in previous works [2], [5], this makes the amplitude and the integrated charge of the energy signal two features allowing for clear pulse shape discrimination (see Fig. 2 top right).Through a change of coordinates -from (Amplitude, Int.Charge) to (Energy in EJ232, Energy in BGO)) -the energy deposited in the two materials was determined for each event.The pure BGO events were selected as those with 0 keV deposited in EJ232 and vice versa for the pure EJ232 events.The photopeak events were selected as those with total reconstructed energy (sum of the energy deposited in BGO and EJ232) between 440 and 665 keV.This classification is shown in Fig. 3, while more detailed information about the coordinates transformation can be found in the Appendix.
The scintillation kinetics was evaluated for the different categories of events separately: pure BGO, pure EJ232, all 511 keV (photopeak), and shared 511 keV with a specific amount of energy deposited in EJ232 (see Fig. 3).In all cases, the scintillation time profile was fit with the sum of exponential functions convolved with the system IRF.In (2), τ d i are the decay components and w i the corresponding normalized weights.We wanted to verify that the scintillation kinetics of the heterostructure is given by the linear combination of the scintillation kinetics of the two constituent materials, weighted by the energy deposited in each of them, i.e.,

S(t) H
with S(t) H/P/B being the intrinsic scintillation functions of the heterostructure, plastic, and BGO, respectively, as defined in (2).Therefore, we first performed the fit of the pure BGO and EJ232 events leaving all parameters (exponential decay components and corresponding weights) free.Next, the fit of the shared events at 511 keV was performed by fixing the decay components according to the results obtained from the fit of the pure events and leaving only the corresponding weights as free parameters.The 511 keV shared events were then divided according to the amount of energy deposited in EJ232: from a minimum of 50 keV to a maximum of 300 keV deposited in EJ232, we selected five intervals of 50 keV amplitude each.The distribution of events within each of them was fairly uniform, so we considered the mean energy deposited in EJ232 to be the average of the interval, i.e., for events with plastic reconstructed deposited energy between 50 and 100 keV, the mean energy deposited here is 75 keV, corresponding to 15% of the entire energy deposited (511 keV) in the heterostructure.
For each class of events, the effective decay time, defined as the weighted harmonic mean of all the decay components of the scintillator was evaluated.The effective decay time is the proper figure of merit for timing in the case of materials with multiple exponential decay components, and the CTR has been shown to be proportional to its square root [10].
By studying the dependency of the weights w i and of the effective decay time on the fraction of energy deposited in plastic E P (over the total of 511 keV), we showed experimentally the validity of (3).
For the analysis of CTR measurements, an equivalent events classification was performed and for each class of events, the FWHM of the coincidences peak was measured after applying time walk correction [31].The CTR of the heterostructure was obtained by subtracting in quadrature the contribution of the reference detector (2 × 2 × 3 mm 3 LSO:Ce:0.4%Ca,58 ps CTR).By correlating the measured CTR for the different intervals of energy deposited in plastic with the effective decay time of corresponding categories of events, we verified that (1) holds also for the different classes of events in heterostructured scintillators.

III. RESULTS
The fit results of the scintillation time profile of pure BGO, pure EJ232, and all 511 keV events are summarized in Table I.
The results of pure BGO and pure EJ232 events are in good agreement with the results published in [20], [32], and [33] obtained from the scintillation kinetics measurement of the two corresponding bulk materials.The fast component of BGO include both Cherenkov emission and the fast scintillation process that has already been observed in the aforementioned studies, both under X-ray and 511 keV irradiation.According to the literature, a possible explanation could be the result of quenching process due to the high density of elementary excitations at the end of the recoil electron track [34].
The agreement between the results found for the pure class of events and the bulk material confirms the accuracy of the  events selection: although contamination in the categories we have classified as pure cannot be excluded with certainty, it is in any case negligible and does not affect the results.
As the fastest decay components of BGO and EJ232 are comparable, they were merged into a single one (1.5 ns) for the fit of the shared events, allowing us to reduce the number of the decay components from five to four.
In all cases, the reduced chi-squared (χ 2 ) resulted ≈1, confirming the goodness of fit for the shared events when using the decay components of pure materials and constituting a first validation of (3).
The inset in Fig. 3 shows the fit function obtained for the five classes of shared 511 keV with a given energy deposition in EJ232 (together with the one of pure BGO and EJ232).The effect of a higher fraction of energy deposited in the plastic on the decay kinetics of the heterostructure can be appreciated.
To prove the validity of (3), we first solved the following system: S(t) P and S(t) B are the scintillation time profiles of plastic and BGO, respectively, namely By matching the terms with the same exponential decay component, a linear relation between the weights w 1,2,3,4 and the fraction of energy deposited in plastic was obtained.This analytical result found confirmation in the experimental data, as shown in Fig. 4. It is worth to mention that the reported values are not corrected by the photon detection efficiency of the IDQ (about 15% and 30% at the peak emission of EJ232 and BGO, respectively [26]), which explains the high contribution from BGO even when we are considering events with 55% of energy deposited in plastic.
After obtaining an analytical expression of w 1,2,3,4 as a function of E P , we substituted it into the definition of the effective decay time (4), and it resulted to have a linear dependency on E P τ d,eff ∝ 1/E P .
This relation was experimentally validated, as it can be seen from the linear fit of the measured decay time as a function of the mean fraction of energy deposited in plastic [see Fig. 5(a)].Finally, we measured the CTR of the heterostructure and, by performing an analogous events classifications as for the scintillation kinetics measurements, the CTR of photopeak events with a specific fraction of energy deposited in plastic was correlated with the effective decay time.By approximating BGO and EJ232 to have the same light output, we expect Also, this result was confirmed experimentally, as shown in Fig. 5(b).

A. Mathematical Approach: Assumptions and Approximations
The model presented in (3) was already introduced in [3].In this work, we give an experimental validation of it using BGO and EJ232 plastic scintillator as simplified proof-ofconcept.The advantage of using these two materials comes from their similar LY which allows us to develop a simple but effective model which can be extended to the more general case (see the next section IV-B).
As explained Section II-B, the similar LY together with the difference in the decay kinetics allows for a clear events separation.Moreover, it simplifies the analytic expression of the CTR (1) when dealing with two or more materials.In first approximation, we can consider only the number of produced photons (and neglect the factors which contribute to reduce the effective light output) given by the LY of the material multiplied by the energy deposited in there.When energy sharing occurs, the contributions from the materials involved must be summed up.Therefore, for a BGO + EJ232 heterostructure, the CTR for photopeak events (total deposited energy 511 keV) with a fraction E P of energy deposition in plastic is where LY P and LY B are the LY of plastic and BGO, respectively.However, in this specific case, we can approximate LY P ≈ LY B , and we are released from the dependency on the LY of the two materials, obtaining (8).The linear dependence of the CTR on the effective decay time (under square root) is thus limited to heterostructured scintillators consisting of two materials with equivalent light output.
A further clarification to be made is that the weighting factor that directly determines the contribution of the two materials to the overall scintillation kinetics is not simply the deposited energy but the actual number of photons detected by each (i.e., the light output).It can be described as The system where the PDE is the photodetection efficiency of the photodetector, and the LTE is the light transfer efficiency of the material.
Finally, it should be pointed out that for simplicity, we focused on the decay kinetics (2), while a complete description of scintillation time profile is given by the sum of biexponential functions, taking into account also the rising part.In any case, BGO and EJ232 have ultrafast scintillation rise time, below 50 ps [10], [20], and we can assume that what has been demonstrated for the scintillation decay apply to the rising edge.

B. Mathematical Generalization
So far, we proved, using the simplified case of BGO and EJ232 plastic scintillator, that the scintillation kinetics of a heterostructure can be as a linear combination of the scintillation kinetics of the of which it is by the energy deposited in each of them.Furthermore, have shown that the effective decay time is the appropriate figure of merit to express the of the CTR on the decay kinetics even in the case of heterostructured scintillators.Now, we want to extend this model to any two materials, i.e., with different LY, decay components, PDE, and LTE.To do this, we consider a heterostructure composed of two generic scintillators A and B. Scintillator A has N A decay time components, LY LY A , light transfer efficiency LTE A , light output N ph,A , and scintillation time profile 2) and ( 6)].The same holds for scintillator B. According to what we said in Section IV-A, their contribution to the overall scintillation is given by the number of photons extracted from each of them.(3) can therefore be generalized into with As before, E A (E B ) is the fraction of energy deposited in material A (B), and we are imposing the constraint E A + E B = 1 (i.e., considering only the events with full energy deposition, in this case 511 keV).However, while S(t) H as defined in (3) was automatically normalized due to the approximations made, here it is not the case.By integrating (11) Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
over time ([0, +∞]), we obtained as normalization factor The system in ( 5) can be generalized into By solving the system as before, we found the following expression for the weights w i where n ph,A (E A ) and n ph,B (E A ) are the fraction of photons (over the total N ph,TOT ) produced by the scintillator A and B, respectively, as function of the fraction of energy deposited in A.
The dependency of the weights on the fraction of energy deposited in one material (E A ) is not linear anymore, but if we approximate (as done for the specific case BGO + EJ232), the linearity is restored.
We can now generalize the expression of the effective decay time.By replacing the expressions just found for w i,A and w i,B in the definition of the effective decay time (4), we obtain The CTR of a heterostructured scintillator-for events with a fraction E A of energy deposited in material A-what can be generalized with N ph,A (E A ), N ph,A (E A ), and τ d,eff defined as in ( 12) and ( 15), respectively.

C. Outlook
For heterostructures to become an effective alternative to the current state-of-the-art PET detector, in parallel with the research for the optimal materials to combine, one should also focus on how to properly exploit the fast events.The ultimate goal in PET imaging is the maximization of the signal-to-noise ratio (SNR) of the reconstructed image.Several works have been published proposing a multikernel approach in case of different types of events [35], [36], [37], [38].This approach was first applied to scintillators which are also Cherenkov emitter, e.g., BGO [35], [36], [38] and more recently also to heterostructured scintillators [37], [38], as each event with different energy deposition in the two materials will give rise to a different TOF-kernel.
In [37], an model to the one presented in ( 9) to describe the of heterostructures proposed and used to reconstruct the simulated data from a heterostructure-based PET ring detector.A complementary work was carried out by [38].They first show how to analytically compute the SNR of a single TOF event considering different scenario (Gaussian and non-Gaussian TOF-kernel).Next, they show that in the presence of several event types (as in heterostructure case), the SNR of an event distribution is given by the squared sum of each event type, weighted by its relative abundance.Our work fits into this research line as it provides an experimental validation of the analytic model of the CTR for heterostructured scintillators, which is at the basis of the aforementioned studies.Though the experimental validation was performed on a simplified proof of concept, its generalization is straightforward as shown in Section IV-B.
Being able to correctly model the CTR of heterostructures will therefore simplify and speed up the search for the optimal materials to combine.Combining the results of this work with those of the aforementioned studies will also make it possible to choose materials not only based on the best CTR but also on the overall performance of the final system.

V. CONCLUSION
In this work, we used a simplified heterostructured scintillator (alternating layers of BGO and plastic scintillator EJ232) to demonstrate that the scintillation kinetics of this type of scintillator can be modeled as a linear combination of the scintillation kinetics of its constituent materials, weighted by the energy deposited in them.Furthermore, we confirmed experimentally that the effective decay time is the appropriate figure of merit to describe the overall scintillation decay of heterostructures when studying the dependence of the CTR on it.
The experimental validation was performed using a TCSPC setup, based on irradiation at 511 keV, capable of simultaneously recording the TCSPC signal and pulse shape on an event-by-event basis.This allowed us to classify events according to both the total amount of energy deposited (thus selecting photopeak events) and the material in which the energy was deposited.The scintillation kinetics were then analyzed separately for each class of events.
This study contributes to a cohesive understanding of heterostructured scintillators and at the same time provides an analytical model capable of predicting the temporal performance of heterostructures simply by knowing the scintillation properties and stopping power of the two materials to be combined, allowing Research and Development on this technology to be accelerated.

APPENDIX
As mentioned in Section II-B and discussed in literature [2], [5], the correlation between integrated charge and amplitude allows for a good discrimination of BGO and EJ232 events.The density scatter plot in Fig. 6(a) shows this correlation for the measured heterostructure.The pure events lie along a straight line because while the integrated charge and amplitude depend on the amount of energy deposited, their ratio depends only on the scintillation kinetics and light output, therefore is fixed for a given material.First, we evaluated the ratio between integrated charge and amplitude for each event, which allowed us to compute the angle between the x-axis and each line passing through the origin and a point on the scatter plot θ = 180 • /π • arctan Int.Charge Amplitude .
θ B and θ P , the angles of BGO and plastic events, respectively, correspond to the two maxima of the resulting distribution [Fig.6(b)].Through a change of coordinates, we can pass from the (Amplitude, Int.Charge) to the (EJ232 events, BGO events) coordinates system Assuming the SiPM saturation to be negligible, which is reasonable considering the relatively low LY of the two materials and the light loss due to not optimal light transport in heterostructure, to perform the energy calibration, we need at least one point per axis for which we know the energy deposited.One is clearly the photopeak of BGO events, while for the energy calibration along the x-axis (EJ232 events), we can estimate the hypothetical photopeak in plastic.In Fig. 6(a), we can observe an accumulation region extending from the photopeak in BGO to higher amplitude values: these are photopeak shared events with increasing energy deposited in plastic.Following this line, if photoelectric interaction for 511 keV γ -ray in plastic would be possible, we would find the photopeak in EJ232.The hypothetical photopeak in EJ232 was therefore estimated from the intersection of shared photopeak events and EJ232 events.To find the line outlining the shared photopeak events, we evaluated the angle of each line passing through the photopeak in BGO and any point in the scatter plot.First, we obtained the coordinates of the BGO photopeak (X BGO ph.peak , Y BGO ph.peak ) from its fit in the integrated charge and amplitude distribution [Fig.6 The main peak of the resulting distribution corresponds to BGO events, while the second one to the shared photopeak events.At this point, evaluating the intersection between the latter and the EJ232 line, we estimated the hypothetical photopeak in EJ232 and could perform the energy calibration along both axes.The resulting density scatter plot is shown in Fig. 6(f).
The selection of the pure events was done by considering the distribution of the energy deposited in BGO and the energy deposited in plastic (i.e., the projection of Fig. 6(f) along the yand x-axis, respectively), and taking the events under the peak at 0 keV (in BGO for pure EJ232 events and vice versa) within the FWHM of the same.

ACKNOWLEDGMENT
This work was carried out in the frame of the Crystal Clear Collaboration.The authors are thankful to Dominique Deyrail for his help in the heterostructure preparation.They express their gratitude to Stefan Gundacker for the valuable discussions that helped them in the improvement of the TCSPC setup allowing them to carry out this study.Author Contribution: Fiammetta Pagano: methodology, investigation, validation, data curation, software, formal analysis, writing-original draft, and visualization.Nicolaus Kratochwil: conceptualization, methodology, validation, investigation, writing-review-editing, and supervision.Loris Martinazzoli: methodology, software, validation, and writing-review-editing.Carsten Lowis: software and writing-review-editing.Marco Paganoni: supervision, project management and coordination, and funding.Marco Pizzichemi: methodology, supervision, writing-reviewediting, project management and coordination, and funding.Etiennette Auffray: conceptualization, methodology, supervision, writing-review-editing, project management and coordination, and funding.

Fig. 1 .
Fig. 1.Illustration of the concept of heterostructure and of the mechanism of energy sharing by the recoil photoelectron.Figure taken from [5].

Fig. 2 .
Fig. 2. Scheme of the experimental setup (left) and examples of the output of the energy signal for different types of events (top right) and of the scintillation time profile resulting from TCSPC measurement (bottom right).The blue points represent the data, the green line is the smoothing of the histogram (performed via moving average), and the red one is the fit function.

Fig. 3 .
Fig. 3. Density scatter plot showing the distribution of events depositing energy in BGO and in EJ232 and their classification.In the top right of the figure, the decay curves of different categories of events-pure BGO, pure EJ232, and shared 511 keV with increasing energy in plastic-are also shown.The corresponding categories of events are highlighted in the scatter plot with boxes of same colors.

Fig. 4 .
Fig. 4. Weights of the three decay components of the 511 keV shared events are represented as a function of the fraction of the mean energy deposited in plastic (E P ), for each group of events.The linear fit (dotted line) confirms the expected trend.

Fig. 5 .
Fig. 5. Effective decay time as a function of the inverse of the fraction of the mean energy deposited in plastic (1/E P ).The linear trend confirms the analytical result obtained in 7 (a).CTR as a function of the square root of the effective decay time.The obtained linear trend confirms the effective decay time an appropriate figure of merit to express the dependency of the CTR on the decay kinetics of heterostructure (b).

Fig. 6 .
Fig. 6.Illustration of the steps of the coordinates transformation.(a) Density scatter plot of integrated charge vs amplitude with BGO, EJ232, and 511 keV events outlined by red straight lines.(b) Angular distribution of the lines passing through the origin of the axis and each point in the density scatter plot.(c) Integrated charge distribution of all the events.(d) Amplitude distribution of all the events.(e) Angular distribution of the lines passing through the position of the photopeak in BGO and each point density scatter plot.(f) Density scatter plot resulting from the coordinates transformation.

TABLE I RESULTS
FROM THE FIT OF DECAY SCINTILLATION OF PURE BGO, PURE EJ232, AND ALL 511 keV EVENTS.BOTH THE DECAY TIME CONSTANTS (τ d ) AND THE CORRESPONDING WEIGHTS (w) ARE REPORTED