Kinetic Monte‐Carlo Simulation of Exciton Hopping: Urbach Tails in Gas‐Molecule Decorated MoSe2

Disorder parameters of gas‐molecule decorated monolayer MoSe 2 are quantitatively investigated. This material system is interesting because disorder may be introduced and removed at will by regulating the number of adsorbed gas molecules through laser annealing. These molecules electrostatically trap excitons leading to localized defect states, which are exponentially distributed in energy. Herein, experiments are described by kinetic Monte‐Carlo simulations, in summary enabling richer studies than within crystalline materials with a fixed degree of disorder. It is found that the surface coverage of the MoSe2 may reach up to one molecule per 2 nm2 and that the density of adsorbed molecules depends on the laser power by a power law.

Decoration of MoSe 2 by air molecules at the surface, for example, leads to local, electrostatic trapping of excitons and subsequent hopping relaxation dynamics. [8][9][10] Selenium is more electronegative than sulfur, and therefore MoSe 2 attracts gas molecules, which adsorb at the surface of the material. [11] One of the most interesting facts about this system is that this kind of disorder can be tuned, e.g., by laser illumination, which leads to desorption of the gas molecules. This can be proved by photoluminescence (PL) spectra: we explicitly showed that this effect is reversible by the (dis-)appearance of the gas-molecule localized exciton peak upon laser annealing. [8] Furthermore, partial annealing reduces the density of hopping sites N 0 and impacts the hopping dynamics, as described in our results. The quantitative description of these disorder effects is the main focus of this publication, which reveals interesting aspects of the physical origin of disorder.

