Time of flight positron emission tomography towards 100ps resolution with L(Y)SO: an experimental and theoretical analysis

Scintillation crystals have a wide range of applications in detectors for high energy and medical physics. They are recquired to have not only good energy resolution, but also excellent time resolution. In medical applications, L(Y)SO crystals are commonly used for time of flight positron emission tomography (TOF-PET). This study aims at determining the experimental and theoretical limits of timing using L(Y)SO based scintillators coupled to silicon photomultipliers (SiPMs). Measurements are based on the time-over-threshold method in a coincidence setup utilizing the ultra-fast amplifier-discriminator NINO and a fast oscilloscope. Using a 2 × 2 × 3 mm3 LSO:Ce codoped 0.4% Ca crystal coupled to a commercially available SiPM (Hamamatsu S10931-050P MPPC), we achieve a coincidence time resolution (CTR) of 108±5ps FWHM measured at E=511keV. We determine the influence of the data acquisition system to 27±2ps FWHM and thus negligible as compared to the CTR. This shows that L(Y)SO scintillators coupled to SiPM photodetectors are capable of achieving very good time resolution close to the desired 100ps FWHM for TOF-PET systems. To fully understand the measured values, we developed a simulation tool in MATLAB that incorporates the timing properties of the photodetector, the scintillation properties of the crystal and the light transfer within the crystal simulated by SLITRANI. The simulations are compared with measured data in order to determine their predictive power. Finally we use this model to discuss the influence of several important parameters on the time resolution like scintillation rise- and fall time and light yield, as well as single photon time resolution (SPTR) and the detection efficiency of the SiPM. In addition we find the influence of photon travel time spread in the crystal not negligible on the CTR, even for the used 2 × 2 × 3 mm3 geometry.

1 Introduction Time-of-Flight (TOF) techniques are widely used in positron emission tomography (PET) to improve the quality of the reconstructed image [1]. They refine the initial information about the exact localization of the positron emission point and thus contribute to the rejection of background events outside the region of interest, reducing the noise in the reconstructed image. Commercial full-body PET currently achieves a coincidence time resolution (CTR) of around 500ps FWHM [2][3][4]. More advanced research solutions aim at a CTR around 200ps FWHM [5], corresponding to a zone of -1 -∼3cm around the point of emission, sufficient to remove coincidence events outside the organ of interest. An improvement of the CTR towards ∼100ps would provide even more precise determination of the point of emission. It is thus interesting to study how an improvement of the full detector chain comprising the scintillating crystal, the photodetector and the electronics contributes to achieving this goal. Silicon photomultipliers (SiPMs) or multipixel photon counters (MPPCs) are promising candidates towards achieving excellent time resolution [6,7]. With 10mm long crystals, these detectors can reach CTRs of about 200ps FWHM [8]. These values are comparable to and even better than the best values obtained with PMTs [9]. The photon travel jitter within the crystal is not a negligible factor and has to be taken into account in order to understand the system's time resolution [10,11]. With dimensions of 3 × 3 × 5 mm 3 it is for example possible to achieve CTRs of 147ps FWHM with LSO:Ce doped crystals. This value degrades to 186ps FWHM for 3 × 3 × 20 mm 3 [12]. The research performed so far highlights the strong advantages of SiPMs in terms of timing. On the other hand, their drawbacks are the high dark count rate (DCR), high temperature dependence and the need for fast and low noise electronics to amplify and read out the signals generated by one pixel. This paper discusses the limits of commercial multi pixel photon counters (MPPCs) from the producer Hamamatsu in terms of timing. We used only LSO:Ce codoped 0.4% Ca crystals [13] with a length of 3mm from the producer Agile with properties similar to ones used in PET systems, i.e. being non hygroscopic and with high stopping power. There were several attempts to model the time resolution of a photomultiplier based scintillation detector, e.g. [14,15]. To our knowledge however there is still a missing analysis on how to incorporate the crystal geometry into the statistical model, because even for short crystals, e.g. 2 × 2 × 3 mm 3 , one cannot omit this influence. This is due to the light transfer inefficiency and the photon travel time spread in the crystal in a specific configuration (Teflon wrapping, grease, surface state, . . . ). The light transfer efficiency is the ratio between the photons arriving at the photodetector and the intrinsic light yield, i.e. the total amount of scintillation photons produced. The photon travel time jitter is a combination of the light transfer time jitter of the scintillation photons and the gamma ray interaction point fluctuation in the crystal. This paper is organized in two main sections. First we present experimental results on coincidence time measurements using the so called NINO chip [16] and newly developed electronics. In the second part we present a complete Monte Carlo simulation tool modeling: (1) the scintillation statistics, (2) photon propagation and loss in the crystal (simulated by SLITRANI), (3) photodetector properties like single photon time resolution (SPTR) and detection efficiency, (4) signal pile-up considering the single cell signal of the SiPM, (5) dark count rate and optical crosstalk of the SiPM and (6) the effect of the electronics, i.e. bandwidth limitations and noise. A comparison of the simulation with experimental data shows good agreement for the influence of the bias voltage and NINO threshold on the CTR. We will use this simulation tool to discuss the nature of several important parameters to the time resolution, like scintillation rise time and fall time and light yield, as well as single photon time resolution (SPTR) and the photo detection efficiency of the SiPM.

