Edge states and spin-valley edge photocurrent in transition metal dichalcogenide monolayers

# Edge states and spin-valley edge photocurrent in transition metal dichalcogenide monolayers

V.V. Enaldiev Kotel’nikov Institute of Radio-engineering and Electronics of the Russian Academy of Sciences, 11-7 Mokhovaya St, Moscow, 125009 Russia
July 2, 2019
###### Abstract

We develop an analytical theory for edge states in monolayers of transition metal dichalcogenides based on a general boundary condition for a two-band -Hamiltonian in case of uncoupled valleys. Taking into account edge spin-orbit interaction we reveal that edge states, in general, have linear dispersion that is determined by three real phenomenological parameters in the boundary condition. In absence of the edge spin-orbit interaction, edge states are described by a single real parameter whose sign determines whether their spectra intersect the bulk gap or not. In the former case we show that illumination by circularly polarised light results in spin and valley polarised photocurrent along the edge. Flow direction, spin and valley polarisation of the edge photocurrent are determined by the direction of circular polarisation of the illuminated light.

## I Introduction

The optical properties of monolayer crystals of transition metal dichalcogenides (TMDs) (like MoS, MoSe, MoTe, WS, and WSe) have recently attracted a considerable interest due to possible optoelectronic applicationsSun et al. (2016); Pospischil et al. (2014). This is due to the direct band gap of the TMD monolayers whose value corresponds to the visible and infrared light frequenciesKormanyos et al. (2015). The optical response of the bulk materials at absorption edge is dominated by excitonsMoody et al. (2015). However, recent experiments on second harmonic generation at sub-band gap frequencies Yin et al. (2014), scanning tunneling microscopy and spectroscopyZhang et al. (2014); Bollinger et al. (2001) as well as microwave impedance microscopyWu et al. (2016) in MoS-monolayers have also exhibited edge state (ES) signs.

From theoretical point of view properties of ESs in atomically thin MoS are extensively studied in frames of density functional theoryBollinger et al. (2003, 2001); Ataca et al. (2011); Erdogan et al. (2012); Vojvodic et al. (2009); Li et al. (2008) as well as tight-binding approximationKhoeini et al. (2016); Rostami et al. (2016) . However description of the ESs in the TMD monolayers within -approach allows one to describe the ESs without going into details of edge microscopic structure and to take into account effects of external fields. This enables to construct an analytic theory for the ESs in the whole class of materials in a unified way. Such a general theory relies on a boundary condition (BC) that describes the edge structure by means of several phenomenological parameters. The values of these parameters can be obtained by fitting with experimental data or other calculations based on density-functional or tight-binding approximations. Authors of Ref.[Péterfalvi et al., 2015] derive a general boundary condition taking into account valley coupling at the edge and neglecting by edge spin-orbit interaction which may, in general, exist (see Ref.[Volkov and Enaldiev, 2016] and references therein). Recently, ES spectra in the TMD monolayer nanoribbonSegarra et al. (2016) and optical absorption in TMD nanoflakes involving transitions between the bulk and edge statesTrushin et al. (2016) have been studied in the -approximation. However, these studies were restricted by some certain values of the phenomenological boundary parameters.

The aim of this paper is twofold. First, we develop an analytical theory for the ESs in the TMD monolayers in the -approach taking into account edge spin-orbit interaction (ESOI) in the case of uncoupled valleys. Second, we consider optoelectronic properties of the ESs and demonstrate an emergence of spin and valley polarised edge photocurrents due to illumination of the monolayer by circularly polarised light. Origin of the effect concerns with selection rules for optical transitions from bulk states to edge states caused by circularly polarised light in the two valleys. The selection rules are essentially different from those for interband optical transitions in the bulk that give rise, for example, to the valley Hall effectMak et al. (2014). The paper is organized as follows: in Sec.II we derive a general boundary condition and resulting ES spectra for two-band continuum model, Sec.III is devoted to derivation of the edge photocurrent.

## Ii Edge states in two-band kp-approximation

In the TMD monolayers conduction and valence band edges are located in and valleys of the honeycomb lattice Brillouin zone. Within a two-band -approach, dynamics of electrons in the valley is described by the HamiltonianKormanyos et al. (2015); Xiao et al. (2012)

 Hτ=⎛⎜ ⎜ ⎜ ⎜⎝m+τΔcvp−,τ00vp+,τ−m+τΔv0000m−τΔcvp−,τ00vp+,τ−m−τΔv⎞⎟ ⎟ ⎟ ⎟⎠ (1)

