A brief overview of existence results and decay time estimates for a mathematical modeling of scintillating crystals

Inorganic scintillating crystals can be modelled as continua with microstructure. For rigid and isothermal crystals the evolution of charge carriers becomes in this way described by a reaction-diffusion-drift equation coupled with the Poisson equation of electrostatic. Here we give a survey of the available existence and asymptotic decays results for the resulting boundary value problem, the latter being a direct estimate of the scintillation decay time. We also show how to recover various approximated models which encompass also the two most used phenomenological models for scintillators, namely the Kinetic and Diffusive ones. Also for these cases we show, whenever it is possible, which existence and asymptotic decays estimate results are known to date.


Introduction
A scintillator crystal is a material which converts ionizing radiations into photons in the frequency range of visible light, hence its name. It acts as a true "wavelength shifter" and in such a role is used as radiation sensor into high-energy physics, in medical imaging and in security applications [1]. The physics of scintillation, which is a complex multi-scale phenomenon (see e.g. [2]) can be described within a continuum approach at three scales: at a Microscopic scale the incoming energy E generates a population of charged energy carriers which moves in straight directions for few nanometer [3] and whose density N = N(E) can be found by the means of approximated solutions of the Bethe-Bloch equation [4]- [7]. 2 These energy carriers wander and migrate within a greater region either generating other energy carriers or recombining with emission of photons hν. In the process some energy is lost and a scintillator is a material in which such a loss reduces the frequency of the incoming ionizing energy to that of visible light. We call this scale the Mesoscopic scale: for P the region occupied by the crystal we denote Ω ⊂ P the mesoscopic volume in which the recombination of charge carriers into photons takes place (Fig. 1). Finally the light rays propagate within the crystal at a Macroscopic scale according to the laws of classical optics. Scintillation is a fast and dissipative phenomenon and there are two major physical parameters which are to be improved in a scintillating material: 2 In [8] we show how the density of excitation carriers N induced by an ionizing energy E which hits the crystal at a given point x * , can be obtained by the means of a suitable rescaling to the mesoscopic scale of the approximate solutions of the Bethe-Bloch equation along an elementary cylindrical track: in such a way N = N (E) maintains informations on both the initial energy E and the material properties of the crystal.
(i) the decay time τ d , which is the time required for scintillation emission to decrease to e −1 of its maximum and is a measure of the scintillator resolution; (ii) the light yield LY , which is the ratio between the collected light energy and the energy of the incoming ionizing radiation and which is a measure of the scintillator efficiency. Scintillators can be modelled as Continua with microstructure [8]- [11] to arrive at a Reaction-Diffusion-Drift (RDD) equation for the energy carriers descriptors, coupled with the Poisson equation of electrostatic and both with Neumann-type boundary conditions. Further in [12], by following the results obtained in [13] and [14], we showed that for these equations it is possible to proof the global existence of renormalized and weak solutions and also how the decay time can be estimated explicitly in terms of these equations constitutive parameters. Here, by using the results obtained into [15], we expand the results of [12] to show the existence of weak-strong renormalized solutions: we further show also (whenever it is possible) existence and asymptotic decay results for the phenomenological models widely used in the literature, which can be obtained from our model by introducing suitable approximations as we did into [9]. We remark that scintillation is strongly affected by temperature and the crystal are linearly elastic deformable bodies, but here we limit our analysis to the isothermal case (the temperature being at most a parameter in the constitutive quantities) and to rigid crystals as in [9]. The effects of temperature are dealt with into [8] and [10] whereas an insight into the deformable case is provided in [11]. A complete treatment of these electromagnetical, thermal and mechanical interactions will be provided in a forthcoming book [16]. The paper is organized as follows: in §.2 we give an overview of the model proposed into [8]- [11], which leads to a RDD system: then we study the properties of the stationary solutions, which are important when we deal with the solutions asymptotic properties. In §. 3 We extend with the help of [15] results presented into [12] and we give results on existence of weak-strong renormalized solutions and on asymptotic decay, this last result being a direct estimate of the decay time. This result allows, for the first time, to estimate the decay time in terms of the parameters describing the physical properties of a given scintillating crystal and which appears explicitly in the equations we obtained. Finally in §.4 first of all we put the coupled RDD and Poisson equations in an adimensional form by introducing characteristic length and time related to the scale of scintillation, as we did into [9]. The resulting boundary value problem depends on a set of three adimensional parameter related to diffusion, drift and recombination respectively. Then by choosing a suitable measure of smallness we obtain, to within higher-order terms in such a measure, different approximated models which represents different physically meaningful regimes for scintillators. In such a way we not only recover the most used phenomenological model in use, namely the Kinetic and Diffusive ones, but also show under which hypotheses scintillation could be described either by a Reaction-Diffusion or a Diffusion-Drift equation. Also for these approximated boundary value problems we give a survey of the existence and asymptotic decay results, either by taking them straight from the available literature or by adapting them to the specific context of scintillation as described by our model.
2. The continuum model for isothermal and rigid scintillators: the boundary value problem

