DOI estimation through signal arrival time distribution: a theoretical description including proof of concept measurements

The challenge to reach 10 ps coincidence time resolution (CTR) in time-of-flight positron emission tomography (TOF-PET) is triggering major efforts worldwide, but timing improvements of scintillation detectors will remain elusive without depth-of-interaction (DOI) correction in long crystals. Nonetheless, this momentum opportunely brings up the prospect of a fully time-based DOI estimation since fast timing signals intrinsically carry DOI information, even with a traditional single-ended readout. Consequently, extracting features of the detected signal time distribution could uncover the spatial origin of the interaction and in return, provide enhancement on the timing precision of detectors. We demonstrate the validity of a time-based DOI estimation concept in two steps. First, experimental measurements were carried out with current LSO:Ce:Ca crystals coupled to FBK NUV-HD SiPMs read out by fast high-frequency electronics to provide new evidence of a distinct DOI effect on CTR not observable before with slower electronics. Using this detector, a DOI discrimination using a double-threshold scheme on the analog timing signal together with the signal intensity information was also developed without any complex readout or detector modification. As a second step, we explored by simulation the anticipated performance requirements of future detectors to efficiently estimate the DOI and we proposed four estimators that exploit either more generic or more precise features of the DOI-dependent timestamp distribution. A simple estimator using the time difference between two timestamps provided enhanced CTR. Additional improvements were achieved with estimators using multiple timestamps (e.g. kernel density estimation and neural network) converging to the Cramér–Rao lower bound developed in this work for a time-based DOI estimation. This two-step study provides insights on current and future possibilities in exploiting the timing signal features for DOI estimation aiming at ultra-fast CTR while maintaining detection efficiency for TOF PET.


