Fundamental limits of scintillation detector timing precision

In this paper we review the primary factors that affect the timing precision of a scintillation detector. Monte Carlo calculations were performed to explore the dependence of the timing precision on the number of photoelectrons, the scintillator decay and rise times, the depth of interaction uncertainty, the time dispersion of the optical photons (modeled as an exponential decay), the photodetector rise time and transit time jitter, the leading-edge trigger level, and electronic noise. The Monte Carlo code was used to estimate the practical limits on the timing precision for an energy deposition of 511 keV in 3 mm × 3 mm × 30 mm Lu2SiO5:Ce and LaBr3:Ce crystals. The calculated timing precisions are consistent with the best experimental literature values. We then calculated the timing precision for 820 cases that sampled scintillator rise times from 0 to 1.0 ns, photon dispersion times from 0 to 0.2 ns, photodetector time jitters from 0 to 0.5 ns fwhm, and A from 10 to 10 000 photoelectrons per ns decay time. Since the timing precision R was found to depend on A−1/2 more than any other factor, we tabulated the parameter B, where R = BA−1/2. An empirical analytical formula was found that fit the tabulated values of B with an rms deviation of 2.2% of the value of B. The theoretical lower bound of the timing precision was calculated for the example of 0.5 ns rise time, 0.1 ns photon dispersion, and 0.2 ns fwhm photodetector time jitter. The lower bound was at most 15% lower than leading-edge timing discrimination for A from 10 to 10 000 photoelectrons ns−1. A timing precision of 8 ps fwhm should be possible for an energy deposition of 511 keV using currently available photodetectors if a theoretically possible scintillator were developed that could produce 10 000 photoelectrons ns−1.

In this paper we review the primary factors that affect the timing precision of a scintillation detector. Monte Carlo calculations were performed to explore the dependence of the timing precision on the number of photoelectrons, the scintillator decay and rise times, the depth of interaction uncertainty, the time dispersion of the optical photons (modeled as an exponential decay), the photodetector rise time and transit time jitter, the leading-edge trigger level, and electronic noise. The Monte Carlo code was used to estimate the practical limits on the timing precision for an energy deposition of 511 keV in 3 mm × 3 mm × 30 mm Lu 2 SiO 5 :Ce and LaBr 3 :Ce crystals. The calculated timing precisions are consistent with the best experimental literature values. We then calculated the timing precision for 820 cases that sampled scintillator rise times from 0 to 1.0 ns, photon dispersion times from 0 to 0.2 ns, photodetector time jitters from 0 to 0.5 ns fwhm, and A from 10 to 10 000 photoelectrons per ns decay time. Since the timing precision R was found to depend on A −1/2 more than any other factor, we tabulated the parameter B, where R = BA −1/2 . An empirical analytical formula was found that fit the tabulated values of B with an rms deviation of 2.2% of the value of B. The theoretical lower bound of the timing precision was calculated for the example of 0.5 ns rise time, 0.1 ns photon dispersion, and 0.2 ns fwhm photodetector time jitter. The lower bound was at most 15% lower than leading-edge timing discrimination for A from 10 to 10 000 photoelectrons ns −1 . A timing precision of 8 ps fwhm should be possible for an energy deposition of 511 keV using currently available photodetectors if Content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