where is the band gap without spin splitting, is the value of spin splitting in the conduction and valence band correspondingly, the index denotes the valley, ( are components of in-plane momentum), is the velocity matrix element between the band extrema. The Hamiltonian (1) possesses diagonal form in the spin subspace, with the upper-left (lower-right) block acting on two-component wave functions of spin up (down) states. To describe edge of the TMD monolayer one should supplement the Hamiltonian with a BC for envelope wave functions. In present work we consider a translation invariant edge for which projections of the valley centers onto the edge direction are well distant from each other (like at zigzag or reconstructed zigzag edges). Therefore, we will neglect by the valley coupling at the edge.

However, even at the translation invariant edge, additional ESOI can mix the spins. Let us derive the most general BC that describes this entanglement in our model. Since (1) is of the first order in momentum, a general BC has a form of a linear combination of two-component wave functions belonging to spin up (down) states:

 [Ψ↑τ−MτΨ↓τ]edge=0 (2)

here is the second order square matrix consisting of phenomenological parameters that characterize the edge structure. Explicit form of the matrix is determined by the Hermiticity of the Hamiltonian in a confined region:

 [M+τnστMτ+nστ]edge=0 (3)

where is vector of Pauli matrix, is a outer unit normal to the edge. Eq.(3) can also be obtained from requirement of vanishing normal component of probability current at the edgeAkhmerov and Beenakker (2008). Time reversal symmetry relates the matrices in BC (2) from the two valleys: . These conditions lead us to the four-parametric form for matrix :

 Mτ=ieiχ[σ0sinhξ+τσzcoshξcoshηcosν−iτnστcoshξ(sinhη+σzcoshηsinν)] (4)

where are identity and the third Pauli matrices respectively, (, ) are real phenomenological parameters characterizing edge properties.

As matrix depends on the parameter through a common phase factor that we assume is constant as a function of coordinate along translation invariant edge, we eliminate it by means of a unitary transformation where . Thus, only three boundary parameters characterize translation invariant edge in the TMD monolayer in absence of intervalley interaction. Their physical meaning may be obtained by considering limiting cases of the BC (2). In the limit , two-component wave functions with opposite spins become decoupled and satisfy BC (for derivation, see Appendix A):

 [ψ↑(↓)c,τ−ias,τe−iτφψ↑(↓)v,τ]edge=0 (5)

where , is an angle characterizing the unit normal , for spin up (down) states. However, disentanglement of the wave functions belonging opposite spin projections in BC (5) does not mean vanishing ESOI, as in general case. It is knownZagorodnev et al. (2015) that BC (5) is equivalent to insertion of diagonal in spin subspace potential in the Hamiltonian (1), which is a combination of electrostatic () and pseudo-electrostatic () potentials , where tends to infinity outside the 2D material and is zero inside of it. The sign of is determined by the sign of outside of TMD monolayer. In addition, at in BC (5), values of the boundary parameters coincide and equal to , which is determined only by the parameter . Therefore, we can conclude that and describe different types of ESOI, like Rashba and Dresselhaus interface spin-orbit parameters for Schrodinger’s electrons in GaAs/AlGaAs quantum wellsDevizorova et al. (2014), but is responsible for coupling of bands with the same spin.

Now we are able to calculate spectra of ESs with the general BC (2),(4). Suppose that the 2D crystal fills a half-plane . Then, ES wave function is of the form:

 (6)

where is a normalization constant, is a decay length of ESs. Substituting wave function (6) in the BC (2) we obtain a general dispersion equation for ESs:

 ⎡⎢ ⎢⎣1+a+1,τℏv(τkx−κ↑τ)ε+m−τΔv⎤⎥ ⎥⎦⎡⎢ ⎢⎣1+a−1,τℏv(τkx−κ↓τ)ε+m+τΔv⎤⎥ ⎥⎦+τv(tanhξ−1)sinhη+coshηsinν⎡⎢ ⎢⎣ℏv(τkx−κ↑τ)ε+m−τΔv−ℏv(τkx−κ↓τ)ε+m+τΔv⎤⎥ ⎥⎦=0. (7)

On Fig.1 we reveal typical dispersion of ESs in the and valleys given by the previous equation. In general case, ESs possess linear dispersion and exist for those longitudinal momenta when their energies are not overlapped with projections of bulk bands. The second term in the left-hand side of Eq.(7) is responsible for coupling of spins due to ESOI as it goes to zero at , which is the discussed above limit of spin disentanglement. In this limit, ES spectra are determined by its own parameter for each spin projection:

 ε↑(↓)e,τ(px)=−τ˜vs,τpx+εs,τ, (8)