Experimental methods and results
We use a time coincidence setup (see figure 1) to test the time resolution of the system. The radioactive β + decay of a 22 Na source produces a positron which emits two back-to-back 511keV γs -2 -  by annihilation with an electron. The 22 Na source has an activity of 3MBq and a spherical geometry with 1mm diameter. In the scintillating crystal these 511keV gammas are converted into visible photons with wavelength around 420nm. These photons are then detected by the photodetector (SiPM). The SiPM signals are fed into the CERN-developed NINO chip [16]. NINO is an ultrafast discriminator-amplifier employing the time over threshold technique. It produces a squared output pulse with the leading edge determining the time information (leading edge discrimination) and the width being a function of the pulse charge and thus delivering energy information (see figure 2). For the measurements in this paper however we will use the output of the linear voltage amplifier to get the energy information. NINO uses a common mode configuration with differential input signals coming from the SiPM [17] and hence can actively suppress common ground noise and pick-up. Despite the advantages of NINO being able to encode time and energy in one single pulse the method has an inherent drawback. Let us consider the SiPM pulse consisting of a single exponential with decay constant τ and amplitude A, The term V noise represents the electronic noise as well as noise generated by SiPM dark count events. It is modeled as a Gaussian distribution around the baseline. The NINO output pulse width (PW ) as a function of the threshold voltage V th can be calculated by setting equation (2.1) equal to V th : V th = A · exp(−PW /τ) +V noise . Considering that the amplitude of the SiPM signal is proportional to the gain G: A = a · G, this gives: The first term in equation (2.2) is related to the unit gain output a of the SiPM signal. It logarithmically encodes the energy resolution through the signal a and fluctuations of V noise . The gain of a SiPM can be increased by operating the device at higher overvoltages. This however also increases the DCR and thus leads to a higher noise component (V noise ). If a linear amplifier was used, the additional noise would be compensated by an increase in SiPM pulse height. On the other hand, equation (2.2) shows that in the logarithmic encoding used in NINO, an increase of the gain only leads to an additional constant factor in the pulse width (right term of (2.2)). Since the fluctuation of V noise also increases with higher overvoltage and is not offset by an increase in gain (left term of formula (2.2)), the energy resolution degrades. This behavior represents a disadvantage in the photopeak selection process when using the logarithmic energy measurement as compared to a linear one. This statment is supported by the findings in [18]. There we find an optimum in CTR around 1.7V overvoltage for the Hamamatsu 50µm MPPC. At higher bias overvoltages the CTR value degrades from 140ps (at 1.7V overvoltage) to 200ps (at 2.2V overvoltage). This degradation is almost entirely caused by a poor photopeak selection from the NINO pulse width energy measurement.

Electronics for SiPM readout
We developed a board with a parallel analog readout of the SiPM signal to overcome the problem of the NINO pulse width energy selection and to exploit the full capabilities of the SiPM. Figures 3  and 4 respectively show a picture of the circuit board and the simplified schematics of the electronics. The voltage amplifier has a high resistive input impedance to avoid additional noise at the NINO input. The board is designed to read out two SiPMs independently. The CTR can thus be measured with a single board that hosts two independent signal chains for a SiPM, NINO and a voltage amplifier.    Figure 5 shows the NINO and amplifier output signals of SiPM dark count events. Figure 5(a) shows the left arm of the coincidence setup. A cascade of afterpulse events can be seen at the analog amplifier output, together with the corresponding NINO signals. In figure 5(b) the right branch shows a typical dark count event at 0ns and a crosstalk event at 1000ns. The long fall time observed for the amplifier signals is generated by the 200Ω bias resistors. Due to the coupling capacitor between the SiPM and NINO, leading to a differential input impedance of 40Ω [16], this long fall time is not observed by NINO. It is only sensitive to signal changes in the nanosecond range. This behavior is intended as it corrects baseline shifts.
We measured the NINO output pulse slew rate to dV /dt = 179V /µs. Together with the rms noise floor of the oscilloscope of σ noise = 1.45 mV we can calculate the influence of our acquisition system on the time resolution according to formula (2.3). The factor 2 * 2 * ln(2) * √ 2 = 3.33 converts to coincidence time resolution FWHM. Table 1 summarizes the parameters for the SiPMs from Hamamatsu on this experimental board, including the measured breakdown voltage at a temperature of 20 • C.

NINO threshold calibration
The calibration of the NINO threshold is necessary to compare measurements with simulations. It relates the NINO threshold value to the SiPM output amplitude, expressed in multiples of single SPAD pulse height. These quantities can be derived from so called staircase plots, produced by measuring the number of triggers in counts per second versus the NINO threshold for different -5 -   SiPM overvoltages (see figure 6). The measurements were performed in a temperature stabilized dark box at a temperature of 20 • C. No crystal was attached to the SiPM.
Each step in the staircase plot represents the single SPAD amplitude incremented by an integer number [19]. NINO shows a non-linear behavior between its threshold and the equivalent cell amplitude N (see figure 7). Superimposed on the expected linear behavior one observes an exponential component in the threshold-to-amplitude response.
To account for this non-linearity we calibrate the NINO threshold in terms of the equivalent cell (SPAD) amplitude N with equations (2.4)-(2.7), by fitting exponentials to the measurements    signal can be described by a parabolic function, in a first order approximation and for low signal values. Thus the slew rate (dV /dt) is getting smaller for lower NINO thresholds, and consequently the electronic noise starts to dominate the CTR value. On the other hand the smearing of the scintillation signal by the photodetector's SPTR causes a significant change to the influence of the photon statistics. Triggering on the first photon emitted by the scintillation will not give the best CTR values anymore. This is mainly caused by the long tails of the SPTR well described by a Gaussian distribution [20]. These two effects are responsible for the rise of the CTR values for low NINO thresholds. We want to mention that these two effects are not independent, for a higher SPTR will as well influence the slew rate of the SiPM output signal.
The improvement of the CTR values with bias overvoltage (figure 10) can be mainly described by an enhancement of the photon detection efficiency (PDE) of the SiPM. The PDE includes the geometrical fill factor, quantum efficiency of the silicon and the probability of triggering an avalanche [19]. With increasing bias voltage the avalanche triggering probability becomes higher, allowing to detect more photons from the scintillation process and thus to improve the timing by higher statistics. The dark count rate and optical cross talk are increasing with bias overvoltage as well. At a certain bias overvoltage, however, this increase in DCR and optical cross talk will begin to deteriorate the time resolution, giving rise to the observed minimum.
3 Simulation of the system time resolution

