The influence of Coulomb interaction screening on the excitons in disordered two-dimensional insulators

We study the joint effect of disorder and Coulomb interaction screening on the exciton spectra in two-dimensional (2D) structures. These can be van der Waals structures or heterostructures of organic (polymeric) semiconductors as well as inorganic substances like transition metal dichalcogenides. We consider 2D screened hydrogenic problem with Rytova–Keldysh interaction by means of so-called fractional Scrödinger equation. Our main finding is that above synergy between screening and disorder either destroys the exciton (strong screening) or promote the creation of a bound state, leading to its collapse in the extreme case. Our second finding is energy levels crossing, i.e. the degeneracy (with respect to index \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu $$\end{document}μ) of the exciton eigenenergies at certain discrete value of screening radius. Latter effects may also be related to the quantum manifestations of chaotic exciton behavior in above 2D semiconductor structures. Hence, they should be considered in device applications, where the interplay between dielectric screening and disorder is important.


Disorder and heavy-tailed Lévy distributions
The most important question of the present consideration is how fractional derivatives, generating heavy-tailed (i.e. decaying slower than Gaussian, typically in a power-law fashion x −1−µ , where x is any physical quantity like energy or potential barrier height) Lévy distributions, are related to the disorder. The common wisdom is that disorder is a lack of regularity. In this case, the above physical quantities are not under precise control so that the properties of such disordered systems are described in terms of distribution functions. The goal now is to predict global properties shared by almost all such systems, i.e. to acquire knowledge of universal features independent of the precise realization of the disorder. If the disorder is weak (i.e. a crystal has a small number of noninteracting defects and/or impurities), its properties are well described by the Gaussian distribution function. As this distribution falls off rapidly, its width is usually modest so that uncommon, "highly disordered" configurations have (very) small statistical weights and do not contribute to the observable properties of such (ordered or weakly disordered) systems. In a highly disordered material, atoms are not arranged in crystalline periodic patterns but appear in random positions. This means that the actual statistical weight of the above "highly disordered" configurations grows, sometimes enormously. This large statistical weight is indeed described by the above non-Gaussian, heavy-tailed distributions. As we discussed above, in an amorphous substance, the electronic states are no more Bloch functions due to severe violation of the translational symmetry. The simplest model that describes the electronic states in the highly disordered matter was introduced by Anderson 32 and leads to the famous Anderson localization phenomenon.
Namely, the main theorem of Anderson 32 states that if the breadth of energy distribution W is around mean value V of the interaction potential between disorder constituents (i.e. the particles like electrons, excitons, spins) located at the sites of some lattice, the transport in such system is absent. Note that to this end, the specific form of the above energy distribution is not given, the only condition is that it should be sufficiently broad. In other words, in a system with the aforementioned strong disorder, having a wide enough distribution of its physical characteristics, all states are strongly localized. In our context, the above "disordered lattice" means considerable substance amorphization or other types of strong disorder, when defects or impurities, having large concentration, start to strongly interact with each other. Such a strong disorder implies that the above distribution function Scientific Reports | (2021) 11:11956 | https://doi.org/10.1038/s41598-021-91414-w www.nature.com/scientificreports/ of the energy or any other physical characteristic is strongly non-Gaussian. At this point, it is reasonable to assume that this distribution is Lévy one. The next observation is that the Eq. (1) of seminal Anderson's paper 32 is very similar (with respect to substitution t → it , i is an imaginary unit) to Langevin equation, governing the stochastic dynamics of the (in our case two-dimensional) coordinate r(t) of a particle Here m is a particle mass and γ is a friction coefficient. Also, z(t) is a noise function, responsible for stochastic behavior. Once more, we assume that z(t) obeys Lévy statistics with probability distribution, most conveniently defined by its Fourier image or characteristic function f (k) = exp(−σ µ k µ /µ) , where k ≡ |k| and 0 < µ < 2 is Lévy index, see 22,23,33 and references therein. The case µ = 2 corresponds to Gaussian distribution with variance σ . This characteristic function is valid for any space dimensionality. In 2D it generates following explicit expression where J 0 (x) is Bessel function 34 . Note that the distribution (2) at µ < 2 have an infinite variance as the corresponding integral becomes divergent. For such heavy-tail distribution to have higher moments, we need some external potential (like in Fokker-Planck or Schrödinger equation), which "tames" the corresponding Lévy flight (common name for the processes, described by distribution (2)), see 35 and references therein. Both Eq. (1) of the paper 32 and Langevin equation (1) are stochastic differential equations. The probability density p(r, t) for the Eq.
(1) is given by the (fractional for µ < 2 ) Fokker-Planck equation Here |�| µ/2 is the fractional Laplacian, which in two dimensions reads Here Ŵ(x) is Ŵ-function 34 , r is a two-dimensional vector. For µ = 2 the operator (4) gives ordinary 2D Laplacian 36 . The details of the derivation of fractional Fokker-Planck equation (2) from Langevin equation with Lévy noise (1) can be found, for instance, in Ref. 37 . Such derivation reduces to the "extraction" of probability density from the noise characteristics. This "extraction" procedure is also present both in the derivation of ordinary Schrödinger equation from Feynman path integral 38 and in the derivation of fractional one by Laskin using the same construction but with Lévy measure 26 . Moreover, the free versions (at potential V = 0 ) of both equations can be reduced to each other by obvious transformations t → it and D = 2 /(2m) . This means that we can safely assume that the substitution of the underlying Gaussian distribution by slowly decaying, heavytailed one in the Schrödinger equation is equivalent to phenomenological account for (strong) disorder. In this case, the Lévy index µ serves as an indicator of the degree of disorder. Note that such suppositions have been made in Ref. 39 in the context of spectral narrowing of nuclear magnetic resonance lineshape. The seminal paper of Sher and Montroll 40 is dealing actually with the same situation. Here we note one more pioneering work 41 , considering the mesoscopic phenomenon of conductance of 1D quantum wire with the strong disorder. Namely, the randomness in the scatterers ensemble had Lévy distribution. It had been shown in Ref. 41 that the conductance statistics in this case is governed primarily by the Lévy index (the scaling exponent in the notations of 41 ), which also signifies the degree of disorder.
The pictorial demonstration of this fact is shown in Fig. 2, where the radial part ψ 00 (r) of the ground state wave function of the 2D screened fractional hydrogen atom is reported. The main panel of Fig. 2 shows the progressive localization of ψ 00 (r) at Lévy index diminishing. At the same time, the inset shows the width W of initial Lévy distribution (2). The comparison of the inset and main panel of Fig. 2 demonstrates that as initial width W increases, the wave function becomes more localized, tending to Dirac δ function as µ → 1 . This shows that for our system the conditions of Anderson theorem 32 is fulfilled, which supports our assumption that the introduction of fractional Laplacian in the Schrödinger equation effectively describes disorder. Note, that the breadth of 1D Lévy distribution, realizing for scalar quantities like energy, behaves qualitatively similar to the inset in Fig. 2. Moreover, our preliminary studies of continuous spectrum in the problem under consideration show that the corresponding wave functions start to localize at µ ≈ 1.9 , which is also related to the disorder. This means that as µ deviates from 2, the lower bound of the continuous spectrum starts to shift so that at µ → 1 only the small portion of continuous spectrum (if any at all) may remain. The answer to this question is so far unclear as we are dealing with the 2D system. It is well known (see, e.g. 42 ) that in the 1D system at strong disorder all states are localized. This problem, however, is out of the frame of present consideration and will be published elsewhere. (2)

