Pushing Cherenkov PET with BGO via coincidence time resolution classification and correction

Bismuth germanate (BGO) shows good properties for positron emission tomography (PET) applications, but was substituted by the development of faster crystals like lutetium oxyorthosilicate (LSO) for time-of-flight PET (TOF-PET). Recent improvements in silicon photomultipliers (SiPMs) and fast readout electronics make it possible to access the Cherenkov photon signal produced upon 511 keV interaction, which makes BGO a cost-effective candidate for TOF-PET. Tails in the time-delay distribution, however, remain a challenge. These are mainly caused by the high statistical fluctuation on the Cherenkov photons detected. To select fast events with a high detected Cherenkov photon number, the signal rise time of the SiPM was used for discrimination. The charge, time delay and signal rise time was measured for two different lengths of BGO crystals coupled to FBK NUV-HD SiPMs and high frequency readout in a coincidence time resolution setup. The recorded events were divided into 5 × 5 categories based on the signal rise time, and time resolutions of 200 ± 3 ps for 2 × 2 × 20 mm3 and 117 ± 3 ps for 2 × 2 × 3 mm3 were measured for the fastest 20% of the events (4% in coincidence). These good timing events can provide additional information for the image reconstruction in order to increase the SNR significantly, without spoiling the detector sensitivity. Putting all photopeak events together and correcting for the time bias introduced by different numbers of Cherenkov photons detected, time resolutions of 259 ± 3 ps for 20 mm long and 151 ± 3 ps for 3 mm long crystals were measured. For a small fraction of events sub-100 ps coincidence time resolution with BGO was reached for a 3 mm short pixel.


