Decay Time Estimates by a Continuum Model for Inorganic Scintillators

We use the phenomenological continuum model for inorganic scintillators proposed by the author to give decay time estimates for four scintillators previously studied, namely NaI:Tl, CaF2, Gd2SiO5:Ce (GSO:Ce), and LaCl3:Ce. We show that, in order to obtain a good estimate of the decay time, we need to know (besides other well-known parameters) either the excitation carriers’ mobility or the structure and the parameters of the recombination mechanism. For these four materials, we know the data for the recombination term, whereas we have very scarce information about mobilities. However, we show that also in absence of experimentally-measured mobilities, with reasonable assumptions about them, we can obtain a good estimate for the slow component of the decay time. We show also when it is appropriate to model scintillation with one of the two most-used phenomenological models, the kinetic and the diffusive. The main point of the present approach is that it requires a limited set of experimentally-measured data and can be hopefully used in conjunction with more sophisticated and detailed models to design faster inorganic scintillators.


Introduction
Phenomenological models are widely used for scintillators, starting from the first model proposed by [1,2] to more recent contributions like, (e.g., [3]), because of their capability to describe the main features of scintillation with a limited number of parameters.In [4], we showed how these models can be encompassed within the general theory of the continua with a microstructure endowed with suitable thermodynamics, as was done for semiconductors, (e.g., in [5]).The resulting evolution equation for the excitation carriers is a reaction-diffusion-drift equation, which was the object of great attention from the mathematical point of view, (see e.g., [6]).Among the many results associated with the study of this equation, a result about the asymptotic decay of the solution allows for an estimate of the scintillator decay time.Such an estimate has an explicit dependence on the equation parameters, namely the carriers' mobility, the recombination rate, the dielectric constant, and the scintillation volume, i.e., a reference volume within which the scintillation associated with incoming ionizing energy is contained.In this paper, we use these results to show how a decay time estimate can be obtained by starting from the material parameters measured for specific scintillators and how it compares with the experimentally-measured decay times.
The paper is organized as follows: in Section 2, in order to make the paper self-contained, we briefly describe the continuum model obtained in [4] (many mathematical details we omit here are fully illustrated in [4]).We arrive at a system of partial differential equations, which describe the non-isothermal excitation carriers' evolution within the mesoscopic volume, and we show the asymptotic decay estimate for the excitation carriers.In Section 3, we put these equations in a dimensionless form to show how, for the isothermal case, the kinetic and diffusive models can be obtained as special cases.Then, the decay time estimates are evaluated for the four materials studied in [7], provided we use the data from the kinetic model the authors measured.Moreover, we show which model, between the kinetic and the diffusive, is the most appropriate to describe scintillation in these materials in relation to the amount of ionization energy.
As far as the notation is concerned, we put all the relevant symbols used in the Appendix A.

