Complex diffusion-based kinetics of photoluminescence in semiconductor nanoplatelets

We present a diffusion-based simulation and theoretical models for explanation of photoluminescence (PL) emission intensity in semiconductor nanoplatelets. It is shown that the shape of PL intensity curves can be reproduced by the interplay of recombination, diffusion and trapping of excitons. The emission intensity at short times is purely exponential and is defined by recombination. At long times it is governed by the release of excitons from surface traps and is characterized by a power-law tail. We show that the crossover from one limit to another is controlled by diffusion properties. This intermediate region exhibits a rich behaviour depending on the value of diffusivity. Proposed approach reproduces all the features of experimental curves measured for different nanoplatelet systems.

The hunt for materials and systems with better optical properties has always been one of the focal points of semiconductor research [1]. The first classical optoelectronical applications were based mostly on bulk properties of semiconductors. Recent progress in material synthesis [2][3][4] led to wider research efforts concentrated on low-dimensional structures. In particular, significant attention is drawn by the semiconductor nanocrystals, such as quantum dots (QDs) and nanoplatelets (NPLs) [5,6]. QDs and NPLs show large exciton binding energy [7][8][9] and possess such remarkable properties as narrow emission lines even at room temperature, tunable emission wavelength, short radiative lifetimes, giant oscillator strength, high quantum yield and vanishing inhomogeneous broadening [10][11][12][13][14]. These properties make them excellent prospective building blocks for future optoelectronic devices such as bright and flexible light emitters [15], colloidal lasers [8,16] or for biomedical labeling applications [17].
Modern techniques of precise colloidal synthesis allow to control the shape [18], size [19] and crystal structure of the nanocrystals [20]. Among these a growing popularity belongs to colloidal QDs and NPLs of CdSe with atomically defined thickness [5,21,22] and various shapes such as pure QDs and NPLs (core) as well as composite coreshell QDs and NPLs or core-crown NPLs. The composites are made of two different semiconductor materials, for instance, CdSe and CdS.
The low-dimensionality of all these structures effects in a decisive contribution of a crystal surface to control of optical and electronic properties [23]. Specifically, the nanocrystal surfaces contain numerous dangling bonds. The dangling bonds which correspond to the undercoordinated surface atoms rearrange themselves by surface reconstruction or adsorb surfactant ligands [24]. The lat-ter are utilised during the colloidal synthesis as they increase the stability of nanocrystals, prevent crystal stacking and screen them from the environment as well as, most importantly, stabilise the surface by saturating dangling bonds [25]. As a result ligands influence surface trap states and, consequently, control photoluminescence (PL) properties [26,27]. Since the PL signal occurs due to the electron-hole recombination the optical properties of colloidal nanocrystals are directly controlled by different kinds of surface traps (reversible and irreversible).
An example of direct experimental manifestation of the surface trap role in the PL signals is the phenomenon of QDs blinking [28,29]. Most probable mechanisms of this phenomenon were discussed in details in [30,31] with time scale-free properties of the blinking process shown to be a signature of complex surface-dependent behaviour of the system.
The NPLs represent a quasi-2D-systems with a large surface to volume ratio and thicknesses down to few atomic layers [27,32,33]. In [32] the authors experimentally demonstrate a major effect of temporary charge carrier trapping on the decay dynamics of pure, core-shell and core-crown CdSe NPLs with results matching the previous measurements for NPLs [34,35]. The powerlaw decay of PL intensity is observed at long times and authors even discerned transitional power laws in some of the cases (for the description of this observation in Ref. [32] the term multiexponential was used). This type of decay was associated with reversible trapping without non-radiative losses. The power laws changed with temperature rather weakly (the authors drew parallel lines on a log-log plot, i.e. stated that the effect is temperature independent, however, looking carefully one is able to identify some variation with temperature). For blinking from single CdSe QDs the power-law statistics has been found to be independent on the temperature [30,31] which points out that blinking and delayed emission may share the same physical origin. Moreover, the model of irreversible charge trapping which leads to the nonradiative recombination was proposed in [36] to explain the discrepancy between PL and transient absorption measurements on NPLs. Among other structures with nontrivial power-law behavior of PL two-dimensional films of GaTe or GaSe are worth mentioning [37].
Currently for description of PL intensity curves most of the theoretical/simulation approaches use kinetic rate equations [27,32,35,37]. The difficulty of this paradigm applied for spatially distributed systems lies in a challenge of interpreting the meaning of the rates from a physical point of view. That makes the experimental estimation and verification of the parameters often impossible. Also simple kinetic rate equations can not produce an explanation for power-law statistics observed in experiments [32]. In the latter reference an ad-hoc term in their kinetic explanation is used to obtain the powerlaw. However, the origin or meaning of this term is not elucidated. A few attempts were made to include diffusion of excitons explicitly. For instance, in Refs. [38,39] CdSe/CdS core/crown nanoplatelets were analysed and it was experimentally demonstrated that exciton localisation efficiency is independent of crown size and increases with photon energy above the band edge, while the localisation time increases with the crown size. To analyse their experimental results authors of Refs. [38,39] also proposed a theoretical 2D model which considers the competition of in-plane exciton diffusion and selective hole trapping at the core/crown interface. In practice the approach was in a reduction of description to the system of kinetic equations and is not capable of prediction of any non-trivial power-law kinetics. The complexity of exciton dynamics was also found for colloidal halide perovskite nanoplateletes [40], where the presence of both dark and bright exciton populations was responsible for the complex excitons dynamics.
The experimental observation of exciton diffusion was found in monolayer semiconductors for both freestanding and SiO 2 supported W S 2 monolayers [41,42]. The supporting theoretical work considered non-equilibrium exciton transport in monolayer transition metal dichalcogenides where the interactions between excitons and nonequilibrium phonons were taken into account [43]. In [34] the influence of defects on the spontaneous and stimulated emission performances of solution-processed atomically flat quasi-2D nanoplatelets (NPLs) as a function of their lateral size using colloidal CdSe core NPLs was systematically investigated. It was revealed that the photoluminescence quantum efficiency of these NPLs decreases with increasing lateral size of the colloidal particles while their photoluminescence decay rate accelerates. This strongly suggests that nonradiative channels prevail in the NPL ensembles with platelets of extended lateral size. The effect clearly happens due to the bigger defected NPL subpopulations in NPLs with larger surfaces. Exci- ton diffusion in 2D metal-halide perovskites was studied in [44]. In this experiment the transition from diffusive to subdiffusive regime is clear which points out to the existence of traps with a power-law distribution of trapping times [45].
In this paper we draw attention to the important role of diffusive nature of excitons in their free state and a power-law statistics of escape from the traps. Our diffusion-based simulations and theoretical analytically solvable model with rates connected to diffusion/trapping properties explain non-trivial kinetics of carriers in the semiconductor nanoplatelets observed in experiments [32] and allow clear physical interpretation of the rates.
In the following sections we present our simulation and theoretical models and compare the results with experimental data.