Introduction
In the past, bismuth germanate (Bi 4 Ge 3 O 12 ) used to be a good choice for positron emission tomography (PET) in terms of its high density, stopping power Moses (2007), high photo fraction and low production cost. However, due to its slow decay time and low light yield it was substituted by lutetium -based (LSO, LYSO, LGSO) scintillating crystals for time-of-flight PET (TOF-PET) to ensure a fast time resolution for noise equivalent count -rate improvement and, hence, better image quality. Nowadays, the appearance of fast silicon photomultipliers (SiPMs) with good photon detection efficiency (PDE) in the near ultraviolet and low single photon time resolution (SPTR) values combined with fast readout electronics allows to make the best use of the Cherenkov signal produced in BGO for a 511 keV interaction Kwon et al (2016), Piemonte et al (2016), Brunner and Schaart (2017), Cates et al (2018), Gundacker et al (2019), Cates and Levin (2019).
The detection of few prompt photons gives a significant improvement in time resolution, and is one possible pathway towards a time resolution of 10 ps for TOF-PET Lecoq (2017), Gundacker et al (2016).
BGO has a density of 7.13 g cm −3 Moses (2007), intrinsic light yield of 10.7 photons per keV and a decay time of around 360 ns Gundacker et al (2020), which makes it much slower than LSO:Ce codoped with Ca. The mean Cherenkov photon yield upon 511 keV interaction for BGO is about 17 ± 3 photons between 310-850 nm Gundacker et al (2020). Considering the light transfer efficiency (LTE) from the crystal to the photo detector (LTE = 0.58 for 2 × 2 × 3 mm 3 pixel, Teflon wrapped and optically coupled with high  Gundacker et al (2020). Due to non-isotropic emission from an event to event basis, the number of detected Cherenkov photons underlies high fluctuations Gundacker et al (2020). These fluctuations are responsible for the long tails in the 511 keV coincidence delay histogram, when events with very few or no Cherenkov photons are detected, hence only the BGO scintillation is used to determine the time stamp. To suppress these tails and to better determine the timestamps, a discrimination between events with high numbers of Cherenkov photons and low numbers of Cherenkov photons is necessary. One possible way to carry out the event separation is by employing a double-sided readout, either with two SiPMs Kwon et al (2019) or by combining a SiPM and a micro channel plate photomultiplier to select only fast events. One drawback of this method is the need for an additional readout channel to exploit the presence of Cherenkov photons and, therefore, distinguish between events giving a good or worse time resolution.
In this contribution, we introduce a single -sided readout approach to separate events by the use of the signal rise time of the SiPM, which is correlated with the photon density and, therefore, connected to the time resolution of the events. It was already demonstrated that, for a sampling pixel composed of BGO and plastic, the signal rise time allows a separation between BGO and plastic events Turtos et al (2019).

Coincidence time resolution (CTR) setup
A 22 Na source with an activity of 260 kBq (Nov. 2019) emits two gammas with 511 keV back to back. The gammas are detected in coincidence by two BGO crystals wrapped in Teflon and glued with Meltmount (refractive index n = 1.582) to 4 × 4 mm 2 FBK NUV-HD SiPMs with 40 µm 2 single photon avalanche diode (SPAD) size. The SiPMs are read out by the high frequency (HF) electronics for the time signal and an analog operational amplifier for the energy signal. The signals are digitized by a LeCroy DDA735Zi oscilloscope with 3.5 GHz bandwidth and a sampling rate of 20 Gs s −1 for four channels. The crystals, SiPMs and the HF electronics are in a light tight box which is kept stable at 16 C • . An illustration of the setup is shown on the left of figure 1. A detailed description of the HF readout is discussed in Gundacker et al (2019).

Data recording
The energy signal of both channels is integrated over a width of 160 ns and stored for offline photopeak selection. For the time signal the leading edge threshold is set on the oscilloscope calculating the signal crossing via interpolation. To find the best configuration, SiPM bias voltage and leading edge threshold scans were performed. The single SPAD amplitude is 44 mV for the best performing SiPM bias of 39 V (10.8 V overvoltage). In addition, the signal rise time of each channel between a 10 mV and 50 mV threshold was recorded and stored. The lower value of the rise time threshold was set to be well above the electronic noise floor, while the upper value was set to be above the signal amplitude of a single SPAD. In total, the charge and signal rise time per channel and the coincidence time delay were stored for each event and analyzed offline.

Energy selection
The charge-histogram was fitted with a Gaussian function around the photopeak. All events that had an energy lower than 440 keV and higher than 665 keV were rejected. This energy range was taken to be consistent with clinical TOF-PET scanners Kolthammer et al (2014). Using a narrower photopeak selection mainly decreases the number of events without changing the time resolution but for a few picoseconds. The impact on the time resolution for different energy windows is discussed in section 4.1.

Event identification based on signal rise time
The motivation for measuring the signal rise time is that for events with higher numbers of Cherenkov photons the photon density at the very beginning of the scintillation signal is higher, resulting in better time resolution. We call this type of events fast events. In contrast, for pure BGO scintillation, the photon density is lower resulting in a relatively flatter SiPM signal (slow event) and worse time resolution. This means, for BGO, which has a mixture of Cherenkov photons with high photon density and scintillation with low photon density, we can correlate the signal rise time with the photon density and thus estimate whether high or low numbers of Cherenkov photons are detected for an individual event, as illustrated on the right -hand side of figure 1.

Time-delay correction based on signal rise time
By using a leading edge threshold, the SiPM signal passes the threshold earlier for fast events than for slow events. This results in a time bias, or time walk, which needs to be corrected. For normal scintillation, the classical approach is a time -walk correction based on the amplitude of the signal. For the special case of BGO, the crucial Cherenkov photons only contribute to a tiny fraction of the amplitude, while the BGO scintillation dominates the height of the pulse, such that the classical approach cannot be used and we use the signal rise time instead.
To correct for this time bias, the events were split into five parts based on the signal rise time (from fast to slow) per channel, giving 25 different event categories, representing all (100%) of the coincidence events detected. For each event category j a histogram of the coincidence time delays is drawn and the data are fitted with a double Gaussian, as stated in equation (1) . (1) The first part is used to model the time delay coming from Cherenkov photons, while the second part is used to model the scintillation, such that σ 1 < σ 2 . The position of the maximum of the first Gaussian µ j 1 gives the time bias for each category. The corrected time for each delay time event i inside category j is stated in equation (2), where t j i (x source ) is a function of the annihilation source position and µ j 1 (rise times) is a function of the signal rise timet After the correction all events are merged together to calculate the time resolution with correction, using again a double Gaussian fit. For a TOF-PET detector, µ j 1 , which depends on the rise time of the two detectors, needs to be measured once as the calibration and must not be modified afterward to preserve the translation linearity of the positron emission point.

Results
The results are divided into different aspects of the analysis. In section 3.1 the impact on the event separation is presented by showing the extreme cases for fast and slow events and the CTR distribution for different rise times. In section 3.2 the rise-time-induced time delay is discussed, and for all the photopeak events a time -delay correction is performed.

Time resolution on different event categories
At the top left of figure 2 the signal rise time distribution for a 2 × 2 × 20 mm 3 BGO crystal is shown. Each zone between the red lines represents 20% of the total number of selected photopeak events. The absence of events with a signal rise time faster than 100 ps is due to the SiPM and electronic response, while the tail at higher rise times is due to the absence of Cherenkov photons. Carrying out the splitting on both channels leads to a grid with 25 bins, where each bin contains 4% of the events in coincidence. On the top right of figure 2 the different mixtures of fast and slow events per channel are illustrated. In this context fast-slow means, that the signal rise time is fast (low value) on the left channel, while it is slow (high value) on the right channel. On the bottom left of figure 2 the delay histogram for the categories fast-fast and slow-slow is plotted. By selecting the fastest 20% of the events per channel (4% in coincidence) the time resolution is 200 ± 4 ps full width at half maximum (FWHM) compared to 333 ± 6 ps FWHM for the slowest 20% of the events. On the bottom right of figure 2 the delay histogram for the non-symmetric events (fast-slow and slow-fast) is shown for the same crystal. The center of the distribution moves, as discussed in section 2.5 and 3.2, while the time resolution is for both categories around 250 ps FWHM. One can also see, that for the fast-slow category the tails are more pronounced to the left, while for slow-fast they are on the opposite side.
For crystals showing a non-trivial (not single Gaussian) distribution of the time delay not only the FWHM, but also the full width at tenth maximum (FWTM) is of importance. The individual time resolutions for all 25 bins are shown in figure 3 for 2 × 2 × 20 mm 3 (left) and 2 × 2 × 3 mm 3 (right) pixels in units of FWHM (top) and FWTM (bottom), respectively. It is to highlight that the best 4% of the events for the 20 mm long crystal have a time resolution of 200 ± 4 ps FWHM, which is only a factor of two worse than the best time resolution achieved with the same configuration and geometry of LSO:Ce codoped with Ca, as reported in Gundacker et al (2019).

CTR with time-delay correction
As illustrated in figure 1, different signal rise times introduce an offset (or timewalk) in the time delays, which is identified by the event selection. At the top of figure 4 the time delay, which is the position of the centroid of the first part of the double Gaussian fit from the delay histogram, is plotted for different rise time categories for 20 mm and 3 mm long crystals. For the longer crystal a total time difference of 172 ps between the two extreme antisymmetric event categories is measured, while for the 3 mm long crystal it is only 130 ps. This difference could point to a sensitivity in the depth-of-interaction (DOI), since for 3 mm long crystals the photon time spread is smaller while the LTE is larger, in contrast to 20 mm long crystals.
After correcting for the individual time delays all events are put together in a delay histogram and the time resolution after time-delay correction is calculated. The histograms for both tested lengths are shown at the bottom of figure 4. We tried to use a single Gaussian fit, but this neither matches the peak height of the distribution, nor is it in agreement with the pronounced tails of the distribution. The blue line represents the applied double Gaussian fit with the individual contributions in dotted blue and green. The horizontal lines represent the 50% height (FWHM), 10% height (FWTM) and 1% height (FW100M). For 20 mm long crystals we calculated a time resolution of 259 ± 3 ps FWHM and for 3 mm of 151 ± 3 ps FWHM, compared to 288 ± 3 ps (20 mm) and 169 ± 3 ps (3 mm) without time -delay correction.
The time delay after correction is centered around 0, since for the correction the centroid of the first Gaussian was subtracted. However, the delay time displacement in accordance with a source position  The dotted lines (blue + green) represent the individual contributions of the double Gaussian fit (blue) compared to a single Gaussian fit (red) and the measured data (black). movement or different electron positron annihilation positions in PET is preserved, since once the correction values are determined they are not modified. The whole correction procedure can be understood to be equal to a standard time walk or amplitude walk correction. The only difference in this case is that we use the signal rise time instead of the signal amplitude, as discussed on the right of figure 1.

Impact of the energy range
The photopeak window was set to be similar to a commercial PET scanner using a lower threshold 440 keV. Generally it can be noted that for lower energy thresholds, fewer Cherenkov photons are produced and detected, leading to an overall deterioration of the time resolution. Looking only at the time resolution for a set of events, it would even be beneficial to take all events without cutting on the energy but to cut on the rise time instead. However, for a PET scanner one would sacrifice on spatial resolution by accepting more scattered gammas in the body. The optimum trade off between spatial resolution and sensitivity might also depend on the used reconstruction algorithms and on how different timing kernels are incorporated, which is a subject for future studies. Narrowing the energy window to 500 keV the time resolution slightly improves, but at the same time the number of coincidence events drops by a factor of three compared to a cut at 440 keV. The time resolution results for different energy ranges are shown in table 1.

Modifying the bin size
So far we have always kept the number of category bins to 5 × 5, hence each bin represents 4% of the total number of events in coincidence. For the corrected time resolution, such a splitting is sufficient, since the maximum time walk is in the range of 130 ps (172 ps) for 3 mm (20 mm) long crystals, and applying five bins is already sufficiently lower than the achieved CTR per bin. Already, by splitting and correcting into 2 × 2 bins, most of the time bias is corrected leading to a time resolution of 157 ± 3 ps FWHM for 3 mm and 267 ± 3 ps FWHM for 20 mm. It was tested to slice into 10 × 10 bins, hence each bin represents 1% of the events in coincidence, observing that the improvement in the CTR after delay correction is within the error of the measurement. However, for the fastest 1% of the events in coincidence time resolutions of 107 ± 4 ps FWHM for 3 mm and 183 ± 4 ps FWHM for 20 mm were found. Further splitting into 14 × 14 bins was tested for the 3 mm crystal, leading to the best time resolution of 99 ± 3 ps FWHM for the fastest category (around 0.5% of all photopeak events in coincidence). Having this splitting in many TOF-kernels needs to be seen as additional information to the overall timing performance of BGO. Although only a subsection of events have a good time resolution, the very good sensitivity of BGO is preserved, as the remaining events with moderate time resolution are used as well.
A test was performed to modify the rise time window instead of 10-50 mV to 10-120 mV, leading to a negligible deterioration of the timing performance compared to the default setting (10-50 mV). This is interesting for integrating the proposed method into systems, because the exact rise time thresholds applied seem to be less critical.

Tails in the distribution
It is important to point out that the constant tails in the delay histogram at the bottom of figure 4 have their origin in the configuration of the setup. Since the leading edge discrimination was put at very low thresholds it can mean that a dark count event fires on one channel before the real signal is detected within the allowed coincidence window. In this case the signal rise time of the dark count event was slow, such that these tails only occur in slow events and could be removed by removing some events. Also, if using a second validation trigger put at a higher threshold these tails would disappear. Since the dark count event can happen on the left or on the right channel, this events are present on both sides of the CTR distribution. Zooming into the CTR histogram we see that the height of the tails is a bit higher on one side than on the other one, indicating that the dark count rate (DCR) of the SiPMs is not identical.
We also tested this method for standard scintillators of 2 × 2 × 3 mm 3 LSO:Ce codoped with Ca and could not find a significant CTR improvement, since this method depends on the fluctuation in photon density, which for LSO is smaller compared to BGO, due to its high light yield and fast scintillation response. In contrast, for other Cherenkov radiators such as PWO, PbF 2 (no scintillation) or LuAG:Pr an improvement is expected, which is a subject for further studies.

Conclusion
We applied an event selection based on the electronic SiPM signal rise time for BGO to classify and separate events based on their timing performance. We measured the timing performance for two lengths of BGO and classified the events into 25 different categories. For the best 20% of the events per channel, which corresponds to 4% in coincidence we found time resolutions of 200 ps FWHM for 20 mm and 117 ps for 3 mm long crystals. We have shown that the high fluctuation of detected Cherenkov photons leads to a time -walk effect, which can be corrected by the signal rise time, leading to time resolutions of 259 ps FWHM for 20 mm and 151 ps FWHM for 3 mm long crystals.
Up until now it is difficult to estimate the full impact of event TOF-kernel classification on image quality after reconstruction, since algorithms for testing are still under development. Nevertheless, we can conclude that event classification and time -walk correction can improve the best achievable BGO timing significantly, which opens new doors for high sensitivity and cost -effective TOF-PET machines.