Introduction
Time-of-flight positron emission tomography (TOF-PET) with detectors having ultra-high coincidence time resolution (CTR) would bring many imaging advantages due to a lower uncertainty of the annihilation localization along the line of response (Lecoq 2017, Conti andBendriem 2019). TOF information, provided by excellent CTR, increases the effective sensitivity and enhances the reconstructed image quality by reducing the statistical noise (Vandenberghe et al 2016). In particular, the recently launched 10 ps TOF-PET challenge (The 10 ps challenge 2020) opens new perspectives , but also imposes more constraints on the exact determination of the annihilation and optical photons travel path and time. CTR below 100 ps FWHM has recently been achieved using 2 × 2 × 20 mm 3 LSO:Ce:Ca crystals and FBK NUV-HD SiPMs read out by fast high-frequency electronics (Cates et al 2018. The achievable CTR with such long crystals using standard single-ended readout is however approaching a limit imposed by the blurring effect of depth-ofinteraction (DOI) on the timing precision. Indeed, even with instantaneous emission and zero-jitter timing photodetection, reaching 10 ps CTR remains inaccessible because of the DOI variability between coincidence events. A DOI difference between two coincident detectors was shown to induce a bias that degrades the timing accuracy, hence the overall CTR for all DOI combinations (Toussaint et al 2019, Loignon-Houle et al 2020. Also, a well-known problematic effect of DOI variability is the parallax effect which degrades the radial spatial resolution uniformity (MacDonald and Dahlbom 1998).
Multiple approaches to retrieve DOI information in PET detectors were explored over the years (Ito et al 2011). DOI encoding methods making use of the temporal information contained in the detection signal shape are common and typically require detectors made from assemblies of crystals (phoswich concept) with distinct temporal kinematics (Saoudi et al 1999, Bergeron et al 2015. Nonetheless, using the time distribution of the signal was proposed in many previous works with detector concepts using a single crystal type. Spatio-temporal localization of an event in a monolithic scintillator using both the signal time and spatial distribution on a segmented photosensor was previously shown (van Dam et al 2013, Iltis andSnoussi 2015, Tabacchini et al 2015). Detectors with double-sided readout typically use the relative signal intensity at both crystal ends to estimate the DOI (Ren et al 2014, Seifert and Schaart 2015, Selfridge et al 2018, but can also exploit the temporal information of the signal to further refine the interaction positioning resolution (Kang et al 2015, Han et al 2019. A detector concept with single-sided readout estimating the DOI through a light sharing scheme in a scintillator array was recently improved by using the signal timing information to enhance the DOI correction efficiency on the measured CTR (Pizzichemi et al 2019). Most of the detector concepts using a temporal scheme rely on signal dispersion among multiple (2) photosensors. A combined estimation of the DOI, energy and timing information in a detector concept based on single-ended readout phosphor-coated crystals was also proposed using the signal waveform obtained with a PMT, achieving CTR values in the range of ∼250-300 ps (Berg et al 2016).

Motivation of a time-based DOI estimation
With the advent of 20 mm long scintillation detectors capable of reaching <100 ps FWHM CTR with a single new fast SiPM readout (Gundacker et al 2019), a pure time-based DOI estimation approach could become feasible. The scintillation light propagation time, which depends on the DOI, becomes increasingly influential on the signal shape. With a single-ended readout of the crystal, the light propagation time profile contains two principal modes. The earliest main wave of photons is going directly towards the photodetector, while a second main wave traveling in the opposite direction is reflected back before reaching the photodetector later in time and thus follows a longer travel path. In addition to these two main modes are the photons undergoing total internal reflections on the crystal side faces or being reflected on an external reflector. These photons typically arrive at the photosensor either between the two main waves or after the second main wave. Scintillation transport time profiles as a function of the DOI were presented in many studies , Cates et al 2015, ter Weele et al 2015, Toussaint et al 2019.
A key element having the potential to be exploited is that DOI information is naturally encoded in the time difference between the two main photon waves. The backreflected photons can arrive closely in time with the direct photons for events with DOI near the crystal entrance face, but will arrive later for events with deep DOI since they have to propagate almost twice the full crystal length. However, there are other processes such as the non-instantaneous scintillation emission and photodetection that disperse the photons in time, therefore possibly squandering the DOI information contained in the waves. The continuous progress on timing at all levels of the detection readout chain could steadily strengthen time-based DOI estimation approaches in a near future due to the reduced time jitter contribution of the detector components. In return, an enhanced DOI estimation will enable progress on the timing precision of TOF-PET detectors. Guided by this joint connection between timing and DOI estimation, the present paper is divided in two main steps.
First, taking advantage of recent scintillation detectors (LSO:Ce:Ca crystals read out by FBK NUV-HD SiPMs) and fast electronics that now provide sufficient time resolution, this paper demonstrates experimentally the possibility of measuring a CTR variation along the DOI induced by the DOI-dependent signal shape. It will be shown that this cannot be attributed only to the commonly known effect of the modulation of the light collection efficiency with respect to DOI and crystal length. In addition, as a proof of concept for a time-based DOI estimation, we confirm in this work the ability to coarsely estimate the DOI using a double-threshold scheme on the analog timing signal. The precision of time-based DOI estimation is however limited by the current detector technology. As a second step, simulation studies were carried out to further explore the ability of future detectors to reach high-resolution time-based DOI estimation. Four time-based estimators to exploit either generic or more precise features of the DOI-dependent timestamp distribution detected by ultra-fast photodetectors were investigated, in anticipation of the potential of next-generation detectors with improved performance to achieve enhanced CTR from a DOI time bias correction.

Material and methods
2.1. Experimental proof of concept with state-of-the-art detectors This first part of the study was to evaluate the DOI influence and time-based estimation in currently available detectors. The CTR behavior as a function of DOI using state-of-the-art scintillation detectors was experimentally measured to show an observable DOI effect on CTR and we investigated the feasibility of extracting the DOI from the timing signal.

Effect of the DOI on CTR
It was demonstrated that a CTR of 98 ps FWHM @511 keV was achievable with two 20 mm long LSO:Ce:Ca coincident crystals using 4 × 4 mm 2 FBK NUV-HD SiPMs having 38 ps SPTR (sigma) read out by fast highfrequency (HF) electronics (Cates et al 2018. This constitutes the state-of-the-art CTR for LSO-type crystals at the moment of writing. The HF-readout has the distinct advantage to preserve time structures in the signal rising edge to a greater extent as compared to more standard readout with lower bandwidth electronics. In the present work, using the same SiPMs and HF-electronic readout, a DOI collimation along a 2 × 2 × 20 mm 3 LSO:Ce:0.2%Ca crystal was performed. This crystal was in coincidence with a 2 × 2 × 3 mm 3 LSO:Ce:0.4%Ca reference crystal. Both crystals were wrapped in Teflon and coupled to the SiPMs with Meltmount glue (n = 1.582). Measurements at different DOIs were performed using electronic collimation. We placed a 22 Na positron source (1 mm diameter) close to the long crystal and the reference detector far away, about five times the distance of the source and long crystal. The DOI positioning was then made with a micrometer manual translation stage, moving the long crystal along its axis, at six positions with ∼4 mm steps from 0 to 20 mm (extreme edges still reaching into the crystal, respectively far and near of the SiPM). The setup had a DOI collimation of 2-3 mm. The time signal was provided by the HF-electronics and the energy signal by an analog operational amplifier. The electronic signals were digitized by a LeCroy DDA735Zi oscilloscope (3.5 GHz bandwidth) using a sampling rate of 40 Gs s −1 reduced to 20 Gs s −1 with the four channels (energy and time in both detectors). A Gaussian fit was applied on the photopeak in the energy spectrum of both crystals, and only photoelectric events between −1.5σ and +2σ were kept (approximate energy window of 490-540 keV), where σ is the standard deviation of the fit. The CTR at each DOI irradiation position was evaluated using a Gaussian fit on the coincidence timestamp distribution obtained from leading edge discrimination with a threshold at 20 mV (the single SPAD signal was 44 mV). The leading edge threshold was set on the oscilloscope calculating the signal crossing time via linear interpolation. This 20 mV threshold provided the best CTR for all measurements, as also seen in Gundacker et al (2019).

DOI discrimination with a double-threshold approach
A second higher threshold in the rising edge signal (350 mV) was applied to extract a second value of the timing signal crossing time in the DOI-collimated detector. We define t ref as the crossing time of the signal at the firstlower threshold in the reference detector. Two delays t 1 and t 2 were then evaluated: t 1 corresponding to the time difference between the crossing time at the first threshold in the collimated detector and t ref , and t 2 corresponding to the time difference between the crossing time at the second threshold in the collimated detector and t ref . A time difference Δt 2−1 = t 2 − t 1 was then evaluated and since this is a differential measurement, the radioactive source position between the two detectors does not affect this time difference. If the DOI effect on the signal is sufficiently strong, the hypothesis is that Δt 2−1 should increase with increasing DOI position, since the two main waves should be more separated in time, thus enabling a possible DOI estimation using this time difference. The second threshold at 350 mV was found sufficient to capture nearly all time differences larger than the theoretical largest possible time difference between the onsets of the two main waves given by scintillation propagating at speed c/n twice the crystal length for an deep DOI near the photosensor. The second main wave contribution to the signal would therefore be included at this threshold level. Throughout this paper, the DOI is defined as the interaction position of the annihilation photon along the length of the crystal with larger DOI values corresponding to positions closer to the photodetector side.

Simulation study with time-based DOI estimators
The second part of the study was aimed towards providing, from simulations assuming digital-like SiPM readout, general guidelines on the main critical characteristics of future detectors required to provide efficient time-based DOI estimation. In the following, we describe the simulation model and parameters, present four DOI estimation methods and show the correction method of CTR using the DOI estimation. We also detail the methodology to define the theoretical limit on time-based DOI estimation with single-ended readout.

Simulation model and parameters
Simulations were performed using the analytical timing model described in Toussaint et al (2019) which deals separately with the light propagation time profiles as a function of DOI to keep separated the variance-related effects (emission, transport and photodetection) and the bias induced by DOI, then combined in a root mean squared error to describe the global timing error. The influence of DOI bias with detectors affected by a strong DOI impact was shown experimentally in Loignon-Houle et al (2020), validating the model presented in Toussaint et al (2019).
Probability density functions (PDFs) of the physical phenomena governing the detection chain were modeled. First, the scintillation emission PDF was modeled as a single bi-exponential function with a rise time constant and decay time constant, denoted τ rise and τ decay respectively. Combinations of a rise time of 1 or 10 ps and a decay time of 100 ps, 1 ns, 10 ns, or 30 ns were studied. Second, the transport time PDFs at each 1 mm DOI step along a polished 20 mm crystal (n = 1.82) with 98% reflectivity specular reflectors were obtained with Geant4 (Agostinelli et al 2003, Allison et al 2006. Third, the SiPM response PDF was modeled as a Gaussian function with standard deviation given by the SiPM SPTR. SPTR values of 10 ps and 30 ps (sigma) were evaluated. The emission and SiPM response PDFs were assumed the same for all DOIs. Finally, the PDFs of three processes were numerically convolved for each DOI step. Scintillation photons were then extracted from these PDFs with a quantity based on the light output given by the product of the light yield (LY), the light transfer efficiency (LTE) and the SiPM photon detection efficiency (PDE). The photons were then ordered in time to provide final timestamps where a multi-digital SiPM was assumed so that every trigger comes with a digital timestamp. LY values of 10 000, 20 000 or 40 000 photons MeV −1 were simulated. The LTE was obtained for each DOI from the Geant4 simulations and the PDE was kept fixed at 55%, an achievable value by SiPMs (Otte et al 2017).
The effective initial photon time density, denoted ρ ph in this paper and sometimes referred more simply as photon time density, was used to incorporate the rise time, decay time and LY into a metric on which the CTR typically strongly depends . It was defined as LY/(τ decay · τ rise ) and thus has the units of photons · MeV −1 · ps −2 , considering the decay time in picoseconds instead of the more conventional nanoseconds. Four different DOI estimation methods were tested and are detailed in the following.

Single time difference method
A first DOI estimation method based on the time difference between two timestamps was explored. The DOI estimation procedure for this single time difference method is illustrated on the left of figure 1. Since the detection signal contains two main photon waves in a single-ended readout as discussed in section 1, using a detected photon from the first wave and comparing its detection time to a photon from the second wave could in principle uncover the DOI. An interaction deep in the crystal (near the SiPM) should have a larger time difference between the two timestamps than an event at a shallow DOI since in the latter case, the two main waves are mostly merged in time.
We consider a crystal virtually divided in z DOI sections. Following an annihilation photon interaction in the crystal, the n first ordered scintillation photon timestamps of all N detected scintillation photons are collected. Two timestamps are used from the n timestamps dataset. The first timestamp t early is chosen among the first detected photons since they typically have the lowest statistical time variability (Gundacker et al 2015, Toussaint et al 2019. The second timestamp t late is chosen so that there is ideally the strongest DOI dependence on the time difference between the two timestamps. For this simulation study, one can assess all possible timestamp order combinations and compute the resulting DOI correction efficiency which enables the best DOI-corrected timing. During a calibration step where the DOI is known, the average of t late − t early of M detected annihilation photons for each DOI d of the z DOI sections is stored in a reference look-up table (LUT): where the brackets represents the average over the M events. The calibration procedure in this study used 10 000 photoelectric events per DOI (20 DOI sections, 200 000 events total) for each studied detector type defined as a combination of ρ ph and SPTR. Then, the DOIs of coincidence events is estimated. For every detector type, 40 000 photoelectric events were simulated (20 000 coincidences) having DOIs randomly chosen according to the 511 keV attenuation probability. The time differences t late − t early were evaluated for these events and the DOI of each event was estimated from the minimum absolute value of the time difference between the measurement and the reference LUT:ˆ( The term DOI estimation is loosely used instead of the term DOI classification although no smooth DOI function was obtained during the calibration which used 1 mm DOI steps.

Multiple time difference method
Another DOI estimation method, the multiple time difference method, was tested (see middle of figure 1). As a calibration step, the DOI d in the crystal is fixed (among the z DOI sections) and the n first timestamps for M detected annihilation photons are measured. For this fixed DOI, a starting photon order s 2 and a photon order step δ 1 are chosen, and every time difference of timestamps separated by an order step δ is calculated. For example, if s = 3 and δ = 1, then starting from the 3rd photon, the time differences t 4 − t 3 , t 5 − t 4 , t 6 − t 5 , ... are calculated until photon order n is reached. The averages of these time differences over the M measurements are then calculated. The standard deviation of these time differences over the M measurements, a vector denoted This calibration step creates a matrix of the average time differences for all DOIs in the crystal: Then follows the measurements without DOI information. For a given interaction, the n first timestamps are collected and a vector of the time differences using the starting order s and step δ is computed: The minimum mean squared difference between the measurement and the reference matrix is evaluated to estimate the DOI:ˆ( where the time differences are weighted by their average standard deviation to give more weight to differences with lower statistical fluctuations.

Kernel density estimation (KDE) method
A third DOI estimation method used the KDE technique. The estimation procedure is illustrated on the right of figure 1. The KDE is a non-parametric method used to estimate the underlying PDF of an observed finite dataset (Parzen 1962). A kernel density estimator of the PDF is also consistent, that is the estimator converges in probability towards the PDF for a large number of observations. Since the underlying PDF can be smooth, the KDE approach, compared to a binned histogram of the data, has the advantage of producing a continuous estimate of the PDF. In our case, the KDE method was used to estimate the PDF of a given ordered set of timestamps and consequently to estimate the DOI since there is a specific PDF for each DOI. The estimation procedure starts by constructing reference KDEs. The analytical PDFs, which the KDEs should approximate, would not be easily available in practice. To circumvent this, at each DOI, the average of KDE curves obtained from multiple photoelectric events is used. Using this approach presents the advantage of choosing a minimum trigger order other than the first photon. It is not easily possible to cut the analytical PDFs to simulate low-order photons discarding, whereas reference KDEs can be built from any subsets of triggers.
All the method parameters are adapted for each studied detector. First, the number of triggers used from the complete timestamps dataset of an event is chosen. Typically, a larger number of triggers is required for high-ρ ph detectors since the second main photon wave comes at higher trigger orders because of the high initial photon time density. The dataset is rescaled in time so that the first chosen trigger (not necessarily the first detected photon) is at t = 0. The fifth trigger was chosen since it provided the best performance. This rescaling ensures that the time-of-flight information is removed from the simulation and does not influence the DOI estimation. Then, a Gaussian kernel is used and its bandwidth (the standard deviation in the case of a Gaussian kernel) is chosen. A too small bandwidth might overfit statistical fluctuations whereas a too large bandwidth might miss relevant DOI-dependent features of the timestamp distribution. The optimal bandwidth might vary as a function of DOI but since in practice the DOI information is not available, a single value of bandwidth throughout all DOI events is used. The Gaussian kernel is applied on each kept timestamp and a KDE is calculated. The KDE parameters for each detector type were chosen by optimizing the DOI-corrected timing using small test datasets.
For each DOI d, a thousand KDE curves are averaged together, giving 20 reference KDE curves (one for each 1 mm DOI section): where M = 1000 in the simulations. A higher value of M did not give significantly more stable curves. Then follows the DOI estimation using 20 000 new sets of timestamps for each detector type. The same parameters, depending on the detector type, used for the reference KDEs are applied for the calculation of the KDE curves of the new datasets. For each timestamp set, a KDE curve is computed: , , , 7 n meas 1 2 Q = ¼ and a mean squared difference with all the DOI-wise reference KDE curves is calculated and the DOI is estimated from the lowest obtained value:ˆ( . A thorough mathematical description of the KDE technique is found in Fan and Gijbels (1996).

Neural network estimation method
A fourth estimation method based on a neural network was also tested, although not as extensively as the previous three methods because of the vast range of possible network architectures that could constitute a separate study. The goal was to use a fairly simple architecture without excessive fine-tuning to assess a possible DOI estimation capability with a classification neural network. The DOI-classification network was implemented using the PyTorch v1.3 library (Paszke et al 2019) with the following steps. For each DOI, 10 3 timestamp sets were used (2 × 10 5 in total) and 80% of the data was randomly chosen as the training dataset while the remaining 20% was used as the validation dataset. The timestamps of all sets were first rescaled in time to have the fifth timestamp at t = 0, then rescaled to have values closer to zero to make the training easier (LeCun et al 2012). The first rescaling ensures that the DOI learning is made only from the scintillation signal time distribution and that the time offset related to the radioactive source position is excluded. Since the timestamps of a set are ordered in time, the second rescaling was done by subtracting the set mean from every timestamp and dividing by the last timestamp of the set. The rescaled timestamps were then fed into an input layer followed by a single hidden layer with a rectifying linear unit activation function. The output layer of the network provided a DOI classification with 1 mm stepping (size of 20). True DOI values were joined with their corresponding timestamp set during the training. During training, mini-batches of 200-500 sets of timestamps were used depending on the detector type. We used a stochastic gradient descent optimization algorithm with a crossentropy loss function. The learning rate was fixed at 10 −3 with momentum of 0.9 and the training was done on 10 000 epochs. The number of detected photons in each set and the input and hidden layer sizes were adapted for each detector type. They were chosen by roughly optimizing the DOI-corrected timing using small test datasets before performing the complete training. High-ρ ph detectors typically required more input timestamps and a larger hidden layer size.
2.2.6. Timing correction from DOI estimation At this point, the DOI in both detectors for all interactions was estimated using one of the methods presented above. The next step is to perform a timing correction using these DOI estimations to improve the CTR. Before performing any correction, a DOI time bias LUT must be filled since we work with DOI classification. For each DOI combination between two coincident detectors, the coincidence time spectrum was assessed and its mean value (time bias) was stored in the LUT. This step can be done simultaneously with the reference DOI estimator assessment of the estimation methods. The timestamps forming the coincidence time spectrum were obtained from an average timing estimator (Gundacker et al 2015) with the best value using the first n primary triggers with n ranging from 1 to 30. The optimal depth order of a timing estimator varies depending on the detector characteristics such as ρ ph and SPTR (Gundacker et al 2016a, Toussaint et al 2019 as well as the DOI estimation method. Since the coincidence events are DOI-symmetric, it is not necessary to assess all m × m DOI combinations, but rather only m(m − 1)/2 combinations. For an even more practical LUT filling, a small-length crystal in coincidence with a DOI-collimated crystal could be used to measure only m time biases and the time biases of all DOI combinations of two long crystals could then be retrieved by assuming symmetric (negative) time biases that would have been measured in another coincident DOI-collimated long crystal. The time bias associated with the estimated DOI in detector A and detector B is retrieved from the LUT and the corresponding value is subtracted from the coincidence time: The precision on the estimated time bias is therefore directly dependent on the precision of the DOI estimation. The same estimation procedure was repeated for all M measured coincidences giving a complete DOI-corrected dataset of coincidence timestamps with which the FWHM CTR, obtained from the average timing estimator, was evaluated from the standard deviation of a Gaussian fit on the dataset distribution: corr.
s =´D 2.2.7. DOI positioning error Assuming a dataset with the same number of events per DOI, a positioning error (ˆ) d e was calculated as the absolute value of the distance (in mm) between the true DOI d * of an event and its estimated DOId, weighted by the true DOI probability P(d * ): where M is the number of events per DOI and z is the total number of DOIs (each DOI is indexed by i). The mean DOI error was weighted by the attenuation probability considering facing crystals since its the most probable and used case for timing optimization in scintillation detectors. The DOI positioning error is directly indicative of the efficiency of reducing the parallax effect, but the main aim of the current paper is to evaluate the effect of DOI correction methods on the CTR. Therefore, the parameters of the DOI estimation methods were optimized to achieve the best DOI-corrected CTR instead of focusing on minimizing the DOI positioning error.

Cramér-Rao lower bound (CRLB) on DOI positioning
A CRLB on DOI positioning was developed to assess the theoretical limit of a time-based DOI estimation. We define ( | ) p t d , tot q as the scintillation signal PDF, obtained from the convolution of the scintillation emission, transport at DOI d and photodetector response. The travel time of the annihilation photon from its source until the crystal entrance face is defined as θ. A DOI-wise offset is applied on ( | ) p t d , tot q so that its onset starts from the earliest possible detection time of an instantaneously-emitted scintillation photon having the shortest path time to the photosensor. This ensures that the DOI-wise PDFs all start at the same time to correctly represent the information available in practice where the first trigger alone should provide no information on d. The PDF ( | ) p t d , tot q is therefore dependent on d and θ, but the observed data reside purely on the time axis. Also, θ applies the same translation on ( | ) p t d , tot q for all DOIs and thus has no affect on the Fisher information. We can thus arbitrarily fix θ and define The Fisher information about d is defined as: assuming that the function giving the correspondence between d and ( | ) p t d tot is differentiable (regular variation along d), which is accurate with the small DOI steps used for this CRLB study. The CRLB on DOI positioning of a single detector is given by: withd designing the best unbiased estimator of d for ( | ) p t d tot and N d being the expected number of scintillation photons detected from an annihilation photon interacting at DOI d. This approach is different from the standard usage of the Fisher information to evaluate the lower bound on timing resolution of scintillation detectors (Seifert et al 2012, Cates et al 2015. Instead of performing the derivative along θ, it was performed along the DOI dimension for every time step (equation (12)). This formulation enables the calculation of the lower bound on variance (in units of distance, not of time) on DOI positioning. The light transport profile as a function of DOI, giving the DOI dependence of ( | ) p t d , tot q , was obtained with the analytical light propagation model detailed in Cates et al (2015), which was shown to have very good correspondence with Monte Carlo simulations. It also enables faster computation time, which becomes more crucial as the DOI-dependent PDF must be extracted with fine precision to better compute its Fisher information. DOI steps of 0.05-0.1 mm (depending on the detector type) were used.
The CTR obtained after a DOI correction having an error based on the CRLB limit can then be evaluated with the following steps. For every DOI combination between two coincident detectors, 10 000 coincidences were assessed. The DOI in both detectors for every coincidence was simulated using a random number extracted from a normal distribution having a variance based on equation (13) and a mean centered on the true DOI. Finally, for each coincidence, the DOI-induced time bias corresponding to the combination of estimated DOIs in both detectors was removed from the coincidence time estimate, as detailed in section 2.2.6. The CTR of the DOI-corrected coincidence time estimates distribution was then evaluated (with DOI events weighted by the annihilation interaction probability), giving the CTR achievable with a DOI correction method mimicking the CRLB in performance.

Results
3.1. Experimental proof of concept with state-of-the-art detectors 3.1.1. Effect of the DOI on CTR Figure 2 shows the integrated energy signal (related to the light output) and the measured CTR as a function of the DOI using the experimental setup and detectors described in section 2.1. If only the light output (signal integral) effect on CTR is considered, the DOI-wise CTR should scale inversely with the square root of the light output N(d), i.e. CTR ( ) N d 1 , where d is the DOI. This expected CTR (based on the CTR and light output at the first DOI) as a function of the DOI is also plotted in figure 2. Higher measured CTR values are observed for DOI positions closer to the crystal center, although a light output increase is observed with respect to the DOI. The measured CTR values are therefore inconsistent with the expected behavior typically strongly governed by light output. The contribution of the light collection effects on timing was established long ago (Cocchi and Rota 1967, Bengtson and Moszyński 1970, Yeom et al 2013. Here, the observed non-monotonic DOI dependence of CTR performance, not apparent before with slower electronics (Brown et al 2014), now becomes observable with the newly improved readout electronics and motivates investigation of a time-based DOI estimation.

DOI discrimination with a double-threshold approach
The procedure described in section 2.1 was then used to extract the time differences, Δt 2−1 = t 2 − t 1 , from the crossing times at the two aforementioned thresholds on the DOI-collimated detector. Figure 3(a) displays a scatter plot of Δt 2−1 as a function of the integrated charge of the signal. Different blobs of points correspond to different DOIs, from small to large DOIs going from left to right on the horizontal axis (integrated charge). Interestingly, the Δt 2−1 distribution shifts towards higher values for the first DOIs (from 0 to 12 mm) even though the light output was increasing as a function of DOI, indicating a DOI-related light transport effect disturbing the timing signals shape. However, the Δt 2−1 distribution shifts towards lower values for the larger DOIs (from 12 to 20 mm), indicating that the light output then dominates the signal shape and thus the measured Δt 2−1 . Indeed, a higher light output generally causes a steeper signal, reducing the delay time. The overall DOI trend is thus non-monotonic, squandering the DOI estimation potential using Δt 2−1 due to possible intermixing of small and large DOI positioning. A time-amplitude correction was performed to retrieve a monotonic evolution of Δt 2−1 . A linear function was fitted through the scatter plot of Δt 2−1 against the integrated charge of the first and last DOIs, giving a negative slope. The resulting fit was then mirrored with higher slope (changing the sign of the slope and multiplying by a scaling factor, see correction curve in figure 3(a)) for better separation of the different DOI regions by exploiting the signal intensity information. Correcting each event with the parameters of the correction curve yielded the distributions seen on the upper part of figure 3(a) leading to a monotonic increase of Δt 2−1 . The choice of the slope extension was done by sweeping the scaling factor until there was minimal overlap between the distributions. A higher slope gave no significant modification of DOI resolution. Indeed, as soon as a monotonic behavior is achieved, separating even more the regions is not useful as intra-region separation of the events occurs at the same time, so the overlap between the adjacent regions remains of the same order. Even for the first shallowest DOIs, extending the slope is less useful since there is already a small monotonic behavior. We also tested some non-linear fits which gave no significant difference against a simple linear fit. Figure 3(b) displays the measured and corrected Δt 2−1 distribution for all DOIs, showing a passage from non-monotonic to monotonic evolution. An approximate  DOI resolution was computed using the corrected histograms. The average FWHM of the histograms was divided by the slope of the centroids calculated as the ratio of the time difference between the extreme centroids (at DOI 0 and 20 mm) and the crystal length, giving a DOI resolution of 8.5 ± 0.4 mm.
Since there are still overlaps between the corrected histograms, a coarser discretization grouping all first three DOIs together and the last three in a second set (dividing the crystal in two 10 mm regions) was analyzed to evaluate the DOI-positioning accuracy in two DOI regions instead of six. The same time-amplitude correction presented in figure 3(a) was performed. The resulting Δt 2−1 distributions for the two regions are shown in figure 4. The overlap between the two distributions is ∼10% of the total counts, representing a classification accuracy of ∼90%. Good positioning accuracy using the timing signal shape and intensity can therefore now be achieved due to faster electronics and excellent SPTR of FBK-NUV SiPMs.

Simulation study with time-based DOI estimators
Having now evaluated the DOI effect on CTR and possibilities of time-based estimation with current detectors in section 3.1, the following simulation study is aimed towards exploring future detectors with better photon time density and SPTR in order to give a guideline for future research efforts. Figure 5 displays the DOI-wise PDFs of a detector having nearly ideal properties and a detector having the current state-of-the-art properties of an LSO:Ce:Ca and SiPM pair. A clear DOI dependence is observed for the nearly ideal detector and to a lesser extent for the current LSO:Ce:Ca/SiPM detector. For both detectors, the two main photon waves are mostly merged in time for smaller DOIs and the time separation between the waves increases as a function of DOI. Depending on the sampling (light output) of these PDFs, one may find DOI information from the timestamp distribution.   Figure 6(a) shows the average of the time difference (see the single time difference method in section 2.2.2) between a range of late triggers and an early trigger depending on the DOI for two detector types. The early trigger was chosen to be the fifth timestamp since it led to the best performance. The figure also displays the time difference distributions obtained at the trigger order where maximal DOI separation is measured (see rotated distributions). Overlaps between adjacent time difference distributions are more pronounced for the slower detector (lower photon time density ρ ph ) and a clearer DOI distinction is achieved with the higher-ρ ph detector, illustrating that ρ ph is correlated with better DOI classification. High-ρ ph detectors however require a higher order for the later trigger since a possible large quantity of detected photons comes from the first main wave. Figure 6. (a) Average (full line) and FWHM error (shaded regions) of the time difference evolution between the late trigger (abscissa) and the early trigger chosen to be the fifth here. The distributions of the time difference at the trigger order where maximal DOI separation is measured are also displayed on the right. Two detectors, one having a ρ ph of 0.4 ph MeV −1 ps −2 and 30 ps SPTR, and the other a ρ ph of 10 ph MeV −1 ps −2 and 10 ps SPTR are displayed. (b) Average time difference as a function of DOI for the two detector types at their best trigger order. These curves are used as a DOI estimator reference (single time difference method) for the corresponding detector. The standard deviation of the time differences is also shown. Figure 6(b) shows the average time difference between the later trigger order which gives the best DOI separation (see figure 6(a)) against the fifth trigger order as a function of the DOI for two detector types. The faster detector shows a steeper DOI dependence with smaller errors and should therefore provide better DOI classification. The curves of figure 6(b) directly represent the DOI estimator reference used with the single time difference method with values to be used in equation (1). The same examples for the multiple time difference method are not as straightforward to display due to the matrix nature of the estimator parameters. With this method, choosing few of well-separated timestamps to calculate the time differences provided enhanced positioning precision instead of using a lot of photons close in time. For instance, with a 1 ns τ decay detector, using only four timestamps with a large step of 200 trigger orders is already sufficient to give slightly lower DOI error compared to using many timestamps with a finer step of ten trigger orders. Examples of KDE curves are shown in figure 7. At first glance, it seems that the KDE can efficiently retrieve the true underlying PDF. Fluctuations of the KDE curves are present due to the finite number of timestamps, similar to histogram bins discrete fluctuations, and to the chosen kernel bandwidth. Figure 8 displays the estimated versus true DOI intensity map obtained using the four DOI estimation methods for detectors having a 10 ps SPTR and a τ decay either of 1 or 10 ns. The same intensity map for an estimator having a DOI error based on the CRLB limit is also displayed. The mean DOI error (weighted by the DOI probability) is displayed in the top-left corner of each heatmap. A higher intensity near the diagonal is observed for all methods showing their ability to provide unbiased DOI positioning. The mean DOI positioning error appears to be rather similar between the four methods, but the multiple time difference and KDE methods provide the best estimation for the 1 ns τ decay detector, and the KDE method is the most effective for the 10 ns τ decay detector. The intensity profile in the CRLB heatmaps is more closely concentrated around the diagonal compared to the four methods which have a little more than twice higher mean DOI errors. Nonetheless, this higher error is only in the range of 0.4 to 1.6 mm compared to the CRLB case. The heatmaps show an intensity Figure 7. Kernel density estimation (orange curves) of the underlying PDF (blue distributions) using a timestamp dataset example (green dots) for a detector with ρ ph of 0.4 ph MeV −1 ps −2 and 10 ps SPTR is shown for three DOIs (1, 10, 19 mm). excess at 0 and 19 mm estimated DOIs with the 10 ns τ decay detector. This could probably be explained in part in light of figure 6(a) (right-hand side rotated distributions) on which the estimation is subjected to misclassification due to the overlap between the estimator value distributions. When a time-based DOI estimation is calculated with either of the four methods and the DOI estimator closest value is determined, boundary conditions imply that the closest DOI class for the estimation either smaller or larger than the DOI estimator boundary values (at 0 and 19 mm) are always classified to the closest DOI candidate which are either one of the edge DOIs. These events, estimated outside the physical crystal, could also be discarded. Figure 9 shows the DOI-corrected CTR as a function of ρ ph using the four DOI estimation methods and a DOI estimator having an error determined by the CRLB, for two SiPM SPTR values. The CTR of the uncorrected DOI and known DOI (using only the signal variance with zero bias) is also displayed. First, as seen from the uncorrected curves, even with high ρ ph the CTR saturates around 60-70 ps because of the effect of DOI variability on timing. On the contrary, the known DOI curves exhibit strong CTR improvements as a function of ρ ph with a nearly linear behavior on a log-log scale (fitted slope near −0.5). Thus, the CTR is approximately improving, as expected, with the inverse square root of ρ ph , that is

CTR improvements
. The DOI-corrected CTR curves are close to the known DOI curves, with enhanced CTR improvements for higher ρ ph . The photon time density is of utmost importance to enhance the time-based DOI correction efficiency on CTR. While varying the detector parameters of the ρ ph metric individually, decreasing the decay time had a strong positive impact on DOI correction. At later times in the scintillation process, a mixing of two statistics occur since the photons that have been reflected are indistinguishable from those that went directly to the SiPM. Therefore, the second main wave, intrinsically with a good time resolution, is now mixed with photons going directly towards the SiPM but that were emitted later in the scintillation emission process with larger time jitter. This conclusion is relevant for classic scintillators having τ decay ? τ rise , but is not as valid with faster emitters of prompt photons. By comparing the left/right plots of figure 9, the 10 ps SPTR case achieves better DOIcorrected CTR, reaching ∼10 ps with the DOI estimation methods, whereas the 30 ps SPTR case would require even higher ρ ph to fully reach 10 ps CTR.
The four DOI estimation methods converge to the known DOI case and follow a similar trend as ρ ph increases. The simple single time difference method provides effective DOI correction and the multiple time difference method offers some further DOI-corrected CTR improvement. The KDE method, which harvests more information from the timestamp distribution, overall achieves the best DOI-corrected CTR. This method can still suffers from a low-sampling of the PDF it aims to estimate since the number of detected photons remains limited, explaining in part why the method does not achieve the optimal performance predicted by the CRLB limit. The single and multiple time difference methods use few but optimally chosen timestamps. The KDE method uses a larger set of timestamps even though some orders suffer from higher variance, especially in the decay part of the scintillation kinematics. The DOI estimation method based on a neural network does not provide the best CTR improvement, but follows the same trend as the other methods even though no thorough optimization on the network architecture was done. We could expect better performance with deeper networks able to catch more features of the timestamp distributions. Finally, the CRLB curves in figure 9 establish what the best time-based DOI estimator could provide in CTR improvement. The CRLB curve for both SPTR values Figure 9. CTR as a function of ρ ph using a DOI correction from the four estimation methods for SPTR of 10 ps (left) and 30 ps (right). The rise time value was fixed at 10 ps. The starting trigger is the fifth for all methods. The CTR curve obtained from an estimator having an error determined by the CRLB on DOI positioning is also displayed. The CTR shown is the best obtained with the average timing estimator of order one to thirty. The CTR curves of the uncorrected DOI and known DOI (using only the signal variance with zero bias) are also displayed. Note the log-log scale (base 10). appears below all estimation methods and the KDE method approaches the CRLB curve the most closely. The CRLB curve remains slightly above the known DOI curve since there always remain a non-perfect separability between the DOI-wise PDFs and the sampling is finite (N d in equation (13)).
As it might be challenging to precisely extract a particular timestamp order among a detection dataset, we can assess how the choice of the two triggers for the single time difference method influences its efficiency. Figure 10 shows heatmaps of the DOI-corrected CTR of two detectors as a function of the triggers t early and t late . The heatmaps show fairly stable CTR for a broad range of the two triggers. Hence, the exact choice of the triggers orders is not absolutely paramount to have the best possible DOI-corrected CTR, allowing some leeway in the choice without highly degrading the DOI correction efficiency.
4. Discussion 4.1. Experimental proof of concept with state-of-the-art detectors A DOI-dependence of the CTR measured with LSO:Ce:Ca crystals, FBK NUV-HD SiPMs and HF electronics was experimentally shown in this work. A poorer CTR for DOIs closer to the center of the crystal was observed (figure 2). One explanation of this observation might be related to the HF electronics that provides higher bandwidth which keeps the two main photon waves distinguishable and separated. For small DOIs, the two main waves are mostly merged in time bringing a high photon time density at the photosensor resulting in a good CTR. For large DOIs, the waves are well separated in time and the light output was significantly higher (+23% for the larger DOI compared to the smaller DOI), providing a high photon time density within the first main wave. For central DOIs, there might be a mixing of two statistics when the second main wave arrives at the SiPM. When the time difference between the direct and reflected waves is of the same order as the SPTR, then the SPTR smearing leads to a mixing of the photons of the two waves, which is the cause for the measured higher CTR values in the middle of the crystal.
A DOI estimation using the time difference between two thresholds on the analog signal was also investigated. Some DOI dependence was observed for the first few DOIs (Figure 3) although squandered by the higher light output at the deeper DOIs (figure 2) which dominated the signal shape. Using a signal amplitude correction provided a monotonic DOI dependence, enabling a more valuable DOI estimation. The same procedure applied to two coarser DOI regions (top versus bottom half of the crystal) showed a good DOI classification of 90% (figure 4). It could therefore be possible to artificially divide the single crystal in two regions like a phoswich with a single-ended readout. Using co-doped LSO crystals with smaller calcium content for lower light self-absorption as shown in Gundacker et al (2016b) could be used to mitigate the strong light output DOI dependence. Crystals with >20 mm length could also be tested to provide better resolvability of the light transport features. With the detector used in this work, the second main wave has an intrinsic good time precision (the second main peak of the photon travel spread (PTS) profile has a rather high intensity within a short lapse of time) but can be spoiled with the photons from the first main wave, which is already later in its scintillation emission where the photons have higher jitter. Overall, this experimental study exhibited the increasing influence of the DOI effect on CTR with fast scintillation detectors, and motivated the measurement of the DOI by exploiting the very fast timing signal since it naturally encodes DOI information. A comprehensive optimization work could therefore follow the current proof of concept study by using different crystal types and geometries. A more thorough assessment of the threshold requirement and necessity of a combined time-energy Figure 10. Heatmap of the DOI-corrected CTR as a function of the triggers choices t early and t late of the single time difference method for detectors having a ρ ph of 4 (left) and 0.4 ph MeV −1 ps −2 (right). The SPTR, τ rise and light yield were respectively fixed at 10 ps, 10 ps and 40 000 ph MeV −1 . Contour lines are shown to better appreciate the CTR evolution as a function of both trigger choices. Note the logarithmic color scale.
(integrated charge) metric to gain a monotonic DOI estimator behavior could also be done. Further exploitation of time-based concepts could also be applied in a digital approach with next-generation 3D digital SiPMs (Nolet et al 2016, Roy et al 2017.

Simulation study with time-based DOI estimators
It is already possible to reach few tens of picoseconds CTR with short crystals (3 mm length) coupled to a FBK NUV-HD SiPM, e.g. 58 ps with LSO:Ce:0.4%Ca, 51 ps with BaF 2 , and 35 ps with BC422 (Gundacker et al 2020). The DOI variability limits longer crystals, essential to provide sensitivity for PET, to achieve such values. The DOI estimation methods tested by simulations in this work present inherent advantages and difficulties.

DOI estimation methods
In spite of the simplicity of the single time difference method, needing only two triggers among all the acquired timestamps, this method provided a steady CTR improvement with ρ ph (figure 9). It was also shown to be fairly robust with only minor CTR degradation related to a variation in the choice of the two triggers (figure 10). As it might be difficult to precisely catch a given trigger order in practice, the observed stability strengthens the method practicality. Also, averaging around the chosen orders, e.g. with three triggers around the late and early triggers 〈t late±3 〉 − 〈t early±3 〉, produced slightly better DOI correction. Test cases with ρ ph of 0.4 and 4 ph MeV −1 ps −2 and 10 ps or 30 ps SPTR had 1%-4% improvement in DOI-corrected CTR when incorporating this averaging.
The multiple time difference method requires only a few more timestamps to improve the CTR. Interestingly, using coarser steps achieves similar or sometimes better DOI positioning than using a fine stepping. With high-ρ ph detectors, the neighbor timestamps (of unity-step) are close to each other in time and the method loses its efficiency since it must deal with many tiny time differences that hide the timestamp distribution pattern. For example, with a detector having ρ ph of 4 ph MeV −1 ps −2 and 10 ps SPTR, using <5 separated triggers was already enough to provide better DOI estimation than the single time difference method, whereas using more triggers provided no significant gain. The multiple time difference is nonetheless more complex to implement as it requires the choice of a starting trigger order and a stepping value between the triggers, complicating the calibration phase.
The KDE method is a powerful approach based on a solid mathematical framework able to catch more information of the timestamp distribution. It typically provided the best DOI estimation efficiency and DOIcorrected CTR among the tested methods. The method was however suffering from low sampling (low light output) and slow decay time since the underlying PDF containing the DOI information is difficult to retrieve when too few timestamps dispersed in time are available. The choice of the kernel bandwidth is also known to be a complex task (Chiu 1991). The method could be refined by using the information about the timing behavior as a function of the trigger order since the first triggers typically have a lower timing variability than later triggers. Therefore, more than one kernel bandwidth could be used to incorporate this behavior in the estimation, using a smaller bandwidth for earlier triggers and a larger bandwidth for later triggers. The effect of multiple kernel bandwidths could thus be worth investigating in a study fully dedicated to the KDE method.
The fourth DOI estimation method based on a neural network was able to decode the DOI rather efficiently using the timestamp distribution and to gain a CTR improvement (figures 8 and 9). Machine learning involvement in the field of (TOF-)PET imaging and detectors is growing . DOI positioning in monolithic crystals using an algorithm based on gradient tree boosting was recently demonstrated (Müller et al 2019). Artificial neural networks were used as a timing estimator for (multi-)digital SiPMs in presence of dark counts (Lemaire et al 2020). Convolutional neural networks were also explored for TOF estimation with PMTbased scintillation detectors waveforms and demonstrated a CTR improvement by 20% compared to a leading edge discrimination (Berg and Cherry 2018). The authors mentioned a possible DOI influence on the waveform shape, which was observed in the current work using faster detectors and electronics. The present paper only focused on a proof of concept to show the possibilities of using a basic neural network for a time-based DOI estimation. The results motivate further studies to explore deeper network architectures.
Finally, the CRLB provided, as expected, the best DOI-corrected CTR values (figure 9). Finding a better DOI estimator would provide only a minor gain for the considered detectors since the tested methods are rather close to the CRLB limit, notably the KDE method. Detectors presenting a strongly structured time profile as a function of the DOI would yield better CRLB on DOI positioning through an increased Fisher information. In light of the strong synergy existing between time-based DOI estimation and intrinsic time resolution of the detectors, the standard relations found for CTR with the main characteristics were also found in the present study. The CRLB on time-based DOI estimation has a typical dependence on LY decay t . Then, an interplay exists between the PTS, SPTR and τ rise . Particularly, for low SPTR values, the CRLB has a rise t dependence and saturates at the PTS variance. These dependencies were notably highlighted in Vinogradov (2018) for CTR. Time resolution improvements consequently lead to enhanced time-based DOI estimation efficiency.

Limitations of the study and future investigations
The present work focused on achieving the best DOI-corrected CTR by only keeping photoelectric events and discarding scattered events, as it is typically done in experiments aiming at the best timing. Scattered events lead to lower rate of scintillation production therefore to a lower rate of scintillation detection which degrades the time resolution. The deposited energy in the crystal could be evaluated under similar principles using the timestamp temporal information, as shown in Therrien et al (2018). It would unfortunately appear arduous to both have a combined DOI/energy estimation based uniquely on the timestamp distribution since both the DOI and energy event-by-event variations spread the triggers in time. A low-energy/large-DOI could have similar timestamp distribution as a high-energy/small-DOI event, therefore prompting event intermixing. Methods to compensate for signal time walk error when a larger energy window is used were recently proposed (Xie et al 2020) and time-energy relations of the signal were detailed with leading-edge discriminator triggering.
Integrating a DOI component in those relations could be worth future investigation. Multi-interaction events (e.g. Compton interaction followed by a photoelectric interaction in the same crystal) also constitute a complex problem worthy of future investigation, especially for monolithic crystals. In our exploration of timing limits, noise sources such as dark counts, which can induce a less precise time estimation, were not included in the simulations and no associated skipping effect (Venialgo et al 2015) due to possible sub-one-to-one ratio of time-to-digital converter (TDC) to single-photon avalanche diode (SPAD) in the SiPM architecture was considered. Every photon measured by a SPAD was considered associated with an available TDC, but other SPAD-TDC ratios, as explored in Tétrault et al (2017), could also be considered in further studies. With the proof of concept measurements done in this work using a double-threshold method with an analog SiPM, an additional TDC channel for the recording of the second crossing time would be required. Saturation effects can occur in SiPMs depending on the number of available SPAD cells following an annihilation photon interaction. This was considered out of the scope of this paper since we focused on exploring the timing improvement limits with DOI estimation for detectors varying by a single effective metric, but saturation effects would reduce the efficiency of the methods for high-ρ ph detectors. The photon time density metric may also not be suitable for detectors producing scintillation with an additional burst of prompt photons (multiple τ rise , τ decay and LY values), thus limiting to some extent the classical interpretation of the CTR relation with the metric.
The present study considered polished crystals with highly-reflective surfaces, but surface and wrapping conditions can alter the light collection efficiency (Loignon-Houle et al 2017a, 2017b) and the light transport time profile (ter Weele et al 2015). The DOI resolution and CTR of double-readout LYSO crystals of different lengths was studied by Kang et al (2015), showing that etched crystals can provide better DOI resolution than polished crystals, but at the expense of degraded CTR. Berg et al (2015) also demonstrated that a 20 mm long crystal etched on all sides had degraded CTR compared to a polished crystal, but partially etching a fraction of one lateral face of the crystal yielded improved CTR. Pizzichemi et al (2019) have recently shown a better DOI estimation with depolished crystals although accompanied by a CTR degradation, but they were able to recover similar CTR as the one obtained with polished crystals by correcting the CTR using the DOI estimation. Future dual CTR-DOI studies in the context of ultra-fast timing could therefore include considerations on a variety of surface conditions of the detector.
Detectors with 20 mm length and 2 × 2 or 3 × 3 mm 2 cross-section were used in this work, but the study could be extended to other geometries, even to monolithic-like detectors. DOI estimation in monolithic detectors was proposed in several studies, relying on photodetector(s) attached to the block crystal to estimate the interaction position from the signal spatial distribution (Li et al 2010, Borghi et al 2016. For high aspect ratio crystals, longer lengths could enhance the DOI estimation efficiency with a higher DOI-wise time separation between the two main waves. Coarser DOI discretization could also be considered depending on the desired aimed CTR, that is considering >1 mm DOI regions for the estimation, which could also simplify the rather cumbersome DOI bias LUT acquisition.
Developments on materials and photodetectors to fully benefit from a time-based estimation approach are still needed. The simulation study gave general trends and guidelines on the required main detector characteristics, highlighting a connection between pure timing improvements and efficiency of time-based DOI estimation. The tested DOI estimation methods do not require detector alterations typically necessary for other DOI detectors, such as phoswich, reflector pattern and double-sided readout. Indeed, the signal is not altered to create a DOI dependence, but rather the intrinsic DOI-dependent timing signal shape is exploited with detectors that provide a high photon time density. Combining two or more DOI estimation methods could be possible, e.g. making use of both the time and spatial distributions of photons to refine the DOI estimation. A spatial distribution method could guide the time-based DOI estimation and restrict the range of possible DOIs.
Finally, the results showed the necessity of a high photon time density to benefit from pure-time based DOI estimators. New assemblies (metamaterials) of dense crystals with fast nanostructured scintillators are developed to strongly enhance this factor with the production of prompt photons towards the goal of 10 ps CTR (Turtos et al 2019a(Turtos et al , 2019b, a value out of reach from classical scintillators. In addition to improving the timing, progresses in such metamaterials, alongside the necessary parallel progresses on the rest of the detection chain, will have the potential to create strongly structured timing signals along the DOI, motivating the extraction of spatial positioning of the interaction from time-based approaches.

Conclusion
A fully time-based approach for DOI estimation with fast TOF-PET scintillation detectors was investigated. State-of-the-art detectors consisting of LSO:Ce:Ca crystals coupled to FBK NUV-HD SiPMs and read out by fast high-frequency electronics provided evidence of a DOI time spread effect on CTR not observable before with slower detectors. DOI discrimination using a two-threshold arrangement on the rising slope of the signal was explored experimentally in combination with the signal intensity information, which gave a monotonic DOI dependence of the estimation metric as well as a DOI classification accuracy of 90% with a 20 mm long crystal virtually divided in two 10 mm regions. A conventional single-ended readout of crystals without any need for DOI-based detector alteration can be used with this method, which takes advantage of the DOI-dependent time profile of the signal. Simulations were also undertaken to explore four pure time-based DOI estimation methods, based either on the time difference between only two timestamps, on the time difference between multiple timestamps, on the timestamp distribution, and on a neural network. The four methods were compared with the CRLB on DOI positioning using the DOI-dependent timing signal profile. Investigation of the scintillation detector design criteria for efficient DOI measurement with the estimation methods highlighted the strong importance of a higher initial photon time density and revealed that the timing response of recent SiPMs is close to be suitable to achieve efficient DOI correction with those methods. The more simple twotimestamp time difference method already provided a CTR enhancement, and incremental improvements were achieved with the more complex methods using more timestamps, even approaching the CRLB limit. With the steady progress on low time jitter detectors for TOF-PET, DOI estimation using pure time-based methods is becoming an achievable realistic feature to push the timing towards record values with long detectors necessary to preserve sensitivity in PET imaging.