Model description
This section gives an overview of the photon statistics to model the scintillation process. It provides a mathematical description of the signal pile-up as well as an overview (flow chart) of a dedicated Monte-Carlo simulation tool.

Scintillation mechanism and photon statistics
The photon emission rate of a scintillation process can be described by bi-exponential functions with one or several rise and fall times [21]. We use this phenomenological model to describe the photon generation process with the following probability density function (PDF): The parameter ∆t describes the start of the scintillation process, i.e. the absorption of a gamma in the crystal, and Θ denotes the Heaviside step function. Each contribution in the sum models different energy transfer mechanisms that populate and de-populate the luminescence centers by convoluting two exponential functions for each mechanism. Y l denotes the yield of every single energy transfer mode with ∑ l i Y l = 1. The cumulative distribution function (CDF) is then defined in equation (3.2) as: The area below f (t) is normalized to one. By multiplying f (t) with the intrinsic light yield n we obtain the average photon rate of the scintillation process f (t) = n · f (t). In the special case of L(Y)SO:Ce crystals we assume equation (3.1) and (3.2) to be simple bi-exponentials (l = 1) [22,23].
We are able to describe the k th photon's generation time PDF using the average photon rate and a modified binomial distribution [22,24,25], with f (t) and F(t) being the propability density function and cumulative density function of the photon generation process, respectively. Equation 3.3 describes the probability of the k th photon being emitted in the time interval t and t+dt. To better understand formula (3.3) we ask: "What is the probability of a photon y out of n to be the k th one and to be emitted in the time interval t and t+dt? " It is the probability that k-1 photons out of n-1 have already been emitted as described by the binomial distribution p = n−1 k−1 F(t) k−1 (1 − F(t)) n−k , multiplied by the probability of photon emission in the specific time interval t and t+dt, i.e. f (t)dt. Since one can repeat this for every photon of n, we apply the factor n in equation (3.3). This treatment is valid only if the rate of photon emission is time independent and identically distributed for each event. For n· F(t) reaching a constant value and in the limit of n → ∞, i.e. k being small compared to the total number of photons n, the binomial expression (3.3) can be reduced to the Poisson distribution: With F(t) = n· F(t) and f (t) = n· f (t). Equation (3.4) was already used by Post and Schiff [26] to describe the time response of a scintillator coupled to a PMT.

SLITRANI for light ray tracing
LITRANI [27] is a Monte Carlo program based on ROOT [28]. In this work, the newest version of LITRANI, called SLITRANI for "Super LIght TRansmission in ANIsotropic media", was used. Any three dimensional setup that can be described by the TGeo [28] class of ROOT can be used as a setup in SLITRANI. Each volume can be of a different material, and each material can have different optical properties. The optical characteristics of materials can be modeled with a dielectric constant, absorption length and diffusion length. Rayleigh and Compton scattering can also be simulated. Wavelength dependency can be added to most of the parameters. In contrast to other Monte Carlo light ray tracing programs (e.g. GEANT4), LITRANI can handle an anisotropic dielectric constant and absorption length. One has various possibilities to produce photons in SLITRANI: spontaneous emission of photons, photons coming from an optical fibre, photons generated by ionizing particles, photons generated by gamma rays of energy between 0.1 and 1MeV or photons generated by a high energy electromagnetic shower. One can define any volume and material inside the setup as the detector. All photons are tracked until they are absorbed or detected. Most of the transition, reflection and absorption processes can be recorded during the simulation for further analysis. The standard histograms are showing the amount of detected photons, the arrival time, the different wavelengths, materials and volumes where photon absorption or reflection appears. To handle unpolished surfaces and imperfect polishing, we added a self-built function where the roughness can be controlled by a single parameter.
The SLITRANI simulations can be summarized by two important parameters, i.e. the light transfer time jitter ξ and the light transfer efficiency η of the scintillaton light in the crystal. Both parameters are dependent on the gamma interaction point χ within the crystal. The light transfer time jitter describes the variation of the time between the generation of a scintillation photon and the moment it reaches the photodetector. We assume an isotropic emission angle of the photons. The transfer efficiency describes the fraction of photons reaching the photodetector compared to the total number of generated photons. It is dependent on the photon absorption and scattering in the crystal as well as on Teflon wrapping and also on the position of the gamma interaction point. We can define an average transfer efficiencyη, averaging the transfer efficiency η over all possible gamma interaction points.

Organisation of the Monte Carlo simulation tool
A multi pixel photon counter (MPPC) or SiPM consists of many micro cells (e.g. 3600 for the Hamamatsu 50µm type). If one of the micro cells or SPADs is fired, they always give rise to the same output signal, no matter how many photons initially triggered the avalanche. The single cell signal is defined by the SiPM's equivalent circuit and the circuitry around the SiPM [29]. In a first order approximation it arises from RC-filters and thus should only have exponential character. The time constants are dependent on the overvoltage, i.e. on the cell capacitance changing with the overvoltage as the depletion zone will change in dimension. However these changes are assumed to be negligible, and to good approximation, we can further assume that the time constants determining the single cell signal are independent of the overvoltage. We model the single cell signal with a bi-exponential function (3.5) with τ Mr and τ Md being the rise time and fall time of the signal, respectively. The parameter A denotes the maximum amplitude of the single SPAD signal, and ε -10 -