A reaction-diffusion-drift equation for scintillators
The charge carriers with density N generated at the point x * and time t * and whose dimension is (length) −3 , represent a population of carriers which can differ by their sign (e.g. negative electrons, positive holes, neutral excitons and so on) and by their recombination mechanism. Accordingly we may differentiate these carriers by introducing an ordered array n, the charge carrier densities vector: clearly, the greater is k, the finer will be our description of the phenomena: the simplest non-trivial choice is k = 2 as in [17], where n 1 represents the electrons population (which is equal to the holes population) and n 2 represents the population of excitons, which are bounded electron-hole pairs evolving together. In the various phenomenological models for scintillation thus far proposed we may have m = 3 as in [18] and [19], m = 7 as in [20] whereas in [21] we have m ≥ 11. The charge carrier densities are related to N by: the set {α 1 , α 2 , . . . , α k } depends on the specific scintillator, on the initial energy E and on the initial data (cf. the discussion in [22]). We identify the charge carriers densities at (x * , t * ) with fields defined on the whole Ω × [0 , τ ) and whose regularity will shall made precise when needed: without loss of generality we may assume that the mesoscopic volume Ω is a ball of a still unprescribed radius R, centered at x * . We assume accordingly the k−dimensional field, the excitation carrier densities vector as the main state variable in our description of scintillation. In view of (1), since n j ≥ 0, j = 1 , 2 , . . . , k, we can view the extension of N to a field on Ω as the L 1 (Ω) norm of n: Let e be the elementary charge, then the excitation carrier vector induces a free charge density within the scintillation volume Ω with z = (z 1 , z 2 , . . . z k ), z j ∈ Z, j = 1 , 2 , . . . k the charge vector. By the Maxwell-Lorentz equations in absence of magnetic fields [23] these free charges induce a local electric potential (x , t) → ϕ(x , t) which, for a given time t ∈ [0 , τ ), is the solution of the Poisson equation of electrostatic with associated Neumann boundary conditions: here ǫ is the permittivity of the crystal (which at this stage we assume isotropic or at most cubic), 3 and m is the outward unit normal to ∂Ω. We notice that in (7) we do not take into account either bound charges or external charges since we are mainly interested into X− or γ−rays which have zero charge, whereas for α− and β−rays, which have respectively positive and negative charges, an external charge contribution q * should be added to (7). We further remark that an implicit way to select the radius of Ω s indeed that at its boundary the Neumann condition holds, i.e. there is no electric field outflow trough ∂Ω. The boundary-value problem (7) admits an unique (up to a constant) weak solution ϕ ∈ H 1 (Ω) provided the total charge is conserved: annd the constant can be conveniently determined if we setφ = 0, 4 thus making the solution unique. Following [24], we assume that the electrostatic free-energy associated to the electric potential has, besides the classical conservative term which depends on n by the means of (7), a dissipative term F (n) of entropic nature which also depends on the excitation carrier vector: where θ > 0 is the (fixed) absolute temperature; in [10] we showed how to 3 For anisotropic crystals equations (7) become: where ǫ o is the vacuum permittivity and K is the symmetric and positive definite permittivity tensor ; relations (11) changes accordingly with ǫ o K[∇ϕ] · ∇ϕ in place of ǫ ∇ϕ 2 4 For any given integrable function f we shall denote withf its mean value: such an energy it can be associated an electrostatic self-power: where the elements of the array s ≡ (s 1 , s 2 , . . . s k ) represent the scintillation potentials associated to the various charge carriers: We shall call s, with an abuse of terminology, the scintillation potentials vector : and it is easy to show that: where D denotes the Frechet derivative.
In [8] and the related papers [9], [10] we showed how, by using a continuum with microstructure approach and the classical dissipation inequality, the evolution equation for the charge carriers in scintillators can be represented as: where S(n) and H(n) are two definite-positive k × k matrices and The boundary-value problem (16) has an equivalent variational formulation which leads to a gradient-flow type problem (vid. e.g. [25] and for more recent results [26]): where the conjugate dissipation functional Ψ(n , s) is related to the thermodynamical dissipation D by The gradient-flow classical structure is and U is the driving functional: starting with [27], provided there exists K(n) = G −1 (n), the following formulation was proposed and successfully used into e.g.
[28]- [35]: where the Onsager structure K ca be splitted additively into different contributions. For instance, in our case we may set: Furthermore, in [8] it is show that, provided we identify the entropic term with the Gibbs entropy, with c j , j = 1 , 2 , . . . , k, normalizing constants and k B the Boltzmann constant then, by using (23) into (14) and (17), we arrive at a Reaction-Diffusion-Drift equation for the evolution of charge carriers in scintillators with Neumann boundary conditions and initial data: and r(n) is the k−dimensional recombination array.
The relation between the matrix S in (16) and the matrix M in (24) where N −1 (n) denotes the diagonal k × k matrix whose entries are n −1 j when n j = 0 and 0 when n j = 0. The semi-definite positiveness of M accounts for excitation carriers whose mobilities (and hence by (25) the associated diffusivities) are either zero or negligible with respect to those of other charge carriers. However in §. 3 we shall see that in order to get decay time estimates we shall require the stronger requirement of positive-definiteness for M. Finally, the k × k matrix H(n) and the k−dimensional array r(n) are related by: where, by following [27], we assume that In (28) h = 1 , 2 , . . . s denote the number of recombination processes with rates k h with the two k−dimensional arrays describing the h th recombination mechanism with rate k h : and the function ℓ(x , y) is the logarithmic mean: Whenever we assume for F (n) the Gibbs entropy (23), then from (27) and (28) we arrive at a polynomial expression for the recombination term r(n) 5 We assume that the mobility is independent on n, as pointed out into [38]: in such a case S(n) must be restricted to the form where S o is a k × k matrix with constant components.
(vid. e.g. [28] for the details): where we used the identity log v a = a · log v. In order to arrive to (31) we implicitly assumed that all the recombination mechanisms are detailed balanced [27], that is there exists a steady recombination state which in our case coincides with c; further, let and then from (31) we have If we consider now the differential form of the charge conservation (10), namely: d dt then from (24) we obtain which in turn, by (31), implies the electrical neutrality of each recombination mechanism: and hence z ∈ S ⊥ in such a way that the local orthogonality condition holds: In [27] it is remarked that the term r(n) can be represented as a polynomial relation in n only if we identify F (n) with the Gibbs entropy (23), whereas for a different entropic term like e.g. the Fermi-Dirac potentials this is not possible. Equation (24) 1 , which represents the conservation of electric current normalized with respect to e, was first proposed for scintillators in [21] by following [36] and [37] and with r(n) a polynomial, at most cubic, function of n; it was used for scintillators into [17], [18], [38]- [41]. Moreover, special cases of this equations were widely used to model scintillation at a phenomenological level in many theoretical and experimental paper, as we shall describe in details in §.4. We also remark that the boundary value problem (24) is the same obtained, by starting from a different approach and with a different reaction term r(n), in [24] for semiconductors (vid. also [27], [29]). In the available phenomenological models for scintillation the recombination term r(n) is generally assumed as a cubic expression in n: (39) The terms r o i describes the excitation carriers creation rate in the crystal under irradiation; the other terms in (39) are further splitted in order to represent different recombination mechanism. The terms A ij , which accounts for linear recombination, can be further decomposed into three terms which represent respectively radiative recombination with photon emission, non-radiative recombination without photon emission and exchange between excitation carriers. The quadratic recombination is represented by the terms B ijh which can also be splitted into B ijh = B r ijh + B nr ijh with the same meaning of superscript as in (40), whereas the third-order or Auger totally nonradiative recombination is described by the terms C ijhm = C nr ijhm . In order to reconcile (39) with (31) we consider here the following example for k = 3 where n 1 = n e and n 2 = n h represent respectively the electrons and holes densities and n 3 = n ex represents the excitons density. By following, for instance, the description provided in §.3.3 of [42] then we may have these recombination mechanisms: and then from (31) we get: which can be trivially put into the form (39), with the appropriate identification of radiative, non-radiative or exchange terms: we notice that the mechanisms 3, 4 and 6 are Auger recombination, whereas 2, 5 and 7 represents scattering and the mechanism 1 is the simple electron-hole recombination (cf. e.g. the models proposed either into [22] or [43] where n e = n h = n eh ).
As we already remarked more complex expressions for (39) can be proposed: for instance into [20], for k = 7, we had: r 1 (n) = r 0 1 + A 14 n 4 + B 113 n 1 n 3 + B 115 n 1 n 5 + B 117 n 1 n 7 , r 2 (n) = r 0 2 + A 23 n 3 + A 25 n 5 + B 223 n 2 n 3 + B 224 n 2 n 4 + B 227 n 2 n 7 , r 3 (n) = A 33 n 3 + B 323 n 2 n 3 + B 336 n 3 n 6 , r 4 (n) = A 44 n 4 + B 417 n 1 n 7 + B 424 n 2 n 4 , (42) r 5 (n) = A 55 n 5 + B 515 n 1 n 5 + B 527 n 2 n 7 , r 6 (n) = r 0 6 + A 66 n 6 + B 613 n 1 n 3 , r 7 (n) = A 77 n 7 + B 715 n 1 n 5 + B 724 n 2 n 4 , with n 1 and n 2 the electron and holes densities, n 3 and n 6 the self-trapped electrons and excitons respectively, n 4 and n 5 the electron and holes captured by the activation centers and finally n 7 denotes the excited activator centers: in [20] a detailed description of the physical motivation of these relations is provided. We notice that in this model the Auger mechanism is missing and the recombination terms are at most quadratic.