Introduction
The fundamental limits of timing precision in scintillation detectors have been studied for many years, and the physical effects that influence timing have been long understood. But because there are at least ten different contributing factors, most analyses have chosen to concentrate on a few factors and assume that the rest are held constant at 'typical' values (Choong 2009, Vinke et al 2009, Seifert et al 2012b, Lecoq 2012, Gundacker et al 2013, Auffray et al 2013, Lecoq et al 2013. Other treatments use advanced statistical analysis to estimate the lower bound on the timing precision, where the timing from some or all the photoelectrons is used (Seifert et al 2012a). In contrast this paper presents a brute-force approach that uses hundreds of Monte Carlo calculations of the timing precision to sample the most useful region of the parameter space and then presents a closed-form expression that allows the timing precision to be estimated at any point in the space. The advantages of this approach over previous work are that it (1) models the rise time and optical photon time dispersion separately, (2) makes numerical values for comparing different designs more accessible, and (3) can be applied to currently available scintillators as well as 'next-generation' scintillators that approach fundamental limits in terms of luminosity and speed. It is intended for preliminary design studies and is not a substitute for more accurate calculations that capture the details of the system being designed such as knowledge of the depth of interaction, a more accurate treatment of the time dispersion of the scintillation photons, and the properties of the photodetector used.
The paper is organized as follows. In section 2, we describe the temporal sequence of events from the arrival of a nuclear particle to the generation of a timing pulse and identify six parameters that have the most effect on the timing precision, specifically (1) the depth of interaction, (2) the scintillator rise time and (3) decay time, (4) the optical photon time dispersion in the scintillator due to variations in path length, (5) the number of photoelectrons detected, and (6) the photodetector time jitter. In section 3 we describe the Monte Carlo calculation that models these factors as well as pulse height variability and electronic noise. In section 4 we explore the time jitter of the individual photoelectron pulses. In section 5 we tabulate the calculated timing precision for Cherenkov radiation produced over time periods from 0 to 0.2 ns and detected by a photodetector with Gaussian time jitter. In section 6 we compute the timing precision for six commonly used scintillators and three photodetectors (18 combinations). In section 7 we show the trigger level optimization for several cases. In section 8 we show that the rise time of the photodetector single photoelectron output pulse is not an important factor in the timing precision. In section 9 we tabulate the results of Monte Carlo calculations for 656 different cases in the parameter space of scintillator rise time 1.0 ns, photon dispersion 0.2 ns, photodetector time jitter 0.5 ns fwhm, and 10 to 10 000 photoelectrons ns −1 decay time. In section 10 we fit a closed-form expression to these using the Eureqa system (Nutonian, Inc.). In section 11 we describe some of the issues involved in combining the information from two photodetectors, one at each end of the scintillator. In section 12 we show the effect of electronic noise on the timing precision. In section 13 we show an example comparing the lower bound on the timing precision with leading-edge timing. Section 14 discusses the concepts developed in the paper and section 15 summarizes the conclusions. Note: since the zero of the time scale depends on signal paths and only the variance of the time measurements is of concern in this work, the term precision is used rather than accuracy or resolution.

Physical effects in the detection processes
Many sequential processes are involved in the conversion of a nuclear particle to electronic signals that can be used to estimate the time when the particle enters the scintillator. These are discussed in temporal order below.

Depth of interaction
Gamma rays, neutrons, and scintillation photons travel at different speeds in the scintillator. As a result, the time between the entrance of the gamma ray or neutron and the detection of the scintillation light depends on the depth of interaction. For example, in a scintillator with an index of refraction n and gamma rays entering the surface opposite a photodetector, moving the depth of interaction by z toward the photodetector decreases the arrival time of the scintillation photons by about t = z (n − 1) 0.033 ns cm −1 . By reading out both ends of a 20 mm long Lu2SiO5:Ce (LSO) crystal, a depth of interaction can be measured with an uncertainty of 3-4 mm fwhm (Yang et al 2006). This corresponds to a timing uncertainty of about 0.010 ns fwhm.
In most cases thick scintillators are needed for good detection efficiency and singleended readout severely limits the ability to measure the depth of interaction, resulting in a corresponding uncertainty in the timing. Thus, double-ended readout is assumed throughout this paper. Amplitude information and time stamps are used to estimate both the depth of interaction and arrival time (figure 1). In section 11 we describe the issues involved. Depth of interaction information is also important for reducing the effects of parallax error when imaging radiation from oblique angles, such as in positron emission tomography.
In positron tomography the two annihilation photons are emitted simultaneously and the difference between the arrival times at the coincident detectors depends only on the position of the annihilation point along the line between them.

The ionization process
When charged particles (electrons, protons, alphas, other ions) enter the scintillator, they create an ionization track as they lose their energy. Gamma rays and x-rays produce ionization by interacting with electrons in the scintillator. Neutrons produce ionization by interacting with atomic nuclei in a variety of elastic and inelastic processes. These processes can be very fast. For example, a 1 MeV electron in a NaI crystal has a path length of about 2 mm and for energies above 80 keV the velocity is greater than 50% of the speed of light. A 5 MeV recoil proton in plastic has a path length of about 0.4 mm. In both cases the charged particles stop in a time scale of 0.010 ns. The result of this process is the creation of energetic electrons and missing electrons (holes) in the scintillator that thermalize to the bottom of the conduction band and the top of the valence band, respectively over a time scale less than 0.001 ns (Rodnyi 1997). The number of electron-hole pairs initially created is given by E/(βE g ), where E is the energy deposited and βE g is the ionization energy required to make an electron-hole pair, expressed as a multiple of the band gap E g (Lempicki et al 1993).

Transport of ionization energy to radiative and non-radiative centers
There are several mechanisms that reduce the transport efficiency S of the electron-hole excitation energy to the radiative centers. One is when intermediate excitonic states are lost through the Auger process. This causes quenching at high ionization density and is responsible for the lower light output of ions relative to electrons, and of electrons at the end of their tracks. The latter is one of the causes of gamma-ray non-proportionality. Another mechanism that reduces transport efficiency is electron and hole trapping on non-radiative centers. These trapping centers can be on crystal structure defects or impurity atoms.
In some cases (e.g., Bi 4 Ge 3 O 12 , BaF 2 , and PbWO 4 ) radiative centers are intrinsic to the host crystal. In other cases activator ions (e.g., Ce 3+ , Eu 2+ , Tl + ) are added. The scintillation efficiency in photons per MeV is determined by how many electron and hole pairs are created, how efficiently their energy is transported to the radiative centers, and the ability of the radiative centers to produce useful optical photons when returning to the ground state. The best timing precision will be achieved if the radiative centers are excited rapidly and efficiently and then radiate rapidly and efficiently. Examples of rapid excitation include CaF 2 :Eu, Lu 2 SiO 5 :Ce, CeF 3 , BaF 2 , ZnO:Ga and the plastic scintillator BC422 (Derenzo et al 2000), which have rise times 0.030 ns. Examples of delayed excitation include Gd 2 SiO 5 :Ce, LaBr 3 :Ce, CsI:Tl, and most plastic scintillators, which have rise times on the order of 0.5 ns. Energy transport to the radiative centers that occurs in a time shorter than the decay time appears as a rise time of the scintillation light.

Thermal quenching
Excited radiative centers can be quenched by a number of non-radiative processes. Configurational quenching results when the excitation of the radiative center causes the lattice to distort to a configuration that is the same as that of a high-temperature ground state. This provides a rapid de-excitation path that is mediated by lattice phonons. An equivalent picture is that the Stokes shift (the difference between the optical absorption and emission energies) is so large that it leaves little energy for emission. Thermal phonons rather than optical photons can then return the system to the ground state. Another quenching process occurs in luminescent ions (e.g., Ce 3+ ) when the excited electron in the 5d orbital are thermally activated into the conduction band and are lost to electron traps in the host. None of these processes reduces the initial number of radiative centers. They shorten the decay time but have little effect on the rise time or the timing precision. Thermal quenching is reduced at lower temperatures and can result in essentially complete quenching at higher temperatures.

Luminosity and decay time
As will be seen in the following sections, the number of scintillation photons and the decay time are important factors for the timing precision. Combined with the photon transport efficiency and the quantum efficiency of the photodetector, they determine N pe /τ d , the number of photoelectrons ns −1 , which is one of the most important factors that determine the timing precision. Scintillators used for fast timing typically have luminosities >20 000 photons MeV −1 and decay times <40 ns. The fundamental limits are 200 000 photons MeV −1 and ∼1 ns, respectively.

Absorption and re-emission of scintillation light
Scintillation photons emitted from a single radiative center will have a spread of energies (the emission band), and the same radiative center can be excited by photons with a range of energies (the absorption band). If these bands overlap, some of the scintillation photons will be absorbed and re-emitted. These photons will be re-emitted with the decay time of the radiative center, and only the earliest will be useful for timing purposes. The primary effect is an increase in the observed decay time, which decreases the factor N pe /τ d and degrades the timing precision (Glodo et al 2010).

Light collection efficiency
Some of the scintillation light is absorbed by the external reflectors and by internal defects, and some may also be trapped by total internal reflection. The number of photons initially created is EL, and the number entering the photodetector is ELF, where E is the energy deposition (MeV), L is the scintillation luminosity (photons MeV −1 ), and F is the fraction of the light entering the photodetector. Ideally, the sum of the F values for the two photodetectors (figure 1) should equal one.

Optical photon dispersion in the scintillator
Optical photons are emitted isotropically from points along the ionization track and reach the photodetector at different times, depending on their paths in the scintillator. The time distribution depends on the specific geometry and surface treatment of the scintillator and can be determined accurately only by experimental measurement or by a Monte Carlo code designed for the purpose. The number of potential geometry and surface treatment combinations is too large to present a comprehensive treatment, but we present results from several representative cases. Figure 2 shows the simulated time dispersion for an emission point at the center of a 3 mm × 3 mm × 30 mm crystal of LSO with double-ended readout, calculated with Geant4 (Agostinelli et al 2003). For this calculation the four side surfaces were simulated as polished and covered with reflective Teflon. The resulting time dispersion can be modeled with the function exp(−t/d) and an optical transport decay time d = 0.042 ns. The values of d for emission points 7 and 23 mm from the photodetector are d = 0.024 ns and 0.059 ns, respectively. As expected, the time dispersion increases with increasing distance to the photodetector. The arrival time of the first photon also depends on the interaction position as well as the index of refraction, as described in section 11. The above values are for excitation on the central axis of the crystal, and distributions for emission points 1.4 mm from the central axis are essentially the same. Similarly, simulations of etched surfaces give essentially the same values for d and the time delay as for polished surfaces. Increasing the size of the crystal to 10 mm × 10 mm × 100 mm, the best-fit value of d for an emission point at the center is 0.127 ns and scales very closely with the size of the crystal, as expected (figure 2). Double-ended readout has the advantage over single-ended readout in that photons emitted away from one photodetector are collected by the opposite photodetector rather than arriving later after many more reflections.

The photodetection process
The photodetector converts scintillation photons into photoelectrons at a photocathode or into electron-hole pairs in a semiconductor. This process is essentially instantaneous and occurs with a probability called the quantum efficiency. These conversion products are then drifted and amplified by electric fields to produce useful pulses. For estimating timing precision, the most important characteristics of the photodetector are the quantum efficiency as a function of photon wavelength, the shape of the single photoelectron output pulse, and the output pulse time jitter. Specific examples are given in later sections.

Electronic noise
Most modern photodetectors have sufficient internal gain so that individual photoelectron pulses are greater than the noise of any downstream electronics. However this noise can still cause variations in timing by varying the time at which the trigger level is reached. If the photodetector output has a slope dV/dt at the trigger level, a voltage variation V will result in a timing variation t = V/(dV/dt). Increasing the slope improves the timing precision. It will be seen in section 12 that electronic noise is not an important factor for systems with good timing precision.
The proportionality constant (ns 1/2 ) that relates the timing precision R (fwhm) to A −1/2 (R = BA −1/2 ) βE g Energy required to create an electron-hole (e − -h + ) pair, expressed as a multiple of the band gap E g (eV) F Fraction of light entering the photodetector (i.e., light collection efficiency) d Time dispersion (ns) of optical photons in the scintillator, modeled as exp Distance between the ends of the scintillator (i.e., the spacing between the two photodetectors) C Duration in time (ns) of the Cherenkov light pulse reaching the photodetector E Energy (MeV) deposited in the scintillator in the form of ionization L Scintillator luminosity (L = SQ S /βE g ) (photons per MeV)

H
Maximum amplitude of the photodetector output pulse (in units of the single photoelectron response) J Jitter (Gaussian fwhm in ns) in the time between the creation of a photoelectron at the input of the photodetector and the corresponding output pulse n Index of refraction of the scintillator at the wavelength of the scintillation light N pe Number of photoelectrons produced in the photodetector, the product of the energy deposited E (MeV), the scintillator luminosity L (photons MeV −1 ), the light collection efficiency F, and the quantum efficiency Q P of the photodetector (N pe = ELFQ P ) N g The number of ionization events used in the Monte Carlo calculation Q S Radiative efficiency of the radiative centers in the scintillator Q P Quantum efficiency of the photodetector at the emission wavelength. This is the fraction of incident optical photons that produce a useful output pulse.

R
Timing precision, the jitter in the time from a localized ionization event in the scintillator to the timing discriminator trigger (fwhm ns). Does not include variations in the depth of interaction.

S
Transfer efficiency of ionization to the radiative centers (i.e. excitations per e − -h + pair) P r The rise time (ns) of the photodetector single-photoelectron response (SER), if modeled as a bi-exponential (i.e. equation (1)) P d The decay time (ns) of the SER, if modeled as a bi-exponential Trigger level of the timing discriminator in units of the peak of the single photoelectron output pulse τ d Decay time of the scintillator (ns) τ r Rise time of the scintillator (ns)

Introduction
In this section we describe the Monte Carlo code used to compute the fwhm of the distribution of time differences from a localized ionization event to the discriminator trigger. It assumes that the photoelectrons created in the photodetector are amplified and that the resulting pulses are summed to produce a single analogue pulse at the input of a single timing discriminator. In section 14 we compare this approach to the lower bound where each photoelectron is amplified in a separate channel and the photodetection time is estimated by a maximum likelihood fit of the theoretical probability density function (PDF) to the distribution of time stamps of all the amplified pulses. Table 1 describes the factors used in the calculations. The scintillator luminosity L (photons MeV −1 ) is given by L = 10 6 Q S S/(βE g ), where βE g is the energy in eV required to create an ionization electron-hole pair, S is the efficiency by which that energy is transferred to radiative centers, and Q S is the radiative efficiency of those centers (Lempicki et al 1993). For the most luminous scintillators (e.g., SrI 2 :Eu, CsBa 2 I 5 :Eu) S and Q are close to one and their yield is 100 000 photons MeV −1 . Unfortunately the Eu 2+ decay time is so long that these scintillators do not provide good timing precision. Faster scintillators (e.g., Lu 2 SiO 5 :Ce, LaBr 3 :Ce) have the shorter decay time of Ce 3+ and are those currently used for timing gamma rays and annihilation photons. Core-valence scintillators (e.g., BaF 2 ) are ultra-fast but do not provide a significant improvement in timing precision due to their low scintillation photon yield.
The number of photoelectrons N pe is the product of the deposited energy E, the scintillator luminosity L, the efficiency of transport to the photodetector F, and the quantum efficiency Q P of the photodetector (N pe = ELFQ P ).
For a scintillator with a single rise time τ r and decay time τ d , the rate of emission of photons that will produce photoelectrons in the photodetector is given by The cumulative integral of f (t) is given by The optical photon time dispersion in the scintillator is modeled by u(t) The photodetector introduces a Gaussian time jitter described by v(t), where the fwhm of the distribution is converted to the standard deviation using σ = J/2.355 To generate a random distribution of N pe emission times distributed according to equation (1), the times were found for which g(t) (equation (2)) was equal to a set of random numbers uniformly distributed between 0 and N pe . The logarithm of a set of uniformly distributed random numbers is then added to these times to include the exponentially distributed optical time dispersion in the scintillator (equation (3)). Finally, a set of Gaussian distributed times is added to include the time jitter in the photodetector (equation (4)). These processes are described in detail in the following section.

Steps in the Monte Carlo calculation
(1) Compute the g(t k ) (equation (2)) on a time grid t k . To provide a fine grid at early times and a large dynamic range, we used t 0 = 0 and t k = t k-1 + (0.0001 ns)(1.05) k-1 (2) Generate N pe random numbers R n uniformly distributed from 0 to N pe (each represents a photoelectron). Any statistical variations in the pulse heights about their average are introduced here. (3) Use linear interpolation of the g(t k ) array to solve R n = g(t n ) and generate a list of photoelectron times t n . These will be randomly distributed with an average distribution f (t).
(4) For a scintillator, convolve with the exponential optical photon dispersion in the scintillator (equation (3)) by subtracting d ln(R k ) from each photoelectron time, where R k is a set of random numbers uniformly distributed between 0 and 1. For a Cherenkov counter, convolve with the temporal distribution of the light by adding a random number CR k , where R k is a set of random numbers uniformly distributed between 0 and 1. (5) To convolve with the transit time jitter of the photodetector (equation (4)), add a random Gaussian time jitter with mean zero and fwhm = J to each photoelectron time. (6) Add the single electron response function for each photoelectron into a histogram array with 0.001 ns time bins. After summing over all photoelectrons this histogram contains the output pulse of the photodetector for a single ionization event. (7) Find the earliest time when the output pulse histogram reaches T (the leading-edge discriminator trigger level) and histogram those times into an array with 0.001 ns time bins. (This describes the times for individual ionization events.) To simulate electronic noise, do a random Gaussian variation of T (fwhm ns). (8) Loop back to step 2 for each ionization event (9) Compute the fwhm = 2.355 rms of the discriminator time distribution accumulated in steps 2-7.
We found that for all values τ r and J computed, the optimal trigger level was proportional to A. An example of this will be shown in section 11. These proportionality constants were used to compute the starting value for the trigger level and four values above and below. A bi-quadratic was fit to the nine values to obtain an accurate value for the optimal trigger level and the corresponding timing precision.
In almost all cases 100 000 ionization events were used in the calculations. The standard deviation in R is R(2N g -1) −1/2 , where N g is the number of ionization events. For N g = 100 000 one standard deviation is R/447.

Individual photoelectron timing statistics
In this section we explore the time jitter distributions of the individual photoelectron pulses. These distributions are generated at step (5) of the Monte Carlo procedure listed in section 3.2.
Post and Schiff calculated that for τ r = 0 and J = 0 the first photoelectron has a timing jitter rms = τ d /N pe (the same as the photon spacing in time) and that the jitter rms of the nth following photoelectron is larger by the square root of n (Post and Schiff 1950). This is shown in the lowest curves of figure 3 for the first ten photoelectrons at N pe /τ d = 100 photoelectrons ns -1 and in figure 4 for the first 100 photoelectrons at N pe /τ d = 1000 photoelectrons ns -1 . Increasing the rise time τ r from 0 to 0.1 ns increases the average time between the first photoelectrons and so increases all their timing jitters. Increasing the photodetector timing jitter J from 0 to 0.1 ns fwhm makes a qualitative difference by introducing a shallow leading-edge to the probability distribution. This makes the timing jitter of the first photoelectron pulse larger than those that follow. The photoelectron pulses that contain the best timing information lie in a broad minimum as seen in figures 3 and 4. Figure 5 shows the PDF for eight combinations: τ r = 0.0 and 0.1 ns, d = 0.0 and 0.1 ns, and J = 0.0 and 0.1 ns fwhm. The scintillator decay time was chosen to be 10 ns, but this has little effect on the shape of the curves shown.

Cherenkov counters
In this section we present the timing precision that can be achieved by detecting the Cherenkov light produced when a charged particle moves though a medium faster than the speed of light in the material. This is a simpler situation than scintillation detectors whose timing accuracy depends on depth of interaction, rise time and decay time.
To estimate the timing precision in this case, we model the light reaching the photodetector as a rectangular pulse in time with width C. If the track length is short, the light pulse is essentially instantaneous. The photodetector single photoelectron output pulses are modeled with a bi-exponential having a rise time of 0.2 ns, a decay time of 2 ns, and a Gaussian time jitter J (fwhm ns). Figure 6 plots the timing precision as a function of N pe for several cases.
For the case C = 0, the Cherenkov pulse is a delta function. The N pe delta functions from a single event will have a Gaussian distribution of width J (fwhm) and their average will have a Gaussian distribution with fwhm J N pe −1/2 . Table 2 shows that for C = 0, RN pe 1/2 is only slightly larger than J (the lowest possible value) and that the optimized leading-edge discriminator does almost as well as the best statistical estimator.

Analysis of six scintillators and three photodetectors
In this section we present calculations of the timing precision for three example photodetectors and six example scintillators. The three photodetectors are a fast phototube    (R-9800 photomultiplier tube (PMT)), a microchannel plate (Planacon MCP) and a silicon photomultiplier (SiPM). Their shapes are shown in figure 7. For the PMT we use the single photoelectron response function derived in ref (Choong 2009), namely V(t) = (t/P t ) Pn exp(−t/P t ) where P n = 4, P t = 0.27 ns, and the transit time jitter J = 0.26 ns fwhm. For the MCP we use the same expression with P n = 6 and P t = 0.12 ns, and J = 0.1 ns fwhm.  Since the rise and decay times of the SiPM arise from different processes (the Geiger discharge and the voltage recovery through an external resistor, respectively) we use the biexponential V(t) = [exp(−t/P d ) − exp(−t/P r )]/(P d − P r ), where P r = 0.2 ns, P d = 20 ns, and J = 0.1 ns fwhm. It will be shown in section 8 that varying P r from 0 to 0.3 ns has very little effect on the timing precision. Table 3 lists the six example scintillators: five commonly available and an 'ideal' scintillator that approaches fundamental limits in terms of luminosity and decay time. The 18 combinations of three photodetectors and six scintillators are shown in tables 4-9. The photon time dispersion d = 0.1 ns was used in the tables, which corresponds to scintillators between 30 and 100 mm deep (figure 2).
The best literature measurements for scintillator timing precision at 511 keV are 0.099 ns fwhm for Lu 1.8 Y 0.2 SiO 5 :Ce and 0.068 ns fwhm for LaBr 3 :30%Ce, using a SiPM with 0.28 ns fwhm time jitter (Seifert et al 2012a); 0.073 ns fwhm for LaBr 3 :30%Ce using a Hamamatsu H5321 photomultiplier with 0.16 ns fwhm time jitter (Glodo et al 2005); and 0.10 ns fwhm for BaF2 using a XP2020Q photomultiplier with 0.6 ns fwhm transit time jitter (Laval et al 1983). These compare favorably with the Monte Carlo calculated values 0.094 ns fwhm for LSO, 0.057 ns fwhm for LaBr 3 :30%Ce, and 0.074 ns fwhm for BaF 2 , listed in tables 5, 6, and 7, respectively. The scintillators currently used in time-of-flight applications (e.g., LSO, LaBr 3 :Ce) are far from fundamental limits in terms of luminosity, decay time, and stopping power. The fundamental upper limit in luminosity is about 200 000 photons MeV −1 , set by the lowest band gap (about 2.5 eV) in which a radiative transition can produce photons compatible with available photodetectors. A high transfer and radiative efficiency is also necessary, such as occurs in SrI 2 :Eu and CsBa 2 I 5 :Eu that have luminosities of 100 000 photons MeV −1 , limited by their band gaps. One of the most luminous scintillators is ZnS, activated by the donor Cl − and the acceptor Ag + , but the decay is very slow. The fundamental limit in the decay time of a scintillator is ∼1 ns, set by the oscillator strength of a fully allowed radiative transition. Examples include organic scintillators and donor-acceptor scintillators such as ZnO:Ga and CdS:In. Core-valence scintillators (e.g., CsF, BaF 2 ) also have short decay times, but their luminosities are low. The fundamental limit in stopping power would be achieved by highdensity compounds of elements of the highest possible atomic number. One possible nextgeneration high-performance scintillator that might approach these limits could be a dense Tl, Hg, Pb, or Bi semiconductor activated by donor and acceptor impurities. Previous work along these lines showed that crystals of PbI 2 could have a luminosity of 40 000 photons MeV −1 and an initial intensity over five times that of LSO, but the donors and acceptors were native defects and not optimized (Derenzo et al 2013). For 511 keV annihilation photons, a scintillator with an optical photon time dispersion of 0.05 ns, a luminosity of 140 000 photons MeV −1 , a decay time of 1 ns, and producing 10 000 photoelectrons in a photodetector with a timing jitter of 0.1 ns fwhm, the timing precision would be 0.008 ns fwhm.

Timing precision versus trigger level
In this section we present plots of timing precision versus the leading-edge trigger level. The Monte Carlo code was used to compute the timing precision as a function of leading-edge trigger level for several cases (figure 8). As described in section 3 the code convolves the photoelectron time distribution with the single-photoelectron photodetector output pulse. This acts to average the individual photoelectron time distributions seen in figures 3 and 4 over the width of the output pulse. When the photodetector timing jitter J is zero, triggering on the first photoelectron provides the best timing precision, even for non-zero values of the rise time τ r .
Increasing J to 0.1 and 0.5 ns changes the leading edge to the gradual rise of the Gaussian function as seen in figure 5. This increases the variance of the timing of the first photoelectron pulse and the spacing of the first few photoelectron pulses. Increasing the leading-edge trigger

Timing precision versus trigger level for different photodetector rise times
In this section we compute the timing precision versus the leading-edge trigger for output pulse rise times P r = 0.   , table 10). Other parameters are N pe /τ d = 1000, τ r = 0, and d = 0. While the optimal trigger level depends strongly on the value of P r , there is much less effect on the timing precision.