JINST 8 P07014
the moment when each microcell fires. Typically the rise time τ Mr is in the order of hundreds of ps and the fall time τ Md in the order of tens of ns. A single cell amplitude jitter caused by gain fluctuations can be modelled by a Gaussian distribution , withĀ being the mean value and σ A the standard deviation of the fluctuation. The value σ A /Ā gives the percental gain fluctuation, and the inverse,Ā/σ A , will be called SNR gain fluctuation .
Coupling the SiPM to a scintillator with an intrinsic light yield n and a light transfer efficiency η it will give rise to n · PDE · η avalanches. Because not all photons generated by the scintillation n will actually produce an avalanche due to absorption in the crystal and limited photo detection efficiency of the SiPM. The probability of one photon emitted by the scintillation to produce an avalanche is n · PDE · η. This process of losing photons from the scintillation to the production of an avalanche is denoted as random deletion. The output signal of the device can be described as the sum of the single cell signals (3.6). Because of the long decay time of a single cell signal this process is often referred to as signal pile-up.
The function "random(1)" generates a random number between 0 and 1 for each k. Thus the Heaviside function Θ [PDE · η(χ) − random(1)] represents the random deletion due to the limited detection efficiency of the SiPM (PDE) and the transfer efficiency of the crystal (η). The transfer efficiency η(χ) depends on the scintillation origin in the crystal (χ) and thus changes with each gamma interaction.
In our simulation ε(k) is a crucial parameter as it models the entire timing behaviour of the system. We identify four contributions to ε(k), i.e. the time ∆t from the emission of the 511keV gamma until its absorption in the crystal, the scintillation statistics (see section 3.1.1), the light transfer time jitter ξ and the detection transit time spread of the SiPM. Supposing that the T k (p) operator generates a single random time stamp with the PDF p, then ε(k) can be expressed as (3.7).
The function T k ( f ) gives one time stamp per k th photon of the scintillation process, according to the method described in section 3.1.1. The function ξ k (χ) gives the random time spread caused by the light transfer in the crystal for the k th photon simulated by SLITRANI (see section 3.1.2). It should again be noted that ξ k (χ) is dependent on the position of the scintillation origin χ. T k (g) models the photodetector's transit time spread (see formula (3.8)) and gives a random time stamp generated by a Gaussian distribution. The sigma of this Gaussian is the single photon time resolution (SPTR) of the SiPM.
The parameter ∆ M denotes a possible electronic delay and σ SPTR is the single photon time resolution of the SiPM. The Monte Carlo simulation is organized in the following way: • Generate a 511keV gamma and calculate the time ∆t until it is absorbed in the crystal at the coordinates χ according to the exponential gamma ray absorption law with an absorption length of 11.5mm [30]. The point χ will be the scintillation point of origin.
• At point χ calculate from SLITRANI simulations n random photon transfer time jitters ξ 1 . . . ξ n . Add to every photon transfer time jitter ξ k a random time stamp t k , with k being the photon rank k ∈ n. t k is subject to the scintillation statistics described in section 3.1.1.
• Perform random deletion of photons according to the product of PDE and the photon transfer efficiency of the scintillator (η) simulated by SLITRANI. The transfer efficiency depends on the scintillation origin χ.
• Generate a random noise floor taking into account the DCR and the optical crosstalk of the SiPM.
• With (3.6) generate signal pile-up for those photons that have "survived" the random deletion process. This also includes signals from crosstalk which are assumed to be simply smeared by the SPTR because at this time no model of time evolution of crosstalk events is known to the authors.
• Filter the obtained signal with a first order Butterworth low pass filter with a cut off frequency of 1GHz, matching the bandwidth of NINO.
-12 - We wish to note that we did not include any saturation effects of the SiPM in the simulations. Within the first nanosecond approximately 100 avalanches are created at the photodetector and therefore we assume that the used SiPM with a total number of 3600 SPADs will show no signs of saturation. Therefore our method should be a good representation of the SiPM signal in the first few nanoseconds of the scintillation pulse. We then apply leading edge discrimination producing the time stamp of the simulated output pulse. Electronic noise is added by Gaussian distributed white noise on the threshold value. To obtain sufficient statistics we repeat this procedure about 10000 times. A flow chart of the program and processing steps can be seen in figure 11.
As a last step we quadratically add CT R acquisition = 27ps FWHM to the simulated CTR value, to account for the time jitter generated by the oscilloscope (see formula (2.3)). In the simulation we did not include any fluctuations of the intrinsic light yield, i.e. the intrinsic energy resolution of the crystal.

Model parameters
In this section we describe the different measurements and methods to obtain values for all input parameters necessary for the simulation.