Theoretical approach
To model the problem of an exciton in 2D disordered dielectric, here we consider 2D fractional (mimicking disorder) hydrogenic problem, where the exciton "resides" in the xy plane. Geometry of our problem is shown in Fig. 1. The question about the form of interaction potential U(r) in the structure, shown in Fig. 1, is usually thought of as a solution of Poisson equation (see, e.g. 43 ) �U(r) = −4πρ(r) , where ρ(r) is a total charge density at the point r , which is responsible for screening of the Coulomb interaction. The explicit form of the screened interaction potential U(r) had been found by Rytova 28 and subsequently by Keldysh 29 . It reads where r = |r| , β = e 2 /ε (e is electronic charge, ε is the dielectric permittivity of middle layer in Fig. 1), www.nature.com/scientificreports/ is the screening radius, see Fig. 1. Here H 0 (z) and Y 0 (z) are Struve and Neumann functions respectively 34 . At zero screening r 0 → 0 , the argument z = r/r 0 of the functions H 0 (z) and Y 0 (z) tends to infinity so that 34 and which corresponds to the ordinary 3D Coulomb interaction. As demonstrated in 43 , in 2D layers of finite thickness, the solution of the Poisson equation has the above 1/r form (rather than well-known logarithmic law, inherent in 2D systems 44 ) due to the presence of a hole (which, say, is our exciton "component") at the origin. The case of fractional 2D exciton in the layers of finite thickness with unscreened interaction (9) has been considered in Ref. 45 . At r 0 → ∞ , the argument of the functions in (6) z = r/r 0 → 0 . The asymptotics reads 34 where γ ≈ 0.5772 is the Euler's constant 34 . It is easy to see that this asymptotics corresponds to ordinary 2D Coulomb interaction in the form ln r , signifying the Green's function of 2D Poisson equation 44 . In other words, the Rytova-Keldysh interaction (6) gives the transition from ordinary 2D Coulomb interaction ln r in the case of infinite screening r 0 → ∞ to the case of 3D one ("ordinary" Coulomb ∼ 1/r ) at zero screening r 0 → 0.
Having potential (6), we can write down the fractional Scrödinger equation for our 2D fractional hydrogenic problem. It reads where U(r) is defined by Eq. (6) and A µ by (5). The Eq. (11) is a (singular as the integral in it exists in the sense of Cauchy principal value only) two-dimensional integral equation with fractional Laplacian having Lévy index µ . Latter index, in a sense, defines how far the fractional Laplacian deviates from ordinary one so that the eigenvalue problem (11) deflects from conventional 2D hydrogenic one 46 . Note that if we define the fractional Laplacian (4) in 3D (1D) space, the Eq. (11) with properly tailored potential (that is to say, 3D or 1D Coulomb interaction) will give 3D (1D) fractional hydrogenic problem.
Indices n and m denote, respectively, the principal and orbital quantum numbers, which are different for any specific µ value. Below we shall see that for µ < 2 , the so-called Coulomb degeneracy is lifted and that is the reason why the eigenenergy E has now two subscripts. To be more specific, it can be shown that in the fractional case, the "accidental" Coulomb degeneracy (stemming from Runge-Lenz vector conservation for µ = 2 47,48 ) is lifted so that the eigenenergy starts to depend on the orbital quantum number. Note, that this feature does not preclude the rotational (around z or k z axis in our case) invariance of the Schrödinger equation so that its solutions can still be separated on the radial and angular parts.
Here we use modified (for the fractional case µ < 2 ) Rydberg units 26 , i.e. we measure the energy E and coordinates r in the units respectively. Parameter β = e 2 /ε , while D µ is a mass term 26 . At µ = 2 units (12) convert into standard Rydberg units. Note that at µ = 1 both E 01 and r 01 in Eq. (12) are divergent. This is one more manifestation of the fact, that discrete spectrum of 2D quantum fractional hydrogenic problem exists at µ > 1 only 45 .
The Eq. (11) is indeed an integral equation, sometimes called pseudo-differential one 33 . Moreover, the integral (4) exists only in the sense of Cauchy principal value 33,36 , which is a source of additional difficulties in such equations solutions. To bypass this, here we transit to k = (k x , k y ) space as the operator (4) becomes −k µ , where k ≡ k 2 x + k 2 y . Although the second term in (11) converts to the integral in momentum space, it proves to be much easier to handle then initial one (4). This becomes clear if we consider the Fourier image of the potential U(r), which is much simpler, then initial expression (6): The particularly simple expression (13) permits to do the angular integration exactly, i.e. without expansion over spherical harmonics, see 49 and references therein. Latter fact permits us to solve the corresponding fractional Schrödinger equation by the method similar to the case without screening 45 . where ϕ and ϕ ′ are asimuth angles of the vectors k and k ′ respectively. Here, generalizing the case of ordinary ( µ = 2 ) 2D hydrogen atoms 49 , we denote It can be seen that at µ = 2 and r 0 = 0 , Eq. (15) gives the well-known momentum space representation of 2D hydrogenic problem with µ = 2 and unscreened "sheet" Coulomb interaction (9) 49 .
For any central (i.e. spherically symmetric) force potential (which is also the case for potential (13)), the angular and radial variables (in k space in our case) can be separated so that we can look for a solution of the Eq. (15) in the form As usual in hydrogenic problems, radial functions ψ nmµ (k) are real. The normalization condition for these functions read