Stationary solutions
It is trivial to show that the stationary solutions, that is equilibrium solutions n ∞ for the boundary value problem (24) withṅ = 0, can be obtained by setting s = 0, which by (27) leads to the equilibrium condition: When F (n) is identified with the Gibbs entropy (23), then from s = 0 we have where the stationary electric field ϕ ∞ is the solution of: with Neumann boundary conditions on ∂Ω. From (44) then we have and (45), (46) together leads to the semilinear Poisson-Boltzmann equation for the stationary electric field: into [44] an uniqueness result in H 1 (Ω) for (47) with Dirichlet boundary conditions was obtained. It is trivial to show that the equilibrium solution n ∞ is also a steady state c for the detailed balance condition and indeed in (27) we have: by (46) and (37), and hence the equilibrium condition (43) is satisfied. We finally deal with the problem of the determination of c = (c 1 , c 2 , . . . , c k ) in (46): if we define the k × k diagonal matrix L(x) then (46) can be written as: and then, by (10) L(x)c ∈ S.
To obtain the explicit value of c we may use (43), whereas the uniqueness of c follows instead, as pointed out in [14], by monotonicity, (10) and the uniqueness of ϕ ∞ ∈ H 1 (Ω) withφ ∞ = 0.

Existence and asymptotic decay
In this section we shall give an account of the existing results concerning the existence of solutions and the asymptotic estimates for the boundary value problem (24): we remark that the latter results are important in order to get a scintillation decay time estimate. We shall not enter into the mathematical details which can be found in the references we quote, rather we shall adapt if necessary these results to the specific cases of our boundary value problems. The problem of finding existence, asymptotic estimates and qualitative bounds for the solutions for the coupled boundary value problem (24) has received a strong attention in the recent years, vid. e.g. [13], [15], [81], [45]- [54] and the many references quoted therein: to this regard it is important to remark that most of these results deal with semiconductors or chemical reactions which differ from scintillators by the structure of the reaction term r(n). Most of these results are based on the so-called Entropy method, whose importance is explained in full in these words taken from [55]: The entropy method refers to the general idea of a functional inequality relationship between an entropy functional of a system and its monotone change in time, usually called the entropy dissipation. Such an entropy-entropy dissipation inequality entails convergence to an entropy minimizing equilibrium state, at first in entropy and further in L 1 using Cziszár-Kullback-Pinsker-type inequalities. The entropy approach is per se a nonlinear method avoiding any kind of linearization and capable of providing explicitly computable convergence rates. Moreover, being based on functional inequalities rather than particular differential equations, it has the advantage of being quite robust with respect to model variations.
In the next subsection we shall show how some of these results can be extended to the RDD equations for scintillators. As far as the decay time is concerned, the available experimental data (vid. e.g. the recent analysis in [56]) and the numerical solution of phenomenological models as in [57], show that the excitation carriers decay exponentially in time to an asymptotic value n ∞ , namely: where the indeces f and s denotes the so-called fast and slow components of the excitation, respectively. Accordingly, since by definition the Decay time is the time required for scintillation emission to decrease to e −1 of its maximum, then we get a Fast Decay Time τ f and a Slow Decay Time τ s . In many cases one of the components is negligible and the decay obeys a simple exponential law, which can be also used to describe an average decay time.