where is an effective speed of ESs, is energy of ESs at center of the corresponding valley, and . Valley index determines chirality of ESs (see Fig. 2) which exist for those wave vectors while their decay length is positive:

 κ↑(↓)τ=τkx˜ms,τ(εs,τ−sτΔv+Δc2)+˜ms,τ˜vs,τℏv2>0. (9)

We note that condition (9) allows spin-polarised ESs to exist even when their spectra are overlapped with projection of bulk band characterized by opposite spin. In particular, inequality (9) specifies that at ES spectra intersect the bulk gap, but at they are out of the gap (see Fig. 2). Wave functions of ESs with spectrum (8) read as follows:

 (10)

where is the normalization factor.

Here we point out that in absence of ESOI (), ES spectra for two spins are parallel to each other with constant energy difference for every momenta: .

## Iii Spin-valley edge photocurrents

In this section we reveal that illumination of semi-infinite 2D TMD crystal by circularly polarised light induces spin and valley polarised edge photocurrents. As we mentioned in Introduction the effect is due to difference in transition probabilities between bulk and edge bands in the two valleys caused by the light.

For definiteness and simplicity we consider absence of spin-orbit interaction at the edge. Therefore ESs are described by a single parameter for both spins and valleys (we suppressed spin and valley indexes). However, obtained results are also valid when the two spin-polarised ES branches are described by different parameters (8), while one can populate only one of them in each valley. We will regard that , so that ES spectra intersect the gap as shown on Fig.3. This situation agrees with tight-binding calculations of ES spectra at zigzag edge of MoSRostami et al. (2016). Further throughout this section we suppress the valley index everywhere except where it is needed.

Now we turn to the bulk states. Below we are only interested in valence band states with energies around band extremum. In this limit () spectra of the spin up (down) states in valence band of the valley are expressed as follows:

 εv=εt−(vp)22˜m, (11)

where is energy of the upper valence band top in each valley, is 2D momentum modulus. Under assumption about type of the edge mentioned above, bulk states should satisfy the BC (5) with the same parameter for both spins. This BC is satisfied by superposition of incident plane wave and plane wave scattered off the edge with a common longitudinal momentum :

 Ψpx,εv=1√2[ψpx,−py+Rεv,kxψpx,py] (12)

where are plane wave solutions of the Hamiltonian (1) in the limit , is a reflection coefficient that is determined by the BC. Plane wave states around extremum of valence band are of the form (see Appendix B):

 ψpx,py=1√LxLy⎛⎜ ⎜ ⎜⎝−vp2˜meiϑp⎞⎟ ⎟ ⎟⎠eikr, (13)

where is the system area, , is 2D wave vector. In fact in semi-infinite sample is not a good quantum number and should be treated as a function of energy and from Eq. (11): . For valence band states with energies around band extremum reflection coefficient can be approximated as follows .

We suppose that the semi-infinite 2D crystal is illuminated in negative direction of -axis by a clockwise polarised light with frequency (in case of counter clockwise polarisation one should exchange in final formulae (19),(20)). Due to spin splitting of the valence bands (order of eV) one can tune frequency of illuminated light and Fermi-energy in the monolayer, to induce electrical dipole transitions only from the upper valence band in the () valley to ESs with the corresponding spin (see Fig.3). Our aim is to calculate induced edge photocurrent owing to these transitions. To this end we derive kinetic equation for distribution function in frames of Keldysh formalismArseev (2015) (see Appendix C):

 ∂fe(ε)∂t=Windε−ℏω0,ε[fv(ε−ℏω0)−fe(ε)]−fe(ε)∫+∞0Wspε−ℏω,ε[1−fv(ε−ℏω)]dω−fe(ε)−feq(ε)τR, (14)

where is Fermi-Dirac distribution function for states in edge/valence band respectively, is equilibrium Fermi-Dirac function, is a rate of induced transitions between ES, characterized by longitudinal momentum and energy , and the valence band state with the same longitudinal momentum and energy , is a rate of spontaneous transitions caused by interaction with ground state of electromagnetic field, is a phenomenological relaxation time that describes other relaxation processes between edge and valence band states caused by, for example, phonon or electron-electron scattering (below we discuss range of the relaxation time). The rate of the induced transitions reads as follows

 Windε−ℏω0,ε=2πℏIΩcℏω0∣∣V(+,τ)ε−ℏω0,ε(q0)∣∣2ρv,e(ε−ℏω0) (15)