Theory and Kinetic Monte-Carlo Simulation
Kinetic Monte-Carlo (KMC) simulation of exciton hopping by means of Miller-Abrahams transition rates [12,13] has been applied to many semiconductor materials with disorder, e.g., (GaIn)(NAs), [14] (Zn,Cd)Se/Zn(S,Se) quantum wells, [15] Ga(AsBi), [16] and for gas-molecule decorated MoSe 2 . [8] Using the approach by Baranovskii et al., [17] potential fluctuations are treated by generic sites, where each site is localized in a potential minimum. Charge carriers or excitons (or any stable, fermionic quasiparticles), which are located on a site i, may undergo a thermally assisted transition to another site j by the rate where ν 0 is the Debye frequency, α is the exciton localization radius, r ij is the intersite distance, E i,j are the eigenenergies, k B is the Boltzmann constant, and T is the temperature. The sites are distributed randomly in space according to their spatial density N 0 , which corresponds to the density of adsorbed molecules in our case. Apart from hopping, excitons can also recombine radiatively with an average radiative exciton lifetime τ 0 . At this point, radiative decay is considered, only. In addition, a nonradiative decay time τ nr can be considered in the simulation, but does not change the shape of the time-integrated PL spectrum. This effect just quenches the PL spectra and alters the PL dynamics, i.e., the net decay time τ ¼ ðτ À1 0 þ τ À1 nr Þ À1 describes the decay dynamics. As we investigate the steady-state, only, nonradiative effects can be neglected at this point.
The aforementioned theory contains two important assumptions: first, the exciton wave function does not change significantly between the initial and the final state. There is negligible electronic DOI: 10.1002/pssb.202100186 Disorder parameters of gas-molecule decorated monolayer MoSe 2 are quantitatively investigated. This material system is interesting because disorder may be introduced and removed at will by regulating the number of adsorbed gas molecules through laser annealing. These molecules electrostatically trap excitons leading to localized defect states, which are exponentially distributed in energy. Herein, experiments are described by kinetic Monte-Carlo simulations, in summary enabling richer studies than within crystalline materials with a fixed degree of disorder. It is found that the surface coverage of the MoSe 2 may reach up to one molecule per 2 nm 2 and that the density of adsorbed molecules depends on the laser power by a power law.
relaxation-in contrast to, e.g., organic molecules-where there is a significant amount of electronic relaxation due to different shapes of the orbitals. The Marcus theory accounts for this relaxation by the reorganization energy. [18] In the limit of vanishing reorganization energy, both theories, Marcus and Miller Abrahams theory, coincide. Second, the spatial dependence of the wave function far from its center is exponential and the overlap of the initial state and the final wave function is small enough such that the carrier-carrier transition can be treated in a perturbative manner. This assumption is introduced in the theory of Miller and Abrahams [12,13] In our case, both assumptions are presumably fulfilled, as the A-excitons at sites i and j show comparable character, both being in their excitonic ground state. Second, the most dominant optical response stems from A-excitons, whose Bohr radius is in the range of 1 nm, whereas the site distances are in the range of a few nm. Thus, the overlap between the wave functions of neighboring sites is small enough to be treated perturbatively and the decay of the wave function is exponential.
The energy distribution of the localized sites cannot directly be obtained from experiments, but only indirectly. In the past, energy distribution functions using a single characteristic (mean) energy have often described the system very well, an exception being Ga(AsBi). [16] The distribution function for lowdimensional systems [8,14] is usually exponential and for bulk systems and molecules, [19][20][21] it is Gaussian. In some rare cases, one can observe different disorder mechanisms resulting in more than one independent distribution functions. [16] For an exponential density of states (DOS) gðEÞ, the characteristic localization energy is identical to the Urbach energy E U [22] gðEÞ where E A is the energy of the A-exciton and the density is usually located energetically below the optical transition without disorder. This corresponds to localized states below the conduction band or above the valence band edge. Equation (1) constitutes a large system of master equations, which we solve by a KMC simulation. An exciton is set to an arbitrary starting site i and may either recombine or relax into another state j. The sites are distributed locally homogeneous and energetically according to a given distribution function (here: exponential). At each time step, hopping rates and the radiative decay rate are calculated and normalized. Pseudorandom numbers then select the respective process until the exciton recombines. After that, the simulation continues on a new lattice realization.
The hopping rates are directly obtained from Equation (1). However, radiative decay itself is another stochastic process, where the distribution of the exciton lifetime is usually an exponential function with the expectation value τ 0 . This means that discrete, random decay times t k are generated at each time step k from uniformly distributed random numbers ξ ∈ ½0, 1 according to the transformation law Its inverse, the decay rate at each time step, ν k ¼ t À1 k , is a possible decay channel at each hopping step. Another random number χ ∈ ½0, 1 then selects the respective process. This method is identical to the method of Bortz, Kalos, and Lebowitz (BKL). [23] One could have also implemented the first reaction method (FRM) [24] -both methods perform similarly and end up in the same results after averaging statistically. Thus, the BKL and the FRM method are equally useful in terms of speed as well as physical results.
The results of the theory scale with three dimensionless disorder parameters: ν 0 τ 0 is related to average number of hopping processes per exciton before recombination at a given temperature and site density. One can justify this argument by averaging The disorder parameters, ν 0 τ 0 , N 0 α d , and E U k B T , can be extracted directly from temperature-dependent PL spectra. For example, the temperature dependence of the PL maximum exhibits a typical S-shape, [16,17,25] which means a low-temperature red and higher temperature blueshift of the PL maximum with increasing temperature. From its shape, one can deduce the kind of disorder (Gaussian or exponential) and the value of E 0 assuming that E 0 is temperature independent. Figure 1 shows temperature-dependent PL spectra of MoSe 2 decorated with gas molecules (wet oxygen); schematically depicted in the inset. The PL spectra are obtained from continuous wave (cw) illumination. The A-exciton (A) and the trion (T) are indicated in this figure. The broad feature at lower energies is the disorder- Data: Temperature-dependent optical spectra from photoluminescence with cw excitation, normalized by to the intensity of the 6 K spectrum: the Stokes shift and the spectral width increase with temperature. A denotes the A-exciton, whereas T is the trion energy. The Urbach energy E U is depicted for the spectrum at 6 K.

