# Resonant Enhancement of Second Harmonic Generation by Edge States in Transition Metal Dichalcogenide Monolayers

###### Abstract

We derive a low-energy theory for edge states in transition metal dichalcogenide monolayers for a two-band -Hamiltonian in case of uncoupled valleys. In the absence of spin-orbit interaction at the edge, these states possess a linear dispersion described by a single phenomenological parameter characterizing the edge structure. Depending on the sign of the parameter, the edge state spectrum can either cross the band gap or lie outside of it. In the first case, the presence of edge states leads to resonant enhancement of the second harmonic generation at frequencies about the half of the band gap, in agreement with recent experiments. The value of the phenomenological boundary parameter is extracted from the resonance frequency position.

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 the possible optoelectronic applicationsbib:Sun (); bib:Posp (). This is due to the direct band gap of the monolayer TMDs whose value corresponds to the visible and infrared light frequenciesbib:Falko (). The optical response of the bulk materials at the absorption edge is dominated by excitonsbib:Moody (). However, recent experimentbib:Yin () have also demonstrated resonant enhancement of non-linear response near the edge of MoS monolayer membranes at frequencies about the half of the band gap. This resonance was attributed with the existence of the edge states (ESs).

The metallic ESs in the band gap have been lately observed in a MoS monolayer on graphite using scanning tunnelling microscopy and spectroscopybib:Zhang (). The dependence of ES properties in MoS monolayer nanoribbons on the type of edge termination, passivation, and reconstruction were studied using density functional theorybib:Boll (); bib:Ataca (); bib:Erdogan (); bib:Vojv (); bib:Li () (DFT) as well as tight-binding approximationbib:Khoeini (); bib:Liu () (TBA). However, description of the ESs in the TMD monolayers within the -approach allow one to describe the ESs without going into details of microscopic structure of the edge 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 parametersbib:Korman (). The values of these parameters can be obtained by fitting with the experimental data or other calculations based on DFT or TBA. Recently, the ES spectra in the TMD monolayer nanoribbonbib:Segarra () and optical absorption in TMD nanoflakes involving transitions between the bulk and edge statesbib:Trushin () have been studied in the -approximation. However, these studies were restricted by some certain values of the phenomenological parameters.

The aim of this communication is to construct an analytical theory for the ESs in TMD monolayers in the -approach and to reveal the effects of ESs on the non-linear optical response. We demonstrate that the ESs possess linear spectra which are described by a single real phenomenological parameter (for each spin value) in the absence of valley coupling and spin-orbit interaction at the edge. Sign of the parameter determines whether the ES spectra intersect the band gap or not. We show that the second harmonic generation is resonantly enhanced by the ESs crossing the band gap, and extract the value of the parameter from the experimentally observedbib:Yin () resonance frequency.

In the TMD monolayers the conductance and valence band edges are located in the and valleys of the honeycomb lattice. Within a two-band -approach, dynamics of electrons with spin in the valley is described by the Hamiltonianbib:Falko (); bib:Xiao ()