Single cell signal
We model the single cell signal by a bi-exponential function with a single rise time τ Mr and fall time τ Md (equation (3.5)). Figure 12 shows the measurement of a single cell signal coming from the voltage amplifier. The NINO discriminator, although not used here, was nonetheless coupled to the SiPM in order to preserve the same circuit conditions as for real timing measurements.
The measured signal is the direct output voltage of the SiPM. We can observe two fall time components. A slow one with 260ns and a fast component with 11ns. The slow 260ns component -13 -is generated by the electronics around the SiPM, i.e. the 200Ω resistors (as can be seen in figure 4). The fast fall time component is generated by the 1nF coupling capacitor at the front end of NINO and its additional input impedance of 40Ω. It corresponds to the fall time τ Md in equation (3.5) which is 11ns. It should be noted that a very small fraction of the slow 260ns-component can still couple into NINO. This fraction however is small enough and will be neglected in our further discussions. In addition we measured the single cell signal with 50Ω bias resistors to determine the SiPM capacitance. The decay time was measured to be 50ns hence corresponding to a SiPM capacitance of 500pF. Compared with the values measured in [29] we suspect a high passive capacitance in our system generated by rather long leads and jumpers used to guide the SiPM signal to NINO. The rise time can be determined to τ Mr = (C D +C q )R D [31]. It is difficult to find the proper values for the "micro plasma" resistance R D , the SPAD capacitance C D , and the quenching capacitance C q . Together with the values given in [32] we deduce τ Mr = 400ps using a value of R D = 3kΩ and C D +C q = 130pF. This value is most probably too high. However, because of the low pass filtering in our simulations, the single cell rise rime t Mr is not a critical parameter. It should be noted that also the exact value of the fall time t Md is not very critical since we are only interested in the pile-up process in the first nanoseconds.

Electronic noise and gain fluctuation
NINO operates as a current mode amplifier and has a low differential input impedance of 40Ω [16]. A schematic of the SiPM bias circuit can be seen in figure 4. It is important to determine the exact electronic noise of the system, in particular when measuring the SPTR. Gain fluctuations of the SiPM have a similar influence on the time resolution as electronic noise. These gain fluctuations can be directly modeled by an amplitude fluctuation of the single cell signal, supposing that the rise and fall times of the single cell signal are independent of SiPM overvoltage and identical within every cell (section 3.1.3). We developed a technique to estimate the electronic noise and the gain fluctuation by means of measurements. With the oscilloscope we measure, at the same time, both the analog SiPM signal from the amplifier and the output pulse of NINO, being operated at a specific threshold. For each dark count event in the SiPM we record the SiPM pulse amplitude and if there was a parallel output signal from NINO or not. In the offline analysis, we select single dark count events using the recorded SiPM amplitude histogram to exclude crosstalk and afterpulsing. We plot the number of events where NINO does not generate a signal versus the threshold (see figure 13). For low thresholds, NINO triggers on every signal, whereas for high thresholds the input signal from a single cell is too low to give rise to a trigger, and thus no signal at the NINO output is seen. For a perfect input signal with no noise and gain fluctuations this transition should be a sharp step function. In a noisy system, however, this transition is smooth and blurred by electronic noise and amplitude fluctuations. We are able to deduce the combined influence of electronic noise and gain fluctuations by calculating the derivative of this transition and applying a Gaussian fit. The ratio of sigma and mean of the Gaussian fit denotes the measured signal to noise ratio SNR meas = µ σ . It should be noted that throughout this work the signal in the SNR term is always equal to the cell amplitude. This equivalent cell amplitude is changing with bias overvoltage and thus SNR meas as well. We can model the measured SNR as shown in equation (3.9).  ((14*U ov ) -2 +26 -2 ) -0.5 measured SNR Figure 13. Number of no NINO-output-signalseen versus the NINO threshold. On the derivative we apply a Gaussian fit which gives the measured SNR. Figure 14.
Measured SNR versus bias overvoltage.
The fit is parametrized as SNR electronic noise = 14 * U ov and SNR gain fluctuation = 26.
Because the electronic noise is independent of overvoltage in contrast to the single cell amplitude changing linearly with gain and hence with bias overvoltage, SNR electronic noise is directly proportional to the overvoltage U ov . We further assume that the gain fluctuations are mainly generated by non-uniformities of the 3600 SPADs of the SiPM. Thus, the relative fluctuation of the gain (SNR gain fluctuation ) should be independent of the bias overvoltage. We model the two contributions in equation (3.9) as SNR electronic noise = a ·U ov and SNR gain fluctuation = b with a and b being two constant parameters. In figure 14 we show the measured SNR versus the bias overvoltage. Applying a fit according to formula (3.9), we can deduce the parameters as: a = 14/V and b = 26. Thus the electronic signal to noise ratio is SNR electronic noise = 14 · U ov , and the gain (single cell amplitude) fluctuation is 3.8%.

Dark count and crosstalk
DCR and optical crosstalk are dependent on the bias overvoltage and can be derived from the staircase plots in section 2.2. The first plateau is at thresholds lower than a single SPAD signal and hence corresponds to the dark count rate. The second plateau is generated from a single crosstalk, the third from two simultaneous crosstalk events and so on. In figure 15 we plot measured dark count and crosstalk values for four different bias overvoltages. We can fit the DCR and crosstalk rates with exponential functions depending on the overvoltage. The functions are summarized in equations (3.10) to (3.13). The parameters deduced from these fits are not directly related to physical processes. They rather provide a phenomenological tool to describe DCR and crosstalk in an easy and integral form for our simulations. As such they are look-up tables in a mathematical description.