Global existence
The first results concerning existence theory for the boundary value problems like (24) was obtained in [50], [51]: such a result, which relies on the notion of global renormalised solutions leaves still open, as pointed out in [51], the problem the existence of weak or even smooth global solutions in time: the reason, as pointed out in detail into [51] are the growth condition on the reaction/recombination terms. These results were the extended into [13] and [15] and here we shall follow the latter. The main strongpoint of these results is that no growth condition are imposed on the reaction/recombination term. In our case the diffusion is linear: unfortunately we cannot extend their results to the non-linear case they study, because for (25) this would imply also a non linear mobility, whereas in all these paper the mobility is implicitly assumed M = I with I the k × k identity. First of all we recall the notion of renormalised solution, first introduced into [58]- [60] for the Boltzmann and transport equations, as it was given into [50]. An excitation density vector n is a renormalised solutions for (24) if for all functions ξ : M → R with compactly supported derivative ∇ n ξ , the function ξ(n) must satisfy the equation derived from (24) by a formal application of the chain rule in a weak sense. As it is pointed out into [61], the function ξ must belongs to a well-choosen class of admissible solutions. The physical interpretation of these renormalised solution is that they gives a distributional sense to the boundary value problem (71).

Asymptotic estimate for the decay time
In [14] an explicit estimate of the asymptotic convergence was obtained for the Rosbroeck model for semiconductors with Shockley-Read-Hall potential and k = 2: here we shall show how the results obtained there can be adapted to the case of scintillators in order to obtain an explicit estimate for the decay time. Once again we shall leave out the technicalities and details and we shall give only the main results using our language and notations.
We have already seen that we can safely assume uniqueness for c: given this in [14] some preliminary bounds are necessary, namely the following for c = (c 1 , c 2 ) and n ∞ = (n ∞ 1 , n ∞ 2 ): Moreover, for the recombination r(n) given by (39) there exists a constant k o such that: If we consider the free-energy U(n) defined by (11) with the choice of the Gibbs entropy (23) for the term F (n), then the main results of [14] are based on the derived notion of Relative Gibbs free-energy (or relative entropy according to [14]): by an explicit calculation it can be show that and then, by an easy calculation, it can be shown that the Dissipation D defined by (19) is given by In [14] by starting from (63) and by the means of a repeated use of Csiszár-Kullback-Pinsker type inequalities, provided (58) hold, the two following estimates were obtained: where ϕ o is the unique solution of (7) for the initial data n o . The most remarkable feature of this results is that both C 1 and C 2 have an explicit dependence on the parameters of (71): in (65) L(Ω) > 0 is the constant in the Poincaré inequality: which in L 2 (Ω) reduces to L(Ω) = λ −1 1 with λ 1 the first eigenvalue of: We notice that, for Ω a sphere of radius R it is L(Ω) ≤ (2R/π) 2 [62]: for L = 2R then we have L(Ω) ≤ 0.1. Such a value is consistent with the value L(Ω) ≈ 0.07 which was obtained into [9] where the λ 1 was calculated as the the square root of the reciprocal of the first zero of the Bessel function J ′ o . The expression for the decay time τ = C −1 1 depends, by (65) 1 , in an explicit manner on the mobility parameter m, the reaction term r(0), the initial data n o , and the scintillation volume Ω: as far as we know it is the first explicit estimate of scintillator decay time which depend on the (measurable) constitutive parameters of the model, albeit limited to the case k = 2. In [9] we showed how, for four different scintillators these results gave a very good estimate for the experimentally measured fast decay time.