Experiment
www.advancedsciencenews.com www.pss-b.com related peak, which emerges due to the adsorbed gas molecules on the MoSe 2 surface. [8,9] The spectra exhibit an incomplete S-shape, which means that E U k B is much larger than 100 K (cf., Imhof et al. [16] ). The physical reason for this is that the PL is quenched for temperatures above 100 K due to nonradiative decay channels. The Urbach energy E U can further be obtained from the exponential slope of the low-energy spectrum, which is visualized exemplarily for the 6 K data. The reason for this is that the PL is equivalent to the DOS in the low-energy part. Furthermore, black double arrows indicate the spectral width at 3/4 of the PL maximum, which we simply refer to as "width" or Δ 3=4 . We do not use the full width at half maximum (FWHM, Δ 1=2 ), as it cannot be obtained reliably for temperatures above 60 K-exemplarily indicated by the blue-dashed arrow at the spectrum at 100 K.

Extracting Material Parameters
To fit the experimental data with our model, we perform a numerical parameter study, where we trace the Stokes shift and the FWHM as a function of the hopping parameters at a temperature of 6 K in Figure 2. The axes are scaled in units of E U , as this is a universal parameter at a constant temperature; for the realization of the data, we chose an arbitrary E U of 11 meV. The value of E U , however, depends on the amount of disorder, which may increase with temperature for surfacesensitive materials as considered here. We observe an increase in the Stokes shift with both, ν 0 τ 0 and N 0 α 2 , where N 0 represents the density of adsorbed molecules in our situation. The FWHM approaches a peak value of about Δ 1=2 % 2.6E U (Δ 3=4 % 1.75E U ) for values of N 0 α 2 ≳ 0.1 (see Figure 2). This value is nearly temperature independent for temperatures below half of the critical temperature T < 0.5E U =k B . [16,17] The experimental spectra show that the ratio of Δ 3=4 =E U % 1.75: Δ 3=4 ranges up to 50 meV (at 80 K), while the corresponding energy E U % 30 meV, extracted from the low energy derivative of the same spectrum (see Figure 3).
Usually, E U is independent of the temperature, as the disorder is often related to lattice defects. Defect energies are rather robust at temperatures far below the melting point. [14][15][16] Here, we observe a significant temperature dependence: by keeping the parameters ν 0 τ 0 ¼ 10 6 and N 0 α 2 ¼ 0.1 constant and taking E U from the spectra, we obtain a nearly perfect agreement with experimental data (see Figure 3). The figure shows the temperature-dependent PL spectra (colored dots) and the simulation results (black dots). Figure 3c shows the extracted data, Stokes shift, and spectral width Δ 3=4 as a function of temperature and also indicates the Urbach energy. The relation between Urbach energy and the PL width becomes apparent, as both quantities are nearly linearly interdependent. The strong redshift of the experimental spectra with increasing temperature can be explained by a superposition of two effects: the redshift due to hopping toward the global energetic minimum at increased temperature becomes even more pronounced by increasing the Urbach energy E U . If temperature is increased above the critical temperature, the exponent E i ÀE j k B T becomes smaller than one and the relative energy landscape becomes very flat. This results in the observation that the emitted PL represents the localized DOS in that temperature regime and this would correspond to an effective blueshift toward the absorption edge. In our case, this critical temperature is about room temperature, which means that nonradiative effects are so strong that PL from localized states will not be visible anymore and this regime cannot be investigated in our experiment. The large size of the Urbach energy by adsorbed molecules, however, is not far-fetched: a recent estimation of the impact of local screening contrasts on the position of the A-exciton yields relative energetic shifts of about 20⋯30 meV in 2D materials. [9] Interestingly, lattice Figure 2. a,c) Stokes shift and b,d) spectral width as a function of the hopping parameters ν 0 τ 0 and N 0 α 2 . E U is fixed at 11 meV, but the results are independent of this specific value. The temperature is kept at 6 K.
www.advancedsciencenews.com www.pss-b.com disorder in a related system, Mo 0.3 W 0.7 S 2 , generates a value of 7 meV disorder potential, which is small compared with our aforementioned value. [26] This underlines that in 2D materials, surface effects can even have larger impact on the material properties than intrinsic defects. If one ignored the temperature dependence of the Urbach energy, the Stokes shift would be reproducible by assuming an extremely large Debye frequency ν 0 and the spectral width of the experiment could not be reproduced within simulations. This is shown in Figure S1, Supporting Information.
The reason for the increase in the Urbach energy E U with temperature could be that the disorder is induced by the electrostatic influence of the gas molecules at the surface, whose kinetic energy increases upon rising temperature. This effect of an increased Brownian motion on the disorder, however, needs to be confirmed by complementary, atomistic investigations.