Timing precision for the general case
In this section we present Monte Carlo calculations of the timing precision R (fwhm ns) for many different values of the A factor (A = N pe /τ d ), the optical photon dispersion exponential coefficient d, the scintillator rise time τ r , and the photodetector Gaussian jitter J. The photodetector single electron response was modeled with a bi-exponential function (i.e. equation (1)) with rise time P r = 0.2 ns and decay time P d = 2.0 ns. Figure 10 shows the timing precision versus A = N pe /τ d for three cases: (1) τ r = 0.0 ns, d = 0.2 ns, (2) τ r = 0.2 ns, d = 0.0 ns, and (3) τ r = 0.1 ns, d = 0.1 ns for J = 0 ns ( figure 10). This shows that τ r and d are symmetric and that the timing precision is exactly the same when they are interchanged. However, when both are non-zero the leading edge of the probability distribution has a gradual rise, as can be seen in figure 5. Figure 10 also shows the timing precision for τ r = 0.0 and 0.5 ns when d = 0 and J = 0.5 ns fwhm. Figure 11 shows B versus A for the same cases (R = BA −1/2 ). The fact that all the curves are nearly parallel (on a log-log plot) implies that they all follow the same power law relationship, up to a multiplicative constant, which in this case is the constant B that relates the timing precision to A.
In this section we use the Monte Carlo code described in section 3 to compute the timing precision for a range of values of τ r , d, J, and A. Tables 11 to 15 list the B factors that relate the timing precision R (fwhm ns) to the A factor (A = N pe /τ d ), the optical photon dispersion exponential coefficient d, the scintillator rise time τ r , and the photodetector Gaussian jitter J.
The precision is ± 1 in the last decimal place. These tables support an important conclusion. It has long been known in the ideal case that the timing precision is proportional to the inverse square of the initial photoelectron rate A, that is, R ∼ A −1/2 . In this work we find that when virtually all the relevant effects are considered, to a very good approximation R = BA −1/2 . This implies that the consequences of all these effects can be incorporated into a single constant B, and although the magnitude of B depends on the magnitude of these effects, B is near one (between 0.5 and 2.0) for combinations of parameters that span virtually all known scintillation materials and photodetectors.

