Strain engineering of electronic properties of transition metal dichalcogenide monolayers

Author’s Accepted Manuscript Strain engineering of electronic properties of transition metal dichalcogenide monolayers Aristea E. Maniadaki, Georgios Kopidakis, Ioannis N. Remediakis www.elsevier.com/locate/ssc

PII: DOI: Reference:

S0038-1098(15)00410-X http://dx.doi.org/10.1016/j.ssc.2015.11.017 SSC12818

To appear in: Solid State Communications Received date: 3 October 2015 Accepted date: 20 November 2015 Cite this article as: Aristea E. Maniadaki, Georgios Kopidakis and Ioannis N. Remediakis, Strain engineering of electronic properties of transition metal dichalcogenide monolayers, Solid State Communications, http://dx.doi.org/10.1016/j.ssc.2015.11.017 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting galley proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Strain engineering of electronic properties of transition metal dichalcogenide monolayers Aristea E. Maniadaki a, Georgios Kopidakis a,b,∗, Ioannis N. Remediakis a,b b Institute

a Department of Materials Science and Technology, University of Crete, GR-71003 Heraklion, Crete, Greece of Electronic Structure and Laser (IESL), Foundation for Research and Technology Hellas (FORTH), GR-71110 Heraklion, Crete, Greece

Abstract We present Density Functional Theory (DFT) results for the electronic and dielectric properties of single-layer (2D) semiconducting transition metal dichalcogenides MX 2 (M = Mo, W; X = S, Se, Te) under isotropic, uniaxial (along the zigzag and armchair directions), and shear strain. Electronic band gaps decrease while dielectric constants increase for heavier chalcogens X. The direct gaps of equilibrium structures often become indirect under certain types of strain, depending on the material. The effects of strain and of broken symmetry on the band structure are discussed. Gaps reach maximum values at small compressive strains or in equilibrium, and decrease with larger strains. In-plane dielectric constants generally increase with strain, reaching a minimum value at small compressive strains. The out-of-plane constants exhibit a similar behavior under shear strain but under isotropic and uniaxial strain they increase with compression and decrease with tension, thus exhibiting a monotonic beahavior. These DFT results are theoretically explained using only structural parameters and equilibrium dielectric constants. Our findings are consistent with available experimental data. Keywords: molybdenum disulfide, monolayer, transition metal dichalcogenides, strain, band gap, dielectric constant

The isolation of atomically thin sheets from layered materials has generated enormous interest in two-dimensional (2D) crystals [1]. During the last decade, single atomic layers cleaved from materials such as boron nitride, graphite, molybdenum disulfide, etc, exhibit novel electronic and optoelectronic properties, different from their bulk, three-dimensional (3D) counterparts [2, 3]. Layered materials such as molybdenum disulfide and other Transition Metal Dichalcogenides (TMDs) had been extensively studied over the past few decades as catalysts [4, 5, 6], lubricants [7] and materials for solar cells [8, 9, 10]. Semiconducting TMDs consist of weakly coupled MX2 hexagonal atomic layers (M = Mo, W and X =S, Se, Te), the M atom being sandwiched between X atoms. Similar to graphite, their layered 3D structure allows for the extraction of single or few layers. Unlike gapless, semimetallic graphene though [11], these new 2D materials possess a band gap and their unique electronic and optical properties are currently under intensive investigation. These properties depend on dimensionality and nanostructuring so that TMDs may be tailored for specific applications. When combined in layered 2D structures, with graphene and other materials, they are expected to revolutionize nanoscale devices [12]. It has been shown relatively recently that the indirect band gap of bulk MX2 increases when they are reduced to a few layers and becomes direct for the monolayer, resulting in a dramatic ∗ Corresponding author. Address: Department of Materials Science and Technology, University of Crete, P.O. Box 2208, 71003 Heraklion, Crete, Greece; Tel.: +30 2810 394218; Fax: +30 2810 394273 Email address: [email protected] (Georgios Kopidakis) URL: http://theory.materials.uoc.gr (Georgios Kopidakis)