where is intensity of the incident light, is the speed of light, is a volume for quantization of electromagnetic field, is the matrix element of interaction between valence and edge band states (34), density of the valence band states with definite longitudinal momentum near valence band extrema is expressed by the formula:

 ρv,e(ε)=∑pyδ(ε−εpe,py)=Ly√2˜m2πℏvΘ(εt−v2p2e2˜m−ε)√εt−v2p2e2˜m−ε, (16)

where is the Heaviside step function. The probability of spontaneous transitions is expressed as follows:

 Wspε−ℏω,ε=Ω(2π)2ℏc3∑λ=±∫doq∣∣V(λ,τ)ε−ℏω,ε(q)∣∣2ρv,e(ε−ℏω), (17)

here integration goes over solid angle of light wave vector , summation runs over clockwise () and counter clockwise () polarisation of the electromagnetic field.

It is known that intraband energy relaxation processes have the shortest times in MoS monolayersShi et al. (2013); Wang et al. (2012) (order of picoseconds). As electron-electron scattering is very efficient in one dimensionHirtschulz et al. (2008) and also due to possibility of relaxation in edge band via scattering by bulk phonons, we suppose that relaxation time within the edge band is of the same order. This allows us to solve the kinetic equation (14) in quasi-equilibrium approximation, which is justified when intra-band relaxation times is shorter than edge-bulk energy relaxation times. Therefore, we look for distribution function of the edge (valence) states in the form of Fermi-Dirac function with its own quasi Fermi-level (). We also suppose that additional edge-to-valence band energy relaxation processes, characterized by the relaxation time order of nanoseconds. This is mainly due to suppression of phase-space for these processes caused by great difference in density of bulk and edge states. This allows edge-to-valence scattering only with certain longitudinal momentum mismatch (as valence band quasi Fermi-level locates around the top of valence band).

We are interested in stationary solution that equates the right-hand side of Eq.(14) to zero. After integration the kinetic equation over energy in the limit of zero temperature we arrive to the first implicit relation for and (D). Another relation between quasi Fermi-levels is imposed by a particle conservation rule (i.e. the number of holes in the valence band equals the number of photo-excited electrons in the edge band):

 (εt−μv,τ)ρv=(μe,τ−εF)ρe (18)

where is density of states in valence band, is density of edge states. By virtue of ratio between the densities one has inequality , which leads us to a simpler equation for (in comparison with Eq.(D)):

 II0[δτ,1+δτ,−1(ℏvκF2˜ma)2]G(μe,τ−ℏω0)−μe,τ−εFℏω20˜τR=0, (19)

here , ( is refractive index of environment), and

 G(ε)=arctan(ky0(ε)κF)−κFky0(ε)κ2F+k2y0(ε),

where . Deriving Eq.(19) we neglected the dependence of decay length of ESs on the energy (i.e. ). We also disregarded dependence of momentum component on , since we consider the valence band states around the band extremum, i.e. with longitudinal momenta . It is the expression in the square brackets of Eq.(19) that describes distinction in probabilities of transitions from valence band states (12) to ESs (6) in the and valleys, which significantly differs from selection rules for interband transitions in the bulkXiao et al. (2012). This results in uncompensated edge photocurrent with a specific spin in a particular valley. In the limit we obtain explicit expression for electron quasi Fermi-level in the edge band:

 (20)

With known quasi Fermi-levels of the ESs in the two valleys, we use a standard formula for one-dimensional current along the sample edge:

 j=eh(μe,1−μe,−1). (21)

Dependence of the photocurrent on the light intensity is plotted on Fig.4. Solid lines show photocurrent obtained by solving Eq.(D) for electron quasi Fermi-levels, dashed lines represent approximated solution (20). At small intensities, the current (21) is a linear function of the intensity with a tilt that are determined by the edge-valence band transition probability difference in the and valleys. However, in the limit of high intensities the current goes to zero as edge bands in both valleys tend to equal population. It gives rise to a maximum photocurrent in middle region of intensities. Here we notice that light intensity as high as hundreds mW/m in continous wave regime for near-infrared range of wavelengths can be realized by means of fiber lasersTer-Mikirtychev (2014). Using solutions (20) one can obtain an approximated expression for maximal current in Fig.4:

 jmax=eh(εt−εF+ω0)1−(ℏvκF2˜ma)21+(ℏvκF2˜ma)2, (22)