Photon detection efficiency
For the photon detection efficiency we used the results published in [33]. In this article measurements of the PDE versus the bias overvoltage at different wavelengths, ranging from 465nm to 870nm, are reported. We are only interested in the measurements for 465nm because of the emission peak at 420nm for LSO:Ce codoped 0.4% Ca.
We can parametrize the PDE as a function of overvoltage for 465nm with equation (3.14) (see figure 16).
PDE(U OV ) = 0.591 * 1 − e −0.379 * U OV (3.14) Figure 16. Photon detection efficiency versus bias overvoltage for the Hamamatsu 50µm type, experimental data and fit for 465nm. Data was taken from [33]. -16 -These parameters only provide a phenomenological description of the data without a direct physical explanation. However, as we can see in figure 16, this parametrization of the PDE is a good representation of the experimental data. If we couple the photodetector to a scintillator, we also have to account for the emission spectra of the scintillator. In figure 17 we plot the emission spectrum of the used LSO:Ce codoped 0.4%Ca scintillator together with the PDE of the photodetector versus the wavelength. If we normalize the PDE at a wavelength of 465nm we can calculate the weighted integral of the PDE with the emission spectrum which results in a correction factor of 0.96 applied to equation (3.14):

Single photon time resolution
The single photon time resolution (SPTR) describes the time jitter of the photodetector when detecting single photons and is commonly quoted in sigma, i.e. not in FWHM. The SPTR is influenced by the characteristics of one single pixel, i.e. the avalanche spread in the SPAD, the uniformity of the single cells and the uniformity of the signal transmission to the output. As the signals from the SiPM are small, the electronics plays an important role for the SPTR value. It is therefore important to measure the SPTR with the same electronics used in the CTR measurements. We tested the time response of the SiPM using a femtosecond laser operating at 400nm wavelength and with a pulse width of 200fs [20]. The measured SPTR values can be seen in figure 18.
To get an idea of the influence of the electronic noise and gain fluctuations, we used the values derived in section 3.2.2 and simulated the time jitter of the noise contributions taking into account the NINO threshold and its shift in terms of equivalent cell amplitude when changing the bias voltage. Subtracting the noise contribution from the measured SPTR values we obtain an almost constant value of σ SPTR = 66ps. Although this result is still under investigation, it seems that the SPTR is not strongly dependent on the bias overvoltage. This could be explained by the fact that the SPTR in our MPPC with 3600 SPADs is dominated by uniformity jitter and signal transmission delays. Thus, a device with smaller geometrical dimensions could lead to a better SPTR.

Scintillation rise and fall time
For the scintillation decay time measurement we used the time correlated single photon counting technique originated by Bollinger and Thomas [34]. This way we determined the single fall time of our LSO:Ce codoped 0.4%Ca scintillator to 30.3 ± 1ns (see figure 19), in good agreement with the values stated in [35]. To measure the scintillation rise time is more difficult as one needs to deconvolute the impulse response of the measurement apparatus from the intrinsic rise time of the scintillation signal [36]. The rise time values stated in the literature for LSO:Ce range from 0ps to 70ps [23,36]. In this work we choose a rise time of 70ps, comparable to the value stated in the most recent work [23]. Using equation (3.1), we derive the parameters for the PDF of the photon generation process as a single bi-exponential (l = 1) with a rise time τ r = 70ps and a fall time

Light yield measurements for the adaptation of SLITRANI crystal model parameters
The different parts of the SLITRANI setup (crystal, photodetector, wrapping, etc.) are modelled as three dimensional objects that are provided by the ROOT geometry package. In figure 20 a schematic of the main parts of the light yield measurement setup can be seen, as well as the different regions of a scintillating crystal which were modeled with different surface roughnesses. The photons in our setup are produced with a Cs137 gamma source. Other parameters necessary for photon generation by gamma excitation, like the photo-electric cross section, the gammaabsorption length, and the emission spectrum, were taken from the work of [39].
There are a number of parameters that are well known like the index of refraction of the crystal or of the glue, but there are nonetheless others that change from crystal to crystal and are difficult to measure (e.g. the roughness of surfaces and edges, the diffusion length inside the crystal). Therefore these parameters were fitted to our model using the model-matching technique [37]. In the work of [38] the parameters, like the crystal diffusion and absorption length as well as the surface roughness, were fitted making use of light yield measurements with crystals of differently shaped -18 -  Table 2. Parameters used in our SLITRANI simulation setup. Parameters were taken from the work of [38] and [39]. walls and different contact agents. The parameters for the crystal used in this work were partly imported from the work of [38], but since wrapping was not considered there, some parameters had to be estimated by running another set of simulations as part of a second model-matching technique for this work. Similar to the work of [37,38] the unknown parameters of the setup were adjusted to fit the results of independent measurements with crystals of different size. The unknown parameters for our model were: (1) airgap between wrapping and crystal, (2) roughness of walls and edges (3) absorption and diffusion of the Teflon wrapping and (4) index of refraction of the wrapping. To fit those parameters we made light yield measurements with a crystal with a common 2 × 2 mm 2 face that is in contact with the photodetector but different lengths of 3, 5, 10 and 20mm. The crystals were wrapped in Teflon, and optical glue was used for the coupling to the photodetector. The fitted parameters together with parameters taken from other references can be seen in table 2 and  table 3.  In table 4 we summarize the measured light yield and the simulated light transfer efficiencies for different crystal configurations. The light yield was measured with a Photonis XP2020Q PMT. The values are corrected for the quantum efficiency of the PMT and given in photons per MeV. We measured the light yield in different configurations, i.e. with 5 (5FP) or 6 polished faces (6FP) and different crystal lengths. The configuration "5 faces polished" means that only the side opposite to the SiPM is unpolished. The crystal was always wrapped in Teflon, 0.075mm thick and with at least 5 layers, and coupled to the photodetector with optical grease: Rhodorsil 47V.
The measured light yield uncertainty was estimated to be about 5%, including uncertainties in quantum efficiency and collection efficiency of the PMT. This number also includes handling uncertainties caused by the mounting of the samples onto the PMT with possible misalignment and -19 -