Preprint submitted to Solid State Communications

increase in photoluminescence and in novel excitonic effects [13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23]. Breaking of inversion symmetry and strong spin-orbit coupling gives rise to valley selective circular dichroism [24]. Quasi-one-dimensional (1D) and zero-dimensional (0D) structures such as nanoribbons and flakes (dots), respectively, exhibit robust metallic edge states [25, 26, 27, 28, 29]. It was predicted [30] and verified experimentally [31] that single-layer TMDs exhibit piezoelectric behaviour [32]. Thus, TMD 2D semiconductors of atomic thickness show great promise for optoelectronic applications and devices have already been fabricated, e.g., transistors [33, 34], solar cells [35, 36, 37], light-emitting diodes [38]. Tunability of optoelectronic properties, which is essential for applications, can be achieved by controlling dimensionality, number of layers, size and shape of nanostructures, chemical substitutions, and strain. Strain is often intentionally or unintentionally present in TMDs monolayers, depending on synthesis method, substrate, stacking, and several other conditions. These materials are flexible and large strains can be applied without damage [39], allowing for electronic stucture engineering applications without any other physical or chemical modification [40]. A number of theoretical and experimental works have examined the role of strain on the electronic properties of several TMDs [41, 42, 43, 44, 45, 46]. Different types of strain on different TMDs have several interesting effects. In many cases, the direct gap becomes indirect and reduced in value, to the point that the material becomes metallic under extreme strain. Detailed knowledge of the variation of electronic band structure under strain is important. Moreover, screening effects determine electron and exciton dynamics, as well as optoelecNovember 26, 2015

ulation suite [51]. The dielectric function is calculated using linear-response theory optimized for the GPAW code [52]. Electron-electron interactions are taken into account using the random phase approximation (RPA). For the in-plane polarization the electric field is parallel to the MX 2 layer (E ⊥ c) and for the out-of-plane polarization it is perpendicular (E  c). The dielectric constant is obtained in the limit ω → 0 of ε(ω). For 2D systems, we found that a vacuum of 12 Å between monolayers gives ε values in reasonable agreement with available experimental data. We are interested in the variation of ε with strain and exact values are beyond the scope of this work. Standard DFT calculations, such as the ones we perform here, agree with photoluminesce experiments for the lowest excitation energy of single-layer TMDs. However, DFT is well known to underestimate electronic band gaps, with actual values approximately one and a half times the DFT values. Dielectric screening effects and excitons should be included. Interestingly, for MX2 , these two contributions appear to almost cancel each other. For instance, for MoS 2 the DFT band gaps are found between 1.5 and 2 eV. Inclusion of screening employing a GW scheme opens the band gap up to about 3 eV. Adding excitons by means of the Bethe-Salpeter equation brings the minimum excitation below 2 eV [18, 19, 20, 36], in agreement with experiments and close to the values found in standard DFT. Apart from this fortuitous agreement between DFT and many-body calculations or experiments, standard DFT can be trusted for band structures and the calculation of the static relative permittivity, εr . As an average over the whole band structure, it is rather insensitive to details of the computational method and the simple DFT-RPA method employed here gives values that are within 10% off values of much more detailed calculations and experiments [52, 53].

Figure 1: a) Top and side view of the hexagonal lattice structure of MX2 . M and X atoms are represented by blue and yellow spheres, respectively. The shaded area corresponds to a unit cell. b) The corresponding Brillouin Zone. c) Schematic representation of the different types of strain applied on the unit cell.