An analytical formula for the B factor
In this section we present an analytical formula that can be used to estimate the timing precision R of a scintillation detector for any value of τ r from 0 to 1.0 ns, d from 0 to 0.2 ns, J from 0 to 0.5 ns (fwhm), and A ≥ 10 photoelectrons ns −1 . This formula was found using the Eureqa system (Nutonian, Inc., www.nutonian.com) to fit the 656 Monte Carlo values in tables 11-15 plus 164 additional values for A = 2000 and 5000 photoelectrons ns −1 . This system tries successively more complex analytical formulas and retains those that have the best success Table 11. B factors (R = BA −1/2 ) for scintillators with rise time τ r = 0.0.   in fitting the data. The user decides which formula provides the best compromise between complexity and goodness of fit. The result is equation (5), the square root of a polynomial in the variables τ r , D, J, and A that has 16 terms. The symmetry B(τ r , d, J, A) = B(d, τ r , J, A) was imposed on the solution. The 16 coefficients were varied to minimize the rms deviation between the 820 Monte Carlo values and the formula. The best fit had an rms deviation of 0.020 and a maximum deviation of 0.058. Figure 12 shows that all 820 formula values plotted against the Monte Carlo values Table 13. B factors (R = BA −1/2 ) for scintillators with rise time τ r = 0.2 ns.     To use tables 11-15 and equation (5) it is necessary to know the scintillator rise and decay times, the optical photon dispersion in the scintillator, the number of photoelectrons produced, and the timing jitter of the photodetector. With the exception of the photon time dispersion in the scintillator, reasonable estimates of all these values can generally be obtained from the literature or the photodetector manufacturer. As described in section 2.8, the photon dispersion in the scintillator depends on the scintillator geometry and the depth of interaction, and while less widely known it can be estimated for a specific system using Monte Carlo simulation and can be validated by measurements using a collimated radiation source as described in (Moses and Derenzo 1999).