Substitution of ansats (18) into Eq. (15) yields
Here for clarity, we suppress the subscripts n, µ. Note that in "ordinary" (non-fractional) 3D quantum hydrogenic problem, the corresponding Scrödinger equation is usually solved by stereographic projection method, which is due to Fock 48 (see also Ref. 49 for 2D case and Ref. 50 for multidimensional case). In this method, the problem can be solved exactly using spherical harmonics expansion [48][49][50] . It can be shown that in 2D screened case (13), such stereographic projection is impossible. The same is true for fractional case µ < 2 . It turns out, however, that the simple form of Fourier image of Rytova-Keldysh potential (13) admits the possibility of exact analytical calculation of the integrals I m (k, k ′ ) (21) for each specific m. This permits to easily solve the problem numerically as it now reduces to the effective one-dimensional integral equation 51 . The integral for m = 0 yields where �(m, n) is a complete elliptic integral of the third kind 34 . The asymptotics of the integral I 0 (k, k ′ ) (23) for unscreened case r 0 = 0 has the form 45 where K(m) is a complete elliptic integral of the first kind 34 . Note that at k = k ′ the argument of the complete elliptic integral in (24) equals to 1, which implies that K is divergent. This reflects the spherical harmonic series divergence at k = k ′ for unscreened case. In the screened case (23) the situation is more complicated as function �(m, n) is divergent at both m = 1 and n = 1 . However, as it can be shown, these divergences are well compensated by other parts of corresponding integrands. We note also, that function �(m, n) has different  34 for details. This fact should be taken into account in numerical calculations. The expressions for the integrals I m (k, k ′ ) at arbitrary m can also be obtained in closed form, although the expressions for higher m's become more and more cumbersome. After lengthy calculations, we arrive at the following form for I 1 (k, k ′ ) The set of Eq. (20) define (for each specific m) the spectrum of fractional 2D screened hydrogen atom for that particular m and all n ≥ m . To be specific, the equation for m = 0 (see (20) and (23)) defines the wave functions ψ 00µ (ground state), ψ 10µ , ψ 20µ etc. The equation for m = 1 determines ψ 11µ , ψ 21µ , etc. The same is true for higher m.
Equation (20) is the main theoretical result of the present consideration. As we shall see below, it is especially useful in our case of screened Coulomb interaction as it reduces the problem solution to the 1D integral equation, which numerical solutions behave well, see below.