tronic properties in 2D TMD semiconductors. The effects of strain on the dielectric constant are crucial in this aspect and relevant studies are limited [45]. In this work, we systematically investigate the MX 2 band structure and dielectric constant variations under principal types of strain with Density Functional Theory (DFT). The rich behavior of these systems is summarized and general trends are established. We show that our results are consistent with available experiments and in some cases we provide theoretical explanations which validate our DFT results. The equilibrium atomic structure of an MX 2 monolayer (with M = Mo, W and X =S, Se, Te) is shown in Fig.1(a). The 2D lattice has hexagonal symmetry (top view) and consists of an M atomic layer sandwiched between two X atomic layers (side view). It is characterized by the distance between neighboring M atoms, a, and an internal parameter, u, which describes the relative distance between M and X atoms. The x- and ydirections coincide with the zigzag and armchair directions, respectively. The corresponding Brillouin Zone (BZ) is shown in Fig.1(b). The 3D lattice consists of stacked MX 2 monolayers which are weakly bonded with Van der Waals interactions. The vector in the z-direction perpendicular to the monolayer is c, so that the distance between M and X layers is ( 43 − u) c. Fig.1(c) shows the effect of the different types of strain we apply on the unit cell, i.e., isotropic, uniaxial-x, uniaxial-y, and shear-xy. We perform DFT total energy calculations for all strain-free (equilibrium) MX2 with full atom relaxation. Then, for all strained structures, we only allow for u-parameter relaxation. We perform DFT calculations with the Grid-based Projected Augmented Wave (GPAW) open-source implementation [47, 48, 49]. GPAW uses a real space grid to describe electron densities and wave functions. We use the Generalized Gradient Approximation (GGA) Perdew-Burke-Ernzerhof exchangecorrelation functional (PBE) [50] and, for atomic relaxation, a conjugate-gradient minimization algorithm as implemented in the open-source Atomic Simulation Environment (ASE) sim-

Structure MoS2 WS2 MoSe2 WSe2 MoTe2 WTe2

E g (eV) 1.34 1.29 1.24 1.21 0.87 0.88

E g -exp (eV) 1.29[54], 1.23[55] 1.35[55] 1.10[54], 1.09[55] 1.20[55] 1.00[54]

ε⊥ 12.3 11.0 14.0 12.3 19.2 16.7

ε 7.2 7.2 9.0 8.6 13.0 11.7

ε (exp) 17.0[56] 13.5[57] 18.0[56] 14.0[57] 20.0[56]

Table 1: Bulk (3D) MX2 band gap, theoretical and experimental values, inplane and out-of plane static relative permittivity.

The lattice parameters we obtain from DFT calculations compare very well with experiments. Table 1 includes our results for the electronic band gap, E g and the static relative permittivity, εr , compared to experiment, for 3D materials. For simplicity, we use the notation ε r = ε and we often use the term dielectric constant when referring to relative permittivity. In equilibrium, ε ⊥ = ε xx = εyy and ε = εzz . For 2D MX2 (monolayers), in addition to band gap and dielectric constant, Table 2 shows the structural parameters a, u. The distance, d, between nearest chalcogen atoms X along the perpendicular to the monolayer (z-direction), which is a measure of monolayer thickness, is also shown. Our results are in good agreement with experiments and other DFT calculations. Before examining the effects of strain, it is interesting to point 2

E-EF(eV)

4 3 2 1 0 -1 -2 -3 -4 M

MoS2

Γ

MoSe2

KM

Γ

MoTe2

Γ

KM

WS2

Γ

KM

WSe2

KM

WTe2

Γ

KM

Γ

K

Figure 2: TMDs band structure along the M-Γ-K path in the Brillouin Zone

MX2 MoS2 WS2 MoSe2 WSe2 MoTe2 WTe2

a (Å) 3.20 3.19 3.30 3.31 3.56 3.56

u 0.635 0.630 0.631 0.628 0.628 0.629

d (Å) 3.15 3.18 3.35 3.39 3.59 3.62

E g (eV) 1.69 1.85 1.51 1.63 1.07 1.06

ε⊥ 6.1 5.3 7.0 6.0 9.1 8.2

ε 3.6 3.4 4.1 3.8 4.8 4.6

WS2

MoS2

Eg (eV)

Table 2: Equilibrium MX2 single layer lattice parameters, X-X distance (d), band gap, in-plane and out-of plane static relative permittivity.