Combining the photodetector information from dual-ended readout
In this section we review the pulse height and timing information provided by dual-ended photodetector readout and their relationship to the depth of interaction and the time at which the nuclear particle (e.g., gamma ray, neutron) entered the scintillator. The two photodetectors (A and B) are connected to individual amplitude and time digitizers that provide the output pulse amplitudes N A and N B , and the leading-edge trigger times t A and t B (figure 1). A digital processor converts these signals into the best estimates of energy deposited, the interaction depth, and the time when the radiation entered the scintillator. In the following analyses we place photodetector A at z = −D/2 and photodetector B at z = +D/2, where D is the length of the scintillator (i.e. the spacing between the two photodetectors).
One common method for estimating the depth of interaction uses the ratio of the two photodetector signals N A and N B (Moses and Derenzo 1994). Yang et al (Yang et al 2006) used avalanche photodiodes coupled to opposing ends of 20 mm long LSO crystals. They found linear relationships between the amplitudes N A and N B and the depth of interaction, and they were able to achieve a depth of interaction uncertainty of 0.3 cm fwhm. Their experimental data can be modeled as where a was approximately 0.4. This leads to the estimate The gamma ray or neutron enters the scintillator at t = 0, travels at a speed c and interacts at depth z at time t = (z + D/2)/c. Timing discriminator A triggers at time t A = (z + D/2)/c + (z + D/2)(n/c) + δ A . Timing discriminator B triggers at time t B = (z + D/2)/c + (−z + D/2)(n/c) + δ B . Here n is an effective index of refraction, and δ A and δ B are delays in the photodetector and electronics. These lead to another estimator for the depth of interaction By subtracting the signal delays from the leading-edge discriminator times t A and t B we have two separate estimators (t 1 and t 2 ) for the time of entry: An algorithm for the best estimator for the time of entry is beyond the scope of this paper, but it is important to point out the following.
• A weighted average of t 1 and t 2 will require combining the uncertainties in z, t A , and t B to determine the weights. • For a given uncertainty in z, t 2 will be a better estimator than t 1 because the dependence on z is less. • The dispersion factor d and the number of photoelectrons N pe will depend on the depth of interaction. Their values are needed for each ionizing event to estimate the uncertainties in t A and t B from equation (5). These dependences can be measured by using a collimated radiation source as described in (Moses and Derenzo 1999).