Solution of the spectral problem
Let us pass to numerical solution of the set (20) for each specific m. The same for m = 1 has the form where the parameters ξ , n and m are defined by (27).
It is instructive to consider the known case of "ordinary" (i.e. with normal Laplacian in the Schrödinger equation) 2D hydrogen atom with unscreened Coulomb interaction. Such spectral problem is well studied (see 46 for coordinate space and 49 for momentum space) and its energy spectrum in our Rydberg units reads Figure 3 shows the ground state energy (proportional to exciton binding energy) as a function of Lévy index µ (a) and dimensionless screening radius ξ (b). Similar to the case without screening 45 , the discrete spectrum in our problem exists at 1 < µ < 2 . Also, at µ = 2 and ξ = 0 (case of "ordinary" 2D quantum mechanical hydrogen atom) the exciton ground state energy equals to -4, which follows from Eq. (29). It can be shown that at ξ = 0 and µ = 2 our numerical calculation gives correct values for the entire spectrum (29) within 1% accuracy. This fact can be thought of as a consistency check for our numerical solution.
It is seen from Fig. 3a that for relatively small screenings 0 < ξ < 1.5 the ground state energy goes to minus infinity at µ = 1 . At the same time at ξ = 1.5 , this energy goes to zero which means that strong screening "kills" the bound state between electron and hole in an exciton. Our analysis shows that this is the consequence of the competition between disorder (defined by fractional Laplacian) and screening effects. Namely, while the former deepens the "potential well", where the bound state exists (say, it promotes the bound state), the latter (at large screenings) makes this well become shallow so that the bound state becomes progressively "less bound" www.nature.com/scientificreports/ as ξ increases. It can be shown that at ξ ≤ 1.42 the ground state energy goes to minus infinity at µ = 1 (weak screening, disorder prevails, making the exciton particularly stable), while at ξ > 1.42 (strong screening) the situation is opposite, which is shown by the curve for ξ = 1.5 . At ξ > 1.5 the ground state energy in the entire domain 1 < µ < 2 tends to zero. The second important effect is energy levels crossing 52 , which is seen in Fig. 3b. The crossing occurs around ξ = 1.18 for all admissible 1 < µ < 2 . That is to say, in the point ξ = 1.18 we have the infinitely degenerate (for all continuously varying 1 < µ < 2 ) ground state energy of an exciton in a 2D system. As the ordered case corresponds to the only µ = 2 , we can assert that the energy level crossing in our system occurs as a result of the synergy between disorder (mimicked by fractional Laplacian with Lévy index µ ) and screening. Note that there is no energy levels crossing without screening in the disordered 2D excitons 45 . Figure 4 portrays the energies of the first excited state, corresponding to n = 1 and m = 0, 1 . The main effect here is lifting of orbital degeneracy due to nonconservation of the Runge-Lenz vector 47 . To be specific, if for µ = 2 and ξ = 0 , the energy is independent of orbital quantum number m (see Eq. (29)), this is not true for µ < 2 (disorder) and ξ = 0 (finite screening). In other words, in our problem, the Runge-Lenz vector nonconservation comes from both disorder and finite screening effects. This fact should be taken into account while designing optoelectronic devices based on disordered heterostructures, interfaces and other 2D semiconductor structures. Note, that nonconservation of Runge-Lenz vector does not mean the nonconservation of z -projection of angular momentum, i.e. central symmetry of the problem. Latter conservation law is reflected in the independence of energy of the sign of orbital index m, i.e E n,m = E n,−m ; it is lifted, when time inversion symmetry is broken by, e.g., external magnetic field 47 .
The energy level crossing is also present for the states with n = 1 . They occur both for different ξ at µ ≈ 1.95 (Fig. 4a) and different µ at ξ ≈ 0.3 for energy level E 10 as well as ξ ≈ 1.25 for energy level E 1,±1 , Fig. 4b. The electron-hole bound state destruction by strong screening is also present in the first excited state, although the picture is more diverse (then in the ground state) as now orbital index m comes into play. Namely, if the energies E 10 are less susceptible to screening (they go to zero at µ = 1 for quite strong screening ξ = 1.5 ), the levels E 1,±1 go to zero already at ξ = 0.5 , showing higher sensitivity to screening effects. This feature becomes especially notable at strong disorder around µ = 1 , where (for instance for small screenings ξ = 0.5 ) the energy E 10 goes to zero, while E 1,±1 -to minus infinity. This means that at strong disorder and finite screening, our 2D excitons are more stable for higher values of orbital quantum number m. Our numerical analysis shows that this effect realizes for higher excited states also. Of course, for higher n this effect is even more diversified as now we have all m < n at our disposal. The latter feature is also important for the proper functioning of different optoelectronic devices, based on disordered 2D structures. It can also be shown, that energy level crossing (having again much more crossing points than for n = 1 ) occurs also for higher n > 1. www.nature.com/scientificreports/ One more interesting feature is seen in Fig. 4b. Namely, while E 10 grows monotonously at screening radius ξ increase, the energy E 1,±1 (ξ ) is notably nonmonotonous having minimum around ξ = 0.5 , which depends on Lévy index µ . The higher is µ deviation from 2, defining the degree of disorder, the deeper is minimum. In other words, we see that the branch E 1,±1 has lower energy than that of exciton in the ordered 2D model without screening. It can be shown that this feature occurs also for E 2,±2 , E 3,±3 and possibly higher excited states. This shows that the synergy between disorder and (nonlocal 27 ) screening in two dimensions minimizes the energy (as compared to the ordered unscreened case) for some values of orbital quantum number m. This effect is also very important for the devices, using transitions (like dipole ones between m and m ± 1 ) between levels with different m's.
The remark about the validity of our fractional hydrogenic model is in place here. We recollect Laskin's construction of fractional quantum mechanics 26 , which uses Feynman path integration but with Lévy measure. Since the tails of Lévy distributions are much "heavier" than Gaussian, the corresponding distribution function (which in resulting Schrödinger equation yields a wave function) spans over many lattice constants. In the preceding discussion, we have shown that one of the sources of Lévy (rather than Gaussian) distributions is strong disorder, i.e. the situation when almost all lattice sites are occupied by impurities and the lattice by itself is highly distorted. In this case, Laskin's construction generates the "enveloping function" giving a good average description of such systems. In that sense our model, based on fractional quantum mechanics, presents the result of "smearing" the microscopic disorder in a 2D semiconductor, representing it as some effective medium. The price which we pay for that is the substitution of the ordinary Laplacian by its fractional counterpart in the corresponding Hamiltonian. Note, that the present approach is complimentary to microscopic one 53,54 , where the explicit averaging over disorder had been performed within the so-called fluctuating field method. There, for the case of strong disorder, the distribution function also appeared to be non-Gaussian. However, in the above microscopic treatment, additional suppositions about averaging over impurity configurations had been made, see Refs. 53,54 for details. At the same time, present construction, being phenomenological by its nature, does not make any assumption about specific properties of (strong) disorder constituents.
The smallest distance between the above disorder constituents in a crystal lattice is lattice constant, which is around 5 Å. This implies that our model is valid at the distances, where the lattice discreteness is unimportant (which is true at strong disorder), i.e. around 50 Å, which is 10 lattice constants. At such distances, due to disorder averaging, also the Anderson localization 32 occurs. Likewise, in our model, at µ ≈ 1.2 , the disorder yields so much localization (exciton collapse), that our model breaks down. Note that the Rytova-Keldysh model also does not describe screening correctly at distances around lattice constant. This is because this model utilizes a   Fig. 3 but for the first excited state E 10 and E 1,±1 . Lifting of orbital degeneracy, related to Runge-Lenz vector conservation at µ = 2 , is seen for entire µ and ξ domains (except one point ξ = 0 , µ = 2 , at which E 00 = E 1,±1 = 0.4444 (29)), where E 00 = E 1,±1 . The curves E 00 and E 1,±1 , corresponding to the same ξ (panel (a)) or µ (panel (b)), are coded by the same color (black for ξ = 0 , red for ξ = 0.5 , green for ξ = 1 and blue for ξ = 1.5 in panel (a) as well as black for µ = 2 , red for µ = 1.9 and blue for µ = 1.5 in panel (b)) and also by figures near curves. As dependencies of E 00 and E 1,±1 on screening radius ξ are nonmonotonous (see panel (b)), we show explicitly which curves correspond to E 00 and which to E 1,±1 in both panels. The positions of level crossing points are also shown. www.nature.com/scientificreports/ continuous approach to the charges interaction. In this approach, the notion of dielectric permittivity is used, which is macroscopic (or at least mesoscopic such that lattice discreteness does not enter) in its nature. To be specific, in order to describe a substance in terms of its permittivities (dielectric, magnetic, etc), the effective averaging over underlying lattice discreteness should be performed 55 . In momentum space, this corresponds to small wave vectors as compared to the size of the Brillouin zone. In real space, this length-scale is different for different substances and by order of magnitude corresponds to screening radius 43 . For instance, for layered transition metals dichalcogenides like WS 2 such radius is around 40Å 56 , i.e. once more about 10 lattice constants. Our analysis shows that the qualitative shape of the exciton wave functions is different for strong (when corresponding energy goes to zero at µ = 1 signifying the exciton ionization, i.e. bound state destruction, see Figs. 3,4) and weak screenings, where, on the contrary, exciton collapses as its energy tends to minus infinity. In both screening regimes, the wave functions are sensitive to the Lévy index µ . Namely, for weak screening, as µ diminishes from 2 to 1, the wave functions become progressively less localized in k space. Also, their amplitudes diminish as µ → 1 so that the wave functions become almost delocalized with zero amplitude. This situation corresponds to "extra strong" localization in coordinate space, where ψ(r) degenerates into (almost) Dirac δ -function, reflecting the exciton collapse, i.e. "falling" of the electron on the hole in it. The situation is opposite for strong screening, where in momentum space the wave function acquires Dirac δ function shape as µ → 1 , which in coordinate space correspond to the almost delocalized wave function, corresponding to exciton ionization. Hence, the interplay between screening effects and disorder (mimicked in our approach by fractional Laplacian with Lévy index µ ) may lead either to exciton ionization or its collapse for the strong disorder. That is to say, the synergy between screening and strong disorder (like substance amorphization, which is common for semiconductor interfaces, where high mechanical stresses occur) may destroy (either ionize or collapse them, depending on the relation between screening radius ξ , Lévy index µ and orbital quantum number m, see above) the excitons, which may preclude the functionality of the devices like solar cells and/or light emitting diodes.