Materials and Methods
2.1.The Continuum Model for Scintillation.From the Microscopic to the Mesoscopic Scale: The Vector of Excitation Carrier Densities Let B be a scintillating crystal, which is hit at a point x ∈ B by an incoming ionization energy E. According to [8], in the very first stage of scintillation, the N excitation carriers generated by E propagate along a cylindrical track of radius r and length l from x, the mean free-path; for E exc , the excitation energy, the relation between the excitation carriers' density across the mean free-path section: and energy decay along the track coordinate z ∈ [0 , l] is given by: Equation ( 2) can be used stepwise to evaluate n: first, we identify the derivative with respect to z of the energy with the stopping power S(E), which is expressed in terms of the Bethe-Bloch equation [9][10][11], in order to obtain the dependence of the length l on the initial energy: An exact solution of (3) is out of question because of the correction terms in the Bethe-Bloch equation and of the singularity for E = 0: accordingly, we follow [12], where the approximate solution: with a o , b k and ξ k depending on the parameters of the Bethe-Bloch equation was obtained.This relation can be inverted, by following the procedure presented in [12], to have: from which the energy decay along the track can be arrived at: and then by (2): The excitation carriers recombine both radiatively and non-radiatively in a volume Ω ⊂ B centered on x at a scale that is greater than the mean free-path scale: we call such a scale the mesoscopic scale.In order to pass from the microscopic mean free-path scale to the mesoscopic one, we introduce a magnifier ε ∈ (0 , 1] and write: in such a way that for ε = 1, we get the microscopic mean free-path, whereas for ε → 0, we are at the mesoscopic scale, where the mean free-path reduces to a point.If we put (r ε , l ε ) into (7) in place of (r , l), we notice that (7) blows-up into the limit for ε → 0; accordingly, we define a rescaled density such that: which represents at the mesoscopic scale the excitation carrier density at a given point x ∈ Ω.
The N excitation carriers represent different types of energy carriers, for instance electron, holes, and excitons (bounded electron-hole pairs), which can recombine with different mechanisms like self-trapping, trapping by activator centers, and many others: in order to describe all these carriers and their different interaction mechanisms, we introduce the "vector" (indeed, a 1 × matrix): which shall be the basic state variable for our phenomenological model at the mesoscopic scale.
For convenience in the sequel, we shall also make use of the × matrix N(n) ≡ diag[n 1 , n 2 , . . ., n ] such that n = tr N; we define also the diagonal matrix N(n) −1 whose entries are n −1 j for n j = 0 and zero for n j = 0.
For e, the elementary charge, the excitation carrier densities generate an electric potential ϕ = ϕ(x , t): where is the permittivity and the components of the charge vector q ≡ [q 1 , q 2 , . . ., q ] take values q j ∈ {0 , ±1 , ±2 , . ..}.We assume that the volume Ω can be chosen in such a way that (11) admits the Neumann-type boundary condition, i.e., ∂ m ϕ = 0, with m the outward unit normal to ∂Ω.Moreover, any time-varying excitation density n h (x , t) generates a current density  h , which obeys the continuity equation: div The idea behind our present approach is to write the left-hand side of ( 12), normalized with respect to the elementary charge e, in terms of the mechanics of the continua with the microstructure [13].
where the 3 × microstress matrix T represents the total current density, the interactive microforce k and volume microforce b "vectors" represent respectively the internal rate of change of charge density and the rate of change of externally-supplied charges, and the surface microforce s represents the external supply of electric current.
A further assumption we are going to make is that within the mesoscopic volume Ω, the classical thermodynamics holds, and an absolute temperature θ is well defined.We can define accordingly an internal energy u = u(ϕ , θ): which, by (11), can be rewritten as: Likewise, we can define an internal entropy s = s(n , θ) of the Boltzmann/Gibbs type: where λ > 0 is the specific heat, k B is the Boltzmann constant, and c is a -dimensional normalizing constant vector.The pair (s , b) expends mechanical power for a quantity d we call scintillation potential in full analogy with the electrochemical potential in semiconductors: For h, the heat flux, and r, the heat supply, we have the energy balance and entropy inequality: By (13), we obtain the local form of (18), and then, upon the introduction of the Gibbs free-energy ψ = u − θs, we arrive at the reduced dissipation inequality: We assume, as a constitutive assumption, that ψ = ψ(n , θ , d , ∇d), and by a standard argument, we have that, so that the Gibbs free-energy is consistent with (19) for all processes, it must obey: with the totally-dissipative terms: where κ > 0 is the thermal conductivity and S , H are two semi-positive definite × symmetric matrices, whereas the local form of the energy balance reduces to: From ( 13) and ( 21), we obtain: where without loss of generality, we set b = 0. Equation ( 23) is expressed in terms of the scintillation potential, and in order to arrive at an expression in terms of the excitation density carriers, then by (20) 2 and ( 16), we arrive at: ∇d = eq ⊗ ∇ϕ then, provided we define respectively the mobility and diffusivity × matrices: and with the position: from ( 23) 1 , we arrive at a reaction-diffusion-drift evolution equation for the excitation carriers' density vector n: To complete the set of evolution equations, from the reduced energy balance (22), we obtain the heat equation with a scintillation-dependent source: where once again, we set the external heat contribution r = 0.As far as the boundary conditions are concerned, once again, we assume the mesoscopic volume Ω such that Neumann-type conditions hold, namely ∂ m ϕ = 0, ∂ m θ = 0, and s = 0 into (23) with (24) 2 , together with the initial conditions: where n is given by ( 10) and ϕ o (x) is the solution of (11) for n = n o (x).
The coupled system of partial differential Equations ( 11), (27), and (28) describes the evolution of excitation carriers within Ω: it recovers and generalizes to the non-isothermal case the equations first proposed in [14].Such an evolutionary PDE system is the counterpart for scintillators of the similar one we have in semiconductor theory; from the mathematical point of view, we have results about the global existence of renormalized and weak solutions, as well as estimates of the solutions asymptotic behavior (see [15] and the references quoted therein).
As far as the solutions' asymptotic behavior is concerned, we say (n ∞ , ϕ ∞ , θ ∞ ) is a stationary solution of ( 11), ( 27) and ( 28), a solution with ṅ = 0 and θ = 0: from Equation ( 23), it is easy to see that for a stationary solution, it is d = 0, and accordingly, from (24) 1 , we obtain: where ϕ ∞ is the solution of (11) for n = n ∞ .Following [6], in [15], we showed that in the isothermal case with θ(x , t) = θ ∞ , we have: where the norm is defined in suitable function spaces and C 1 , C 2 depends explicitly on the equations' parameters.In particular, C −1 1 , which is an estimate for the decay time, reads: where L(Ω) is the Poincaré constant of Ω, K = K(0) , M is the greater eigenvalue of M, and Φ = qϕ ∞ .