Effect of electronic noise on the timing precision
In this section we compute the effect of electronic noise on the timing precision for the cases of τ r = 0.5 ns, τ d = 10 ns, d = 0.1 ns, J = 0.2 ns fwhm, and N pe = 100, 1000, 10 000, and 100 000 photoelectrons (table 16). Because the photodetector output pulse has a finite slope dV/dt, a noise variation V in the output or the discriminator will produce a time variation t = V/(dV/dt). In the case of N pe = 100 the optimum trigger level is only 0.9 of the single photoelectron pulse amplitude, and a noise level of 0.25 fwhm was used. The timing precision is dominated by photoelectron statistics, and trigger level noise has a relatively small effect. In the other three cases with larger photoelectron statistics, even a noise level of 1.0 fwhm has a minor effect since it is much smaller than the optimum trigger level. Systems designed for good timing precision should be able to detect individual photoelectrons, and these will have an electronic noise below 1.0 fwhm. We conclude that electronic noise in modern photodetectors and electronics is not a serious limiting factor in scintillation detector timing precision.

Theoretical lower bound of timing precision
In this section we use a particular example (τ r = 0.5 ns, τ d = 10 ns, d = 0.1 ns, J = 0.2 ns fwhm) to compare the timing precision of leading-edge discrimination of the photodetector output (used in sections 5-9) with the theoretical lower bound. Table 17 lists the timing precision R for N pe = 100, 1000, 10 000, and 100 000 photoelectrons at the optimal trigger level T.
The theoretical lower bound was estimated by the following steps.
(1) Generate N pe random photoelectron pulse time stamps distributed according to the theoretical PDF (described below).
(2) Perform a maximum likelihood fit of the PDF to all the time stamps to determine the best-fit time.
(3) Repeat steps (1) and (2) for 100 000 ionization events (4) Compute the standard deviation of the 100 000 times and multiply by 2.355 to determine the lower bound fwhm of the timing precision R(ML) shown in table 17.
In this case the leading-edge trigger timing precision R(T opt ) is within 15% of the lower bound R(ML). We conclude that summing photoelectron pulses at the output of the photodetector and adjusting the trigger level of the timing discriminator for the best timing precision can be a simple and effective means for estimating the photodetection time. Similar comparisons have been done by Seifert (Seifert et al 2012a).
We derived the theoretical PDF of the photoelectron pulses h(t) (equation (6)) by convolving the photon emission PDF (equation (1)) with the time dispersion of the photons in the scintillator (equation (3)) and then with the photodetector time jitter (equation (4)) Note that h(t) is unchanged when d and τ r are interchanged, so in this sense the photon dispersion acts as a rise time, but only if the rise time is zero. When both are non-zero, the timing precision is worse than if either of them were set equal to their sum and the other were set equal to zero, as may be seen from the τ r d product terms in equation (5), in figure 5, and in tables 11-15.