out some characteristic features of equilibrium 2D structures. There is a very small lattice expansion when going from bulk to monolayer, in agreement with previous calculations. Comparing the MX 2 lattice constant for the monolayers, it is clear that a increases for heavier X. Not surprisingly, unlike their indirect gap 3D counterparts, all strain-free MX 2 monolayers are direct gap semiconductors at the K-point of the BZ as the band structure diagrams of Fig.2 show, with increased E g . The K, K , K and M, M , M points (Fig.1(b)) are equivalent, respectively. Moreover, for all MX 2 the static relative permittivity is almost halved from 3D to 2D, as already found for MoS 2 theoretically [29, 58] and experimentally [59]. The gap decreases for heavier X, the conduction band minimum (CBmin) is mildly lowered and the valence band maximum (VBmax) is lifted up (Fig.2), but the dielectric constant increases. We find that strain significantly affects the electronic properties of semiconducting TMDs. Band gap values as a function of isotropic, uniaxial-x, uniaxial-y, and shear-xy %-strain are shown for all MX 2 in Fig.3. The wide strain range, from −12% (compressive) to 12% (tensile), is justified by experiments (mainly, but not exclusively, on MoS 2 ) which demonstrate the high elasticity of these 2D materials and their resilience to mechanical deformation. Smaller strains are relevant to TMD, graphene, BN, etc, heterostructures, as well as free-standing monolayers with small deformations. A general feature for all MX 2 and for all types of strain is that E g is maximum at, or close to, zero strain (equilibrium) and decreases away from it. For all MX 2 and for very large tensile isotropic strains E g = 0. The same is true for very large compressive isotropic and uniaxial-x strains, except for MS 2 . For all other

3 2.5 2 1.5 1 0.5 0 -15 -10 -5

0

3 2.5 2 1.5 1 0.5 0 5 10 15 -15 -10 -5

Eg (eV)

2

2

1.5

1.5

1

1

0.5

0.5 0

0 5 10 15 -15 -10 -5

5 10 15

0

5 10 15

WTe2

MoTe2

Eg (eV)

0 WSe2

MoSe2

0 -15 -10 -5

shear xy

uniaxial y

uniaxial x

isotropic

1.5

1.5

1

1

0.5

0.5

0 0 -15 -10 -5 0 5 10 15 -15 -10 -5 0 5 10 15 strain (%) strain (%) Figure 3: Band gap, Eg , versus strain for all MX2 and types of strain. Results from calculations are shown with black circles for isotropic strain, red for uniaxial-x, green for uniaxial-y, and blue for shear-xy. Lines are a guide to the eye. For MS2 (top two plots) the direct gap at K-point of the BZ is also shown with open circles and dotted lines.

3

Mo

shear xy

Mo e

isotropic

shear xy

M K

M K

CBmin

uniaxial y

VBmax

uniaxial x

CBmin

isotropic

VBmax

uniaxial y

CBmin

uniaxial x

VBmax

isotropic

-15

-10

0

strain (%)

5

uniaxial y

15

-15

-10

isotropic

0

strain (%)

e

15

-15

M K

M K

CBmin

-15

-10

-5

0

strain (%)

5

10

15

-15

-5

0

strain (%)

e

uniaxial x

5

uniaxial y

10

15

shear xy

M K M K

M K

M K

-10

isotropic

shear xy

VBmax

uniaxial y

10

CBmin

uniaxial x

5

VBmax

shear xy

-5

CBmin

uniaxial x

10

shear xy

M K

VBmax

isotropic

-5

uniaxial y

M K

M K

M K

Mo e

uniaxial x

-10

-5

0

strain (%)

5

10

15

-15

-10

-5

0

strain (%)

5

10

15

Figure 4: CBmin-VBmax diagrams for all six TMDs structures. The black (very small) circles denote the CBmin and VBmax positions for isotropic strain, the red (small) for uniaxial-x, the green (large) for uniaxial-y and the blue (very large) for shear-xy strain. The triangle and square symbols denote the CBmin and VBmax positions along the Γ-K -M -Γ and Γ-K -M -Γ, respectively, and are marked when they differ from the ones along the Γ-K-M-Γ path.The dashed lines indicate the positions of the high symmetry points of the BZ.The top area of each diagram shows the CBmin and the bottom the VBmax positions.