Site Density
From an experimental point of view, the site density is an important piece of information, as it is proportional to the amount of adsorbed gas molecules. This information is assessed by our simulations and enables to judge the quantity of adsorbed gas molecules. For technical reasons, the results from Figure 4 and 5 are obtained from another sample than in the prior paragraph. However, the only quantity that changes between the samples is the Urbach energy E U . The other hopping parameters ν 0 τ 0 and N 0 α 2 remain constant as we show in Section C, Supporting Information. In consequence, we observe the same site density on both samples. An increased laser power in the experiment reduces the site density by evaporating the surface molecules. PL is measured at a low laser power of 0.1 μW, afterward. Figure 4 shows that the absolute value of the Stokes shift decreases with increasing laser power, which is due to the reduced site density. Furthermore, the inset of this figure shows that the Urbach energy, extracted from spectral width and the low-energy slope, is independent of the laser irradiation power. This means that the disorder strength itself remains and that there is no significant sample heating after the annealing step. Local heating during laser annealing is indeed possible, but heat dissipates rapidly in crystalline, low-dimensional systems.
The estimation of the density of available hopping sites N 0 is shown in Figure 5. At this point, the hopping parameter ν 0 τ 0 needs to be fixed. The parameter ν 0 is in the range of 1⋯100 THz, which is related to typical phonon frequencies, and τ 0 is in the range of ns for MoSe 2 . [27,28] This means that ν 0 τ 0 is in the range of 10 3 ⋯10 5 . We also included a large value of 10 7 in Figure 5 to be consistent with Figure 2. Combining the data of Figure 2 and 4, we estimate a range of the density of available hopping sites in the order of 0.1 ≤ N 0 α 2 ≤ 0.5.
To obtain a proper value for the site density N 0 , the exciton localization radius of the A-exciton α is required, which is related    www.advancedsciencenews.com www.pss-b.com to the Bohr radius a B (see discussion in Section D, Supporting Information). Theory values for the Bohr radius range from 0.5 nm up to about 1 nm in MoSe 2 in vacuum [28][29][30][31] (see Table 2, Supporting Information), slightly increasing upon environmental screening as, e.g., demonstrated for MoS 2 . [28] We chose the radius α ¼ a exp: B % 1nm to calculate the area per exciton, where a exp: B is the decay constant in the exponent of the wave function. This yields a site density range of about (0.1⋯0.5) hopping sites nm À2 in the low-power regime. One CO 2 , O 2 , and H 2 O molecule covers approximately 0.05 nm 2 of surface (we estimate this by the van der Waals radius of the atoms and their distances), which finally translates to 0.5⋯2.5% surface coverage.
The parameter N 0 α 2 decreases for an increasing laser power, which reveals a power law behavior with an exponent of À0.12. This means that the molecule desorption is a statistical process similar to a phase transition triggered by laser irradiation. The evaporation of gas molecules is due to local heating of the material by the laser resulting in this phase transition.

Conclusion
We have shown that KMC simulations of exciton hopping reproduce experimental data of gas-molecule decorated MoSe 2 quantitatively with a consistent parameter set. In contrast to many inorganic semiconductors, the Urbach energy increases with rising temperature in this material system, which may be due to an increased Brownian motion of the gas molecules at the surface. At this point, atomic studies of gas-molecule decorated MoSe 2 are required to gain a thorough understanding of this effect. We further quantified boundaries for the gas-molecule density absorbed, which is about (0.1⋯0.5) nm À2 . Finally, we find a power law behavior of the gas molecule density as a function of the laser irradiation power, which indicates a statistical phase transition between adsorbed and desorbed phase.

Supporting Information
Supporting Information is available from the Wiley Online Library or from the author.