Discussion
We have used Monte Carlo methods to predict the timing precision of scintillation detectors that includes all the physical parameters known to affect the precision. We find that the timing precision depends primarily on A = N pe /τ d (number of photoelectrons)/(scintillator decay time in ns) and follows the form R = B/A 1/2 , where R is the timing precision. While the other factors that affect timing (scintillator rise time τ r , time dispersion d of optical photons in the scintillator, and transit time jitter J in the photodetector) affect the factor B, this factor has a very narrow range of values (0.5 < B < 2) for virtually all reasonable values of these other factors and we empirically derived a formula that allows others to obtain an accurate estimate of B.
In a typical case (τ r = 0.5 ns, d = 0.1 ns, J = 0.2 ns fwhm), leading-edge discrimination with the optimal threshold provides a timing estimate that is within 15% of the theoretically best possible estimate, and the timing precision is only weakly dependent on the value of this threshold. While simple theory (τ r = d = J = 0) places this optimal threshold at the single photoelectron level, even very small non-zero values for these other effects (which are almost certain to be present in actual detector systems) cause the optimal threshold to be raised to higher levels. In these cases, we found that the optimal trigger level was proportional to the number of photoelectrons. As a result, the leading-edge threshold should be set to a fixed fraction of the signal amplitude.
As the factor that dominates timing precision is A = N pe /τ d , efforts to improve timing precision should concentrate on this factor, which includes scintillator luminosity and decay time, as well as photodetector quantum efficiency. Values of A up to 10 000 photoelectrons ns −1 are theoretically possible and would yield a timing precision of 0.008 ns fwhm for a 511 keV photon.
Finally, we note that all the timing resolution measurements presented herein are for interactions at a specific position in a scintillation crystal. The penetrating nature of gamma rays implies that most systems will have interactions at many positions within the crystals, and while the timing resolution will depend on the interaction position because of the position dependence of optical dispersion, this position dependence can be estimated. As described in sections 2.8 and 11, there will also be a shift in time due to scintillation photon propagation from the interaction point to the photodetector, but this is easily calculable. Thus, the timing resolution of a system can be predicted using the results presented herein, using the expected interaction position distribution, and folding that in with the appropriate position-dependent timing resolution and time delays.

Conclusions
• The optical photon dispersion in a scintillator with dual-ended readout can be modeled as single exponential optical transport decay with time constant d. • The timing precision R of a scintillation detector can be parameterized as R = BA −1/2 , where A = N pe /τ d (number of photoelectrons)/(scintillator decay time in ns), and B is between 0.5 and 2 for virtually all reasonable values of the scintillator rise time τ r , the photon time dispersion d, and the photodetector time jitter J. • An analytical formula for B(A, τ r , d, J) was found that fit 820 Monte Carlo calculations of the timing precision with a root mean square deviation of 2.2% of the value of B. • For virtually all the cases that we explored, the optimal leading-edge trigger was found to be proportional to the number of photoelectrons.