Discussion: relation to experiment
The above solution permits calculation of the observable characteristics of the excitons. One of the important characteristics is the exciton binding energy. This quantity is proportional to the exciton ground state energy and defines the work needed to remove a bound electron to infinity. The latter quantity has been already calculated and is reported in Fig. 3. One more important physical characteristic is the exciton localization radius. This parameter is a mean value of exciton radius-vector r in the ground state n = m = 0 To obtain the expression (30), we have integrated over the angle ϕ . The results of numerical calculations of r are reported in Fig. 5. It is seen that while at screening radii ξ < 1.5 the exciton localization radius decays monotonously as µ (mimicking the degree of disorder in our problem) approaches 1, at ξ = 1.5 it starts to grow signifying the exciton ionization at µ = 1 . This is a reflection of the fact, that the joint action of strong screening and disorder destroys the bound state in an exciton, leading to its ionization at µ = 1 . This effect has already been shown in Fig. 3, where ground state energy at ξ = 1.5 and µ = 1 was zero. Note that at the same time, at weak screening, the disorder strengthens the bound state (this may be thought of as a kind of Anderson localization 32 ), leading to exciton collapse at µ = 1.
To understand better the physical situation with exciton collapse, we recollect that the bound (localized) state in quantum mechanics corresponds to negative energy, E < 0 47 . Our numerical calculations of the exciton (30) r ≡r 00µ = 2π ∞ 0 r 2 ψ 2 00µ (r)dr. www.nature.com/scientificreports/ localization radius in Fig. 5 show that the latter quantity diminishes as µ → 1 , giving exactly zero (say, "infinitely localized state" of zero spatial extension) at µ = 1 and screening radius ξ = 0 . So, the phenomenon of exciton collapse, corresponding to an "infinitely localized state" must have infinitely large negative energy. This is actually seen in Fig. 3. Our analysis for higher excited states with n > 1 shows the same tendency. Most probably, in real 2D structures having different kinds of disorder (say, defects and/or impurities at the interface), the exciton may be trapped and recombined before its actual collapse. The experimental and theoretical (inclusion of other sources of disorder in our model) studies of this interesting question should be carried out. The inset of Fig. 5 portrays the dependence r(ξ ) at two fixed values of Lévy index µ : µ = 2 (ordered case) and µ = 1.5 . In both cases, exciton localization radius grows with ξ in accord with Figs. 3 and 4, where it is shown that screening breaks bound state between electron and hole in an exciton. Also, in disordered case ( µ = 1.5 ) the exciton localization radius is smaller than that for µ = 2 for the same ξ . This shows that the synergy between disorder and nonlocal screening (peculiar to 2D case 27 ) stabilizes the exciton. Our estimations show that for typical (not very strong, when exciton collapses) strength of disorder µ = 1.5 , the screening radius r 0 , corresponding to dimensionless ξ =1.5 is 50 Å, which is in qualitative agreement with values, known from literature sources 30 Å < r 0 < 80 Å 57 , see also below.
The recent discovery of high photovoltaic efficiency in organic -inorganic halide perovskites like methylammonium lead iodide ( CH 3 NH 3 PbI 3 ) 57-59 requires knowledge of how the disorder influences their excitonic properties, which are responsible for solar to electric power conversion in them. The exciton properties of halide perovskites are still the subject of intense research and active debate, see, e.g. refs. 60,61 . Along with aforementioned halide perovskites, other good candidates for optoelectronic applications are 2D transition metal dichalcogenides (TMD) like MoS 2 , WS 2 and WSe 2 57 . In these substances, the excitons have high binding energy, which makes them extremely thermally stable. However, the influence of defects and other disorder on these excitons properties (like their relaxation and decoherence mechanisms) is still unknown to a larger extent. In the physical units, the typical exciton binding energy in 2D dichalcogenide semiconductor WS 2 is around 0.32 eV, which is much more than that (around 0.05 eV) in ordinary 3D materials 56 . This makes 2D TMD materials to be good candidates for optoelectronic and photonic devices at room temperatures. The in-plane exciton localization radius in them is of the order of 5 nm 56,57 . The numerical calculations with potential (6), performed in Ref. 56 , show that such Wannier-Mott like exciton can be well realized in atomically thin WS 2 (and other TMD 62 ) monolayer with screening radius r 0 ≈ 75 Å. The other sources (see 57 and references therein) list r 0 to be in the range from 30 to 80 Å, obtained from ab initio calculations.
Discussed physical properties of an exciton in 2D materials, related to the interplay between 2D nonlocal screening and disorder can play a role in multiexciton configurations. For example, they can be relevant for interaction between distant 2D excitons in the above interfaces or semiconductor-based ultrathin films. Indeed, the exciton radii (around 5 nm, which is of the order of 10 lattice constants) presented in Fig. 5 are characterized by dipole moments er , enhancing intrinsic electric fields of the excitons and their interactions. On the other hand, these fields will be screened nonlocally so that many defects and impurities will fall in the span of exciton radius. As we have mentioned above, in the semiconductor structures with different degrees of disorder (different Levy indices µ in our formalism) such a random screened exciton-exciton interaction may lead either to their ionization (high screenings) or to the collapse (low screening) at µ = 1 . This for sure will have a detrimental effect on the optoelectronic and/or spintronic device functionalities.