Reaction-Diffusion-Drift Equations for Scintillators
The solutions of the evolution Equations ( 11), (27), and (28) depends on different contributions, which are mutually competitive: for instance, in (27), we have diffusion, drift, and recombination, whereas in (28), there is a heat source given by excitation carriers.All these contributions depend on the given material, and accordingly, in order to evaluate their contribution to the solution, we need to put the evolution equations into a dimensionless form.Accordingly, we choose a characteristic length L, a characteristic time T, and a reference temperature Θ and define the adimensional variables: and the adimensional fields of (ζ , τ): From ( 11), (27), and (28), by (32) and (33), we arrive at the adimensional evolution equations: which depend on the set of five-adimensional parameters {d , m , k , c , f }: with D * = D/ D , M * = M/ M , and K * = K/ K .We remark that k is the only parameter that depends on the energy E by means of n.
First of all, we notice that the electrostatic Equation (34) 2 is independent of the parameters (35); as far as the heat Equation (34) 3 is concerned, for the typical values of permittivity and conductivity and on the whole temperature range, so that f ≈ 1, we need the characteristic length and time at the femto-scale.Accordingly, at the scintillation scale, the heat Equation (34) 3 is independent of scintillation, and we can assume that we may arrive at an equilibrium temperature θ * = Θξ * ; the whole process is isothermal, the temperature appearing as a parameter in the scintillation evolution Equation (34) 2 .
The isothermal case (34) 1 , which contains three contributions due to diffusion, drift, and reaction, admits various regimes depending on if one or two parameters are negligible with respect to the others.The two most successful phenomenological models available in the literature can also be recovered from (34) 1 : • k sup{d , m}; in this case, provided: then from (34) 2 , we obtain the kinetic model, cf.[7, [16][17][18][19][20]: • d sup{m , k}, which explicitly reads: we obtain from (34) 2 the diffusive model; (see, e.g., [21,22]): It is important to remark that the diffusive model, however, is rarely used alone; rather, it is coupled with some reactive terms as in [23,24]; moreover, the equations proposed in [3,[25][26][27] can be viewed as particular cases of Equation (34) 2 .

Explicit Results for a Cubic Model of Recombination
In terms of the equation parameters, the decay time estimate (31) normalized with respect to the characteristic length L reads: with φ in place of ϕ in the definition of Φ.
In order to obtain some explicit results from (40), we assume for the recombination term K(n) the cubic form proposed in [7] for a model with = 2 where: represent respectively the electron-hole pairs and the exciton densities: The R * * , K * * , and γ * * terms in (42) represent respectively the radiative, non-radiative, and exchange recombination parameters, which account for linear, quadratic, and Auger (cubic) recombination mechanisms: they were either evaluated experimentally or estimated for four scintillators (NaI:Tl, BaF 2 , Gd 2 SiO 5 :Ce (GSO:Ce), and LaCl 3 :Ce) in [7] (see also [18,28]); in the sequel, we are going to give explicit estimates for the decay time of these four scintillators by means of the estimate (40).
First of all, we notice that, since both the carrier densities are electrically neutral (i.e., q 1 = q 2 = 0), then by (11) with Neumann boundary conditions, we have an uniform electric potential with Φ = 0: as we show in [4], then by the conservation of charge, we have ϕ ≡ 0, and hence, ϕ ∞ = 0 with n ∞ = 0.It follows that the drift term vanishes identically and (34) 1 reduces to a reaction-diffusion system parameterized on d and k.
Upon the further assumption that Ω is the unit ball with L(Ω) = 0.068, then (40) reduces to: In order to evaluate k from (35), we assume that n in (41) can be evaluated by means of Equation ( 14) of [7] with the coefficients f e−h and f exc given therein.
By looking at that the two terms in (43), it is important to notice that the first one shows the dependence of the decay time on excitation carriers' mobilities, on the material dielectric properties, and on the amount of volume of interestin scintillation, whereas the second accounts for the recombination rate.
When we use the model of [7], K depends only on the linear rate recombination parameters, which are the ones that were determined experimentally; as far as M is concerned, instead, the only experimental result concerns the electron mobility in CsI [29].That was pointed out by many authors (also in [7]), however, that for these materials, the electron mobility is lower than 10 cm 2 /s, whereas that of the holes is an order of magnitude smaller; nothing is known for the exciton mobility.Accordingly, we choose, in order to apply (43), to assume negligible the hole mobility (cf., e.g., [30]) and the exciton mobility (if any) and to give for the electron mobility (when it is unknown) a value between 1 and 10 cm 2 /Vsec, to comply with the experimentally-measured decay times (in fact, one can use (43) to extract the mobility parameter from measured decay times).
We have d = 0.0202, whereas in the range of 1-10 4 KeV, the parameter k decreases from 0.6477-0.0061:for E = 1 KeV, we have d/k 1, and the evolution of excitation carriers is well described by the kinetic model alone.For increasing energies, the diffusive term becomes more and more important, and between 10 and 100 KeV, the kinetic and diffusive terms are comparable with the diffusive contribution, which becomes prevalent at E = 10 4 KeV (k/d = 0.33).
BaF 2 has a slow (τ d, f = 620 ns) and a fast (τ d,s = 0.6 ns) decay time [32]: the greater value of (43) then is a good estimate of the slow decay time.We suggest that the smaller value of (43) could be a lower bound for the fast decay time, a hypotheses that deserves further investigation, as we shall comment on in Section 4.
As far as the adimensional parameters are concerned, we have d/k 1 in the whole energy range 1-10 4 KeV: then, due to the great recombination rate, the scintillation can completely be described by the kinetic model.

GSO:Ce
The dielectric constant is = 3.4 [18]: there are no data concerning the mobility, and once more, we assume M = 1.5 cm 2 /Vsec.From [7], we obtain K = 3.3 • 10 7 s −1 and then τ d ≤ max{669 ; 16} ns.Since GSO:Ce has a slow decay time τ d, f = 600 ns and a fast one τ d,s = 56 ns [33], then we have again a good estimate for the slow value and a possible lower bound for the fast one.
In the usual range, d/k 1 and the kinetic model describe the scintillation well, the diffusive term being non-negligible at the higher energies (d/k = 0.11 for E = 10 4 KeV).

LaCl 3 :Ce
This scintillator has, according to the data provided by [7], the same recombination rate as GSO:Ce.The dielectric constant, from [17], is = 3.64.With M = 5 cm 2 /Vsec, the estimated decay time is τ d ≤ max{215 ; 16} ns against a slow decay τ d, f = 213 ns and a fast one in the range of 20 ns [34].
In the same energy range of the other scintillators, the parameter's behavior is purely kinetic at low energies, the diffusive part playing an increasing role for increasing energies (d/k = 0.37 for E = 10 4 KeV).

Discussion
The continuum model proposed in [4] generalizes the various phenomenological models widely used to characterize scintillators, which can be obtained as special cases, like the kinetic and diffusive models.Moreover, with such a model, we can obtain qualitative results on the evolution of excitation carriers.One of these qualitative results allows for the estimate of the decay time by means of a limited number of parameters.The most important of these parameters are the excitation carriers' mobility, the dielectric constant, the rate of recombination, and the scintillation volume Ω.
By using the recombination model and the data given in [7], we estimated the decay time for four well-studied scintillators (NaI:Tl, BaF 2 , GSO:Ce, and LaCL 3 :Ce) and obtained a good upper bound for the slow decay time.There are however some points to ponder.First of all, the relation (40) tells that the decay time depends either on the carriers' mobility (together with the dielectric constant and scintillation volume) or on the recombination rate.To use such a relation in a reliable way, we need accordingly to have both a good model for recombination (with experimentally-measured parameters) and detailed knowledge of the carriers' mobilities.Whereas in the first case, we have the experimental data collected in [7], for the present one, we lack full knowledge of the mobilities.It would be necessary to have, in order to use the decay time estimate at its fullest, the same degree of information we have for instance for semiconductors.
Further, the estimate (40) gives two values: in the case of the four scintillators analyzed, the greater value (the mobility-driven one) yields a good estimate for the slow component of the decay, also in the case of NaI:Tl, the only one among them for which we have a reliable estimate of the electron mobility.It would be a suggestive idea to assume that the smaller value (which in our case depends on the recombination rate) is a lower bound for the fast decay time.We have not yet enough physical insight into the phenomena to justify such a, rather suggestive, assumption, which will deserve further analysis both from the mathematical and the physical point of view.These results are summarized in Table 1.
Another result that follows from the model is that we can write dimensionless parameters d, m, and k, which account for the importance of the carriers' evolution of the diffusion, drift, and recombination.In particular, since in the model proposed in [7], the drift is unimportant, we show, for any given material and in a range of energy comprised between 1 KeV and 10 4 KeV, if the carrier evolution is of the kinetic, diffusive, or mixed type.For all four scintillators, the behavior at low energy is well described by the kinetic model alone.In the case of BaF 2 , the ratio d/k 1 in the whole energy range and the scintillation is of the kinetic type, as pointed out in [7], where it is stated that "diffusion plays no significant role in the kinetics"; on the other hand, in NaI:Tl, the diffusion plays an increasingly important role for increasing energies, becoming dominant on the recombination mechanism at high energy.The other two scintillators are, as again pointed out in [7], "a midway between BaF 2 and NaI:Tl": indeed, in GSO, the diffusive term plays a minor role only at high energies, whereas in La Cl 3 , again at high energies, the behavior is closer to the NaI:Tl one.These results are summarized in Table 2.
Table 1.Material data, measured decay times, and decay time estimates by Formula (43).The boldface entries are measured quantities; those in italics are estimated; the other being the results we obtained.In the first row, M is expressed in cm 2 /V s, K in s −1 , and the decay times in ns.The references for the experimental data are given in the text.

Table 2 .
Energy-dependent regimes for isothermal evolution Equation (34) 1 with no drift.Here, kstands for purely kinetic and d for purely diffusive; the entry "d + 0.α k" means that the evolution is mainly diffusive with an α% of kinetic behavior, whereas "k + 0.β d" means that the evolution is mainly kinetic with a β% of diffusive behavior.This research was funded by the UNIVPM Progetto Strategico di Ateneo 2017: Scintillating crystals: an interdisciplinary, applications-oriented approach aimed to the scientific knowledge and process control for application concerning life-quality improvement, Grant Number I31I18000210005.