I. SIMULATION MODEL AND COMPARISON TO EXPERIMENTS
Our simulation model assumes that after the laser beam produces excitons the latter start diffusing on the crystal lattice of a nanoplatelet. During this diffusion process a recombination with a spontaneous photoluminescence or trapping on a surface defect (trap) can happen. We assume that the recombination of excitons happens only in the free state. Once the particle hits a trap it is considered to be bound (trapped). Eventually excitons leave the traps and again could undergo diffusion, trapping or recombination.
In the model we neglect the possibility of recombination for the trapped excitons. Usually, the trapped PL is well pronounced only for thin nanoplatelets with the width of 1 or 2 monolayers. Its intensity decreases with the growth of number of monolayers. Since we simulated 5 layers which typically corresponds to experiments. Hence, this effect should not be critical for the predictions. Moreover, it is rather weak for the core-shell and core-crown nanoplatelets [34,46]. We also neglect such effects as nonradiative recombination or exciton dissociation. The role of these effects can be controlled experimentally by measuring the quantum yield of the nanoplatelet. The state-of-the-art nanoplatelet synthesis allows to obtain samples with quantum yield varying in the wide range from 10 percent to nearly 100 percent [47]. The role of nonradiative recombination processes and exciton dissociation processes decreases with the growth of the quantum yield. Samples with the quantum yield around 95 percent were synthesised and studied, for example, in Refs. [48][49][50].
One can estimate the concentration of empty dangling bonds on the nanocrystal surface which are considered to be effective traps following the parameters obtained in [27]. For CdSe nanoplatelets without shell with the typical largest area size about 100 nm 2 there exists one paramagnetic center associated with the dangling bond state per 50 Cd surface atoms. In our simulations we choose the values close to that estimate.
In order to simulate the exciton diffusion we used a discrete random walk on the 2D square lattice. The finite N × N lattice with periodic boundary conditions was taken. At every step the simulated particle performs jumps to the nearest neighbours (see Fig. 1) with equal probability 1 4 . The surface defects (traps) were fixed at the left-down corner of the lattice while the initial exciton position was sampled from the discrete uniform random distribution. The time ∆t of the jump along the lattice was modelled by a constant value much smaller than the time scale of photoluminescence and trapping processes. The exciton in a free state can recombine with a probability p at each jump, i.e. we neglect any interactions between the excitons which is perfectly justified for low-to-moderate intensities of an initial laser pulse [51]. The trapping occurs with the unit probability when the exciton reaches the trap position. Importantly, in our model the probability distribution of trapping times is non-exponential, having power-law tail such that the mean trapping time is infinite. Such trapping, or waiting time, distributions appear in the theory of transport in disordered media called Continuous Time Random Walk model [45,52,53], see also [54][55][56] and references therein. This distribution of trapping times occurs, for instance, in systems with an exponential distribution of depths of potential wells [57]. We sample the trapping times from the one-sided α-stable distribution with the Lévy index µ, 0 < µ < 1 following [58]. Such distribution possesses long time asymptotics γ(t) ∼ t −1−µ . The Laplace transform of the distribution reads γ(s) = e −τ µ * s µ . In simulations we use τ * = 1 ns, for more information about α-stable distributions see Appendix A. After a particle escapes from the trap it can be trapped again, i.e. the multiple trapping of excitons is taken into account by the model.
As a result the motion of an exciton in our model is characterised by a normal diffusive (Brownian) motion at short times and a subdiffusion with a characteristic power-law exponent µ at long times. Fig. 2 shows the dependence of an exciton mean-squared displacement as a function of time in the absence of decay. Indeed, at first, the particles diffuse freely but then the trapping and a subsequent release become important and the effective motion becomes subdiffusive.
The photoluminescence intensity is defined by the number of exciton recombinations per second, i.e. by the distribution of recombination times. In simulations we obtained this distribution from the set of N sim independent runs, i.e. we again use a low exciton density assumption. The normalised emission intensity was calculated by aggregation of the single recombination (photoluminescence) events in the bins of the size t bin (corresponding to the time resolution in the experiment [32]) and further normalisation by the number of events in the first bin in  Fig. 1c in the main text of the reference) (above) with the simulation curves (below) for the case of core CdSe NPLs. The simulations were performed on 15 × 15 periodic lattice (i.e. 1 defect per 225 sites), the exponent for trapping times was µ = 0.8, the step duration ∆t =1 ps and the probability of recombination per step p = 10 −3 . Nsim = 10 7 . Two different binnings are shown in the plot with simulation data. The green curve corresponds to t bin = 2500 ps while the red curve is for t bin = 500 ps. We have recalculated the slopes from Fig. 1c in Ref. [32]. For both experiment and simulations we show the range of slopes obtained with least mean squares method applied to different subintervals of a seemingly straight part of dependence.
the spirit of experimental work [32]. The middle of the corresponding bin was used as a time coordinate for the normalised emission intensity. Thus, our t bin corresponds to the exposition time in experiment [32]. Increase of t bin leads to the noise reduction on the one hand and to smoothing of the curve and possible loss of fine effects on the other hand. Another important fact is that the highest intensity of PL occurs at short times. Hence by enlarging the bin one gets disproportionately high counts to the first bin. In the latter case we effectively decrease the rest of the values.
In Figures 3,4 we show the comparison of our simulation model with fitted parameters (bottom plots) with The simulations were performed on 10 × 10 periodic lattice (i.e. 1 defect per 100 sites), the exponent for trapping times was µ = 0.8, the step duration ∆t =1 ps and the probability of recombination per step p = 10 −3 . Nsim = 10 7 . The green curve in the bottom plot corresponds to t bin = 2500 ps while the red curve is for t bin = 500 ps. We assume that defects are concentrated on the outer surface of CdS. This explains a different selection of trap density in our simulation model for this case (cf. Fig. 3). We have recalculated the slopes from Fig. 2g in Ref. [32] and for both experiment and simulations and show the range of slopes obtained with least mean squares method applied to different subintervals of a seemingly straight part of dependence. experimental data for the photoluminescence intensity in Ref. [32] (upper plots). Fig. 3 shows the PL intensity for core CdSe NPLs and Fig. 4 shows the PL intensity for core-shell CdSe-CdS nanoplatelets with CdS forming the outer layers. In order to plot and analyse the experimental data we have used the plots from the Supporting Information for Ref. [32] and digitised the data with the help of GetData Graph Digitizer program [59]. Black dots in the upper plots of Figs. 3, 4 are the averages of the digitised data with red region around representing the statistical error of measurements. One can see that the power laws including the intermediate ones nicely match between the simulations and experiments. For both the simulation and the experimental Ref. [32] results we show the range in exponents in order to show how sensitive it is to the size of the interval which visually looks "straight". By showing the results for two different bin sizes in simulations we illustrate that while the long time power-law asymptotics stays the same, the intermediate one is affected by renormalisation. This pushes us to the conclusion that the intermediate power laws are transitive phenomena rather than true power laws. Moreover the adjustment of the bin in simulations could produce a very good quantitative match with experimental data (see Appendix B). However, as we discovered the experimental data consist of two differently normalised parts (see the Supplemental Information for [32]) and also reveal a range of values for power law asymptotics depending on the time interval chosen for averaging. Hence we focus here on the explanation of the power laws rather than on the attempt to match the data quantitatively.
Note that we were not able to get a good match of the overall curve shapes with experiments for the core-crown case from Ref. [32] (still the long time power-law tail can be easily reproduced). We believe that the reason is that this case is markedly different. In comparison to pure core or core-shell cases it stands out as a system with a surface divided into two chemically distinctive areas with different properties. This obviously strongly affects the intermediate time scales. Hence, this case would require a more complex modelling approach.
The NPLs studied in experiments normally have a thickness of a few atomic layers. Hence the diffusion into inner layers is also possible. If one takes into account these bulk layers the results stay qualitatively the same which is shown in Appendix C.