which is valid in the limit . Therefore, the greater edge-valence band transition probability difference in the two valleys, the higher maximal value of spin-valley polarised current along the edge. Absolute values of the maximal current can attain several microamperes. In case of clockwise polarization of light, uncompensated photocurrent flows in the valley in negative direction along the edge (see Fig.3). As it was mentioned above for counter-clockwise polarisation one should exchange in (20), which leads to uncompensated photocurrent in the valley but in positive direction along the edge. Thus, direction of light polarisation controls not only spin and valley polarisation of the uncompensated edge photocurrent but also its direction of flowing.

## Iv Conclusion

We developed a theory of ESs in monolayers of TMD crystals which takes into account ESOI within intra-valley approximation. The theory relies on a general BC comprising three real phenomenological parameters that characterize microscopic structure of the edge. We revealed that is responsible for spin coupling, accounts for inequivalence of the edge structure for opposite spin projections in case of decoupled spins, and describes interband interaction of states with the same spin. In general case, ESs have linear spectra determined by the all three parameters. However, in the case of decoupled spins (), ESs spectra become chiral in valley index, and are described by a single parameter (which is function of and ) for each spin projection. Sign of the latter parameter determines whether ES spectra are in the bulk gap or outside of it.

We also considered optical pumping from valence band states to ESs in absence of the ESOI and demonstrate possibility for generation of spin and valley polarised edge photocurrents. We revealed that direction, valley and spin polarisation of the photocurrent are determined by direction of circular polarisation of the light. Maximal values of the photocurrent are order of several microamperes.

###### Acknowledgements.
This work was supported by the Russian Foundation for Basic Research (Project No. 16-32-00655).

## Appendix A Derivation of spin-polarised BC (5)

In the limit the general BC (2) with matrix (4) is reduced to the following one:

 [−ie−ξ−iχ1+τcoshηcosν(100ia1,τe−iτφ)Ψ↑τ+(1−ia−1,τe−iτφ1−ia−1,τe−iτφ)Ψ↓τ]edge=0, (23)

where and we used an identity:

 sinhη−coshηsinν1+τcoshηcosν=−1−τcoshηcosνsinhη+coshηsinν.

In the limit under consideration, wave functions accounting for opposite spin projections are decoupled in the BC, as the coefficient under in Eq.(23) is exponentially small. This leads us to a BC for spin-down states:

 [ψ↓c−ia−1,τe−iτφψ↓v]edge=0. (24)

where . Subtracting the second raw of the BC (23) from the first one, we arrive for a BC for spin up states:

 [ψ↑c−ia1,τe−iτφψ↑v]edge=0. (25)

here .

## Appendix B Approximate plane wave solutions of the Hamiltonian (1) around bulk band extrema

