Complex diffusion-based kinetics of photoluminescence in semiconductor nanoplatelets
Aleksandr A. Kurilovich, Vladimir N. Mantsevich, Keith J. Stevenson, Aleksei V. Chechkin, Vladimir V. Palyulin
CComplex diffusion-based kinetics of photoluminescence in semiconductor nanoplatelets
A. A. Kurilovich, V. N. Mantsevich, K. J. Stevenson, A. V. Chechkin,
3, 4 and V. V. Palyulin Center for Energy Science and Technology, Skolkovo Institute of Science and Technology, 121205, Moscow, Russia Chair of Semiconductors and Cryoelectronics & Quantum Technology Center,Physics department, Lomonosov Moscow State University, 119991 Moscow, Russia Institute for Physics & Astronomy, University of Potsdam, D-14476 Potsdam-Golm, Germany Akhiezer Institute for Theoretical Physics National Science Center”Kharkov Institute of Physics and Technology”, 61108, Kharkov, Ukraine Center for Computational and Data-Intensive Science and Engineering,Skolkovo Institute of Science and Technology, 121205, Moscow, Russia (Dated: January 19, 2021)We present a diffusion-based simulation and theoretical models for explanation of photolumines-cence (PL) emission intensity in semiconductor nanoplatelets. It is shown that the shape of PLintensity curves can be reproduced by the interplay of recombination, diffusion and trapping of ex-citons. 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 bya power-law tail. We show that the crossover from one limit to another is controlled by diffusionproperties. This intermediate region exhibits a rich behaviour depending on the value of diffusiv-ity. Proposed approach reproduces all the features of experimental curves measured for differentnanoplatelet systems.
PACS numbers:
The hunt for materials and systems with better opti-cal properties has always been one of the focal points ofsemiconductor research [1]. The first classical optoelec-tronical applications were based mostly on bulk proper-ties of semiconductors. Recent progress in material syn-thesis [2–4] led to wider research efforts concentrated onlow-dimensional structures. In particular, significant at-tention is drawn by the semiconductor nanocrystals, suchas quantum dots (QDs) and nanoplatelets (NPLs) [5, 6].QDs and NPLs show large exciton binding energy [7–9] and possess such remarkable properties as narrowemission lines even at room temperature, tunable emis-sion wavelength, short radiative lifetimes, giant oscilla-tor strength, high quantum yield and vanishing inho-mogeneous broadening [10–14]. These properties makethem excellent prospective building blocks for future op-toelectronic devices such as bright and flexible light emit-ters [15], colloidal lasers [8, 16] or for biomedical labelingapplications [17].Modern techniques of precise colloidal synthesis allowto control the shape [18], size [19] and crystal structure ofthe nanocrystals [20]. Among these a growing popularitybelongs to colloidal QDs and NPLs of CdSe with atomi-cally defined thickness [5, 21, 22] and various shapes suchas pure QDs and NPLs (core) as well as composite core-shell QDs and NPLs or core-crown NPLs. The compos-ites are made of two different semiconductor materials,for instance, CdSe and CdS.The low-dimensionality of all these structures effectsin a decisive contribution of a crystal surface to controlof optical and electronic properties [23]. Specifically, thenanocrystal surfaces contain numerous dangling bonds.The dangling bonds which correspond to the underco-ordinated surface atoms rearrange themselves by surfacereconstruction or adsorb surfactant ligands [24]. The lat- ter are utilised during the colloidal synthesis as they in-crease the stability of nanocrystals, prevent crystal stack-ing and screen them from the environment as well as,most importantly, stabilise the surface by saturating dan-gling bonds [25]. As a result ligands influence surfacetrap states and, consequently, control photoluminescence(PL) properties [26, 27]. Since the PL signal occurs dueto the electron-hole recombination the optical propertiesof colloidal nanocrystals are directly controlled by differ-ent kinds of surface traps (reversible and irreversible).An example of direct experimental manifestation of thesurface trap role in the PL signals is the phenomenon ofQDs blinking [28, 29]. Most probable mechanisms of thisphenomenon were discussed in details in [30, 31] withtime scale-free properties of the blinking process shownto be a signature of complex surface-dependent behaviourof the system.The NPLs represent a quasi-2D-systems with a largesurface to volume ratio and thicknesses down to fewatomic layers [27, 32, 33]. In [32] the authors experimen-tally demonstrate a major effect of temporary charge car-rier trapping on the decay dynamics of pure, core-shelland core-crown CdSe NPLs with results matching theprevious measurements for NPLs [34, 35]. The power-law decay of PL intensity is observed at long times andauthors even discerned transitional power laws in some ofthe cases (for the description of this observation in Ref.[32] the term multiexponential was used). This type ofdecay was associated with reversible trapping withoutnon-radiative losses. The power laws changed with tem-perature rather weakly (the authors drew parallel lineson a log-log plot, i.e. stated that the effect is tempera-ture independent, however, looking carefully one is ableto identify some variation with temperature). For blink-ing from single CdSe QDs the power-law statistics has a r X i v : . [ c ond - m a t . m e s - h a ll ] J a n been found to be independent on the temperature [30, 31]which points out that blinking and delayed emission mayshare the same physical origin. Moreover, the model ofirreversible charge trapping which leads to the nonradia-tive recombination was proposed in [36] to explain thediscrepancy between PL and transient absorption mea-surements on NPLs. Among other structures with non-trivial power-law behavior of PL two-dimensional filmsof GaTe or GaSe are worth mentioning [37].Currently for description of PL intensity curves mostof the theoretical/simulation approaches use kinetic rateequations [27, 32, 35, 37]. The difficulty of this paradigmapplied for spatially distributed systems lies in a chal-lenge of interpreting the meaning of the rates from aphysical point of view. That makes the experimentalestimation and verification of the parameters often im-possible. Also simple kinetic rate equations can not pro-duce an explanation for power-law statistics observed inexperiments [32]. In the latter reference an ad-hoc termin their kinetic explanation is used to obtain the power-law. However, the origin or meaning of this term is notelucidated. A few attempts were made to include diffu-sion of excitons explicitly. For instance, in Refs. [38, 39]CdSe/CdS core/crown nanoplatelets were analysed andit was experimentally demonstrated that exciton localisa-tion efficiency is independent of crown size and increaseswith photon energy above the band edge, while the lo-calisation time increases with the crown size. To anal-yse their experimental results authors of Refs. [38, 39]also proposed a theoretical 2D model which considers thecompetition of in-plane exciton diffusion and selectivehole trapping at the core/crown interface. In practicethe approach was in a reduction of description to the sys-tem of kinetic equations and is not capable of predictionof any non-trivial power-law kinetics. The complexity ofexciton dynamics was also found for colloidal halide per-ovskite nanoplateletes [40], where the presence of bothdark and bright exciton populations was responsible forthe complex excitons dynamics.The experimental observation of exciton diffusion wasfound in monolayer semiconductors for both freestand-ing and SiO supported W S monolayers [41, 42]. Thesupporting theoretical work considered non-equilibriumexciton transport in monolayer transition metal dichalco-genides where the interactions between excitons and non-equilibrium phonons were taken into account [43]. In [34]the influence of defects on the spontaneous and stimu-lated emission performances of solution-processed atom-ically flat quasi-2D nanoplatelets (NPLs) as a functionof their lateral size using colloidal CdSe core NPLs wassystematically investigated. It was revealed that thephotoluminescence quantum efficiency of these NPLs de-creases with increasing lateral size of the colloidal parti-cles while their photoluminescence decay rate accelerates.This strongly suggests that nonradiative channels prevailin the NPL ensembles with platelets of extended lateralsize. The effect clearly happens due to the bigger defectedNPL subpopulations in NPLs with larger surfaces. Exci- FIG. 1: In our simulations we used a model of a randomwalk on a 2D square lattice modelled as N × N set of vertices(atoms) with periodic boundary conditions. The jump prob-abilities are equal in four directions. Some of the vertices arecapable of trapping an exciton with a heavy-tail distributionof trapping times γ ( t → ∞ ) ∼ t µ , < µ <
1. If an excitonjumps on the trap position it gets bound to the cite until itescapes. The recombination is modelled by decay probability p per step. ton diffusion in 2D metal-halide perovskites was studiedin [44]. In this experiment the transition from diffusiveto subdiffusive regime is clear which points out to the ex-istence of traps with a power-law distribution of trappingtimes [45].In this paper we draw attention to the importantrole of diffusive nature of excitons in their free stateand a power-law statistics of escape from the traps.Our diffusion-based simulations and theoretical analyt-ically solvable model with rates connected to diffu-sion/trapping properties explain non-trivial kinetics ofcarriers in the semiconductor nanoplatelets observed inexperiments [32] and allow clear physical interpretationof the rates.In the following sections we present our simulation andtheoretical models and compare the results with experi-mental data. I. SIMULATION MODEL AND COMPARISONTO EXPERIMENTS
Our simulation model assumes that after the laserbeam produces excitons the latter start diffusing on thecrystal lattice of a nanoplatelet. During this diffusionprocess a recombination with a spontaneous photolumi-nescence or trapping on a surface defect (trap) can hap-pen. We assume that the recombination of excitons hap-pens only in the free state. Once the particle hits a trapit is considered to be bound (trapped). Eventually exci-tons leave the traps and again could undergo diffusion,trapping or recombination.In the model we neglect the possibility of recombi-nation for the trapped excitons. Usually, the trappedPL is well pronounced only for thin nanoplatelets withthe width of 1 or 2 monolayers. Its intensity decreaseswith the growth of number of monolayers. Since we sim-ulated 5 layers which typically corresponds to experi-ments. Hence, this effect should not be critical for thepredictions. Moreover, it is rather weak for the core-shelland core-crown nanoplatelets [34, 46]. We also neglectsuch effects as nonradiative recombination or exciton dis-sociation. The role of these effects can be controlledexperimentally by measuring the quantum yield of thenanoplatelet. The state-of-the-art nanoplatelet synthesisallows to obtain samples with quantum yield varying inthe wide range from 10 percent to nearly 100 percent [47].The role of nonradiative recombination processes and ex-citon dissociation processes decreases with the growthof the quantum yield. Samples with the quantum yieldaround 95 percent were synthesised and studied, for ex-ample, in Refs. [48–50].One can estimate the concentration of empty danglingbonds on the nanocrystal surface which are consideredto be effective traps following the parameters obtainedin [27]. For CdSe nanoplatelets without shell with thetypical largest area size about 100 nm there exists oneparamagnetic center associated with the dangling bondstate per 50 Cd surface atoms. In our simulations wechoose the values close to that estimate.In order to simulate the exciton diffusion we used adiscrete random walk on the 2D square lattice. The fi-nite N × N lattice with periodic boundary conditionswas taken. At every step the simulated particle performsjumps to the nearest neighbours (see Fig. 1) with equalprobability . The surface defects (traps) were fixed atthe left-down corner of the lattice while the initial excitonposition was sampled from the discrete uniform randomdistribution. The time ∆ t of the jump along the lat-tice was modelled by a constant value much smaller thanthe time scale of photoluminescence and trapping pro-cesses. The exciton in a free state can recombine witha probability p at each jump, i.e. we neglect any inter-actions between the excitons which is perfectly justifiedfor low-to-moderate intensities of an initial laser pulse[51]. The trapping occurs with the unit probability whenthe exciton reaches the trap position. Importantly, inour model the probability distribution of trapping timesis non-exponential, having power-law tail such that themean trapping time is infinite. Such trapping, or waitingtime, distributions appear in the theory of transport indisordered media called Continuous Time Random Walkmodel [45, 52, 53], see also [54–56] and references therein.This distribution of trapping times occurs, for instance,in systems with an exponential distribution of depths ofpotential wells [57]. We sample the trapping times fromthe one-sided α -stable distribution with the L´evy index FIG. 2: Mean-squared displacement for exciton diffusion asa function of time on double logarithmic scale in the absenceof decay. The transition from a normal diffusive regime tosubdiffusion is observed. The red curve corresponds to theparameters we use for core structures, i.e. the simulationswere performed on 15 ×
15 periodic lattice (i.e. 1 defect per225 sites), the exponent for trapping times was µ = 0 . t =1 ps. The green curve stands for core-shell NPLs and we used 10 ×
10 periodic lattice (i.e. 1 defectper 100 sites), the exponent for trapping times was µ = 0 . t =1 ps. The curves were obtained byensemble averaging over 1000 trajectories. The numbers 1.0and 0.8 show the corresponding slopes. µ, < µ < γ ( t ) ∼ t − − µ . The Laplace trans-form of the distribution reads γ ( s ) = e − τ µ ∗ s µ . In sim-ulations we use τ ∗ = 1 ns, for more information about α -stable distributions see Appendix A. After a particleescapes from the trap it can be trapped again, i.e. themultiple trapping of excitons is taken into account by themodel.As a result the motion of an exciton in our model ischaracterised by a normal diffusive (Brownian) motionat short times and a subdiffusion with a characteristicpower-law exponent µ at long times. Fig. 2 shows thedependence of an exciton mean-squared displacement asa function of time in the absence of decay. Indeed, at first,the particles diffuse freely but then the trapping and asubsequent release become important and the effectivemotion becomes subdiffusive.The photoluminescence intensity is defined by thenumber of exciton recombinations per second, i.e. by thedistribution of recombination times. In simulations weobtained this distribution from the set of N sim indepen-dent runs, i.e. we again use a low exciton density assump-tion. The normalised emission intensity was calculatedby aggregation of the single recombination (photolumi-nescence) events in the bins of the size t bin (correspondingto the time resolution in the experiment [32]) and furthernormalisation by the number of events in the first bin in FIG. 3: Comparison of digitised experimental data taken fromSI of Ref. [32] (correspond to Fig. 1c in the main text of thereference) (above) with the simulation curves (below) for thecase of core CdSe NPLs. The simulations were performedon 15 ×
15 periodic lattice (i.e. 1 defect per 225 sites), theexponent for trapping times was µ = 0 .
8, the step duration∆ t =1 ps and the probability of recombination per step p =10 − . N sim = 10 . Two different binnings are shown in theplot with simulation data. The green curve corresponds to t bin = 2500 ps while the red curve is for t bin = 500 ps. Wehave recalculated the slopes from Fig. 1c in Ref. [32]. Forboth experiment and simulations we show the range of slopesobtained with least mean squares method applied to differentsubintervals of a seemingly straight part of dependence. the spirit of experimental work [32]. The middle of thecorresponding bin was used as a time coordinate for thenormalised emission intensity. Thus, our t bin correspondsto the exposition time in experiment [32]. Increase of t bin leads to the noise reduction on the one hand and tosmoothing of the curve and possible loss of fine effectson the other hand. Another important fact is that thehighest intensity of PL occurs at short times. Hence byenlarging the bin one gets disproportionately high countsto the first bin. In the latter case we effectively decreasethe rest of the values.In Figures 3,4 we show the comparison of our simu-lation model with fitted parameters (bottom plots) with FIG. 4: Comparison of digitised experimental data taken fromSI of Ref. [32] (correspond to Fig. 2g in the main text)(above) with the simulation curves (below) for the case ofcore-shell CdSe-CdS NPLs with CdS forming the outer layers.The simulations were performed on 10 ×
10 periodic lattice(i.e. 1 defect per 100 sites), the exponent for trapping timeswas µ = 0 .
8, the step duration ∆ t =1 ps and the probabil-ity of recombination per step p = 10 − . N sim = 10 . Thegreen curve in the bottom plot corresponds to t bin = 2500ps while the red curve is for t bin = 500 ps. We assume thatdefects are concentrated on the outer surface of CdS. Thisexplains a different selection of trap density in our simulationmodel for this case (cf. Fig. 3). We have recalculated theslopes from Fig. 2g in Ref. [32] and for both experiment andsimulations and show the range of slopes obtained with leastmean squares method applied to different subintervals of aseemingly straight part of dependence. experimental data for the photoluminescence intensity inRef. [32] (upper plots). Fig. 3 shows the PL intensityfor core CdSe NPLs and Fig. 4 shows the PL intensityfor core-shell CdSe-CdS nanoplatelets with CdS formingthe outer layers. In order to plot and analyse the experi-mental data we have used the plots from the SupportingInformation for Ref. [32] and digitised the data with thehelp of GetData Graph Digitizer program [59]. Blackdots in the upper plots of Figs. 3, 4 are the averages ofthe digitised data with red region around representing thestatistical error of measurements. One can see that thepower laws including the intermediate ones nicely matchbetween the simulations and experiments. For both thesimulation and the experimental Ref. [32] results weshow the range in exponents in order to show how sensi-tive it is to the size of the interval which visually looks”straight”. By showing the results for two different binsizes in simulations we illustrate that while the long timepower-law asymptotics stays the same, the intermediateone is affected by renormalisation. This pushes us to theconclusion that the intermediate power laws are transi-tive phenomena rather than true power laws. Moreoverthe adjustment of the bin in simulations could producea very good quantitative match with experimental data(see Appendix B). However, as we discovered the exper-imental data consist of two differently normalised parts(see the Supplemental Information for [32]) and also re-veal a range of values for power law asymptotics depend-ing on the time interval chosen for averaging. Hence wefocus here on the explanation of the power laws ratherthan on the attempt to match the data quantitatively.Note that we were not able to get a good match of theoverall curve shapes with experiments for the core-crowncase from Ref. [32] (still the long time power-law tail canbe easily reproduced). We believe that the reason is thatthis case is markedly different. In comparison to purecore or core-shell cases it stands out as a system with asurface divided into two chemically distinctive areas withdifferent properties. This obviously strongly affects theintermediate time scales. Hence, this case would requirea more complex modelling approach.The NPLs studied in experiments normally have athickness of a few atomic layers. Hence the diffusion intoinner layers is also possible. If one takes into accountthese bulk layers the results stay qualitatively the samewhich is shown in Appendix C. II. THEORETICAL MODEL.NON-MARKOVIAN KINETIC RATEEQUATIONS
Our theoretical model corresponds to the simulationsetup with one simplification. We additionally assumethat at every step the trapping happens at a constantrate β and is proportional to the number of particles infree state βn f ( t ), n f ( t ) being the number of excitons inthe free state as a function of time t . This coefficient β effectively describes diffusion of excitons and can beconnected to properties of the simulation model (see be-low). The decay (recombination) is assumed to be pro-portional to the concentration n f ( t ) as αn f ( t ), where α is a rate of recombination process. Eventually an ex-citon escapes the trapped state and switches back tothe diffusion-decay mode. For the population of exci-tons in the diffusive state the trapping process presentsa temporary loss with a delayed return. This returncan be described by the term with a kernel γ ( t − τ ), FIG. 5: The outcome of the theoretical model vs the simula-tion results in the double log scale for all relevant time scales.The parameters of theory are β = 7 . − , α = 1ns − , µ =0 .
8. In simulations β = 4 . − . Red dots correspond tothe simulation data. Blue points are the numerical inverseLaplace transform of Eq. (3). t bin = 500ps (cid:82) t γ ( t − τ ) βn f ( τ ) dτ , where t is a release time momentas compared to τ , the time of the trapping event. Thekernel has an explicit meaning of trapping probabilitydensity as a function of time, i.e. it abides a normali-sation property, (cid:82) ∞ γ ( t ) dt = 1. Hence, our theoreticalmodel describes the recombination-caused photoluminis-cence process as a two-state system model with only oneof the states allowing the recombination.Summarising the above model in terms of kinetic rateequations for concentrations of free excitons n f ( t ) andthe trapped excitons n t ( t ) one gets dn f ( t ) dt = − αn f ( t ) − βn f ( t ) + (cid:90) t γ ( t − τ ) βn f ( τ ) dτ, (1) dn t ( t ) dt = βn f ( t ) − (cid:90) t γ ( t − τ ) βn f ( τ ) dτ. (2)The initial conditions are n f (0) = N , n t (0) = 0, i.e. thelaser excitation produces N particles in the free state,but no trapped ones. Applying the Laplace transform L{ f ( t ) } = (cid:82) ∞ f ( t ) e − st dt = ˜ f ( s ) one turns the equationsinto algebraic ones. The solution for ˜ n f ( s ) reads˜ n f ( s ) = N s + α + β (1 − ˜ γ ( s )) . (3)Since γ ( t ) is normalised, 1 − ˜ γ ( s ) never becomes negativeand any possible divergencies are avoided. Also it can beproven 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 thememory kernel. The Markovian case can be obtained forthe particular case of exponential distribution of trappingtimes (for more details and the exact solution of thatcase, see Appendices E,F). We assume that the powerlaw tails observed in experiments appear due to the powerlaw asymptotics of trapping times, i.e. one can use an α -stable law for the kernel, ˜ γ ( s ) = e − τ µ ∗ s µ , < µ < t → ∞ (which corresponds to s →
0) the function ˜ γ ( s ) ≈ − τ µ ∗ s µ . Thus, in this limitthe concentration of freely diffusing excitons is˜ n f ( s ) ≈ N α (cid:18) − τ µ ∗ βα s µ (cid:19) . (4)This expression produces a long-time power-law asymp-totics n f ( t (cid:29) τ ∗ ) (cid:39) Cτ µ ∗ t µ with C = µβN α Γ(1 − µ ) (see thederivation in Appendix G). Then the limit expression forthe emission intensity is I ( t (cid:29) τ ∗ ) (cid:39) µβN α Γ(1 − µ ) τ µ ∗ t µ (5)In order to compare the theoretical model with thesimulations we need to find the correspondence betweenthe simulation parameters and the constants used in thetheory. As a test example we use the simulation resultsfrom Fig. 3 with the following parameters: the simula-tion step ∆ t is 1ps, the bin size 500 ps, the time limit ofthe simulation 10000 ns, 2D square lattice with 1 defectper 15x15=225 points, µ = 0 .
8, probability of recombi-nation per step 10 − . Hence, the parameters to be put inEq. (3) can be estimated as follows. For ˜ γ ( s ) = e − τ µ ∗ s µ the coefficients are µ = 0 . τ ∗ = 1 ns. The rate ofrecombination α can be determined as a probability ofrecombination per step divided by step’s duration, i.e. α = 1ns − . 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 ofsecondary trapping in which case the probability of be-ing caught by the defect is higher after the escape froma defect. For the case of the first trapping event thevalue of the probability is based on the fact that origi-nal excitations occur homogeneously in the sample. Forthe first trapping case β firsttrap can be then estimated as β firsttrap = / ∆ t = 4 .
44 ns − (the trapping rate is al-most constant at short times). In order to account forthe undercounting of trapping events at long times in ourtheoretical model we slightly increase the trapping rateand set β = 7 . − . At long times the secondary trap-ping is quite likely since right after the escape from thetrap a particle is more likely to come back then at earliertimes.Fig. 5 shows the excellent correspondence between thetheoretical and simulation results obtained for these pa-rameters. Since the simulation results are normalised onthe first bin the theoretical curve was scaled such thatthe first points roughly coincide.The study of parameter space both in simulation andtheoretical models shows that at the intermediate times arich behaviour can be observed. In Fig. 6 for the chosenparameters the curve weekly oscillates which is a resultof an interplay of trapping, release and recombinationterms both in simulations (the main plot) and in theory(the inset). FIG. 6: The simulation results (the main plot) with fine bin-ning showing the same intermediate fluctuating behaviour asin theoretical model (shown in the inset). Here the bin sizewas chosen to be t bin = 25 ps, µ = 0 . p = 10 − , N sim = 10 ,one trap was put per 100 lattice sites.FIG. 7: Dependence of logarithm of emission intensity as afunction of logarithm of time for different β . The parameter β is defined by the free diffusion coefficient of excitons. Thevariation of β leads to the change in a characteristic shape ofthe emission intensity curves. The shapes of the emission intensity curves are definedby the interplay between free diffusion, trapping and re-combination of excitons. One can see from Fig. 7 thatthe overall shape is controlled by diffusion properties. Forany set of parameters at short times emission is definedby recombination only and the intensity decays exponen-tially. The latter observation was made experimentally,e.g. see Ref. [27]. Exponential decay at short timesalso follows from Eq.(3). At large s which correspond tothe short time domain ˜ n f ( s ) ∼ N / ( s + α + β ), that is n f ( t ) ∼ N exp( − ( α + β ) t ). At short times most of theexcitons are free. In the course of time they get trappedand escape the traps later. Since the trapping times aredescribed by a power-law with exponent µ + 1 at longtimes the emission is defined by the release of excitonsand characterised by the same power-law. The param-eter β defines both the duration of the short time limitand the properties of the transition region. At small β (slow diffusion) the transition from exponential to power-law decay happens without a noticeable transition region(see the black solid curve in Fig. 7). Then with a growthof β the intermediate region appears which is well illus-trated by the pink dash-dotted curve in Fig. 7. This isa regime when intermediate power laws were suggestedin Ref. [32]. Finally, for fast diffusion (green dash dot-dotted curve in Fig. 7) one can observe a wide transitionregion with oscillations of emission intensity. These os-cillations could be associated with multiple trapping andrelease events. III. SUMMARY
We proposed a diffusion-based simulation model andnon-Markovian kinetic rate equations with a delay func-tion, which describe the interplay between free diffusion,trapping and recombination of excitons in 2D semicon-ductor nanostructures. The parameters of both mod-els have clear physical interpretation. The models wereapplied for the interpretation of experimental data onphotoluminescence kinetics in CdSe/CdS core and core-shell nanoplatelets. The characteristic power-law tailsare the result of exciton escape properties from surfacedefect traps. We show that the intermediate power lawsfound in Ref. [32] are transitional effects rather thantrue power-law dependencies. Moreover, we see that dif-fusion properties control the emission intensity in the in-termediate region between the simple exponential decayat short times and a power-law at long times. The du-ration and the features of this intermediate region varysubstantially with the value of diffusivity. In particular,in the case of a fast diffusion complex oscillating depen-dencies are demonstrated. The proposed model can beapplied for the analysis of exciton kinetics in semicon-ductor nanoplatelets (core, core-shell, core-crown) andfor estimation of microscopic exciton parameters. In theAppendix C we show a possible pathway for the modelgeneralisation for the 3D case. Further generalisationsof the proposed model and additional assumptions maybe necessary to use the model for other low-dimensionalsemiconductor structures.
Conflicts of interest
There are no conflicts to declare.
Acknowledgements
The authors would like to acknowledge the support ofRussian Science Foundation project No. 18-72-10002 aswell as fruitful discussions with Igor M. Sokolov.
Appendix A: α -stable L´evy laws In our model we choose a one-sided (completely asym-metric) α -stable L´evy distribution as a model distribu-tion for the statistics of trapping times. The reasons forthat are both mathematical and practical.Mathematically, according to the generalised centrallimit theorem the distribution of the properly normalisedsum of independent identically distributed variables con-verges to α -stable L´evy law if the variables are drawnfrom a distribution with a divergent second moment [60].The α -stable distributions can only be written analyti-cally in terms of often hard-to-grasp Fox H-functions withthree exceptions (Gaussian, Cauchy and L´evy-Smirnovdistributions). However, their Fourier/Laplace trans-forms adapt a simpler stretched Gaussian/exponentialform [61]. In the case of trapping time distributions thequantities adopt the positive values only, i.e. one needsto consider a particular case of one-sided alpha-stabledistributions. For the latter the Laplace transform reads G ( s ) = exp( − τ µ ∗ s µ ) , (6)where 0 < µ ≤ , s ∈ [0 , ∞ ). At long times this distribu-tion exhibits a power-law asymptotics which decays (upto a numerical prefactor) as τ µ ∗ t µ .The practical reason for our choice is a convenienceof working with exponential dependencies both in termsof expansions and numerical calculations. In simulationsthe distribution was sampled according to Ref. [58]. Appendix B: Match between theory and experimentfor the binning value t bin = 2 . ns If one plots the simulation and experimental resultsfrom Fig. 4 for the binning 2.5 ns = 2500 ps one findseven a very good quantitative match (Fig. 8). However,one has to keep in mind that the values depend on thenormalisation which was not always known from experi-mental data in Ref. [32].
Appendix C: Modelling the process in a fewatomistic layers
For the 3D model the half of NPL vs. the symmetryplane parallel to the surface was simulated with the N l layers ( N × N × N l ). The traps are located only in thesurface layer. Periodic boundary conditions are appliedto in-plane coordinates. The initial exciton positions aregenerated from discrete uniform distributions along allthree coordinates. The exciton cannot escape to the so-lution so the possible jumps on the surface layer ( N, N, , , , ( − , , , (0 , +1 , , (0 , − , , (0 , , +1)]with the equal probabilities of . Inthe internal layer the possible jumps are[(+1 , , , ( − , , , (0 , +1 , , (0 , − , , (0 , , +1) , (0 , , − FIG. 8: Quantitative match for photoluminescence intensitiesbetween the simulations and experiment from Ref. [32] forthe binning t bin = 2 . µ = 0 . , N sim = 10 , t bin = 500 ps . with equal probabilities of In order to account forsymmetry of NPL, the crossing of the symmetry plane(0 , , +1) is replaced by (0 , ,
0) for even N l (2 N l layersin NPL) or by (0 , , −
1) for the odd N l (2( N l −
1) + 1layers in NPL).In Fig. 9 we see that for core NPLs 2D model and themodel of diffusion in a few parallel atomistic layers givethe same results. In this case we assume the density ofsurface defects the same as it was in Fig. 3, i.e. 1 defectper 225 surface lattice sites. Fig. 10 shows comparisonof three different models. Namely, the 2D model, the3D model which assumes the same properties and jumpprobabilities in all 5 layers and the 3D model, where jumpprobabilities differ when the jumps occur from the sur-face layers into the bulk. The latter case is the easiestgeneralisation of a 3D model. It assumes that the jumps
FIG. 10: Comparison of curves for photoluminescence inten-sities for 2D (red), 3D homogeneous (blue) and 3D distinctdomains (cyan) models. The core-shell case is shown with adensity of defects 1 per 100 surface lattice sites. µ = 0 . N sim = 10 , t bin = 500 ps . The kink at short times occurs dueto binning. In the first two cases the jumps are consideredto be equally likely for all directions. For the last case (cyandata) the layers are split between two regions. The outerregion is one material and the inner region is another one.Thus, it corresponds to the core-shell structure more realisti-cally. The jumps from the surface layer into the bulk layersand back are 5 times smaller than the probabilities of jumpswithin the same region. into the bulk and back to the surface are 5 times lesslikely than the jumps within the same layer. The densityof defects corresponds to Fig. 4, i.e. 1 per 100 surface lat-tice sites. All three curves show the same tail behaviouras one would expect. The transitional region changes itsshape and shifts, i.e. more complex models allow for finerdescription. We see that the available experimental datacan be fitted well with a relatively simple 2D approachwhich tells us that 2D model is sufficient as a minimalmodel. Appendix D: Proof of the non-negativity of Eq. (1)solutions
The concentration of excitons must be non-negative.Hence, according to the Bernstein theorem [62], itsLaplace transform ˜ n f ( s ) should be completely monotonic(c.m.). The function ˜ g ( s ) is called c.m. if ( − n ˜ g ( n ) ( s ) ≥ s > n = 0 , , ... . If m ( s ) and m ( s ) are com-pletely monotonic then their product m ( s ) m ( s ) is alsoc.m. Also if ˜ f ( s ) is c.m. and ˜ h ( s ) is Bernstein function(i.e. ˜ h ≥ s and ˜ h (cid:48) ( s ) is c.m.) then ˜ f (˜ h ( s ))is c.m. [62].In order to show the complete monotonicity of ˜ n f ( s )we present it as ˜ n f ( s ) = ˜ f (˜ h ( s )), where ˜ f = 1 /s and˜ h ( s ) = s + α + β (1 − ˜ γ ( s )). The function ˜ h ( s ) is positivesince ˜ γ ( s ) is the Laplace transform of the probability dis-tribution function and, therefore, ˜ γ ( s ) ≤
1. The function˜ γ ( s ) is c.m. as the Laplace transform of the PDF, i.e.( − n ˜ γ ( n ) ( s ) ≥
0. The derivative ˜ h (cid:48) ( s ) = 1 − β ˜ γ (cid:48) ( s ) isc.m. since ˜ γ ( s ) is c.m. Thus, due to the property men-tioned above ˜ n f ( s ) = ˜ f (˜ h ( s )) is c.m. as well and n f ( t ) isnon-negative due to the Bernstein theorem [62]. Appendix E: The derivation of the limit result forthe two state Markovian limit
We show here how one can recover the limit of a twostate Markovian system with constant rates of exchangeand a decay in one of the states (considered, for instance,in Van Kampen’s book [63], Chapter VII, section 5, ex-ercise (5.8)). This corresponds to the Markovian limit ofEqs. (1)-(2).The case of two states with exchange can be writtenas dn f ( t ) dt = − αn f ( t ) − βn f ( t ) + n t ( t ) /τ ∗ , (7) dn t ( t ) dt = βn f ( t ) − n t ( t ) /τ ∗ , (8)where n f ( t ) and n t ( t ) correspond to the free and trappedstates in our model and α , β and τ ∗ are positive con-stants.Let us treat Eq. (8) formally as a non-homogeneousordinary differential equation of the first order. Then itsformal solution reads n t ( t ) = (cid:90) t dt (cid:48) e − ( t − t (cid:48) ) /τ ∗ βn f ( t (cid:48) ) . (9)By substituting the last equation into (7) we get dn f ( t ) dt = − αn f ( t ) − βn f ( t ) + 1 τ ∗ (cid:90) t dt (cid:48) e − ( t − t (cid:48) ) /τ ∗ βn f ( t (cid:48) )(10)and, respectively, dn t ( t ) dt = βn f ( t ) − τ ∗ (cid:90) t dt (cid:48) e − ( t − t (cid:48) ) /τ ∗ βn f ( t (cid:48) ) . (11)We see that the system of equations (10)-(11) is a par-ticular case of our model Eqs. (1)-(2) with exponentialtrapping PDF γ ( t ) = τ ∗ exp ( − t/τ ∗ ). Thus, we have es-tablished a link of the non-Markovian model (1)-(2) andthe classical Markovian two-state equations [63]. Appendix F: Solution of Markovian kinetic rateequations
We consider equations (7), (8) with initial conditions n f ( t = 0) = N , n t ( t = 0) = 0. In the Laplace space s ˜ n f − N = − α ˜ n f − β ˜ n f + ˜ n t /τ ∗ , (12) s ˜ n t = β ˜ n f − ˜ n t /τ ∗ (13) or, ( s + α + β )˜ n f − ˜ n t /τ ∗ = N , (14) − β ˜ n f + ( s + 1 /τ ∗ )˜ n t = 0 . (15)Determinant Det = s + ( α + β + 1 /τ ∗ ) s + α/τ ∗ . (16)Solution reads˜ n f ( s ) = N s + 1 /τ ∗ s + ( α + β + 1 /τ ∗ ) s + α/τ ∗ , (17)˜ n t ( s ) = N βs + ( α + β + 1 /τ ∗ ) s + α/τ ∗ . (18)The denominator can be written as s + ( α + β + 1 /τ ∗ ) s + α/τ ∗ = ( s + | s | )( s + | s | ) , (19)where s and s are the roots of the quadratic equation s + ( α + β + 1 /τ ∗ ) s + α/τ ∗ = 0 , (20)that is s = − α + β + 1 /τ ∗ (cid:114)
14 ( α + β + 1 /τ ∗ ) − α/τ ∗ < ,s = − α + β + 1 /τ ∗ − (cid:114)
14 ( α + β + 1 /τ ∗ ) − α/τ ∗ < . After plugging (19) into (17),(18) we get˜ n f ( s ) = N s + 1 /τ ∗ ( s + | s | )( s + | s | ) = N | s | − | s | (cid:26) /τ ∗ − | s | s + | s | + | s | − /τ ∗ s + | s | (cid:27) , (21)˜ n t ( s ) = N β ( s + | s | )( s + | s | ) = N β | s | − | s | (cid:26) s + | s | − s + | s | (cid:27) . (22)Taking inverse Laplace transform n f ( t ) = N | s | − | s | (cid:110) (1 /τ ∗ − | s | ) e −| s | t +( | s | − /τ ∗ ) e −| s | t (cid:111) , (23) n t ( t ) = N β | s | − | s | (cid:110) e −| s | t − e −| s | t (cid:111) . (24) Remark 1.
It is possible to prove that | s | < /τ ∗ < | s | , thus both terms in braces of (23) are positive and,automatically, the term in braces of (24) is also positive.For example,1 /τ ∗ > | s | = α + β + 1 /τ ∗ − (cid:114)
14 ( α + β + 1 /τ ∗ ) − α/τ ∗ , (25)0 (cid:114)
14 ( α + β + 1 /τ ∗ ) − α/τ ∗ > α + β + 1 /τ ∗ − /τ ∗ . (26)If α + β < /τ ∗ , then the left side >
0, right side < α + β > /τ ∗ , then taking square of both sides,14 ( α + β + 1 /τ ∗ ) − α/τ ∗ >
14 ( α + β + 1 /τ ∗ ) − /τ ∗ ( α + β + 1 /τ ∗ ) + 1 /τ ∗ ⇒ > − β/τ ∗ . In the same way the inequality 1 /τ ∗ < | s | is proved. Remark 2. If α = 0, then s = 0, s = − ( β + 1 /τ ∗ ),then n f ( t ) and n t ( t ) reduce to n f ( t ) = N β + 1 /τ ∗ (cid:110) /τ ∗ + βe − ( β +1 /τ ∗ ) t (cid:111) , (27) n t ( t ) = N ββ + 1 /τ ∗ (cid:110) − e − ( β +1 /τ ∗ ) t (cid:111) . (28) Remark 3. Normalisation. ˜ n f ( s = 0) = N /α = const . This implies normalisation (cid:82) ∞ n f ( t ) dt = N /α what follows directly from Eq. (3). Appendix F: Derivation of the asymptotic power-lawfor n f ( t ) In order to derive Eq. (5) for the emission intensity inthe long-time limit we follow the reasoning from [64] and use an identity N α − ˜ n f ( s ) = (cid:90) ∞ (cid:2) − e − st (cid:3) n f ( t ) dt. (29)From the expansion of ˜ n f ( s ) at small s , Eq. (4), we inferthat n f ( t ) decays at long times as n f ( t (cid:29) τ ∗ ) (cid:39) Cτ µ ∗ t µ .We put the latter expression in the identity and get N α − ˜ n f ( s ) (cid:39) C (cid:90) ∞ T (cid:2) − e − st (cid:3) τ µ ∗ t µ dt. (30)Now we take derivatives from both sides with respect to s and using Tauberian theorem for Laplace transforms[62], we get − ˜ n (cid:48) f ( s ) (cid:39) C Γ(1 − µ ) τ µ ∗ s µ − . (31)After differentiating Eq. (4) one obtains the expressionfor C , C = µβN α Γ(1 − µ ) . (32) [1] E. L. Ivchenko, Optical spectroscopy of semiconductornanostructures (Alpha Science Int., Harrow, UK, 2005).[2] V. L. Colvin, M. C. Schlamp, and A. P. Alivisatos, Na-ture , 354 (1994).[3] V. I. Klimov, A. A. Mikhailovsky, S. Xu, A. Malko, J. A.Hollingsworth, C. A. Leatherdale, H. J. Eisler, and M. G.Bawendi, Science , 314 (2000).[4] A. Y. Nazzal, L. Qu, X. Peng, and M. Xiao, Nano Lett. , 819 (2003).[5] S. Ithurria and B. Dubertret, J. Am. Chem. Soc. ,16504 (2008).[6] M. Xu, T. Liang, M. Shi, and H. Chen, Chem. Reviews , 3766 (2013).[7] R. Benchamekh, N. A. Gippius, J. Even, M. O. Nestok-lon, J.-M. Jancu, S. Ithurria, B. Dubertret, A. L. Efros,and P. Voisin, Phys. Rev. B , 035307 (2014).[8] J. Q. Grim, S. Christodoulou, F. D. Stasio, R. Krahne,R. Cingolani, L. Manna, and I. Moreels, Nat. Nanotech. , 891 (2014).[9] A. W. Achtstein, A. Schliwa, A. Prudnikau, M. Hardzei,M. V. Artemyev, C. Thomsen, and U. Woggon, NanoLett. , 3151 (2012).[10] S. Ithurria, M. D. Tessier, B. Mahler, R. P. Lobo, B. Du-bertret, and A. L. Efros, Nature Mater. , 936 (2011).[11] B. Mahler, B. Nadal, C. Bouet, G. Patriarche, and B. Du-bertret, J. Am. Chem. Soc. , 18591 (2012). [12] M. D. Tessier, P. Spinicelli, D. Dupont, G. Patriarche,S. Ithurria, and B. Dubertret, Nano Lett. , 207 (2014).[13] L. Biadala, F. Liu, M. D. Tessier, D. R. Yakovlev, B. Du-bertret, and M. Bayer, Nano Lett. , 1134 (2014).[14] M. Pelton, S. Ithurria, R. Schaller, D. Dolzhnikov, andD. V. Talapin, Nano Lett. , 6158 (2012).[15] C. She, I. Fedin, D. S. Dolzhnikov, A. Demortiere, R. D.Schaller, M. Pelton, and D. V. Talapin, Nano Lett. ,2772 (2014).[16] B. Guzelturk, Y. Kelestemur, M. Olutas, S. Delikanli,and H. V. Schaller, Demir, ACS Nano , 6599 (2014).[17] M. J. Bruchez, M. Moronne, P. Gin, S. Weiss, and A. P.Alivisatos, Science , 2013 (1998).[18] L. Manna, D. J. Milliron, A. Meisel, E. C. Scher, andA. P. Alivisatos, Nature Materials , 382 (2003).[19] C. B. Murray, D. J. Norris, and M. G. Bawendi, J. Am.Chem. Soc. , 8706 (1993).[20] X. Michalet, F. F. Pinaud, L. A. Bentolila, J. M. Tsay,S. Doose, J. J. Li, G. Sundaresan, A. M. Wu, S. S. Gamb-hir, and S. Weiss, Science , 538 (2005).[21] A. M. Smirnov, A. D. Golinskaya, P. A. Kotin, S. G.Dorofeev, V. V. Palyulin, V. N. Mantsevich, and V. S.Dneprovskii, Journal of Luminescence , 29 (2019).[22] A. M. Smirnov, A. D. Golinskaya, P. A. Kotin, S. G.Dorofeev, E. V. Zharkova, V. V. Palyulin, V. N. Mant-sevich, and V. S. Dneprovskii, J. Phys. Chem. C , , 10934 (2016).[24] J. Owen, Science , 615 (2015).[25] M. A. Boles, D. Ling, T. Hyeon, and D. V. Talapin, NatMater. , 141 (2016).[26] J. Owen and L. Brus, J. Am. Chem. Soc. , 10939(2017).[27] L. Biadala, E. V. Shornikova, A. V. Rodina, D. R.Yakovlev, B. Siebers, T. Aubert, M. Nasilowski, Z. Hens,B. Dubertret, A. L. Efros, et al., Nature Nanotechnology , 569 (2017).[28] M. Nirmal, B. O. Dabbousi, M. G. Bawendi, J. J. Mack-lin, J. K. Trautman, T. D. Harris, and L. E. Brus, Nature , 802 (1996).[29] U. Banin, M. Bruchez, A. P. Alivisatos, T. Ha, S. Weiss,and D. S. Chemla, The Journal of Chemical Physics ,1195 (1999).[30] M. Kuno, D. P. Fromm, H. F. Hamann, A. Gallagher,and D. J. Nesbitt, The Journal of Chemical Physics ,1028 (2001).[31] F. D. Stefani, J. P. Hoogenboom, and E. Barkai, PhysicsToday , 34 (2009).[32] F. T. Rabouw, J. C. van der Bok, P. Spinicelli, B. Mahler,M. Nasilowski, S. Pedetti, B. Dubertret, and A. Van-maekelbergh, Nano Lett. , 2047 (2016).[33] E. V. Shornikova, L. Biadala, D. R. Yakovlev, V. F.Sapega, Y. G. Kusrayev, A. A. Mitioglu, M. V. Ballot-tin, P. C. M. Christianen, V. V. Belykh, M. V. Kochiev,et al., Nanoscale , 646 (2018).[34] M. Olutas, B. Guzelturk, Y. Kelestemur, A. Yeltik, S. De-likanli, and H. V. Demir, ACS Nano , 5041 (2015).[35] F. T. Rabouw, M. Kamp, R. J. A. van Dijk-Moes, D. R.Gamelin, A. F. Koenderink, A. Meijerink, and D. Van-maekelbergh, Nano Lett. , 7718 (2015).[36] L. T. Kunneman, J. M. Schins, S. Pedetti, H. Heuclin,F. C. Grozema, A. J. Houtepen, and B. Dubertret, NanoLett. , 7039 (2014).[37] O. Del Pozo-Zamudio, S. Schwarz, M. Sich, I. A. Aki-mov, M. Bayer, R. C. Schofield, E. A. Chekhovich, B. J.Robinson, N. D. Kay, O. Kolosov, et al., 2D Mater. ,035010 (2015).[38] Q. Li and T. Lian, Nano Lett. , 3152 (2017).[39] Q. Li, K. Wu, J. Chen, Z. Chen, J. R. McBride, andT. Lian, ACS Nano , 3843 (2016).[40] M. C. Weidman, A. J. Goodman, and W. A. Tisdale,Chem. Mater. , 5019 (2017).[41] M. Kulig, J. Zipfel, P. Nagler, S. Blanter, C. Sch¨uller,T. Korn, N. Paradiso, M. M. Glazov, and A. Chernikov,Phys. Rev. Lett. , 207401 (2018).[42] J. Zipfel, M. Kulig, R. Perea-Caus´ın, S. Brem, J. D.Ziegler, R. Rosati, T. Taniguchi, K. Watanabe, M. M.Glazov, E. Malic, et al., Phys. Rev. B , 115430 (2020). [43] M. M. Glazov, Phys. Rev. B , 045426 (2019).[44] M. Seitz, A. J. Magdaleno, N. Alcazar-Cano, M. Melen-dez, T. J. Lubbers, S. W. Walraven, S. Pakdel, E. Prada,R. Delgado-Buscalioni, and F. Prins, Nature Communi-cations , 2035 (2020).[45] H. Scher and E. W. Montroll, Phys. Rev. B , 2455(1975).[46] B. Saidzhonov, V. Kozlovsky, V. Zaytsev, andR. Vasiliev, Journal of Luminescence , 170 (2019).[47] T. Erdem and H. Demir, Nanophotonics , 74 (2016).[48] B. Saidzhonov, V. Zaytsev, M. Berekchiian, andR. Vasiliev, Journal of Luminescence , 117134 (2020).[49] A. Achtstein, O. Marquardt, R. Scott, M. Ibrahim,T. Riedl, A. Prudnikau, A. Antanovich, N. Owschimikow,J. Lindner, M. Artemyev, et al., ACS Nano , 9476(2018).[50] R. Scott, S. Kickh¨ofel, O. Schoeps, A. Antanovich,A. Prudnikau, A. Chuvilin, U. Woggon, M. Artemyev,and A. Achtstein, Phys. Chem. Chem. Phys. , 3197(2016).[51] P. Irkhin and I. Biaggio, Phys. Rev. Lett. , 017402(2011).[52] E. W. Montroll and G. H. Weiss, Journal of Mathemati-cal Physics , 167 (1965).[53] H. Scher and M. Lax, Phys. Rev. B , 4491 (1973).[54] B. D. Hughes, Random Walks and Random Environ-ments, Vol. 1, Random Walks (Clarendon Press, Oxford,UK, 1995).[55] J.-P. Bouchaud and A. Georges, Physics Reports ,127 (1990).[56] R. Metzler and J. Klafter, Physics Reports , 1 (2000).[57] T. Tiedje and A. Rose, Solid State Communications ,49 (1981).[58] J. M. Chambers, C. L. Mallows, and B. W. Stuck,Journal of the American Statistical Association , 340(1976).[59] GetData Graph Digitizer program , URL http://getdata-graph-digitizer.com .[60] B. V. Gnedenko and A. N. Kolmogorov,
Limit dis-tributions for sums of independent random variables (Addison-Wesley Mathematics Series, Cambridge, UK,1954).[61] G. Samorodnitsky and M. S. Taqqu,
Stable Non-Gaussian Random Processes: Stochastic Models with In-finite Variance (CRC Press, 1994).[62] W. Feller,
An Introduction to Probability Theory and ItsApplications, Vol. 2 (Edition, Wiley, Cambridge, UK,1991).[63] N. G. Van Kampen,
Stochastic Processes in Physics andChemistry (Elsevier, 2007).[64] S. Havlin and G. H. Weiss, Journal of Statistical Physics58