II. THEORETICAL MODEL. NON-MARKOVIAN KINETIC RATE EQUATIONS
Our theoretical model corresponds to the simulation setup with one simplification. We additionally assume that at every step the trapping happens at a constant rate β and is proportional to the number of particles in free state βn f (t), n f (t) being the number of excitons in the free state as a function of time t. This coefficient β effectively describes diffusion of excitons and can be connected to properties of the simulation model (see below). The decay (recombination) is assumed to be proportional to the concentration n f (t) as αn f (t), where α is a rate of recombination process. Eventually an exciton escapes the trapped state and switches back to the diffusion-decay mode. For the population of excitons in the diffusive state the trapping process presents a temporary loss with a delayed return. This return can be described by the term with a kernel γ(t − τ ), Summarising the above model in terms of kinetic rate equations for concentrations of free excitons n f (t) and the trapped excitons n t (t) one gets The initial conditions are n f (0) = N 0 , n t (0) = 0, i.e. the laser excitation produces N 0 particles in the free state, but no trapped ones. Applying the Laplace transform L{f (t)} = ∞ 0 f (t)e −st dt =f (s) one turns the equations into algebraic ones. The solution forñ f (s) reads Since γ(t) is normalised, 1−γ(s) never becomes negative and any possible divergencies are avoided. Also it can be proven that the solution n f (t) is non-negative at all times (Appendix D). Clearly, in a general case the system of Eq.
(1) and (2) describes a non-Markovian process due to the memory kernel. The Markovian case can be obtained for the particular case of exponential distribution of trapping times (for more details and the exact solution of that case, see Appendices E,F). We assume that the power law tails observed in experiments appear due to the power law asymptotics of trapping times, i.e. one can use an α-stable law for the kernel,γ(s) = e −τ µ * s µ , 0 < µ < 1. In the limit of long times t → ∞ (which corresponds to s → 0) the functionγ(s) ≈ 1 − τ µ * s µ . Thus, in this limit the concentration of freely diffusing excitons is This expression produces a long-time power-law asymptotics n f (t τ * ) Cτ µ * t 1+µ with C = µβN0 α 2 Γ(1−µ) (see the derivation in Appendix G). Then the limit expression for the emission intensity is In order to compare the theoretical model with the simulations we need to find the correspondence between the simulation parameters and the constants used in the theory. As a test example we use the simulation results from Fig. 3 with the following parameters: the simulation step ∆t is 1ps, the bin size 500 ps, the time limit of the simulation 10000 ns, 2D square lattice with 1 defect per 15x15=225 points, µ = 0.8, probability of recombination per step 10 −3 . Hence, the parameters to be put in Eq. (3) can be estimated as follows. Forγ(s) = e −τ µ * s µ the coefficients are µ = 0.8, τ * = 1 ns. The rate of recombination α can be determined as a probability of recombination per step divided by step's duration, i.e. α = 1ns −1 . The rate of exciton capture by defects (traps) β is roughly the probability of find a defect per step. This probability equals 1/225 if one neglects the effect of secondary trapping in which case the probability of being caught by the defect is higher after the escape from a defect. For the case of the first trapping event the value of the probability is based on the fact that original excitations occur homogeneously in the sample. For the first trapping case β firsttrap can be then estimated as β firsttrap = 1 225 /∆t = 4.44 ns −1 (the trapping rate is almost constant at short times). In order to account for the undercounting of trapping events at long times in our theoretical model we slightly increase the trapping rate and set β = 7.5 ns −1 . At long times the secondary trapping is quite likely since right after the escape from the trap a particle is more likely to come back then at earlier times. Fig. 5 shows the excellent correspondence between the theoretical and simulation results obtained for these parameters. Since the simulation results are normalised on the first bin the theoretical curve was scaled such that the first points roughly coincide.
The study of parameter space both in simulation and theoretical models shows that at the intermediate times a rich behaviour can be observed. In Fig. 6 for the chosen parameters the curve weekly oscillates which is a result of an interplay of trapping, release and recombination terms both in simulations (the main plot) and in theory (the inset). The shapes of the emission intensity curves are defined by the interplay between free diffusion, trapping and recombination of excitons. One can see from Fig. 7 that the overall shape is controlled by diffusion properties. For any set of parameters at short times emission is defined by recombination only and the intensity decays exponentially. The latter observation was made experimentally, e.g. see Ref. [27]. Exponential decay at short times also follows from Eq.(3). At large s which correspond to the short time domainñ f (s) ∼ N 0 /(s + α + β), that is n f (t) ∼ N 0 exp(−(α + β)t). At short times most of the excitons are free. In the course of time they get trapped and escape the traps later. Since the trapping times are described by a power-law with exponent µ + 1 at long times the emission is defined by the release of excitons and characterised by the same power-law. The parameter β defines both the duration of the short time limit and the properties of the transition region. At small β (slow diffusion) the transition from exponential to powerlaw decay happens without a noticeable transition region (see the black solid curve in Fig. 7). Then with a growth of β the intermediate region appears which is well illustrated by the pink dash-dotted curve in Fig. 7. This is a regime when intermediate power laws were suggested in Ref. [32]. Finally, for fast diffusion (green dash dotdotted curve in Fig. 7) one can observe a wide transition region with oscillations of emission intensity. These oscillations could be associated with multiple trapping and release events.