In this section we find plane wave solutions of the Hamiltonian (1) in a system without edge and show that for energies around valence band extrema they are expressed by the formula (13). For simplicity we consider only spin up electrons in the valley (i.e. ) and suppress the spin and valley indexes below in this section. Therefore, components of a plane wave solution satisfy the system:

 {(m+Δc−ε)ψ1+v(τpx−ipy)ψ2=0v(τpx+ipy)ψ1+(−m+Δv−ε)ψ2=0. (26)

Together with normalization condition

 ∫Lx/2−Lx/2dx∫Ly/2−Ly/2dyΨ∗pΨp=1,

solutions for the amplitudes can be represented as follows:

 ψ1=(εc,v+m−Δv)[(εc,v+m−Δv)2+(vp)2]1/2,ψ2=vpeiθp[(εc,v+m−Δv)2+(vp)2]1/2, (27)

where energies for the plane wave solutions read as follows:

 εc,v=(Δv+Δc)/2±√˜m2+(vp)2, (28)

here sign before square root corresponds to bulk states in conduction (valence) band respectively. Finally, expanding energies of valence band states (28) in expressions for the amplitudes (27) around the band extremum (i.e. at ): , we obtain the following expression for the plane wave solutions in vicinity of valence band maximum:

 Ψp=1√LxLy⎛⎜ ⎜ ⎜⎝−vp2˜meiϑp⎞⎟ ⎟ ⎟⎠eipr/ℏ, (29)

which is identical to Eq.(13) used in the main text. For the states around the conduction band minimum , one can obtain the following expressions:

 Ψp=1√LxLy⎛⎜ ⎜ ⎜⎝1vp2˜meiϑp⎞⎟ ⎟ ⎟⎠eipr/ℏ, (30)

## Appendix C Derivation of kinetic equation (14)

In order to derive the kinetic equation (14) we first write down the Hamiltonian of the system under consideration in terms of second-quantization operators:

 H0=∑peεea+peape+∑px,pyεv,px,pya+v,px,pyav,px,py+∑px,pyεc,px,pya+c,px,pyac,px,py+∑q,λ=±ℏω(c+q,λcq,λ+12), (31)

where is annihilation (creation) operator of the edge state (10) with energy , is annihilation (creation) operator of the valence/conduction band state (12) with energy , is annihilation (creation) operator of photon with clockwise () or counter-clockwise () polarisation and energy . In previous equation we suppress spin and valley indexes for brevity. In terms of the creation and annihilation operators of photon field vector-potential reads as follows:

 A(r,t)=∑q,λ=±√2πc2ℏn2rωΩ[cq,λeλei(qr−ωt)+c+q,λe∗λe−i(qr−ωt)], (32)

where is quantization volume, polarisation unit vectors (spherical angles , characterize direction of the photon wave-vector ). Interaction of electrons with electromagnetic field reads as follows:

 Vint=∑pe,px,py,q,λ=±[V(λ,τ)εv,εe(q)a+pecq,λav,px,py+h.c.], (33)

where we take into account only transitions between valence band electrons and edge states that are determined in dipole approximation by the matrix elements:

 V(λ,τ)εv,εe(q)=τδpe,pxvec√2πc2ℏn2rωΩCe√Lx√2Ly2iky(ε,pe)κ2+k2y(ε,pe)×[e−iταq√2(1+τλcosθq)+eiταq√2(1−τλcosθq)ℏvκ+τvpe2˜ma], (34)

From the previous formula it follows that ratio of probabilities for induced transitions (at normal incidence of light ) in the two valleys for definite polarisation has an order of at the boundary parameter values .

Now we introduce Keldysh Green functions of electrons ( means quantum numbers of edge or bulk state) and photons (). Following standard procedureArseev (2015) we obtain a kinetic equation for Green function that determines population distribution of edge state ():

 i∂∂tG

where , . In the Eq. (35) we calculate self-energies in the second order in perturbation (33):

 (36)

where the terms with wave vector describe transitions induced by illuminated light with frequency and clockwise polarisation, is the number of the illuminated light quanta, the term proportional to (where is Fermi-Dirac distribution function in valence band) concerns with spontaneous recombination processes. After substitution of Green functions of zero approximation , in Eq. (35) and accounting for only contributions from poles at integration over time we arrive to a kinetic equation similar to that used in the main text (14):

 ∂fe∂t=IΩcℏω02πℏ∑p∣∣V(+,τ)εv,εe(q0)∣∣2[fv−fe]δ(εv−εe+ℏω0)−fe2πℏ∑p,q,λ=±∣∣V(λ,τ)εv,εe(q)∣∣2[1−fv]δ(εv−εe+ℏω)−fe−feqτR, (37)

where we add the term with phenomenological relaxation time as in the main text, .

## Appendix D General relation for quasi Fermi-levels

In the section we will obtain a general relation for quasi Fermi-levels of electrons in valence and edge bands under illumination of the light. Introducing the density of valence band states with definite momentum (16) we rewrite Eq. (37) in the following form:

 II0[δτ,1+δτ,−1(ℏvκe+τvpe2˜ma)2][fv(ε−ℏω0)−fe(ε)]√εt−v2p2e2˜m−ε+ℏω0[(ℏvκe)22˜m+εt−v2p2e2˜m−ε+ℏω0]2− fe(ε)12π2(ℏω0)2[1+(ℏvκe+τvpe2˜ma)2]∫+∞0d(ℏω)ℏω[1−fv(ε−ℏω)]√εt−v2p2e2˜m−ε+ℏω[(ℏvκe)22˜m+εt−v2p2e2˜m−ε+ℏω]2− (1+a2)n2rc3√2˜m4a2πτR(ve)2ω20ℏvκe[fe(ε)−feq(ε)]=0 (38)

To proceed further analytically, we consider low temperature limit ( to treat Fermi-functions as step-like ones. Then we integrate the above equation over energy and arrive a final relation between quasi Fermi-levels :

 II0[δτ,1+δτ,−1(ℏvκF2˜ma)2]{arctan[κF(ky0(μe,τ−ℏω0)−ky0(μv,τ))κ2F+ky0(μe,τ−ℏω0)ky0(μv,τ)]+κFky0(μv,τ)κ2F+k2y0(μv,τ)}− arctan(ky0(μv,τ)κ0)⎡⎣(