The adimensionalised boundary value problem
In order to put the boundary value problem (24) in an adimensionalised form, we begin with the choice of a characteristic length and time pair (L , T ) related to the space and time scintillation scales, in order to define the dimensionless coordinates (ζ , τ ): and the dimensionless excitation carriers density u = u(ζ , τ ) and electric potential ψ = ψ(ζ , τ ): If we further set: where µ is the greatest eigenvalues of M and and κ the greater component of r or, as in [9], the value of their norm, then equations (24) 1,3 can be rendered dimensionless and dependent on the three adimensional parameters: We remark that the parameter k depends on the incoming energy by the means of r(n), whereas m depends only on constituive or scaling parameters: accordingly, depending on the energy of ionizing radiation we may have different physically meaningful regimes which are described in the following subsections. Further, for the greatest eigenvalue of D and provided we define the diffusion length L D = √ δ T , then we have we notice that an experimental estimate of L D is given e.g. into [22]. The set of dimensionless parameters {d , m , k} which appears into (71) 1 describes three main regimes for the boundary value problem, depending on the crystal constitutive properties and on the initial ionizing energy. In the ext subsections we shall show how many of the most used phenomenological models for scintillators, can be encompassed within (71) 1 by an appropriate choice of these parameters. For these models, which we define approximated because some terms of (71) 1 can be neglected, we shall also briefly describe the results concerning existence and decay time which are available in the literature.