the CBmin at K, so that the gap is indirect. Thus, even small isotropic strain makes 2D WS 2 an indirect gap semiconductor (at Γ-K) for tensile strain. The conclusions are very similar for uniaxial strains in the range −10% to 10%. For shear strain, we find that CBmin and VBmax remain very close to K, so that the material remains direct gap semiconductor. For compressive shear strain, the smallest gap appears at the K  point, although Eg values at K, K , K are very similar. Our results show that gaps change position in the BZ under strain and, very often, the direct K-gap value remains very close to the strain-induced gaps E g . Experimental indications for this behavior have been reported [60]. Apart from details, Figs.3 and 4 show that it is mainly X that determines electronic band structure variations with strain.

cases, large strains induce significant E g reduction but the material remains semiconducting. In the small strain regime (±3%), shear-xy does not seem to significantly affect E g . For very small compressive strains, E g slightly increases, with the exception of MSe2 . From the top two plots of Fig.3, where the direct gap at K is also shown for MoS 2 and WS2 , it is clear that for small strains the gap at K remains very close to E g , something which could be observed in experiments. Fig.4 summarizes the strain-induced changes in the band structure for all MX 2 and types of strain. On these plots, which complement the E g vs strain graphs of Fig.3, the circle points indicate the positions of CBmin and VBmax along the Γ-K-MΓ path of the BZ. The transitions from direct band gaps at the K point for zero strain to indirect for many strained monolayers are evident in the plots. For very small strains, the band structures and, thus, the CBmin and VBmax, are very similar to the one of the paths Γ-K  -M -Γ and Γ-K -M -Γ. However, for shear-xy and large uniaxial strains the lift of degeneracy of the three aforementionned paths becomes evident. This is also suggested for K, K  and M, M from previous calculations [41]. The points corresponding to the different CBmin and VBmax along Γ-K -M -Γ and Γ-K -M -Γ are marked with triangle and square symbols, where K-M are replaced by K  -M and K -M , respectively. Figs.3 and 4 contain all essential information for the band structure of strained MX 2 monolayers. For instance, for WS2 (bottom left plot of Fig.4), even very small compressive isotropic strain brings the CBmin form K to a point between Γ and K. The VBmax remains at K (only for extreme strain it goes to Γ) so that the gap is indirect. Very small tensile isotropic strain shifts the VBmax to the Γ point and leaves

Our results for the dielectric constant are presented in Fig.5 and show a similar qualitative behavior for all TMDs. Since the materials under strain become anisotropic, all three diagonal elements of the dielectric tensor are included in our analysis. Similarly to the electronic band structure, it is the chalcogen X that mostly influences the dielectric constant variation with strain. The in-plane constants, ε xx and εyy , exhibit minima at small tensile strains, with the exception of uniaxial-x which increases almost lineraly with expansion. The out-of-plane ε zz increases for compressive and decreases for tensile strain, with the exception of shear strain, where a relatively weak variation is observed with a minimum close to equilibrium. A more careful analysis reveals that ε zz strongly depends on the distance between X atoms, d, which does not change significantly for shear strain. The dielectric tensor ε depends on strain tensor u 4

Figure 5: In-plane εxx , εyy and out-of-plane εzz static relative permittivities of MX2 as a function of strain.

through [61] ε = ε0 + δε(u)

(1)

where the variation δε(u) from the equilibrium dielectric tensor ε0 is a function of u. It is generally expressed as δ ik = aiklm · ulm + γiklm · ωlm , with aiklm and γiklm denoting electrostriction coefficients; here the rotation tensor ω lm is zero. In the case of out-of-plane ε zz , the variation δε takes the relatively simple form [62, 63] δd δV + a2 (2) d V where δd is the change of the X-X distance d, δV/V the volumetric strain, and δε = a1

2 a1 = − · (ε − 1)2 5 1 2 · (ε − 1)2 . a2 = − · (ε − 1)(ε + 2) + 3 15

