Time resolution and efficiency of SPADs and SiPMs for photons and charged particles
TTime resolution and efficiency of SPADs and SiPMs for photons andcharged particles
W. Riegler a , P. Windischhofer b a CERN b University of Oxford
Abstract
We give an analytic treatment of the time resolution and efficiency of Single Photon Avalanche Diodes(SPADs) and Silicon Photomultipliers (SiPMs). We provide closed-form expressions for structures withuniform electric fields and efficient numerical prescriptions for arbitrary electric field configurations. Wediscuss the sensor performance for single photon detection and also for charged particle detection.
Preprint submitted to Elsevier February 2, 2021 a r X i v : . [ phy s i c s . i n s - d e t ] J a n . Introduction Single Photon Avalanche Diodes (SPADs) have been used for photon detection and photon count-ing since many decades. These semiconductor devices contain a highly doped p-n junction of 0 . − µ mthickness, the so called gain region, that is biased above the breakdown voltage. This means that asingle primary electron or hole entering this region can produce a diverging avalanche through impactionization and therefore lead to a detectable signal. The growth of the diverging avalanche is quenchedby the breakdown of the electric field, which leads to a digital type signal with an amplitude independentof the number of primary charges. The restoration of the electric field is governed by a quench resistoror quench circuit external to the device. n + p + p p-epi E elec. γp + n E |E||E| w=10-100μmd=0.5-2μmd=0.5-2μm c on v e r s i on r eg i on ga i n r eg i onga i n r eg i on a) b)γ elec. γ elec. Figure 1: a) p-in-n SPAD with a gain layer only. b) n-in-p SPAD with a gain layer and a conversion layer.
Two basic SPAD geometries are shown in Fig. 1. For the p-in-n SPAD in Fig. 1a the high field regionis placed close to the surface of the silicon. Photons with a wavelength of <
500 nm have an absorptionlength of less than 1 µ m, so they will be efficiently absorbed in the gain region and trigger the breakdown.For photons with longer wavelengths and therefore larger absorption length, the geometry of Fig. 1b witha so called ’conversion region’ will be more efficient. This region consists of a depleted layer of siliconwith thickness ranging from 10-100 µ m that is adjacent to a gain region. The electrons created in theconversion region will drift to the gain region where they provoke the breakdown. The sensor can beilluminated from the front and the back side. These SPADs are nowadays heavily used for LIDAR [1]applications with wavelengths in the near-infrared above 750 nm.While originally being a single channel device, the advances in silicon industry allowed the arrange-ment of many of these diodes in chessboard structures. Depending on the readout circuitry these arecalled SPAD image sensors or Silicon Photomultipliers (SiPMs). In SPAD image sensors, the SPADpixels are read out individually in order to form an image. In SiPMs the individual pixels are connectedin parallel through series resistors into a single channel. The amplitude of the output signal will thenbe proportional to the total number of fired pixels and therefore proportional to the number photonsthat have hit the pixel matrix, which represents the function of a traditional photomultiplier. Thereare many technological challenges related to the implementation of such pixel structures, specifically theelimination of crosstalk between the channels and the minimization of the dark count rate.2 cope and outline In this report, we discuss the time resolution and efficiency of SPADs and SiPMs of the two types shownin Fig. 1 and evaluate their performance for the detection of photons as well as charged particles. Ourresults are derived from a series of fundamental equations that describe the movement and the interac-tions of the participating charge carriers.We provide analytic expressions for the simplified structures in Fig. 2 with constant electric fields. Fig. 2arepresents the situation where a photon interacts directly inside a gain layer of constant electric field E that is above the breakdown value. Fig. 2b shows a SPAD where we represent the conversion region bya silicon layer of thickness w with constant electric field E below the breakdown field, which is adjacentto a gain layer with constant field E . w hole d E E γ |E| q γ elec. γ elec. elec.elec. b)a) c on v e r s i on l a y e r ga i n l a y e r Figure 2: Two simplified SPAD structures with constant electric fields discussed in this report: a) A photon interactsdirectly in the gain layer and produces an e-h pair that provokes breakdown. b) A photon produces an e-h pair in theconversion layer and the electron moves to the gain layer where it provokes breakdown. A thin gain layer of 0.5-2 µ m is alsohighly efficient for charged particle detection. This report is structured as follows: Section 2 discusses the absorption of photons in silicon, and thecontribution of the conversion layer to the time resolution. Section 3 then discusses the mechanisms ofavalanche breakdown and the average growth of the avalanche. Section 4 describes the contribution tothe efficiency resulting from the avalanche formation in the gain layer. Section 5 discusses avalanchefluctuations and their contribution to the time resolution. This discussion is based to a large degree onthe companion paper [2], which develops the statistics of electron-hole avalanches in great detail. Section6 then discusses the performance of SPADs for the detection of charged particles. We finally drop theassumption of constant electric field in Section 7 and return to the realistic field configurations of Fig. 1.We give efficient numerical prescriptions that extend the analytic results obtained previously.Details of all calculations are given in several appendices. Although we focus specifically on devices basedon silicon, our results are expected to also cover the basic geometries for different types of semiconductors.3 . Conversion layer
We assume a layer of silicon of thickness w extending from x = 0 to x = w as shown in Fig. 3.The probability for a photon to be absorbed between position x and x + dx is given by P ( x ) dx =1 /l a e − x /l a dx , where l a is the photon absorption length from Fig. 4a. The efficiency, i.e. the probabilityfor a photon to convert in the layer, is then given by p = 1 − e − w/l a and the numbers are shown in Fig. 4b.Photons of wavelength <
500 nm are efficiently absorbed in < µ m of silicon while infrared photons of >
750 nm need several tens of µ m of silicon to be absorbed efficiently. x=0 E γ x=wx=x xγ elec. a)b) Figure 3: The conversion layer of a SPAD or SiPM. A photon is absorbed at position x producing an e-h pair. The electronwill drift to the gain layer at x = w . a)
300 400 500 600 700 800 900 1000 1100 12000.0010.10010100010 wavelength ( nm ) a b s o r p t i on l e ng t h ( μ m ) b) μ m100 μ m10 μ m3 μ m1 μ m0.5 μ m
300 400 500 600 700 800 900 1000 1100 12000.1151050100 wavelength ( nm ) a b s o r p t i on e ff i c i e n cy ( % ) Figure 4: Absorption length l a for photons of different wavelengths in silicon [3, 4]. b) Photon absorption efficiency fordifferent values of silicon thickness. First we investigate the case where the photon is arriving from the ’left’ side as shown in Fig. 3a.Normalising the conversion probability to the efficiency, the probability for a photon to be absorbedbetween position x and x + dx is P ( x ) dx = 11 − e − w/l a l a e − x /l a Θ( w − x ) dx (1)The electron is then moving to the edge of the silicon layer at x = w with a velocity v e , where itarrives at time t = ( w − x ) /v e . The velocity of electrons and holes in silicon is shown in Fig. 5a andthe parametrization is given in Appendix A. The arrival time distribution of the electron at x = w istherefore ρ ( t ) = (cid:90) w P ( x ) δ [ t − ( w − x ) /v e ] dx = wl a ( e w/l a −
1) 1
T e twlaT Θ( T − t ) (2)4) v e v h v * - - ( V / cm ) V e l o c i t y ( μ m / p s ) b) nodiffusiondiffusion ( ps ) ρ ( t ) Figure 5: a) Drift velocity of electrons ( v e ) and holes ( v h ) as a function of electric field in silicon. The velocity v ∗ =2 v e v h / ( v e + v h ) that is relevant for the avalanche growth in the gain layer is shown as well. b) Probability for the electronto arrive at x = w between times t and t + dt for w = 10 µ m, for a photon with l a = 1 µ m entering the layer from the leftside. The ’no diffusion’ curve refers Eq. 2 and the ’diffusion’ curve refers to Eq. 6. where we have expressed the velocity v e by the maximum drift time T = w/v e of the electrons inside theconversion layer. An example is shown in Fig. 5b. The variance of the arrival time is then σ t = (cid:90) T t P ( t ) dt − (cid:32)(cid:90) T t P ( t ) dt (cid:33) = T (cid:18) l a w −
14 sinh[ w/ l a ] (cid:19) (3)Including the effect of diffusion we use the fact that an electron deposited at position x at t = 0 andmoving with an average velocity of v e will be found at position x after a time t with a probability of p ( x, x , t ) dx = 1 √ π √ Dt exp (cid:20) − ( x − ( x + v e t )) Dt ) (cid:21) dx (4)The standard deviation of the distribution is given by σ ( t ) = √ Dt . The probability for an electronto arrive at x = w between time t and t + dt is then p ( w, x , t ) v e dt . Since the probability of a photonabsorption at x is given by P ( x ) from Eq. 1 we find the arrival time distribution P ( t ) as ρ ( t ) = (cid:90) w P ( x ) p ( w, x , t ) v e dx = (5) v e l a e w/l a − e ( D + l a v e ) t/l a (cid:32) Erf (cid:34)(cid:32) √ Dtl a + v e (cid:114) tD (cid:33)(cid:35) − Erf (cid:34) √ Dtl a + v e (cid:114) tD − w √ Dt (cid:35)(cid:33) (6)with Erf( z ) = 2 / √ π (cid:82) z e − t dt . The variance evaluates to σ t = T (cid:18) l a w −
14 sinh[ w/ l a ] (cid:19) + T Dv e (cid:18) − e − w/l a − l a w (cid:19) + 8 D v e (7)The first term is the one from Eq. 3 due to the varying position of the photon interaction together withthe average drift time, the second term is due to diffusion. The third term is an artefact of the assumptionthat a charge placed at x = 0 can diffuse into a region of x <
0, so the term does not vanish even for w = 0. In a realistic implementation of a conversion layer, the region of x < l a /w the time resolutionapproximates to σ t = T
12 +
DTv e l a (cid:29) w σ t = T l a w + 2 DTv e l a (cid:28) w (8)We assume D = 35 cm / s for electrons in silicon. For electric fields in excess of 5 × V/cm the electronvelocity is close to saturation and we have D/ ( v esat ) = 0 .
35 ps. For a conversion layer of w = 1 / / µ mthe saturated drift time is T ≈ / / l a (cid:29) w the probability for the photon conversion position becomes uniform across w , diffusion isnegligible and the time resolution is σ t = T / √ ≈ . / . /
289 ps for w = 1 / / µ m. For l a (cid:28) w the photon conversion point is always close to x = 0, diffusion will dominate and the time resolution isequal to σ t = √ DT /v e ≈ . / . / .
46 ps for w = 1 / / µ m.If the illumination takes place from the ’right’ side as indicated in Fig. 3b, the time resolution becomes σ t = T (cid:18) l a w −
14 sinh[ w/ l a ] (cid:19) + T Dv e (cid:18) l a w + 11 − e w/l a (cid:19) + 8 D v e (9)This expression differs from the previous one only by the diffusion term, because the distance betweenthe conversion point and x = w is now different. In the limit of large and small values of l a we have σ t = T
12 +
DTv e l a (cid:29) w σ t = T l a w + 2 DTv e l a w l a (cid:28) w (10)For l a (cid:29) w the expression is equal to the one from above, while for l a (cid:28) w the variance goes to zerobecause the conversion point is close to x = w . The expression does of course not apply if the absorptionlength l a is of the same order or smaller than the gain layer thickness d , because the photons will interactdirectly in the gain layer. 6 . Electron-hole avalanches and breakdown, average signal elec.holes E d x0 x=x n e0 n h0 Figure 6: The primary electrons and holes deposited at x = x are multiplying, which results in a diverging avalanche incase the electric field E is above the breakdown limit. An electron drifting inside the conversion layer will move to the gain layer and trigger an avalanchestarting from x = 0. Alternatively a photon can convert inside the gain layer and the e-h pair at position x = x will trigger the avalanche. To cover both situations, we treat the general case where n e electronsand n h holes are deposited at x = x at time t = 0, as shown in Fig. 6.To derive equations describing the avalanche, we allow for general position-dependent electric fields E ( x ).With the field orientated as shown in Fig. 6, electrons move to the right and holes move to the left withvelocities v e ( x ) , v h ( x ). The probability for an electron to create an e-h pair when travelling a distance dx is α ( x ) dx while the probability for a hole to produce an e-h pair over distance dx is β ( x ) dx , where α ( x ) and β ( x ) are called the impact ionization coefficients. The values for silicon (Fig. 7) are reportedin [7] with the parameters listed in Appendix A. Since 1 /α and 1 /β refer to the average distance thatan electron or a hole has to travel in order to produce one additional e-h pair, we see that only for fieldsin excess of 2 − × V/cm there is an appreciable probability to provoke an avalanche in a few µ m ofsilicon. Fig. 8 shows a Monte Carlo (MC) simulation of a few avalanches starting with a single electron.a) αβ ( V / cm ) α , β ( / μ m ) b) ( μ m ) E ( V / c m ) Figure 7: a) Impact ionization coefficient α for electrons and β for holes as a function of the electric field. b) Minimumelectric field value provoking breakdown for a gain layer with constant electric field across a given thickness d . After some initial fluctuations the avalanche just grows exponentially. There is also a finite probabilitythat no diverging avalanche develops.We denote as n e ( x, t ) dx and n h ( x, t ) dx the average number of electrons and holes between position x and x + dx at time t (note that this is different from the notation used in [2]). These charge densitiesresult in local average current densities of j e ( x, t ) = v e ( x ) n e ( x, t ) and j h ( x, t ) = − v h ( x ) n h ( x, t ). By the7
10 20 30 40 50110100100010 Time ( ps ) N e + N h Average MCAnalyticMC1MC2MC3MC4
Figure 8: Monte Carlo simulation for an electron-hole avalanche starting with a single electron at x = 0 for a 1 µ m diodeat a field of 4 V/ µ m. After some initial fluctuations the number of charge carriers increases exponentially. The dominantterm from Eq. 30 approximates the average signal extremely well for times t > d/v ∗ ≈
11 ps. When the avalanche is stillsmall there is a finite probability that no breakdown occurs, as is the case for the avalanche MC4. continuity equation ∂j/∂x + ∂n/∂t = σ , with σ the generation rate, we therefore have ∂n e ( x, t ) ∂t + ∂v e ( x ) n e ( x, t ) ∂x = α ( x ) v e ( x ) n e ( x, t ) + β ( x ) v h ( x ) n h ( x, t ) (11) ∂n h ( x, t ) ∂t − ∂v h ( x ) n h ( x, t ) ∂x = α ( x ) v e ( x ) n e ( x, t ) + β ( x ) v h ( x ) n h ( x, t )The fact that electrons move to the left and holes move to the right gives the boundary conditions n e (0 , t ) = 0 n h ( d, t ) = 0 (12)Since Eq. 11 represents a set of linear equations we can use the Ansatz n e ( x, t ) = f ( x ) e St and n h ( x, t ) = g ( x ) e St and we find Sf ( x ) + [ v e ( x ) f ( x )] (cid:48) = α ( x ) v e ( x ) f ( x ) + β ( x ) v h ( x ) g ( x ) (13) Sg ( x ) − [ v h ( x ) g ( x )] (cid:48) = α ( x ) v e ( x ) f ( x ) + β ( x ) v h ( x ) g ( x )with f (0) = 0 and g ( d ) = 0. The multiplication of electrons and holes can lead to a finite amount of totalcharge in the avalanche ( S <
0) or it can diverge and cause breakdown (
S > S = 0, so by setting S = 0 in the above equations and solving them with the givenboundary conditions we find the breakdown condition (Appendix B) (cid:90) d α ( x ) exp (cid:20) − (cid:90) x ( α ( x (cid:48) ) − β ( x (cid:48) )) dx (cid:48) (cid:21) dx = 1 (14)The breakdown condition is independent of v e ( x ) and v h ( x ). This relation is usually derived by evaluatingthe point at which the gain for a constant current injected into the gain layer diverges [8]. Evaluatingthe breakdown equation for constant α and β implies that breakdown occurs if d > α − β ln αβ (15)The electric field at which breakdown takes place for a gain layer with a given thickness d is shown inFig. 7b. For general electric field profiles, Eqs. 11 can be efficiently solved with numerical methods. The8olution can be given in analytical form for constant values of α, β, v e , v h and is derived in [2]. It isrepresented as an infinite sum of exponential terms with (generally complex valued) time constants. Atleast one time constant is guaranteed to be real-valued. The largest real-valued time constant definesthe long-term behaviour of the avalanche. Above the breakdown limit, this term determines the rate ofexponential growth of the avalanche. Starting with n e electrons and n h holes at position x at time t = 0,it reads n e ( x, t ) = 1 d a ( x ) (cid:2) n e u e ( x ) + n h u h ( x ) (cid:3) e γv ∗ t n h ( x, t ) = 1 d b ( x ) (cid:2) n e u e ( x ) + n h u h ( x ) (cid:3) e γv ∗ t (16)with u e ( x ) = e − a x sin( k − kx /d ) (17) u h ( x ) = e − a x αd [ k cos ( k − kx /d ) + λ sin( k − kx /d )] (18) a ( x ) = 2 v h ke a x sin (cid:0) k xd (cid:1) ( v e + v h )(1 + λ ) sin k (19) b ( x ) = 2 kv e e a x (cid:2) k cos (cid:0) k xd (cid:1) + λ sin (cid:0) k xd (cid:1)(cid:3) ( v e + v h ) βd (1 + λ ) sin k (20)and the constants v ∗ , a , γ, k are defined by v ∗ = 2 v e v h v e + v h (21) a = αv e − βv h v e + v h + v e − v h v e + v h d λ (22) γ = α + β λ d (23) k = (cid:113) αβd − λ (24)The parameter λ is the largest real solution of the equation λ + (cid:113) αβd − λ cot (cid:113) αβd − λ = 0 (25)It holds that −∞ < λ < d √ αβ . For λ < − d √ αβ the constant k will become imaginary, which will stilllead to real valued expressions for n e ( x, t ) and n h ( x, t ) with sin , cos becoming sinh , cosh. The functionalform of a ( x ) and b ( x ) as well as the equation for γ were already derived in [11]. Fig. 9a shows the functions u e ( x ) , u h ( x ) , u e ( x ) + u h ( x ) that determine how the average growth of the avalanche depends on theposition of a primary electron, hole or e-h pair. They are the mirror images of the functions a ( x ) , b ( x )from Eqs. 19, 20 that determine the distribution of the electrons and holes inside the gain layer [2].The parameter γ defines the exponential growth of the avalanche and is shown in Fig. 9b. It has thefollowing properties: γ = 0 αβd = αβ ( α − β ) ln βα ≤ γ = α + β − (cid:112) αβ αβd = 1 (27) γ = α + β αβd = π ≈ .
47 (28) γ max = α + β (cid:112) αβ αβd → ∞ (29)9) u e ( x ) u h ( x ) u e ( x )+ u h ( x ) x ( μ m ) b) d = μ md = μ md = μ m γ max α + β ( V / cm ) γ ( / μ m ) Figure 9: a) The functions u e ( x ) , u h ( x ) , u e ( x ) + u h ( x ) from Eqs. 17, 18 that determine how the average growth ofthe avalanche depends on the position of a primary electron, hole or e-h pair. The values are for a 1 µ m gain layer at E = 4 . × V/cm. b) γ as a function of electric field in silicon for different values of the gain layer thickness d . At thebreakdown limit we have γ = 0. For higher electric fields γ quickly approach γ max . The total number of electrons and holes is then given by N = (cid:82) d n ( x, t ) dx and evaluates to N e ( t ) = B e (cid:2) n e u e ( x ) + n h u h ( x ) (cid:3) e γv ∗ t N h ( t ) = B h (cid:2) n e u e ( x ) + n h u h ( x ) (cid:3) e γv ∗ t (30)where we have B e = 2 kv h (cid:2) e a d ( a d − k cot k ) + k csc k (cid:3) ( v e + v h )( a d + k )(1 + λ ) (31) B h = 2 kv e (cid:2) e a d ( k + a dλ + k ( a d − λ ) cot k ) + k ( λ − a d ) csc k (cid:3) ( v e + v h ) βd ( a d + k )(1 + λ ) (32)The total induced current becomes I ( t ) = e E w V w [ v e N e ( t ) + v h N h ( t )] (33)= e E w V w ( v e B e + v h B h ) (cid:2) n e u e ( x ) + n h u h ( x ) (cid:3) e γv ∗ t (34)where e is the electron charge. Here we have assumed a constant weighting field E w /V w in the region inwhich the charges are moving. For the single photon detection using the conversion layer we have only asingle electron at x = 0 and therefore n e = 1 , n h = 0 and the expression is I ( t ) = e E w V w [ v e B e + v h B h ] sin k e γv ∗ t (35)Assuming a velocity v e ≈ v h ≈ v sat ≈ . µ m/ps and a weighting field of E w /V w = 1 /d = 1 / (1 µ m), thecurrent corresponding to 10 charges at t = 44 ps in Fig. 8 is 1.6 mA.10 . Efficiency In this section we calculate the probabilities P e ( x ) and P h ( x ) for a single electron or a single holeplaced at position x in the gain layer to cause breakdown. We follow [9] to establish the equations for thesequantities. We start by considering a single electron at position x − dx which moves in positive x -directionin the applied electric field. The probability for it to create a diverging avalanche is P e ( x − dx ). Between x − dx and x two things can happen. 1) With a probability of (1 − αdx ) there is no multiplication of theelectron and then the electron at position x creates a diverging avalanche or 2) the electron is multiplyingover the distance dx and at least one of the two electrons and the hole create breakdown. This can bewritten as P e ( x − dx ) = (1 − αdx ) P e ( x ) + αdx (cid:2) − (1 − P e ( x )) (1 − P h ( x ) (cid:3) (36)Writing the corresponding equation for P h ( x ) and expanding for small dx gives d P e ( x ) d x = − α ( x )[1 − P e ( x )] [ P e ( x ) + P h ( x ) − P e ( x ) P h ( x )] (37) d P h ( x ) d x = β ( x )[1 − P h ( x )] [ P e ( x ) + P h ( x ) − P e ( x ) P h ( x )]Provided α ( x ) and β ( x ) are known, these equations can be integrated with the boundary conditions P e ( d ) = 0 and P h (0) = 0. Following [10] we define P ( x ) = P e ( x ) + P h ( x ) − P e ( x ) P h ( x ) and by differen-tiating this expression and using Eqs. 37 we have d P ( x ) d x = − ( α − β ) P ( x )[1 − P ( x )] (38)We use the boundary condition P (0) = P e (0) = p , with p still to be determined, giving the solution P ( x ) = p p + (1 − p ) exp (cid:2)(cid:82) x ( α ( x (cid:48) ) − β ( x (cid:48) )) dx (cid:48) (cid:3) (39)Since P ( x ) can be written as P ( x ) = 1 − [(1 − P e ( x ))(1 − P h ( x ))] we see that P ( x ) refers to the breakdownefficiency of a single e-h pair, which we will use later for the efficiency calculation for MIPs. Knowing P ( x ) and using P h ( d ) = 0 we can integrate Eqs. 37 and have P e ( x ) = 1 − exp (cid:32) − (cid:90) dx α ( x (cid:48) ) P ( x (cid:48) ) dx (cid:48) (cid:33) P h ( x ) = 1 − exp (cid:18) − (cid:90) x β ( x (cid:48) ) P ( x (cid:48) ) dx (cid:48) (cid:19) (40)And finally P e (0) = p gives the equation that allows us to determine p p = 1 − exp − (cid:90) d p α ( x (cid:48) ) p + (1 − p ) exp (cid:104)(cid:82) x (cid:48) ( α ( x (cid:48)(cid:48) ) − β ( x (cid:48)(cid:48) )) dx (cid:48)(cid:48) (cid:105) dx (cid:48) (41)In general this equation can only be evaluated numerically. This equation also reveals again the breakdowncondition. Close to the threshold of breakdown the value of p will be small. For small values of p theexpression 1 − exp [ − p / ( p + (1 − p ) e v ))] is approximated by p e − v + O ( p ), so the above relation turnsinto the breakdown condition of Eq. 14.For constant α and β the above expressions evaluate to P ( x ) = p p + [1 − p ] e ( α − β ) x (42) P e ( x ) = 1 − e − α ( d − x ) (cid:20) (1 − p ) e ( α − β ) d + p (1 − p ) e ( α − β ) x + p (cid:21) αα − β (43) P h ( x ) = 1 − e − βx (cid:104) (1 − p ) e ( α − β ) x + p (cid:105) βα − β (44)11q. 41 that determines p reads as e − ( α − β ) d = 1 p (cid:104) (1 − p ) − βα − (1 − p ) (cid:105) (45)Fig. 10 shows a few examples.a) ( μ m ) E ff i c i e n cy e 4x10 V / cmh 4x10 V / cme 6x10 V / cmh 6x10 V / cm b) ( V / cm ) E ff i c i e n cy μ m1 μ m0.5 μ m Figure 10: a) Breakdown probability (efficiency) for a single electron and a single hole deposited at position x inside a gainlayer of d = 1 µ m for two values of the electric field. b) Breakdown probability p (efficiency) for a single electron placed at x = 0 for different values of the gain layer thickness d .
5. Time resolution
The statistics of electron-hole avalanches and the resulting contribution to the time resolution arediscussed in detail in [2], where the problem is treated by using the theory of continuous-time Markovprocesses.
For the case of an e-h avalanche in absence of any boundaries and for a constant electric field, analternative approach can be used to derive the avalanche fluctuations [12] that does not require theformalism developed in [2] and which is closely related to the arguments resulting in Eq. 37. We define p e ( n, m, ∆) to be the probability to find n electrons and m holes at time t + ∆ for an avalanche startingwith a single electron at t . For an avalanche starting at t = 0, there are two ways to reach this stateat the later time t + dt . First, the initial electron does not multiply in the first small time interval [0 , dt ](with probability 1 − αv e dt ), but then produces n electrons and m holes during the subsequent timeinterval [ dt, t + dt ]. This happens with a probability p e ( n, m, t ). Second, the electron already multipliesin the interval [0 , dt ] (with probability αv e dt ) and the resulting two electrons and one hole multiply into n electrons and m holes during [ dt, t + dt ]. This is written as p e ( n, m, t + dt ) = (1 − αv e dt ) p e ( n, m, t ) (46)+ αv e dt n (cid:88) i =1 i (cid:88) j =1 m (cid:88) r =1 r (cid:88) s =1 p e ( n − i − j, m − r − s, t ) p e ( i, r, t ) p h ( j, s, t ) (47)where p h ( n, m, ∆) is the probability that an avalanche starting with a single hole at time t produces n electrons and m holes at t + ∆. Writing the corresponding equation for p h ( n, m, t ) and expandingfor small dt results in the equations defining p e ( n, m, t ) and p h ( n, m, t ). The equations have a structuresimilar to Eq. 37. Their solution is given in Appendix C and equal to the one derived in [2]. Having n e electrons and n h holes at t = 0, the probability to have n additionally created e-h pairs at time t is p ( n, t ) = Γ( A + n )Γ( A )Γ(1 + n ) (cid:18) ν ( t ) (cid:19) A (cid:18) − ν ( t ) (cid:19) n ∞ (cid:88) n =0 p ( n, t ) = 1 (48)12 = n e αv e + n h βv h αv e + βv h ν ( t ) = e λ t t λ t = αv e + βv h (49)The average number of e-h pairs and the variance are n ( t ) = [ ν ( t ) − A σ n ( t ) = ν ( t )[ ν ( t ) − A (50)In the continuous approximation for n , which is valid for large n (and therefore for large avalanchesi.e. late times) we have p ( n, t ) = n A − Γ( A ) (cid:18) An ( t ) (cid:19) A e − nA/n ( t ) (cid:90) ∞ p ( n, t ) dn = 1 (51)with n ( t ) = Ae λ t t σ n ( t ) = 1 √ A n ( t ) (52)In order to measure the ’signal time’ we apply a threshold to the signal, which is proportional to thea) Threshold t t t Time ( ps ) N e + N h b) A = = = - λ t t ρ ( t ) Figure 11: a) A threshold is applied to the signal. The fluctuations of the avalanche size result in a fluctuation of thethreshold crossing time, which defines the time resolution. b) Time response function ρ ( t ) for different values of theparameter A . total number of charge carriers. Fig. 11a shows how the avalanche fluctuations lead to fluctuations of thethreshold crossing time, which determine the time resolution. The probability that the signal crosses thethreshold of n e-h pairs between time t and t + dt , the so-called time response function, is given by ρ ( n, t ) dt = λ t Γ(1 + A + n )Γ( A )Γ(1 + n ) (cid:18) ν ( t ) (cid:19) A (cid:18) − ν ( t ) (cid:19) n dt (cid:90) ∞ ρ ( n, t ) dt = 1 (53)The average threshold crossing time and its variance are given by t = 1 λ t [ ψ ( n + A + 1) − ψ ( A )] σ = 1 λ t (cid:112) ψ ( A ) − ψ ( n + A + 1) (54)where ψ ( x ) = d ln Γ( z ) /dz is the digamma function and ψ ( x ) = d ln Γ( z ) /dz is the trigamma function.For large numbers of n the above time response function approximates to ρ ( n, t ) = λ t Γ( A ) exp (cid:2) A ln n − Aλ t t − ne − λ t t (cid:3) (55)13nd the average threshold crossing time and the time resolution approximate to t = 1 λ t [log n − ψ ( A )] σ = 1 λ t (cid:112) ψ ( A ) (56)so the time resolution becomes independent of the threshold. This can be understood when looking atFig. 11a and it is a well-established fact for detectors like Resistive Plate Chambers [13], where avalanchefluctuations dominate the signal characteristics [14, 15]: scaling the threshold by a constant c will justshift the time response function by ∆ t = (ln c ) /λ t without altering its shape. If we are not interested inthe absolute time of the threshold crossing but just the time variations we can arbitrarily set n = 1 andhave the time response function ρ ( t ) = λ t Γ( A ) exp (cid:2) − Aλ t t − e − λ t t (cid:3) (cid:90) ∞−∞ ρ ( t ) dt = 1 (57)An example of the time response function for different values of the parameter A is shown in Fig. 11b.The function (cid:112) ψ ( x ) approaches (cid:112) /x + 1 /x for small and large values x and it is within 10% in theentire range of x >
0, as can be seen in Fig. 12a. The values of (cid:112) ψ ( A ) for a primary electron, a primaryhole and a primary e-h pair are shown in Fig. 12b.We see that when starting with a primary hole, the time resolution is significantly worse when comparedto a primary electron or a primary e-h pair. An electron travels on average a distance of 1 /α and a holetravels on average a distance of 1 /β before creating an additional e-h pair. Since α > β in silicon theholes will travel significantly longer than the electrons before triggering the breakdown and therefore thetime fluctuations are larger. The geometries of SPADs and SiPMs produced from silicon are thereforebuilt such that the electrons trigger the avalanche in the gain region. Charge carrier mobilities are quitedifferent in other semiconductors and the arrangement of doping layers can therefore differ.a) ψ ( x ) x + x b) n e = n h = n e = n h = n e = n h = ( V / cm ) ψ ( A ) Figure 12: a) The trigamma function ψ ( x ). b) Values of (cid:112) ψ ( A ) for different initial conditions. Starting with a single e-hpair we have A = 1 and (cid:112) ψ (1) = π/ √ ≈ . .2. Avalanches in a gain layer of finite thickness Our central result in [2] is the conclusion that the finite thickness of the gain layer will to first ordernot affect the avalanche fluctuations but will only affect the average growth of the avalanche. Thisapproximation works best if the primary charge carrier has the larger impact ionization coefficient, i.e. ifthe primary charge in the gain layer is either an electron or an e-h pair for silicon. Then, the probabilityto have a total number of n electrons and m holes at time t in the gain layer, starting with n e electronsand n h holes at x = x at time t = 0, is p ( n, t ) ≈ [1 − ε ( x )] δ ( n ) + ε ( x ) n A − Γ( A ) (cid:18) AN e ( t ) (cid:19) A e − nA/N e ( t ) m = B h B e n A = n e αv e + n h βv h αv e + βv h (58)with N e ( t ) from Eq. 30 and B e , B h are from Eqs. 31, 32. The number of electrons n and the number ofholes m are considered continuous variables and taken to be fully correlated in this approximation (asshown in [2], this is strictly true only at late times). The efficiency ε ( x ) that the n e electrons and n h holes at x = x trigger a diverging avalanche is ε ( x ) = 1 − [1 − P e ( x )] n e [1 − P h ( x )] n h (59)with P e ( x ) and P h ( x ) from Eqs. 43, 44. The corresponding time response function in reference toEq. 57 is then ρ ( t ) ≈ γv ∗ h A Γ( A ) exp (cid:104) − Aγv ∗ t − he − γv ∗ t (cid:105) (60)with h = A ε ( x ) B e [ n e u e ( x ) + n h u h ( x )] (61)Here we have divided N e ( t ) by the efficiency in order to account for the fact that only diverging avalanchesare crossing the threshold. The corresponding time resolution due to the avalanche fluctuations, whenkeeping n e , n h , x constant, is σ = (cid:112) ψ ( A ) γv ∗ (62)where γ and v ∗ are from Eq. 21 and Eq. 23. The value of 1 /γv ∗ for different values of the thickness d ofthe gain layer is shown in Fig. 13a.a) d = μ md = μ md = μ m1 / γ max v * /( α + β ) v * ( V / cm ) / γ v * ( p s ) b) ( V / cm ) γ v * σ σ av σ pos d = μ m σ pos d = μ m σ pos d = μ m Figure 13: a) The factor 1 /γv ∗ for different values of the gain layer thickness. b) Contributions to the time resolution fora photon interacting in the gain layer. The value of π/ √ For a photon interacting in the conversion layer (Fig. 2b), the electron arriving at the gain layer willtherefore start an avalanche from x = 0 and the time resolution contribution of the gain layer will be σ ≈ (cid:112) ψ ( A ) γv ∗ A = αv e αv e + βv h (63)15s seen in Fig. 12 we have (cid:112) ψ ( A ) ≈ . µ m SPADat 4 × V/cm the contribution to the time resolution will be around 6 ps. For higher fields a timeresolution of 1 ps should theoretically be achievable.In case the photon interacts inside the gain layer and produces an e-h pair, the varying position of thephoton interaction will also contribute to the time resolution. Having one e-h pair at position x as shownin Fig. 2a we have n e = n h = 1 and therefore A = 1, ε ( x ) = P ( x ) from Eq. 42 and the time responsefunction is ρ ( t, x ) = γv ∗ h ( x ) exp (cid:104) − γv ∗ t − h ( x ) e − γv ∗ t (cid:105) h ( x ) = P ( x ) B e [ u e ( x ) + u h ( x )] (64)The probability of conversion at position x is given by p ( x ) dx = 11 − e − d/l a l a e − x /l a dx (65)and the time response function including the fluctuation of the conversion point is therefore defined by ρ ( t ) = (cid:90) d p ( x ) ρ ( t, x ) dx (66)The time resolution σ = (cid:82) t ρ ( t ) dt − ( (cid:82) tρ ( t ) dt ) has then two components, a contribution σ av fromavalanche fluctuations and a contribution σ pos from the varying position of the primary e-h pair in thegain layer ( γv ∗ ) σ av = π γv ∗ ) σ pos = (cid:90) d p ( x ) [ln h ( x )] dx − (cid:32)(cid:90) d p ( x ) ln h ( x ) dx (cid:33) (68)In the case of a very small photon absorption length i.e. l a (cid:28) d the primary e-h pair will always be createdat the very edge of the sensor which is equal to the situation of x = 0 and the time resolution is givenby σ av . In the other extreme of l a (cid:29) d there will be a uniform distribution for the position of the photoninteraction in the gain layer and the contribution to the time resolution is given in Fig. 13b. Whether theavalanche fluctuations or the position fluctuations dominate depends on the gain layer thickness and theelectric field. We can conclude that the time resolution is well approximated by σ ≈ c /γv ∗ , where themain dependence is given by the variation of γ with sensor thickness and electric field (Fig. 13a) while c ≈ − v ∗ is saturated at ≈ . µ m/ps.The term σ pos can also be derived directly from the average avalanche growth. Neglecting avalanchefluctuations the primary e-h pair will simply trigger an average avalanche according to Eq. 30 N ( t ) = B e u e ( x ) + u h ( x ) P ( x ) e γv ∗ t (69)where we have divided by the efficiency to account for the avalanches that do not cross the threshold.Applying a threshold N thr to this signal gives a threshold crossing time of t ( x ) = 1 γv ∗ [ln N thr − ln h ( x )] (70)and the variance t − t is the one from Eq. 68. 16 . Charged particle detection with SPADs elec.holes E x=d xx=0 q Figure 14: A charged particle is leaving clusters of e-h pairs in the gain layer.
SPADs can also be used for detection of charged particles. Charged particles interacting with siliconproduce clusters of e-h pairs along their track, with an average distance of λ ≈ . µ m for MIPs. ASPAD sensor of 1 µ m or 2 µ m thickness will therefore be highly efficient to charged particles and there isno need for a conversion layer. The probability p clu ( n ) for having n > /n distribution, but there are significant deviations at small numbers of n insilicon. Fig. 15a shows the cluster size distribution for silicon as calculated with HEED [16]. Assuming aa) - - - / cluster P r ob a b ili t y b) μ m 1 / n μ m 1 / n ● μ m HEED ■ μ m HEED ● ● ● ● ●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●■ ■ ■ ■ ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ × - × - P ( n ) Figure 15: a) Cluster size distribution i.e. the probability for a single cluster to contain n electron-hole pairs as calculatedwith HEED [16]. b) Probability p ( n, d ) to find n e-h pairs inside a silicon layer of thickness d . The solid lines refer to an1 /n distribution in both plots. sensor thickness d and an average distance between clusters of λ , the average number of clusters in thesensor is n = d/λ . Assuming the cluster size distribution p clu ( n ), the probability p ( n, d ) to find n e-hpairs in the sensor can then be calculated by using the Z-transform as [15] P clu ( z ) = ∞ (cid:88) n =1 p clu ( n ) z n G ( z ) = e d/λP clu ( z ) − e d/λ − p ( n, d ) = 1 n ! (cid:20) d n G (1 /z ) dz n (cid:21) z =0 (71)It is shown in Fig. 15b. 17 .1. Efficiency First we want to calculate the efficiency for a charged particle to cause breakdown in the gain layer.Since P ( x ) from Eq. 39 is the probability for a single e-h pair at position x to trigger a diverging avalanche,the probability Q ( x ) that this single e-h pair does not cause breakdown is Q ( x ) = 1 − P ( x ) = 11 + p − p exp (cid:2) − (cid:82) x ( α ( x (cid:48) ) − β ( x (cid:48) )) dx (cid:48) (cid:3) (72)The probability dq that there is no breakdown caused by the charged particle traversing the interval[ x, x + ∆ x ] is given by a) the probability 1 − ∆ x/λ that there is no interaction in ∆ x and b) theprobability that there is an interaction but it does not lead to breakdown, which we can write as dq = (cid:18) − ∆ xλ (cid:19) + ∆ xλ ∞ (cid:88) n =1 p clu ( n ) Q n ( x ) (73)From this, the probability q that the charged particle does not create a diverging avalanche in any of theslices ∆ x is derived in Appendix D and the efficiency p = 1 − q is given by p = 1 − exp (cid:34) − λ (cid:32) d − ∞ (cid:88) n =1 p clu ( n ) (cid:90) d Q ( x ) n dx (cid:33)(cid:35) (74)If Q ( x ) = 0 i.e. if an e-h pair deposited at x will definitely cause breakdown, we have p = 1 − e − d/λ , whichis the correct probability that there is at least one interaction within d . The evaluation for constant α, β is given in Appendix E and shown in Fig. 16 together with the efficiency for a single electron starting at x = 0 in the gain layer. Above the breakdown field the efficiency rises steeply to the maximum level1 − e − d/λ . For a gain layer thickness of d > µ m the efficiency is larger that 99 % above the breakdownlimit. ( V / cm ) E ff i c i e n cy μ m1 μ m0.5 μ m Figure 16: Efficiency for a MIP to provoke breakdown. Beyond the breakdown field, the efficiency rises sharply from 0 tothe level of 1 − e − d/λ . The thin lines correspond to the efficiency for a single electron at x = 0. .2. Time resolution We first consider the simpler case where the MIP produces a single cluster of m e-h pairs at position x . Then, we have A = m , n e = n h = m and the time response function is ρ ( t, m, x ) ≈ γv ∗ ( m − h ( x ) m exp (cid:104) − mγv ∗ t − h ( x ) e − γv ∗ t (cid:105) (75) h ( x ) = 1 − [1 − P ( x )] m B e [ u e ( x ) + u h ( x )] ≈ B e [ u e ( x ) + u h ( x )] (76)In case the cluster size varies according to p clu ( m ) and the probability to have the cluster at position x varies according to p ( x ), the time response function becomes ρ ( t ) = ∞ (cid:88) m =1 (cid:90) d p clu ( m ) p ( x ) ρ ( t, m, x ) dx (77)and the related time resolution has three contributions( γv ∗ ) σ = ∞ (cid:88) m =1 p clu ( m ) ψ ( m ) (78)+ ∞ (cid:88) m =1 p clu ( m ) ψ ( m ) − (cid:32) ∞ (cid:88) m =1 p clu ( m ) ψ ( m ) (cid:33) (79)+ (cid:90) d p ( x )[ln h ( x )] dx − (cid:32)(cid:90) d p ( x ) ln h ( x ) dx (cid:33) (80)The first term represents the average of the avalanche fluctuations, where ψ ( m ) is decreasing from ψ (1) = π / ≈ /m dependence. For the cluster size distribution in silicon from Fig. 15athis first term evaluates to ≈ . m e-h pairs will on average grow as me γv ∗ t . The term evaluates to ≈ .
39 for the cluster size distribution in silicon, significantly larger thanthe contribution from the avalanche fluctuations.The third term represents the dependence on the position of the primary cluster and for a uniformprobability this term evaluates to the values already shown in Fig. 13b.In conclusion we therefore observe that assuming a single e-h cluster at a random position in the gainlayer, the avalanche fluctuations are negligible and only the average growth of the avalanche as well asthe position dependence play a role. In a regime where the contribution from the position dependenceis negligible ( <
1) only the fluctuation of the total charge is important, which results in a universaldependence of the time resolution on the thickness of the gain layer( γ v ∗ ) σ = ∞ (cid:88) m =1 p ( m, d ) ψ ( m ) − (cid:32) ∞ (cid:88) m =1 p ( m, d ) ψ ( m ) (cid:33) (81)Here, p ( m, d ) is the probability that the passing MIP produces m e-h pairs in the gain layer of thickness d . The resulting time resolution is shown in Fig. 17a. Since γv ∗ σ is close to unity for typical dimensionsof the gain layer, the time resolution of a SPAD for a MIP is essentially defined only by γ and v ∗ andthe values are the ones given in Fig. 13a. Since the cluster size distribution p clu ( m ) has a long tailtowards large values of m , the same is true for the time response function. The standard deviation of thethreshold crossing time is then generally not identical to the parameter σ extracted from a Gaussian fitto the distribution of the threshold crossing time. Both measures are compared in Figure 17a.19) π / r.m.s. HEEDGausfit HEEDr.m.s. 1 / n Gaussfit 1 / n ( μ m ) γ v * σ b) ( V / cm ) γ v * σ μ m0.5 μ m1 μ m1 μ m2 μ m2 μ m Figure 17: a) Time resolution for a MIP when the fluctuations of the cluster position can be neglected, for differentvalues of the gain layer thickness d . The solid lines refer to p clu ( n ) for silicon from Fig. 15a, while the dashed lines assume p clu ( n ) ∼ /n . The time resolution is quantified by the standard deviation of the threshold crossing time (“r.m.s.”) aswell as by the parameter σ extracted from a Gaussian fit. b) Time resolution for a MIP for a gain layer of 0 . , , µ mthickness, taking the fluctuation of the cluster positions into account. The dashed lines refer to the numbers from a). Next, we consider the general case where the MIP produces a variable number of clusters, all of whichfluctuate in size. Since avalanche fluctuations are negligible we perform this calculation by following thediscussion around Eq. 70. We divide the gain layer into N + 1 slices of thickness ∆ x = d/ ( N + 1) andassume that a charged particle leaves m n primary e-h pairs in the n th slice. The efficiency for this caseis very close to unity, and the additional small dependence on the cluster distribution is neglected here.Then, the total average number of charges in the avalanche becomes N tot ( t ) = N (cid:88) n =0 m n h ( n ∆ x ) e γv ∗ t (82)with h ( x ) taken from Eq. 76. Applying a threshold N thr to this signal and shifting it by a constant offsetof ln N thr we find a threshold crossing time of t ( m , m , ..., m N ) = − γv ∗ ln (cid:34) N (cid:88) n =0 m n h ( n ∆ x ) (cid:35) (83)The probability p ( m, ∆ x ) to find m e-h pairs in a slice ∆ x is given by p ( m, ∆ x ) = (cid:18) − ∆ xλ (cid:19) δ m + ∆ xλ p clu ( m ) (84)so the average threshold crossing time and the second moment are t = − lim N →∞ γv ∗ ∞ (cid:88) m =0 ∞ (cid:88) m =0 ... ∞ (cid:88) m N =0 p ( m , ∆ x ) p ( m , ∆ x ) ...p ( m N , ∆ x ) ln (cid:34) N (cid:88) n =0 m n h ( n ∆ x ) (cid:35) (85) t = lim N →∞ γv ∗ ) ∞ (cid:88) m =0 ∞ (cid:88) m =0 ... ∞ (cid:88) m N =0 p ( m , ∆ x ) p ( m , ∆ x ) ...p ( m N , ∆ x ) ln (cid:34) N (cid:88) n =0 m n h ( n ∆ x ) (cid:35) (86)These relations are evaluated in Appendix F, giving a contribution to the time resolution of( γv ∗ ) σ = (cid:34)(cid:90) ∞ w ( y ) ln ydy − (cid:18)(cid:90) ∞ w ( y ) ln ydy (cid:19) (cid:35) (87)with W ( s ) = exp (cid:104) dλ d (cid:82) d P clu (cid:16) s h ( x ) h (0) (cid:17) dx (cid:105) − e d/λ − w ( y ) = L − [ W ( s )] (88)20ere, P clu ( s ) is the Laplace transform of the cluster size distribution and the operator L − denotes theinverse Laplace transform. The evaluation is shown in Fig. 17b. The dashed lines show the time resolutionfrom Fig. 17a where the contribution from position fluctuations is neglected. We see that for fields aroundthe breakdown limit the effect from the positions variations is small and it increases with increasing field.Overall, the time resolution for MIPs stays within σ t = (0 . − . /γv ∗ for the parameters investigated.
7. Realistic field configuration
In this section we finally discuss a realistic field configuration of a SPAD and we apply the insightsfrom all previous sections to assess its performance. Fig. 18a shows an example for the electric field ina SPAD created by a highly doped p-n junction. The specific functional form is defined in AppendixA, Eq. 95. With the impact ionization and drift velocity parameters from Appendix A we obtain α ( x ), β ( x ), v e ( x ) and v h ( x ) as explicit functions throughout the sensor. The impact ionization coefficients α ( x )and β ( x ) are shown in Fig. 18b. For the purposes of our discussion here, we define x = 0 . µ m and x = 1 . µ m as the boundaries of the gain layer, which thus has a thickness of d = 1 . µ m.a) ( μ m ) E ( V / c m ) b) αβ ( μ m ) α , β ( / μ m ) Figure 18: a) Electric field in a realistic SPAD or SiPM. b) Impact ionization coefficients inside the sensor with parametersfrom Appendix A.
Efficiency.
For this geometry, the integral in Eq. 14 evaluates to 1.39 which is larger than unity andtherefore guarantees that breakdown can take place. To find the efficiency of the sensor we solve Eqs. 37numerically, using as boundary conditions P h ( x ) = 0 and P e ( x ) = 0. The solution is shown in Fig. 19atogether with the corresponding efficiencies obtained from a MC simulation of the avalanche development.The efficiency for a MIP passing this sensor can be calculated with Eq. 74 using Q ( x ) = 1 − P eh ( x ). Withthe cluster size distribution from Fig. 15 and λ = 0 . µ m, the efficiency evaluates to p = 1 − . × − .A SPAD of this kind is a highly efficient detector for a MIP. Average signal and contribution to time resolution.
To study the average growth of the avalanche, wesolve Eqs. 11 with n e ( x ) = 0 and n h ( x ) = 0 as boundary conditions. This yields the average chargedensities in the gain layer, n e ( x, t ) and n h ( x, t ). The average total charge present in the gain layer canbe obtained through a numerical integration of these densities. This quantity is proportional to theaverage signal produced by the avalanche. As shown in Section 3, it grows exponentially as e St . Thetime constant S can be directly extracted from the numerical solution and evaluates to S = 0 .
48 ps.In case the position x of the primary charge fluctuates, it generates a contribution to the timeresolution according to Eq. 80. The magnitude of this effect can be estimated from Fig. 19b, which showsthe time at which the average signal crosses the applied threshold. If the position x of the initial chargeis uniformly distributed, the resulting contribution to the time resolution is 2.7/2.1/2.4 ps for an initialelectron, an initial hole, and an initial e-h pair. 21) ● ● ● ● ● ● ●● ● ● ● ● ● ●● ● ● ● ● ● ● x ( μ m ) E ff i c i e n cy P e P h P eh ● MC b) ● ● ● ● ● ● ●● ● ● ● ● ● ●● ● ● ● ● ● ● x ( μ m ) C r o ss i ng t i m e ( p s ) electronholee - h ● MC Figure 19: a) Breakdown probability for a primary electron ( P e ), primary hole ( P h ) and primary e-h pair ( P eh ) placed atposition x in the gain layer. b) Time at which the average total number of charges in the avalanche (proportional to theaverage signal) crosses a threshold of 10 charges. The avalanche is initiated by a primary electron, a primary hole or aprimary e-h pair placed at position x in the gain layer. For both plots, the numerical solutions give rise to the lines, themarkers correspond to the values from MC simulations. Avalanche fluctuations and contribution to time resolution.
According to the discussion in Section 6.2,we expect the contribution of the time resolution from fluctuations in the avalanche development, σ av ,to be of the order of 1 /S ≈ σ av ≈ (cid:112) ψ ( A ) /S . This formula neglects effects due to the finite size of the gain region and was originallyderived for constant impact ionization coefficients, which enter into the computation of the parameter A .For position-dependent electric fields, the largest values of α and β in the gain layer are relevant for theformation of the avalanche fluctuations and can be used to compute A . As the comparison with resultsfrom MC in Fig. 20 shows, this estimates the time resolution for the case of an initial electron to withinat most 20%. As expected, the approximation becomes worse if the initial charge includes holes, whichhave a low impact ionization coefficient. In this case, corrections due to the finite size of the gain layerbecome important. These can also be computed numerically, as shown in [2], but the calculations aremore involved. ● ● ● ● ● ● ●▲ ▲ ▲ ▲ ▲ ▲ ▲■ ■ ■ ■ ■ ■ ■ x ( μ m ) σ av ( p s ) electronholee - h ● MC Figure 20: Time resolution σ av due to avalanche fluctuations, for an avalanche initiated by a primary electron, a primaryhole or a primary e-h pair. The primary charge is placed at position x in the gain layer. The approximation from Eq. 62neglects the x -dependence and gives rise to the horizontal lines. The markers correspond to the values obtained from MCsimulations. . Conclusions We have performed a detailed study of the time resolution and efficiency of SPADs and SiPMs for thedetection of photons and charged particles. Our discussions start from a series of differential equations,which cover the conversion and the drift of charges in the conversion layer as well as the formation of theavalanche in the gain layer. For arbitrary electric field profiles, the equations for the average avalanche de-velopment as well as the breakdown efficiency can be easily solved using numeric solvers. The calculationof the avalanche fluctuations and their impact on the time resolution is more involved in this case and a de-tailed discussion is given in [2]. We have provided analytic solutions for the case of constant electric fields.For the detection of single photons, the contribution of the conversion layer of thickness w to the timeresolution for constant electric field is σ = w/ ( v e √ σ = c γv ∗ (89)with c = 0 . − . . − µ m. This relation holds for single photondetection and MIP detection. It also extends to realistic non-uniform electric fields. Both contribu-tions from avalanche fluctuations as well as the variation of the photon or MIP conversion point in thegain layer are captured. The constant γv ∗ determines the average growth of the avalanche according to N ( t ) ∝ e γv ∗ t , with v ∗ ≈ . µ m/ps and γ saturating at γ max ≈ α + β at high fields. It should be possiblein practice to limit this contribution to the level of a few picoseconds at high fields.The efficiency of a SPAD or SiPM for photons has many contributions, including the photon conver-sion probability, the geometry and fill factor of the sensor as well as the breakdown probability in thegain layer. The contribution from the breakdown probability can be easily calculated by numericallysolving the related equations. For SPADs with a conversion layer or for photons with absorption length < µ m that are absorbed close to the edge of the gain layer this efficiency quickly approaches valuesclose to 100% when biasing the sensor beyond the breakdown field.SPADs or SiPMs with a gain layer of 1 − µ m thickness should be highly efficient for MIP detection.A dedicated conversion layer is not necessary.This report discussed ’one-dimensional’ sensors. For realistic implementations of SPAD pixel sensorsand SiPMs the boundaries of the pixels together with all the elements for limitation of optical crosstalkmake up complex three dimensional electric fields. To study these sensors, the 3D field map togetherwith a full MC simulation with programs like Garfield++ [17] has to be used and our results can serveas benchmarks for these simulations. 23 . Appendix A The velocity of electrons and holes in silicon is shown is parametrized by v e ( E ) = µ e E (cid:20) (cid:16) µ e Ev esat (cid:17) β e (cid:21) /β e v h ( E ) = µ h E (cid:20) (cid:16) µ h Ev hsat (cid:17) β h (cid:21) /β h (90)The parameters from [5] are µ e = 1417 cm /Vs, µ h = 471 cm /Vs, β e = 1 . β h = 1 .
213 and v esat =1 . × cm/s and v hsat = 0 . × cm/s at 300 K . The impact ionization parameters α and β as reported in [7] are given by α ( E ) = α ∞ e − a/E β ( E ) = β ∞ e − b/E (91)with α ∞ = 7 . × cm − a = 1 . × V /cm . × ≤ E ≤ . × V /cm (92) β ∞ = 1 . × cm − b = 2 . × V /cm . × ≤ E ≤ . × V /cm (93) β ∞ = 6 . × cm − b = 1 . × V /cm . × ≤ E ≤ . × V /cm (94)
For a realistic SPAD we assume the following electric field: E ( x ) = E exp (cid:104) − ( x − µ ) /σ − e − ( x − µ ) /σ (cid:105) E = 5 × V/cm µ = 1 µ m σ = 0 . µ m (95)
10. Appendix B
For S = 0 Eq. 13 reads as[ v e ( x ) f ( x )] (cid:48) = α ( x ) v e ( x ) f ( x ) + β ( x ) v h ( x ) g ( x ) (96) − [ v h ( x ) g ( x )] (cid:48) = α ( x ) v e ( x ) f ( x ) + β ( x ) v h ( x ) g ( x ) (97)with boundary conditions f (0) = 0 and g ( d ) = 0. Subtracting the two equations gives[ v e ( x ) f ( x )] (cid:48) + [ v h ( x ) g ( x )] (cid:48) = 0 → v e ( x ) f ( x ) = − v h ( x ) g ( x ) + c (98)Inserting this expression into Eq. 96 we have[ v h ( x ) g ( x )] (cid:48) − v h ( x ) g ( x )[ α ( x ) − β ( x )] = − c α ( x ) (99)with the general solution v h ( x ) g ( x ) = c − c (cid:82) x α ( x (cid:48) ) exp[ − (cid:82) x (cid:48) ( α ( x (cid:48)(cid:48) ) − β ( x (cid:48)(cid:48) )) dx (cid:48)(cid:48) ] dx (cid:48) exp[ − (cid:82) x ( α ( x (cid:48) ) − β ( x (cid:48) )) dx (cid:48) ] (100)The condition f (0) = 0 implies c = c and g ( d ) = 0 then implies (cid:90) d α ( x ) exp (cid:20) − (cid:90) x ( α ( x (cid:48) ) − β ( x (cid:48) )) dx (cid:48) (cid:21) dx = 1 (101)24his is the general breakdown condition which is independent of the electron and hole velocities. Ex-pressing v h g from Eq. 98 and inserting this expression into Eq. 97 we have v h ( x ) g ( x ) = − v e ( x ) f ( x ) + c [ v e ( x ) f ( x )] (cid:48) − v e ( x ) f ( x )[ α ( x ) − β ( x )] = c β ( x ) (102)and therefore v e ( x ) f ( x ) = c (cid:82) x β ( x (cid:48) ) exp[ − (cid:82) x (cid:48) ( α ( x (cid:48)(cid:48) ) − β ( x (cid:48)(cid:48) )) dx (cid:48)(cid:48) ] dx (cid:48) + C exp[ − (cid:82) x ( α ( x (cid:48) ) − β ( x (cid:48) )) dx (cid:48) ] (103)The condition f (0) = 0 gives C = 0 and g ( d ) = 0 gives (cid:90) d β ( x ) exp (cid:20) − (cid:90) x ( α ( x (cid:48) ) − β ( x (cid:48) )) dx (cid:48) (cid:21) dx = exp (cid:34) − (cid:90) d ( α ( x (cid:48) ) − β ( x (cid:48) )) dx (cid:48) (cid:35) (104)which is equal to (cid:90) d β ( x ) exp (cid:34)(cid:90) dx ( α ( x (cid:48) ) − β ( x (cid:48) )) dx (cid:48) (cid:35) dx = 1 (105)As shown in [8] Eq. 101 and 105 are identical because (cid:90) d ( α ( x ) − β ( x )) exp (cid:20) − (cid:90) x ( α ( x (cid:48) ) − β ( x (cid:48) )) dx (cid:48) (cid:21) dx = 1 − exp (cid:34) − (cid:90) d ( α ( x (cid:48) ) − β ( x (cid:48) )) dx (cid:48) (cid:35) (106)
11. Appendix C
The relation that determines the number of electrons n and the number of holes m at a given timewhen starting with a single electron is defined in Eq. 46 as p e ( n, m, t + dt ) = (1 − αv e dt ) p e ( n, m, t ) (107)+ αv e dt n (cid:88) i =1 i (cid:88) j =1 m (cid:88) r =1 r (cid:88) s =1 p e ( n − i − j, m − r − s, t ) p e ( i, r, t ) p h ( j, s, t )Establishing the corresponding equation for p h ( m, n, t ) and expanding for small dx we have1 v e dp e ( n, m, t ) dt = − αp e ( n, m, t ) (108)+ α n (cid:88) i =1 i (cid:88) j =1 m (cid:88) r =1 r (cid:88) s =1 p e ( n − i − j, m − r − s, t ) p e ( i, r, t ) p h ( j, s, t )1 v h dp h ( n, m, t ) dt = − βp h ( n, m, t ) (109)+ β n (cid:88) i =1 i (cid:88) j =1 m (cid:88) r =1 r (cid:88) s =1 p h ( n − i − j, m − r − s, t ) p h ( i, r, t ) p e ( j, s, t )The Z-transform of these equations is1 v e ∂P e ( z , z , t ) ∂t = − αP e ( z , z , t )[1 − P e ( z , z , t ) P h ( z , z , t )] (110)1 v h ∂P h ( z , z , t ) ∂t = − βP h ( z , z , t )[1 − P h ( z , z , t ) P e ( z , z , t )] (111)25he equations have a structure similar to 37 and we therefore form the function P ( z , z , t ) = 1 − P e ( z , z , t ) P h ( z , z , t ) (112)Differentiating this equation and using the above relations gives ∂P∂t = ( αv e + βv h )(1 − P ) P (113)with the solution P ( t ) = e ( αv e + βv h ) t e ( αv e + βv h ) t + c (114)The initial conditions that there is one electron at t = 0 for p e and one hole for p h reads as p e ( n, m, t =0) = δ n, δ m, and p h ( n, m, t = 0) = δ n, δ m, and we have therefore P e ( z , z , t = 0) = 1 z P h ( z , z , t = 0) = 1 z → P ( z , z , t = 0) = 1 − z z (115)and P ( z , z , t ) = e ( αv e + βv h ) t ( z z − e ( αv e + βv h ) t ( z z −
1) (116)We can now write Eqs. 110 and 111 as d ln P e dt = − αv e P d ln P h dt = − βv h P (117)Integrating the equations with the above initial conditions we finally have P e ( z , z , t ) = 1 z (cid:20) z z e ( αv e + βv h ) t ( z z − (cid:21) αveαve + βvh (118) P h ( z , z , t ) = 1 z (cid:20) z z e ( αv e + βv h ) t ( z z − (cid:21) βvhαve + βvh (119)In case the avalanche is initiated by n e electrons and n h holes at time t = 0, we are interested in theprobability p ( n, m, t ) to find n electrons and m holes at time t . In terms of p e ( n, m, t ) and p h ( n, m, t ),it is expressed as an iterated convolution in analogy to the right-hand sides of Eqs. 108 and 109. In the z -domain, P ( z , z , t ) reads P ( z , z , t ) = P e ( z , z , t ) n e P h ( z , z , t ) n h (120)= 1 z n e z n h (cid:20) z z e ( αv e + βv h ) t ( z z − (cid:21) n eαve + n hβvhαve + βvh (121)The inverse Z-Transform of this expression also gives access to p ( n, t ) shown in Eq. 48, which is definedas the probability to find n e-h pairs that are created in addition to the initial n e electrons and n h holes. p ( n, t ) = Γ( A + n )Γ( A )Γ(1 + n ) (cid:18) ν ( t ) (cid:19) A (cid:18) − ν ( t ) (cid:19) n ∞ (cid:88) n =0 p ( n, t ) = 1 (122) A = n e αv e + n h βv h αv e + βv h ν ( t ) = e ( αv e + βv h ) t (123)26
2. Appendix D
We divide the sensor into N + 1 slices of thickness ∆ x = d/ ( N + 1). The probability that there is nobreakdown caused by the particle traversing the slice [ x, x + ∆ x ] is given by the probability 1 − ∆ x/λ that there is no interaction in ∆ x and the probability that there is an interaction but it does not lead tobreakdown. dq = (cid:18) − ∆ xλ (cid:19) + ∆ xλ ∞ (cid:88) n =1 p clu ( n ) Q n ( x ) (124)= 1 − ∆ xλ (cid:32) − ∞ (cid:88) n =1 p clu ( n ) Q n ( x ) (cid:33) (125):= 1 − ∆ xλ f ( x ) (126)The probability q that there is no breakdown in any of the slices of ∆ x throughout the sensor is then q = (cid:20) − ∆ xλ f (0) (cid:21) (cid:20) − ∆ xλ f (∆ x ) (cid:21) (cid:20) − ∆ xλ f (2∆ x ) (cid:21) ... (cid:20) − ∆ xλ f ( N ∆ x ) (cid:21) (127)ln q = N (cid:88) n =0 ln (cid:20) − dN λ f ( nd/N ) (cid:21) (128) ≈ N (cid:88) n =0 − dN λ f ( nd/N ) (129) ≈ − λ (cid:90) ∞ f ( x ) dx (130)= − λ (cid:32) d − ∞ (cid:88) n =1 p clu ( n ) (cid:90) d Q ( x ) n dx (cid:33) (131)and the probability p = 1 − q that the sensor is efficient is therefore p = 1 − exp (cid:34) − λ (cid:32) d − ∞ (cid:88) n =1 p clu ( n ) (cid:90) d Q ( x ) n dx (cid:33)(cid:35) (132)
13. Appendix E Q ( x ) = 1 − P ( x ) = 11 + p − p exp [ − ( α − β ) x ] (133) (cid:90) d Q ( x ) n dx = 1 α − β (cid:20) H (cid:18) n, p exp[ − ( α − β ) d ]1 − p (cid:19) − H (cid:18) n, p − p (cid:19)(cid:21) (134) H ( n, y ) = (cid:90) dyy n (1 − y ) = ln y − y − n − (cid:88) m =1 my m (135)
14. Appendix F
We assume the probability p clu ( n ) to be continuous in n (which we can imagine by expressing it as asum of delta functions centered at integer values of n ) and write the expression in Eq. 84 as p ( n, ∆ x ) = (cid:18) − ∆ xλ (cid:19) δ ( n ) + ∆ xλ p clu ( n ) (136)27e can then replace the sums in Eq. 85 by integrals and have t ( N ) = − γv ∗ (cid:90) dm (cid:90) dm ... (cid:90) dm N p ( m , ∆ x ) p ( m , ∆ x ) ...p ( m N , ∆ x ) ln (cid:34) N (cid:88) n =0 m n h ( n ∆ x ) (cid:35) (137)We change variables according to m = 1 h N (cid:88) n =0 m n h n → m = m − h N (cid:88) n =1 m n h n (138)where we have written h n = h ( n ∆ x ), which gives t ( N ) = − γv ∗ (cid:90) dm (cid:34)(cid:90) dm ... (cid:90) dm N p (cid:32) m − h N (cid:88) n =1 m n h n , ∆ x (cid:33) p ( m , ∆ x ) ...p ( m N , ∆ x ) (cid:35) ln( h m )= − γv ∗ (cid:90) w ( m ) ln( h m ) dm (139)with w ( m ) = (cid:90) dm (cid:90) dm ... (cid:90) dm N p (cid:32) m − h N (cid:88) n =1 m n h n , ∆ x (cid:33) p ( m , ∆ x ) ...p ( m N , ∆ x ) (140)The Laplace transform of this expression is W ( s ) = P ( s, ∆ x ) P (cid:18) h h s, ∆ x (cid:19) P (cid:18) h h s, ∆ x (cid:19) ...P (cid:18) h N h s, ∆ x (cid:19) = exp (cid:34) N (cid:88) n =0 ln P (cid:18) h n h s, ∆ x (cid:19)(cid:35) (141)With P ( s, ∆ x ) being the Laplace transform of Eq. 136 P ( s, ∆ x ) = 1 + ∆ xλ ( P clu ( s ) −
1) (142)we have W ( s ) = exp (cid:32) N (cid:88) n =1 ln (cid:20) xλ ( P clu (cid:18) s h n h (cid:19) − (cid:21)(cid:33) (143) ≈ exp (cid:32) N (cid:88) n =1 ∆ xλ (cid:20) P clu (cid:18) s h n h (cid:19) − (cid:21)(cid:33) (144)= e − d/λ exp (cid:34) λ (cid:90) d P clu (cid:18) s h ( x ) h (0) (cid:19) dx (cid:35) (145)Normalizing to the probability that there is at least one interaction inside the gain layer we have W ( s ) = exp (cid:104) dλ d (cid:82) d P clu (cid:16) s h ( x ) h (0) (cid:17) dx (cid:105) − e d/λ − w ( y ) = L − [ W ( s )] t = − γv ∗ (cid:90) w ( y ) ln[ h (0) y ] dy t = 1( γv ∗ ) (cid:90) w ( y ) ln [ h (0) y ] dy (147)Comparing the expression to Eq. 71 we see that the effect of the position dependence h ( x ) is equivalentto a change of the cluster size distribution according to P clu ( s ) = 1 d (cid:90) d P clu (cid:18) s h ( x ) h (0) (cid:19) dx (148)28
5. BibliographyReferences5. BibliographyReferences