Conclusions
The message of the present paper is that the synergy between the nonlocal screening of 2D Coulomb interaction and disorder in semiconducting (generally speaking dielectric) surfaces, interfaces, thin films, and multilayers has novel properties, which do not occur either in 2D unscreened ordered case or in 3D one. Our main supposition here is that Laskin's construction of path integrals with Lévy measure 26 is equivalent to "extraction" of probability density function from fractional Langevin equation and, in turn, to the assumptions made in seminal Andeson paper 32 . This (along with the fact that as the width of initial distribution grows, the exciton wave function becomes more and more localized, see Fig. 2) permits us to assert that fractional Schrödinger equation accounts for disorder phenomenologically with Lévy index µ being the measure of the degree of disorder. The fact that initial Lévy distributions do not have higher moments of order α > µ 22,23,33 show that such construction describes the systems with broad (wider then Gaussian) distribution of its "Brownian paths" in a generalized Lévy sense. It is almost sure that the physical origin of such broad distribution in solids is disorder. The presence of potential U in the system "tames" the initial Lévy distribution, making it decay faster than that in a corresponding free problem 35,37 . In other words, here once more we have an interplay between the "breadth of disorder distribution" 32 and system potential, which makes the probability distribution (square of modulus of the corresponding wave function) decay faster in space. This makes the problem of a fractional Schrödinger equation resemblant to Anderson localization in disordered systems.
In this context, the very interesting work in mesoscopic physics 41 should be mentioned. This work pioneers the application of Lévy (rather than Gaussian, which are always employed in these kind of problems) distributions for the calculation of conductance through disordered quantum wires.It had been shown in 41 , that the distribution of conductances is entirely due to standard average < ln G > (G stands for conductance) and Lévy index, which is (similar to the present problem) "responsible" for the degree of disorder. All other microscopic characteristics of the defect ensemble were shown to be irrelevant. This means that the theoretical method of 41 is similar to ours in the sense that it gives the effective averaging over disorder without any additional assumptions about its (disorder) microscopic details.
In our problem, one of the physical effects of the interplay of Coulomb interaction screening and disorder is the possibility for exciton to collapse at weak screening or to break up at strong one. As stated above, these www.nature.com/scientificreports/ effects cannot be realized in ordered substance, i.e. at µ = 2 . The next important effect is energy levels crossing 52 , appearing also due to the joint action of Coulomb screening and disorder. The latter effect is related to the levels repulsion and their non-Poissonian statistics, inherent in quantum chaotic systems. This might point to the possible chaotic features of exciton motion in the above 2D dielectric structures, especially in the presence of spin-orbit interaction 63 . The point is that recently such chaotic features have been found in the excitonic spectra of 2D structures with unscreened Coulomb interaction, but with the inclusion of Rashba spin-orbit interaction 64,65 . In context of free electrons and holes, the role of latter interaction in 2D electron gas confined in GaAS quantum well had been studied in 66 . This suggests a generalization of the present problem. Namely, the spin-orbit interaction term can be added to the corresponding fractional Schrödinger equation. In this case, the solution will be more sophisticated as the wave function will be spinor now 65 although the problem can possibly be solved in the momentum space similar to the present consideration. Such a problem turns out to be extremely important for the above classes of substances 67 where chaos can even disrupt the functionality of corresponding optoelectronic and/or photovoltaic devices. It has been shown 65 that the description utilizing the "ordinary" (i.e. that both with conventional Laplacian and unscreened Coulomb interaction) hydrogenic problem does not show strong quantum chaotic features like non-Poissonian energy level statistics, see, e.g., 68 . This may be related to the fact that proper description of such features is possible only within excitonic models, containing fractional Laplacians and screened Coulomb interaction, inherent in the majority of semiconductors. Similar to the above "chaotic models", which stem from averaging over different microscopic disorder realizations, our 2D fractional hydrogenic model is limited to the distances of about 10 lattice constants i.e. around 40-50 Å. This is because at such distances, the description in terms of material constants (like permittivities ε ) is possible and on the other hand, the effective averaging over disorder have been performed resulting in the fractional Laplacian introduction in the Schrödinger equation. The described effects can play an important role in the relaxation of the energies of electron and hole, bound in an exciton. In a disordered 2D substance, instead of a process with well-defined time dependence, the energy relaxation from a highly excited to the ground state may become chaotic. This can obviously hinder useful photovoltaic processes (like conversion of solar energy to electric one) since the disorder may reduce the controllability of the photovoltaic performance.
There is a significant corps of articles, devoted to optical properties of excitons in geometrically confined environments like quantum wells (2D) 69,70 and quantum wires (1D) 71,72 with the disorder, introduced either by fluctuations of the quantum well width 70 or random potential 69 . These kinds of disorder had been treated within the formalism, based on ordinary Schrödinger equation with subsequent averaging over disorder 69 , see Ref. 73 for details of averaging procedure in 70 . Moreover, the randomness of the potential in Ref. 69 was of the white noise type (see expression (12) of 69 ) with ordinary (i.e. non-fractional) Gaussian fluctuations. The main qualitative features of the above considerations coincide with our results. Namely, the more disorder is in the system, the "more localized" in r space the exciton wave function is. Although here we study the spectrum of the 2D fractional hydrogenic problem, this study has been fulfilled in the environment with highly non-Gaussian fluctuations. The results obtained here will in the future be used to calculate the observable characteristics (like spectra of absorption, radiative lifetimes, etc) of the excitons based on our model of fractional Schrödinger equation.
An interesting generalization of our 2D fractional hydrogenic problem is consideration of exciton-phonon bound states-so-called excitonic polarons, see Ref. 74 and references therein. As the coupling between an electron and a hole in such polaron occurs via flexural phonon modes, which are spread between several (usually two) TMD monolayers, this problem may be considered an intermediate between 2D and 3D situations. Once more, the introduction of fractional Laplacian to the corresponding Schrödinger equation (and eventually to the equations, governing the properties of phonons) can be important to consider the effects of disorder, which in such bilayer structure may change the optical absorption and emission spectra drastically, see Ref. 74 and references therein.
One more generalization of the problem considered is to address the 3D fractional quantum-mechanical hydrogenic problem with screened Coulomb interaction. This problem is important for Rydberg excitons (described by the quantum mechanical Kepler problem, see, e.g. 75 ) in amorphous (i.e. highly disordered) bulk semiconductors. This additional dimension may generate certain subtleties (for instance the system may be prone to chaos, see 64,65 ) as compared to present 2D screening, especially in the presence of Rashba spin-orbit coupling 7 .