(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, is the in-plane 2D momentum, is the velocity matrix element between the band extrema. The Hamiltonian (1) acts on the two-component wave function . As it was mentioned above, to describe the edge of the TMD monolayer one should supplement the Hamiltonian (1) with a BC for . Here we consider only the zigzag or reconstructed zigzag types of edges for which projections of the valley centers onto the edge direction are well distant from each other. In this case we can neglect by the valley coupling at the edge. We also suppose that spin-orbit interaction is absent at the edge, so we can independently treat the two spin components in the Hamiltonian (1). Under this conditions the most general BC that entangles the components of the function is the same as in a single valley approximation of graphene bib:VVolkov (); bib:Ostaay (); bib:Zag ():

(2) |

where is the real phenomenological parameter for the valley and the spin , characterizing the edge structure (including passivation, relaxation or reconstraction of the zigzag edge), is an angle between the vector of unit normal to the edge and the -axis. The BC (2) is derived from vanishing of the probability current normal to the edge. The identity between the intravalley BC (2) for TMD monolayers and that for graphene stems from the equality of the current operator for the Hamiltonian (1) and that for graphene. Below we will consider the edges which preserve the time reversal symmetry. This symmetry imposes the following relation between the phenomenological parameters in the two valleys for the opposite spin values: . Here we note that study of Ref.[bib:Trushin, ] was limited by the BC (2) with . As the valleys and spins do not entangle in our consideration, the indexes will be hereafter suppressed everywhere except where they are needed.

Now we consider the monolayer occupying the half-plane with translation invariant edge characterized by constant along the -axis. In this case, the quasimomentum component measured from the projections of valley centers on the edge is a good quantum number. The ES wave function that obeys the equation and the BC (2) is

(3) |

here we have introduced the 2D wave vector ,

(4) |

is the normalization factor,

(5) |

is an inverse decay length of the edge states. The dispersion law of the ESs is expressed as follows (see Fig.1):

(6) |

The condition determines the energy range for the existence of the ESs. In case () the spectrum of ESs (6) intersects the band gap and has the end point in the valence band for () or in the conduction band for () in the valley (). In opposite case () energies of the ESs lies outside the band gap and overlap with the valence band for () or conduction band for () in the valley ().

Now we consider the bulk states. In the infinite monolayer the Hamiltonian (1) is diagonalized in the plane wave basis:

(7) |

here , , is the conduction (valence) band index, the normalization factor is expressed as follows:

(8) |

The last equality in (7) is the definition of the two-component function . Spectra of bulk states is expressed by the formula:

(9) |

where , .

In the half-plane, the wave functions (7) do not satisfy the BC (2). Therefore, we will seek the bulk wave function as a sum of incident and reflected plane waves with identical wave vector component :

(10) |

here . The reflection coefficient is determined by the BC (2) and reads as follows:

(11) |

As a consequence of the probability current conservation we have . The wave functions of the edge (3) and bulk (10) states form the complete set in the system under consideration. Here, we note that for the half-plane problem the mirror symmetry plane may also present. However, existence of the mirror plane does not impose any additional constriction for the boundary parameter . This means that arbitrary observable quantity (for example, the second order conductivity) calculated in the basis (3), (10) will keep the mirror symmetry plane restrictions even if the latter is really not present at the edge on the atomic scale.

Now we turn to the calculation of the second order conductivity which is responsible for the second harmonic generation and show that the ESs lying in the gap (see Fig.1(a)) resonantly enhance the latter. For this aim we solve the quantum kinetic equation

(12) |

for the density matrix . The electric field is determined through the spatial derivative of the potential Below, we are interested in the limit that should be taken with a care as linear terms in to be preserved. The solution of Eq.(12) can be expanded in the series of the electric field amplitude:

(13) |

where is the equilibrium density matrix, , . The second harmonic current is determined by the . In the basis (3), (10) matrix elements of are expressed as follows:

(14) |

here is a composite index running over the bulk state quantum numbers and the edge state quantum number , is the equilibrium density matrix whose diagonal elements are the Fermi-Dirac distribution function .

Using Eq.(Resonant Enhancement of Second Harmonic Generation by Edge States in Transition Metal Dichalcogenide Monolayers) one readily obtains the current density at the double frequency as follows:

(15) |

where is the velocity operator (). Equation (15) allows us to express the non-linear conductivity as , here is a contribution to the conductivity from the electrons with the spin in the valley . The resonant term in the current (15) emerges when the intermediate (generally virtual) states in Eq.(Resonant Enhancement of Second Harmonic Generation by Edge States in Transition Metal Dichalcogenide Monolayers) coincide with some real states. For the frequencies less than the bulk band gap, this can only occur when the state belongs to the conduction band bulk, the state belongs to the valence band bulk, and corresponds to the ESs intersecting the band gap (see Fig.1a). Below, we analyse only the resonance terms in the conductivity . The calculations presenting in detail in Supplemental Materialbib:Support () show that the resonance arises only in . After some algebra one finds

(16) | |||

(17) | |||

where is the Heaviside step function. The resonance frequency in the above equations is determined by the formula:

(18) |

From the Eq.(18) it follows that there are two resonance frequencies characterized by the sign of the product . The time reversal symmetry relates the latter from the different valleys with opposite spins . The resonance contributions (16),(17) as a function of frequency are plotted in Fig.(2). The power-law singularity in Eqs.(16),(17) is provided by the density of bulk states which is proportional to and the product of matrix elements that is proportionalbib:Support () to . Therefore, the main contribution to the resonance stems from the bulk electrons propagating along the edge with . The resonance strength is lessened as tends to unity and it vanishes at (see Fig.(2)). This happens because in numerators of the Eqs.(16),(17) arises at . The system size enters explicitly in the denominators of Eqs.(16),(17) due to matrix elements between the bulk and edge states. This is not unusual as the contribution of transitions between the bulk band states via the intermediate ESs to the non-linear two-dimensional conductivity should tend to zero as the system size tends to infinity. Nevertheless, this contribution leads to the non-vanishing total 2D current (even if the system size tends to infinity) and is readily observable in spatially-resolved experiments.

We also note that the conductivity is proportional to the valley index and, therefore, has opposite sign in the two valleys. Consequently, the total resonance term in the non-linear conductivity is nonzero if only electrons have different energy distributions in the two valleys. It is manifestation of the time reversal symmetry imposing the invariance under the mirror symmetry , which leads to the relation (if there is no imbalance in the valley populations). The previous relation ensures the total conductivity to be zero.

Now we extract the value of the boundary parameter from the recent experiment(bib:Yin, ) reporting the resonant second harmonic generation near the edge of a MoS single-layer membrane. The resonance wavelength nm corresponds to the energy eV. The bulk energy gap of atomically thin MoS layer derived from the absorption spectroscopybib:Mak (); bib:Zande () is eV. The DFT calculationsbib:Falko () provide us with the following spin splitting in the valence band eV and in the conduction band eV. Therefore, we have eV, eV. The solutions of Eq.(18) with respect to unknown exist only if . For the bulk parameters given above it is true only for () in the valley . There are two non-equivalent solutions for the phenomenological parameter and (remind that ). Although we do not have additional information to choose the definite one, we can observe from Fig.(2) that the resonance for is more pronounced than that for . Therefore, it is more likely that the boundary parameter is .

To conclude, we have shown that the ESs of TMD monolayers possess linear spectra describing by the only phenomenological parameter in absence of valley coupling and spin-orbit interaction at the edge. At the ES spectra intersect the band gap in the () valley, which leads to resonantly enhanced second harmonic generation in the TMD monolayers. Based on the experimental results on the second harmonic generation from a single layer membranes of MoS we have found the value of the phenomenological parameter for these structures.

This work was partially supported by the Russian Foundation for Basic Research (project 16-32-00655) and the Russian Ministry of Education and Science (”5 top 100” program).

*

## Appendix A Derivation of the Eqs.(16), (17)

Below throughout this section we suppress the indexes in as we are only interested in the matrix elements within the same valley and spin . First, we give an explicit expression for the matrix element , in the linear in approximation between the bulk (10) and the edge (3) states. They read as follows:

(19) | |||

(20) | |||

One can see from Eqs.(19), (20) that the terms linear in are proportional to , but those that linear in are proportional to . Therefore, we can infer that the resonance in the second harmonic current will be in terms that are proportional to . To derive the resonant terms in the non-linear conductivity we also need the matrix elements of the velocity operator between the bulk states (10):

(21) |

(22) |

where we use the identity which follows from the equality . Summation over in Eq.(Resonant Enhancement of Second Harmonic Generation by Edge States in Transition Metal Dichalcogenide Monolayers),(15) implies the following:

(23) |

We note that at integration over the bulk states we gain an additional factor due to the density of bulk states. Inserting Eqs.(19-22) in Eqs.(15), (Resonant Enhancement of Second Harmonic Generation by Edge States in Transition Metal Dichalcogenide Monolayers) we obtain that ,

(24) | |||

(25) | |||

(26) | |||

(27) | |||

where the integration limits are the following , and

On the final step using the Sokhotskiâ-Plemelj formula in Eqs.(24–27)

(28) |

we arrive to the following formulae for contributions which contain the product of the two Dirac delta-functions in the integrand functions:

(29) | |||

(30) | |||

(31) |

(32) |

The power-law singularity emerges only in those components of the non-linear conductivity where the integrand function in Eqs.(24–27) is proportional to .

## References

- (1) Z. Sun, A. Martinez, F. Wang, Nat. Photonics, 10, 227 (2016).
- (2) A. Pospischil, M.M. Furchi, and T. Mueller, Nat. Nanotech., 9, 257 (2014).
- (3) A. Kormanyos, G. Burkard, M. Gmitra, J. Fabian, V. Zolyomi, N. D. Drummond, and V. Fal’ko, 2D Mater. 2, 022001 (2015).
- (4) G. Moody, C. K. Dass, K. Hao, C.-H. Chen, L.-J. Li, A. Singh, K. Tran, G. Clark, X. Xu, G. Berghuser, E. Malic, A. Knorr, and X. Li, Nat. Commun., 6, 8315 (2015).
- (5) X. Yin, Z. Ye, D.A. Chenet, Y. Ye, K. OâBrien, J.C. Hone, X. Zhang, Science, 344, 488 (2014).
- (6) C. Zhang, A. Johnson, C. Hsu, L. Li, and C. Shih, Nano Lett., 14, 2443 (2014).
- (7) M.V. Bollinger, K.W. Jacobsen, and J.K. Norskov, Phys. Rev. B, 67, 085410 (2003).
- (8) C. Ataca, H. Sahin, E. Akturk, and S. Ciraci, J. Phys. Chem. C, 115, 3934 (2011).
- (9) E. Erdogan, I.H. Popov, A.N. Enyashin, and G. Seifert, Eur. Phys. J. B, 85, 33 (2012).
- (10) A. Vojvodic, B. Hinnemann, and J.K. NÃ¸rskov, Phys. Rev. B, 80, 125416 (2009).
- (11) Y. Li, Z. Zhou, S. Zhang, and Z. Chen, J. Am. Chem. Soc, 130, 16739 (2008).
- (12) F. Khoeini, Kh. Shakouri, and F.M. Peeters, Phys. Rev. B, 94, 125412 (2016).
- (13) G. Liu, W. Shan, Y. Yao, W. Yao, and D. Xiao, Phys. Rev. B, 88, 085433 (2013).
- (14) C.G. PÃ©terfalvi, A. KormÃ¡nyos, G. Burkard, Phys. Rev. B, 92, 245443 (2015).
- (15) C. Segarra, J. Planelles, and S.E. Ulloa, Phys. Rev. B, 93, 085312 (2016).
- (16) M. Trushin, E.J.R. Kelleher, T. Hasan, arXiv:1602.06298v2.
- (17) D. Xiao, G. Liu, W. Feng, X. Xu, and W. Yao, Phys. Rev. Lett., 108, 196802 (2012).
- (18) J.A.M. van Ostaay, A.R. Akhmerov, C.W.J. Beenakker, M. Wimmer, Phys. Rev. B 84, 195434 (2011).
- (19) V.A. Volkov and I.V. Zagorodnev, Low Temp. Phys. 35, 2 (2009).
- (20) I.V. Zagorodnev, Zh.A. Devizorova, and V.V. Enaldiev, Phys. Rev. B 92, 195413 (2015).
- (21) See Supplemental Material at … for derivation of the non-linear conductivity resonant terms (16), (17).
- (22) K.F. Mak, C. Lee, J. Hone, J. Shan, and T.F. Heinz, Phys. Rev. Lett. 105, 136805 (2010)
- (23) A. M. van der Zande, P. Y. Huang, D.A. Chenet, T.C. Berkelbach, Y. You, G.H. Lee, T.F. Heinz, D.R. Reichman, D.A. Muller, and J.C. Hone, Nat. Mater. 12, 554 (2013)