III. SUMMARY
We proposed a diffusion-based simulation model and non-Markovian kinetic rate equations with a delay function, which describe the interplay between free diffusion, trapping and recombination of excitons in 2D semiconductor nanostructures. The parameters of both models have clear physical interpretation. The models were applied for the interpretation of experimental data on photoluminescence kinetics in CdSe/CdS core and coreshell nanoplatelets. The characteristic power-law tails are the result of exciton escape properties from surface defect traps. We show that the intermediate power laws found in Ref. [32] are transitional effects rather than true power-law dependencies. Moreover, we see that diffusion properties control the emission intensity in the intermediate region between the simple exponential decay at short times and a power-law at long times. The duration and the features of this intermediate region vary substantially with the value of diffusivity. In particular, in the case of a fast diffusion complex oscillating dependencies are demonstrated. The proposed model can be applied for the analysis of exciton kinetics in semiconductor nanoplatelets (core, core-shell, core-crown) and for estimation of microscopic exciton parameters. In the Appendix C we show a possible pathway for the model generalisation for the 3D case. Further generalisations of the proposed model and additional assumptions may be necessary to use the model for other low-dimensional semiconductor structures.

Conflicts of interest
There are no conflicts to declare.
In our model we choose a one-sided (completely asymmetric) α-stable Lévy distribution as a model distribution for the statistics of trapping times. The reasons for that are both mathematical and practical.
Mathematically, according to the generalised central limit theorem the distribution of the properly normalised sum of independent identically distributed variables converges to α-stable Lévy law if the variables are drawn from a distribution with a divergent second moment [60]. The α-stable distributions can only be written analytically in terms of often hard-to-grasp Fox H-functions with three exceptions (Gaussian, Cauchy and Lévy-Smirnov distributions). However, their Fourier/Laplace transforms adapt a simpler stretched Gaussian/exponential form [61]. In the case of trapping time distributions the quantities adopt the positive values only, i.e. one needs to consider a particular case of one-sided alpha-stable distributions. For the latter the Laplace transform reads where 0 < µ ≤ 1, s ∈ [0, ∞). At long times this distribution exhibits a power-law asymptotics which decays (up to a numerical prefactor) as τ µ * t 1+µ . The practical reason for our choice is a convenience of working with exponential dependencies both in terms of expansions and numerical calculations. In simulations the distribution was sampled according to Ref. [58]. If one plots the simulation and experimental results from Fig. 4 for the binning 2.5 ns = 2500 ps one finds even a very good quantitative match (Fig. 8). However, one has to keep in mind that the values depend on the normalisation which was not always known from experimental data in Ref. [32]. with equal probabilities of 1 6 In order to account for symmetry of NPL, the crossing of the symmetry plane (0, 0, +1) is replaced by (0, 0, 0) for even N l (2N l layers in NPL) or by (0, 0, −1) for the odd N l (2(N l − 1) + 1 layers in NPL).
In Fig. 9 we see that for core NPLs 2D model and the model of diffusion in a few parallel atomistic layers give the same results. In this case we assume the density of surface defects the same as it was in Fig. 3, i.e. 1 defect per 225 surface lattice sites. Fig. 10 shows comparison of three different models. Namely, the 2D model, the 3D model which assumes the same properties and jump probabilities in all 5 layers and the 3D model, where jump probabilities differ when the jumps occur from the surface layers into the bulk. The latter case is the easiest generalisation of a 3D model. It assumes that the jumps FIG. 10: Comparison of curves for photoluminescence intensities for 2D (red), 3D homogeneous (blue) and 3D distinct domains (cyan) models. The core-shell case is shown with a density of defects 1 per 100 surface lattice sites. µ = 0.8, Nsim = 10 6 , t bin = 500ps. The kink at short times occurs due to binning. In the first two cases the jumps are considered to be equally likely for all directions. For the last case (cyan data) the layers are split between two regions. The outer region is one material and the inner region is another one. Thus, it corresponds to the core-shell structure more realistically. The jumps from the surface layer into the bulk layers and back are 5 times smaller than the probabilities of jumps within the same region.
into the bulk and back to the surface are 5 times less likely than the jumps within the same layer. The density of defects corresponds to Fig. 4, i.e. 1 per 100 surface lattice sites. All three curves show the same tail behaviour as one would expect. The transitional region changes its shape and shifts, i.e. more complex models allow for finer description. We see that the available experimental data can be fitted well with a relatively simple 2D approach which tells us that 2D model is sufficient as a minimal model. The concentration of excitons must be non-negative. Hence, according to the Bernstein theorem [62], its Laplace transformñ f (s) should be completely monotonic (c.m.). The functiong(s) is called c.m. if (−1) ng(n) (s) ≥ 1 for s > 0 and n = 0, 1, 2.... If m 1 (s) and m 2 (s) are completely monotonic then their product m 1 (s)m 2 (s) is also c.m. Also iff (s) is c.m. andh(s) is Bernstein function (i.e.h ≥ 0 for positive s andh (s) is c.m.) thenf (h(s)) is c.m. [62].
Appendix E: The derivation of the limit result for the two state Markovian limit We show here how one can recover the limit of a two state Markovian system with constant rates of exchange and a decay in one of the states (considered, for instance, in Van Kampen's book [63], Chapter VII, section 5, exercise (5.8)). This corresponds to the Markovian limit of Eqs. (1)- (2).
The case of two states with exchange can be written as dn t (t) dt = βn f (t) − n t (t)/τ * , where n f (t) and n t (t) correspond to the free and trapped states in our model and α, β and τ * are positive constants. Let us treat Eq. (8) formally as a non-homogeneous ordinary differential equation of the first order. Then its formal solution reads n t (t) = t 0 dt e −(t−t )/τ * βn f (t ).