Model validation
In this section we compare the simulated CTR values with the measured ones. In table 5 we summarize the input parameters describing the scintillation process and in table 6 those for the SiPM. We emphasize that the simulation parameters were not deduced from any fit of the measured CTR. Thus the simulation outcome is pure prediction without incorporating any prior knowledge of the CTR measurements.

Special case: timing of the SiPM using a femtosecond laser
An important special case is the simulation of the time response of the SiPM plus the accompanying readout electronics alone without the scintillation mechanism and light transfer. Measurements were done at different light levels with a femtosecond laser operating at 400nm wavelength and 200fs pulse width [20]. Figure 21 shows the measured time response of a Hamamatsu 50µm SPAD size MPPC versus the number of fired cells or photo-electrons. The SiPM was operated at an overvoltage of 2.2V. We compare the experimental data with simulations described in section 3.1.3. For that purpose we "switched off" the scintillation and crystal influence. The simulation thus represents the interplay of the signal pile-up with the SPTR jitter and the leading edge discrimination. The solid line shows the simulation with its associated 99.7% confidence interval. These error bars comprise only the statistical Monte Carlo fluctuation since we have not yet incorporated the uncertainties of the input parameters.

CTR as a function of SiPM Bias and NINO threshold
We compare the Monte Carlo simulation results with measured data. In the following plots (figure 22 and 23) the solid lines show the simulation with a 99.7% statistical confidence level (no parameter uncertainties are incorporated) different shaped points are the data points from different CTR measurements, e.g. CTR versus bias overvoltage and CTR versus NINO threshold voltage. Figure 22 shows the CTR versus the SiPM bias overvoltage. Points show the experimental data at a NINO threshold of 80mV. We achieved best CTR values for all SiPM overvoltages around this threshold value. The simulation is in good agreement with the experimental data. The decrease of the plot (which corresponds to an improvement in CTR) as a function of SiPM overvoltage is almost perfectly matched by the corresponding increase in PDE, as we modeled the SPTR to be constant with bias voltage. For higher bias overvoltages the experiment shows a higher saturation or even a degradation in CTR when compared to simulations. One possible explanation can be the neglect of an additional time smearing of the crosstalk events, which we had assumed to be prompt.  We also neglected afterpulsing and its time distribution. However we do not believe that afterpulses influence the time resolution to a high degree as these events happen relatively late compared to the avalanches triggered by the relatively high scintillation photon rate (∼10 avalanches per 100ps).
In figure 23 the CTR versus the NINO threshold for different bias voltages can be seen. For low bias overvoltage (1.1V) the simulation predicts around 5ps lower values of CTR than measured. This could be an indication that the approximation of the SPTR being independent of the bias overvoltage is not absolutely true. By looking more closely at picture 18 in section 3.2.5 we can even see a slight dependence of the SPTR on the bias voltage. Furthermore for 1.1V overvoltage we have no measured SPTR value available, and thus the simulation for this bias overvoltage region is simple extrapolation. Figure 23 shows an interesting behavior at low threshold values. As we had -22 -

JINST 8 P07014
noted in section 2.3 the rise in the CTR value at low threshold values is caused by the SPTR and electronic noise. The fact that, in this low threshold region, the MC simulation predicts consistently lower CTR values than the measured data could be explained by a deficiency in the model of the single cell signal. An additional parabolic component at the beginning of the single cell signal could lead to a lower slew rate of the signal at low threshold values. Because of electronic noise this would lead to a deterioration of the measured CTR values as compared to the simulation, where we assumed a simple bi-exponentially shaped single cell signal.
The deviation of the measurement from the simulation is in the range of the measurement error plus the statistical fluctuations of the Monte Carlo simulation. In view of the fact that no fits were made to the experimental CTR data, the agreement between simulation and measurement is quite good. In a future work we will address the question of the SPTR in more detail and incorporate a full model of the crosstalk and afterpulse events. We will also calculate full covariance matrices of the simulations taking all parameter uncertainties into account.