Figure 6: Out-of-plane static relative permittivity for MS2 as a function of all types of strain. Points are DFT results and continuous curves are theoretical predictions (see text).

(3) gap semiconductors under tensile isotropic and uniaxial moderate strain (up to ∼ 5%) and become indirect in all other cases. MTe2 remain direct for all types of moderate tensile and very small compressive strains. The effects of broken lattice symmetry for uniaxial and shear strain on the band structure were discussed. We also presented results for the in-plane ε xx , εyy , and out-of-plane ε zz diagonal elements of the static relative permittivity tensor. All three constants increase with heavier X and are strain dependent. The in-plane constants, which are equal in unstrained monolayers, generally increase with strain, exhibiting a minimum at small compressive strains. The out-ofplane constant under isotropic and uniaxial strain increases with compressive and decreases with tensile. Under shear strain, it exhibits a similar to the in-plane constants variation. These results are theoretically explained solely based on structural parameters and equilibrium dielectric constant values. Many of our predictions for band structure and dielectric constant are in agreement with available data or could be verified by experi-

(4)

Using the ε values of Table 2, the results from Eqs. (1)-(4) are compared with DFT in Fig.6 for ε zz of MS2 . Theory fits DFT results very well up to moderate strains (∼ ±6%), as shown in the plots. A similar analysis for the other diagonal elements of the dielectric tensor, i.e., the in-plane constants ε xx and εyy is more involved. We have systematically studied the effects of the different types of strain on the electronic band structure and the dielectric constant of 2D semiconducting TMDs. At zero strain, electronic band gaps are direct for all MX 2 and while they generally decrease with strain and for heavier chalcogen X, we find a rich variation of electronic properties depending on X and the type of strain. For MS 2 the direct gap changes to indirect for isotropic and uniaxial, whereas for shear, these structures remain direct gap up to very large strains. MSe 2 remain direct 5

ment.