The Reaction-Diffusion approximation
and let then from (71) 1 we recover, to within higher order terms, the reaction- and (71) 2 is uncoupled, the excitation density vector being a data; this phenomenological model, which traces back its ancestry to the analogous models for chemical reactions [63], [64], is used in many papers dealing with experimental identification of scintillator properties, as in [65]- [71] and did not take into account the drift contribution induced by the electric field ϕ.
Reaction-Diffusion equations were studied with a great interests starting from the kinetics of chemical reactions [64] and there are treatises and textbook devoted to various aspects of them like, e.g., [72]- [74] and many others. However, problems like the existence of solutions and their asymptotic decay have attracted a growing number of studies in recent years, mainly because the difficulties posed by the the lack of control of the recombination terms as pointed out in details into [49]- [51]. Most of the recent results dealing with asymptotic decay deal indeed with reaction terms with quadratic growth for chemical reactions vid.
[75]- [79] and also, for a different point of view, [27], [80]- [81]. A first result concerning the existence of renormalised solutions was presented into [82], whereas the most general result at the present available is [83], where however the nonlinear case of diffusion in porous media is treated. See also the references into [34], where these results were extended to the non-isothermal case.
Here we show how these results can be used in the context of our problem, described by the equation (77). There are, at the best of our knowledge, no results for a general recombination term: however, if we neglect the cubic Auger effect, we can use the existence results of [82] and decay estimates obtained into [79].
We shall give only a brief survey of both these results and how their hypotheses fit within the physics underlying scintillation: we leave out all the technical details contained in the cited works. We remark that a complete existence and decay estimate results for the scintillation models has yet to be done and the task, given the nature of the reactive term, is a far from an easy one. A further result presented into [79] is concerned with the existence of renormalized solutions for (77), whose structure is the same as in (52) when we neglect the drift terms.