Discussion
In the previous section we showed that the developed Monte Carlo simulation tool is in good agreement with experimental data. This gives us confidence to correctly quantify the different model input parameters to the CTR. We change one input parameter of the simulation, e.g. light yield, scintillation rise-, fall time, PDE or SPTR while keeping the others constant. In figure 24 we depict the minimum simulated CTR values versus the normalized parameter value, i.e. the model input parameter under investigation multiplied by a factor ranging from 0.1 to 10.
The CTR is dependent on the square root of the intrinsic light yield (n), the scintillation fall time (τ d ) and the SiPM photon detection efficiency (PDE), i.e. CT R ∝ τ d /(n · PDE). This is in agreement with the findings of [14] and [40]. In the conditions given in table 5 and 6, a lower scintillation rise time (τ r ) than 70ps seems to improve the CTR only marginally. This outcome is also in agreement with [14]. For the minimum CTR values we calculate the number of photoelectrons Because of the electronic noise and the SPTR triggering on the first photon this situation will not necessarilly lead to the best CTR values. This also explains why the scintillation rise time seems to have only little influence on the CTR. If we could reach lower threshold values we might then expect a behavior like CT R ∝ √ τ r as predicted by photon statistics. For the SPTR we observe a similar behavior as for the scintillation rise time: a decrease of the measured SPTR value of 66ps by a factor of 10 only leads to a CTR improvement by ∼ 10ps. The single cell signal rise time has almost no influence on the CTR (see figure 24(b)). This is because of the bandwith limitation of NINO together with the low electronic noise and gain fluctuation in our system.
In figure 25 we show the minimum simulated CTR versus bias overvoltage and the simulated CTR versus threshold expressed in terms of equivalent single cell amplitude at an overvoltage of 2.3V. In these graphs we successively eliminate the contributions of electronic noise and gain fluctuations, the dark count and crosstalk probability, the SPTR and finally the light transfer time jitter and gamma conversion jitter in the crystal. The light transfer time jitter is denoted as ξ k (χ) and the gamma conversion jitter as ∆t in equation (3.7). We always keep the PDE and the photon transfer efficiency of the crystal constant. We can see that electronic noise and gain fluctuations have almost no influence on the final time resolution. We notice an improvement only for very low bias and equivalent cell amplitude values if we set these two terms to zero. If additionally we also set the dark count and crosstalk probabilties to zero, we notice an improvement by ∼5-10ps at 2.3V overvoltage. This improvement is mostly due to the absence of crosstalk. Setting additionally the SPTR to zero leads to an further improvement of around 20ps. In figure 25(b) we see that the CTR changes significantly as a function of the equivalent cell amplitude. Because of the zero SPTR, a minimum in the CTR plot at very low threshold values is observed, meaning that Poisson-photostatistics gain in influence. If we only consider Poisson-photostatistics triggering on the first photon would then lead to the best CTR, which in our system occurs only at very low thresholds. The last test is to set the influence of the light transfer time and gamma conversion jitter -24 -  in the crystal to zero. An additional improvement of approximately 25ps can be seen in figure 25, manifesting the importance of the crystal dimensions on time resolution, and this even for small crystals with dimensions of only 2 × 2 × 3 mm 3 . We want to mention that for these tests we kept the PDE and photon transfer efficiency constant.
In figure 26 we show the improvement of the CTR if we only set the light transfer time plus the gamma conversion jitter to zero and leave all other parameters unchanged. Even in this case we see an improvement of around 20ps. The improvement is smaller because other factors like SPTR, DCR and crosstalk are deteriorating the CTR. We can estimate the influence of the light transfer time plus the gamma conversion jitter to the CTR to (110ps 2 − 90ps 2 ) = 63ps. According to [32] we can determine the gamma conversion jitter for a 2 × 2 × 3 mm 3 size crystal to approximately 6ps FWHM CTR and thus negligible if compared to 63ps.
To understand the rather dramatic influence of the crystal size on the CTR we show the light transfer time jitter in figure 27. For this histogram we simulated random photon emissions in the crystal according to the SLITRANI model also used for the CTR simulations (section 3.2.7). The time from the emission of the photon until it reaches the face of the crystal next to the photodetector is recorded and plotted in the histogram. We can see two initial peaks, the left one correspond to direct photons and the second one to reflected photons. The direct photons are the ones emitted directly towards the photodetector, whereas the reflected ones are emitted in the opposite direction and thus have to undergo at least one reflection at the back end of the crystal in order to reach the front face. In addition we observe a rather long exponential-like tail. This tail is caused by photons that are not in a mode to exit the crystal directly. These photons undergo several reflections in the crystal and eventually exit due to surface roughness and scattering in the crystal. It is still under investigation how this tail influences the CTR, as the total number of photons in this tail is not negligible. From figure 27 we determine the influence of the intrinsic light transfer time -25 -jitter (ξ k (χ)) to 40ps. This value is smaller than the value of 63ps as calculated from the CTR simulations and gives a first hint on how the tail in figure 27 might influence the time resolution.

Conclusion
We developed new electronics to probe the timing capabilities and limitations of a possible TOF-PET system using LSO crystals. An analog voltage amplifier gives the energy information and the discriminator amplifier NINO the leading edge time information. With LSO:Ce codoped 0.4% Ca and crystal dimensions of 2 × 2 × 3 mm 3 coupled to a commercial MPPC from Hamamatsu (S10931-050P) we measured a CTR of 108±5ps FWHM for gammas with an energy of 511keV. To understand the measured CTR values in detail we developed a comprehensive Monte Carlo simulation tool to model the scintillation process and the detector chain, i.e. the light transfer in the crystal, the SiPM properties and the electronics. We took special care of deducing the input parameters for the Monte Carlo simulation without fitting to measured CTR values. Our simulations show very good agreement with the experimental CTRs as a function of both bias overvoltage and NINO threshold. In addition we were able to give a good description of the intrinsic time response of the SiPM when exposed solely to femtosecond laser pulses with different light intensities.
In our system the CTR is strongly dependent on the scintillator's light yield, on the scintillation fall time and the PDE of the SiPM. Due to the interplay of the single cell signal pile-up with the SPTR and the electronic noise we were unable to trigger on the first photon emitted in the scintillation process. In addition, because of the SPTR, triggering on the first photon will not automatically lead to the best CTR values. Consequently the scintillation rise time gives a smaller contribution to the CTR. In addition the assumed rise time of 70ps is low compared to the achieved CTR of 108ps. We used commercial SiPMs from Hamamatsu with a maximum PDE of 33%, a rather modest value. Increasing the PDE of the SiPMs would be a plausible and effective way to improve the time resolution. Also an increase in the scintillation light yield and a faster decay time would both lead to similar improvements.
With the developed model we determined the photon travel jitter in the crystal influencing the time resolution to a high degree. If we artificially set the photon travel jitter for the used 2 × 2 × 3 mm 3 geometry to zero we would obtain a CTR improvement of 20ps. For PET systems, crystals of 15mm length and longer are needed to achieve an acceptable conversion efficiency of the 511 keV gamma rays. Thus the influence of the photon travel jitter will be more significant and consequently a potential problem, if a CTR of 100ps FWHM is to be achieved. Techniques like double-sided readout with travel jitter corrections could be a possible solution in that case.