[32] E.J. Reed, Nat. Nanotechnol. 10 (2015) 106-107. [33] B. Radisavljevic, A. Radenovic, J. Brivio, V. Giacometti, A. Kis, Nat. Nanotechnol. 6 (2011) 147-150. [34] A.K.M. Newaz, D. Prasai, J.I. Ziegler, D. Caudel, S. Robinson, R.F. Haglund Jr., K.I. Bolotin, Solid State Commun. 155 (2013) 49-52. [35] M. Shanmugam, T. Bansal, C.A. Durcan, B. Yu, Appl. Phys. Lett. 100 (2012) 153901-153905. [36] M. Bernardi, M. Palummo, J.C. Grossman, Nano Lett. 13 (2013) 36643670. [37] M.-L. Tsai, S.-H. Su, J.-K. Chang, D.-S. Tsai, C.-H. Chen, C.-I Wu, L.-J. Li, L.-J. Chen, J.-H. He, ACS Nano 8 (2014) 8317. [38] J.S. Ross, P. Klement, A.M. Jones, N.J. Ghimire, J. Yan, D.G. Mandrus, T. Taniguchi, K. Watanabe, K. Kitamura, W. Yao, D.H. Cobden, X. Xu, Nat. Nanotechnol. 9 (2014) 268-272. [39] S. Bertolazzi, J. Brivio, A. Kis, ACS Nano 5 (2011) 9703-9709. [40] S. Manzeli, A. Allain, A. Ghadimi, A. Kis, Nano Lett. 15 (2015) 53305335. [41] P. Johari, V.B. Shenoy, ACS Nano 6 (2012) 5449-5456. [42] E. Scalise, M. Houssa, G. Pourtois, V. Afanasev, A. Stesmans, Nano Research 5 (2012), 43-48. [43] K. He, C. Poole, K.F. Mak, J. Shan, Nano Lett. 13 (2013) 2931. [44] R.C. Cooper, C. Lee, C.A. Marianetti, X. Wei, J. Hone, J.W. Kysar, Phys. Rev. B 87 (2013) 035423. [45] A. Kumar, P.K. Ahluwalia, Physica B 419 (2013) 66-75. [46] D.M. Guzman, A. Strachan, J. Appl. Phys. 115 (2014) 243701. [47] https://wiki.fysik.dtu.dk [48] J.J. Mortensen, L.B. Hansen, K.W. Jacobsen, Phys. Rev. B 71 (2005) 035109. [49] J. Enkovaara, C. Rostgaard, J.J. Mortensen et al., J. Phys.: Condens. Matter 22 (2010) 253202. [50] J.P. Perdew, K. Burke, M. Ernzerhof, Phys. Rev. Lett. 77 (1996) 3865. [51] S.R. Bahn, K.W. Jacobsen, Comput. Sci. Eng. 4 (2002) 56-66. [52] J. Yan, J.J. Mortensen, K.W. Jacobsen, K.S. Thygesen, Phys. Rev. B 83 (2011) 245122. [53] M. Gajdoˇs, K. Hummer, G. Kresse, J. Furthm¨uller, and F. Bechstedt, Phys. Rev. B 73 (2006) 045112. [54] Th. Boker, R. Severin, A. Muller, C. Janowitz, R. Manzke, D. Voß, P. Kruger, A. Mazur, J. Pollmann, Phys. Rev. B 64 (2001) 235305. [55] K.K. Kam , B.A. Parkinson, J. Phys. Chem. 86 (1982) 463-467. [56] R. Beal, H.P. Hughes, J. Phys. C: Solid State Phys. 12 (1979) 881. [57] A.R. Beal, W.Y. Liang, H.P. Hughes, J. Phys. C: Solid State Phys. 9 (1976) 2449. [58] E.J.G. Santos, E. Kaxiras, ACS Nano 7 (2013) 10741-10746. [59] D.O. Sigle, J. Mertens, L.O. Herrmann, R.W. Bowman, S. Ithurria, B. Dubertret, Y. Shi, H.Y. Yang, C. Tserkezis, J. Aizpurua, J.J. Baumberg, ACS Nano (2015) 825-830. [60] Y. Wang, C. Cong, W. Yang, J. Shang, N. Peimyoo, Y. Chen, J. Kang, J. Wang, W. Huang, T. Yu, Nano Research 8 (2015) 2562-2572. [61] L D Landau, J. S. Bell, M. J. Kearsley, L. P. Pitaevskii, E.M. Lifshitz, J. B. Sykes, Electrodynamics of Continuous Media, Vol. 8, 2nd ed., Elsevier, 1984 [62] J.A. Stratton Electromagnetic Theory, McGraw-Hill, New York, 1941. [63] H.Y. Lee, Y. Peng, Y.M. Shkel, J. Appl. Phys. 98 (2005) 074104.