Asymptotic decay estimates
The result obtained into [79] is tailored on chemical reactions, which implies mass conservation and detailed or complex balance of reactive terms, conditions which have no correspondence in the physics of scintillation. In this paper, by using arguments and tools which are unsurprisingly related to those used into [14], it is shown that, for n ∞ a detailed balanced equilibrium solution corresponding to r(n ∞ ) = 0, then: where the relative entropy G is defined as (cf. (62)): C 1 = C −1 CKP is the constant in a Csiszár-Kullback-Pinsker type inequality, whereas the constant C 2 , whose inverse is the estimate for the decay time, is given by: in this relation λ 1 = C LSI min{D i }, with D = diag{D 1 , D 2 , . . . D k } and C LSI is the constant in the logarithmic Sobolev inequality, whereas the two constants K 1,2 depends explicitly on Ω, D, the set k h , a h , b h , n ∞ and on the constant K such that all the details being given in full into [79]. Finally, the term H 1 (t) in (82) is given by:

The Diffusion-Drift approximation
and then we recover, to within higher order terms, the Poisson-Nernst-Planck system, used to model ion fluxes, cell biology and other electrically driven evolution phenomena: with Neumann boundary conditions. This model was used into [17] and [39] to describe the initial stage of the scintillation, before the recombination takes place. The system (87) is well-studied and there are many results, dealing with existence, asymptotic decay, equilibrium solutions and even explicit analytical solutions for k = 2 (e.g. amongst the many [84]- [87] and for more recent advances [88]- [90]). Equation (87) 1 can be put in the equivalent gradient flow formulation (21) with K = K D . Its stationary solutions (n ∞ , ϕ ∞ ) are characterized by s ∞ =const. and therefore (46) is replaced by and the Poisson-Boltzmann equation (47) changes accordingly.

The Reaction-Drift approximation
The boundary value problem (71) reduces in this case to and (7). At the best of our knowledge this equation was never used for scintillators and indeed it can be found in the theory of dopant diffusion in semiconductors (vid. e.g. [91]): however we list it together with those used for scintillators for the sake of completeness. For the boundary value problem (90) and (7) globale existence, uniqueness and asymptotic decay results were obtained into [92] (vid. also [93] and [94]).

The Kinetic approximation
Let ε a small parameter, then for k = O(1) and then (71) 2 reduces, to within higher order terms in ε, to the rate ODE equation: with initial condition n(0) = n o ; we remark that by (92) and (38), then is a first integral of (92). Also in this case (71) 2 is uncoupled, the solution of (92) becoming a charge density supply. With such approximation, which describes phenomena mainly driven by recombination, we recover the so-called Kinetic model, the oldest and most used phenomenological model for scintillators [20], [22], [95]- [100], which is borrowed from the kinetic of chemical reactions, see e.g. [101], [102]. In a series of recent papers, [31]- [35], dealing with the gradient structure of (24) it is show that also (92) admits a gradient structure with K = K R for detailed balanced recombination mechanisms and hence, the energy dissipation methods can be used with success to get asymptotic decay estimates.
In detail, into [35] the well-posedness of (92) was obtained, in the sense that for all initial data n o ∈ M there exists an unique global solution n : [0 , τ ) → M. Moreover it was proved that since the recombination term satisfies the detailed balance conditions for the equilibrium state n ∞ = c, then (92) admits the gradient structure:ṅ = −HDU(n) , with H given by (28).
Finally, into [32] an estimate of the kind of (64) was provided: where the simplified dissipation D * , such that D ≥ k o D * where k o is given by e.g. (60), is defined as: 4.6. The Diffusive approximation For d = O(1) (which, by (74) means that L ≈ L D ), and then from (71) 1 we obtain a classical anisotropic parabolic equation: In this case, which describes the phenomena when the diffusion within the track is the driving mechanism of scintillation, we recover the Diffusive model, which describes the diffusion of excitation carrier within the track vid. [19], and is rarely used alone as in e.g. [38]. Once again the equation of electrostatic (71) 2 is uncoupled. We notice that the ratio k/d is known as the Thiele modulus in the kinetics of chemical reactions. A sharp estimate of the decay time for (98) can be easily found by following [55]: 7 first of all we notice that at the equilibriumṅ = 0 and, by the Neumann boundary conditions, ∇n = 0, hence n ∞ =n ∞ = const and by the charge conservation (10): Further, by multiplying both sides of (98) by n(x , t) − n ∞ , integrating by parts with the Neumann boundary conditions and by the Poincaré inequality then we get: The first and last terms of (100) yield a first order equation which can be integrated to obtain at an expression alike (64): with Accordingly the decay time τ d can be estimated as: Whenever D = diag{D 1 , D 2 , . . . , D k }, equation (98) reduces to the k independent classical diffusion equations: with Neumann boundary conditions, whose a complete mathematical treatment can be found into many books, like e.g. [103].