1. Acknowledgements This research has been co-financed by EU and Greek national funds through the NSRF Heracleitus II and ERC-02 ExEL Programs. 2. References [1] K.S. Novoselov, D. Jiang, F. Schedin, T.J. Booth, V.V. Khotkevich, S.V. Morozov, A. K. Geim, Proc. Natl. Acad. Sci. U.S.A. 102 (2005) 1045110453. [2] S.Z. Butler, S.M. Hollen, L. Cao, Y. Cui, J.A. Gupta, H.R. Guti´errez, T.F. Heinz, S.S. Hong, J. Huang, A.F. Ismach, E. Johnston-Halperin, M. Kuno, V.V. Plashnitsa, R.D. Robinson, R.S. Ruoff, S. Salahuddin, J. Shan, L. Shi, M.G. Spencer, M. Terrones, W. Windl, J.E. Goldberger, ACS Nano 7 (2013) 28982926. [3] F. Xia, H. Wang, D. Xiao, M, Dubey, A. Ramansubramaniam, Nat. Photonics (2014) 899-907. [4] R. Prins, V.H.J. De Beer, G.A. Somorjai, Catal. Rev. 31 (1989) 1-41. [5] B. Hinnemann, P.G.Moses, J. Bonde, K.P. Jørgensen, J.H. Nielsen, S. Horch, I. Chorkendorff, J.K. Nørskov, J. Am. Chem. Soc. 127 (2005) 5308-5309. [6] T.F. Jaramillo, K.P. Jørgensen, J. Bonde,1 J.H. Nielsen, S. Horch, I. Chorkendorff, Science 317 (2007) 100-102. [7] F.J. Clauss, “Solid Lubricants and Self-Lubricating Solids” (Elsevier, New York 1972). [8] J.A. Wilson, A.D. Yoffe, Adv. Phys. 18 (1969) 193-335. [9] H. Tributsch, Solar Energy Mater. 1 (1979) 257-269. [10] G. Kline, K. Kam, D, Canfield, B.A. Parkinson, Solar Energy Mater. 4 (1981) 301-308. [11] A.H. Castro Neto, F. Guinea, N.M.R. Peres, K.S. Novoselov, A.K. Geim, Rev. Mod. Phys. 81 (2009) 109-163. [12] A.K. Geim, I.V. Grigorieva, Nature 499 (2013) 419-425. [13] K.F. Mak, C. Lee, J. Hone, J. Shan, T.F. Heinz, Phys. Rev. Lett. 105 (2010) 136805. [14] A. Splendiani, L. Sun, Y. Zhang, T. Li, J. Kim, C-Y Chim, G. Galli, F. Wang, Nano Lett. 10 (2010) 1271-1275. [15] T. Korn, S. Heydrich, M. Hirmer, J. Schmutzler, C. Sch¨uller., App. Phys. Lett. 99 (2011) 102109. [16] G. Kioseoglou, A.T. Hanbicki, M. Currie, A.L. Friedman, D. Gunlycke, B.T. Jonker, App. Phys. Lett. 101 (2012) 221907. [17] Z.Y. Zhu, Y. C. Cheng, and U. Schwingensch¨olgl, Phys. Rev. B 84 (2011) 153402. [18] T. Cheiwchanchamnanjij, W.R.L. Lambrecht, Phys. Rev. B 85 (2012) 205302. [19] A. Ramasubramaniam, Phys. Rev. B 86 (2012) 115409. [20] K.P. Komsa, A.V. Krasheninnikov, Phys. Rev. B 86 (2012) 241201. [21] E.S. Kadantsev, P. Hawrylak, Solid State Commun. 152 (2012) 909-913. [22] D.Y. Qiu, F.H. da Jornada, S.G. Louie, Phys. Rev. Lett. 111 (2013) 216805. [23] H.R. Guti´errez, N. Perea-L´opez, A.L. El´ıas, A. Berkdemir, B. Wang, R. Lv, F. L´opez-Ur´ıas, V.H. Crespi, H. Terrones, M. Terrones, Nano Lett. 13 (2013) 3447-3454. [24] X. Xu, W. Yao, Di Xiao, T.F. Heinz, Nat. Phys. 10 (2014) 343-350. [25] S. Helveg, J.V. Lauritsen, E. Lægsgaard, I. Stensgaard, J.K. Nørskov, B.S. Clausen, H. Topsøe, F. Besenbacher, Phys. Rev. Lett. 84 (2000) 951-954. [26] M.V. Bollinger, J.V. Lauritsen, K.W. Jacobsen, J.K. Nørskov, S. Helveg, F. Besenbacher, Phys. Rev. Lett. 87 (2001) 196803. [27] M.V. Bollinger, K.W. Jacobsen, J.K. Nørskov, Phys. Rev. B 67 (2003) 085410. [28] T. Li, G. Galli, J. Phys. Chem. C 111 (2007) 16192-16196. [29] D. Davelou, G. Kopidakis, G. Kioseoglou, I.N. Remediakis, Solid State Commun. 192 (2014) 42-46. [30] K.-A. Duerloo, M.T. Ong, E.J. Reed, J. Phys. Chem. Lett, 3 (2012) 28712876. [31] H. Zhu, Y. Wang, J. Xiao, M. Liu, S. Xiong, Z.J. Wong, Z.Ye, Y. Ye, X. Yin, X. Zhang, Nat. Nanotechnol. 10 (2015) 151-155.

6