Conclusions
For the reaction-diffusion-drift equation which describes the evolution and recombination processes of charge carriers in scintillators we gave an overview of the existence and asymptotic decay estimate which are know to date, at the best of our knowledge. Despite the fact that the topics is a well-studied one, as the non-exhaustive list of references shows, there are still many unanswered questions which deserves further investigations: here we shall give a concise list of some them. Such a list is of course far from be exhaustive and its items are those which seem more interesting for the physicists. First of all the mathematical treatment of the RDD system is based on a precise choice of the entropic term F (n), namely the Gibbs entropy: however the physics of scintillation processes suggests that such a choice can be appropriate into describing the behavior over a limited range of energy. The choice of Fermi-Dirac potential in place of the Gibbs entropy will be not only more general but also more related to the true physical nature of the phenomena. Clearly, this leads to a different formulation for the recombination term which requires a different mathematical treatment.
A second point of interest is the decay time estimate obtained for the RDD system: here there are two major directions for further investigations. The first one is to extend the results obtained into [14] for a system of two charge carriers, k = 2 to systems with general k as most of the scintillator phenomenological studies require. The second and more intriguing aspect is the following: the decay time estimate which follows from the result of [14] gives as an upper bound of the decay time the maximum between two values, one which depends on the carrier mobility and the other which depends on the recombination time. When into [9] we applied these results to four scintillating crystal, not only we obtained a very good estimate ad consistent of the fast decay time: we also notice that the minumum between these two values looked like a lower bound for the slow decay time. Clearly the results of [14] tells nothing about this but it would be tempting to prove such and assertion. In such a case we should have two bounds, an upper one on the fast decay time and a lower one on the slow decay time, a result that would be appreciable in terms of material science. A third remark is that the Entropic methods used to get these estimate should be applied also to the approximated phenomenological models, in particular to the kinetic one: we need to remark that a very useful aspects of the results in [14] is the explicit dependence of the estimates on the constitutive parameters of the RDD systems. This should allows for a predictive use of these results, not only for an a-posteriori check with the available experimental data.
As far as the questions related to the existence issues for very general form of recombination term is concerned, the necessity of further generalization doesn't need to be justified, since they are the foundation stone of any numerical procedure we need to implement in order to obtain numerical solutions. Finally there are two side-aspects related to those we looked at in this paper and which we left out. The first one concerns the solution of the Bethe-Bloch equation, which is a necessary requirement to bridge the microscopic world with the phenomenological treatment we want to give to scintillation. To date there are many results which concern mainly the radiation decay length or the Bragg peak: a treatment which conveys in a straightforward way the relevant parameters of the phenomena at the mesoscopic scale would be welcomed. The second and in some sense more important problem we left out is the Ligth yield : such a parameter is indeed a measure of scintillator effectiveness and tells us about the minimum energy we can detect. To date there is not a coherent mathematical definition for this parameter, apart a descriptive one: LY = Number of charge carriers recombined into photons Total number of charge carriers generated , which in turns can be applied either locally or globally, whatever this means.
Once again, a formal definition and its consequent mathematical treatment are, in our opinion, still missing and it would benefit from the huge amount of experimental works which give a precise evaluation of Light Yield.