First principles studies on electronic and transport properties of edge contact graphene-MoS2 heterostructure

Computational Materials Science 133 (2017) 137–144

Contents lists available at ScienceDirect

Computational Materials Science journal homepage: www.elsevier.com/locate/commatsci

First principles studies on electronic and transport properties of edge contact graphene-MoS2 heterostructure Jie Sun a, Na Lin a,⇑, Cheng Tang a, Haoyuan Wang a, Hao Ren b, Xian Zhao a,⇑ a b

State Key Laboratory of Crystal Materials, Shandong University, 250100 Jinan, Shandong, PR China State Key Laboratory of Heavy Oil Processing & Center for Bioengineering and Biotechnology, China University of Petroleum (East China), 266580 Qingdao, PR China

a r t i c l e

i n f o

Article history: Received 22 December 2016 Received in revised form 27 February 2017 Accepted 2 March 2017

Keywords: First principles Lateral heterostructure Projected densities of states Band alignment Transmission spectra

a b s t r a c t Nanodevice based on MoS2 channel lateral connecting with graphene electrode was fabricated in recent experiment. In present paper, first principles calculations are carried out to reveal the relationship between contact geometries and electrical properties of graphene-MoS2 heterostructure. Four different contact edges between graphene and MoS2, namely, Armchair-Armchair, Zigzag-Armchair, ArmchairZigzag, Zigzag-Zigzag, are investigated. Calculations indicate that MoS2 will be metalized as a consequence of junction formation with graphene. The metallic states located at Fermi level are mainly laid at the contact interface and dominated by 4d states of Mo atom as well as 2p states of both S and C atoms. Different contact geometries of graphene-MoS2 result in different charge transfer values in contact interfaces. Investigation on band alignments reveals that n-type Schottky contacts are formed in four graphene-MoS2 lateral heterostructures with barrier heights of 0.45–0.75 eV, which are larger than those of edge contact with Sc and Ti metals. The transmission gap of each configuration obtained using a two-probe system is unexpectedly larger than the intrinsic band gap of MoS2. The discrepancies of current-voltage behavior in two represented configurations demonstrate that contact geometries play an important role in electronic transport properties of graphene-MoS2 junctions. Ó 2017 Elsevier B.V. All rights reserved.

1. Introduction Since the discovery of graphene, two-dimensional (2D) materials have attracted enormous experimental and theoretical interests [1–4]. Among them, the semiconductor transition metal dichalcogenides (SCTMDs), such as MoS2, WS2, are promising materials due to their distinctive electronic and optical properties [5–10]. Unlike zero band gap nature of graphene, the occurrence of moderate band gaps makes them potential candidates for nanodevice applications. Field-effect transistors (FETs) based on MoS2 exhibit large on/off ratios and near theoretical subthreshold swing values [5,11,12]. However, further improvement of the devices performance have been limited by relatively low carrier mobility, which mainly arise from the contact with metal electrode generating large contact resistance owing to the formation of Schottky barrier [13–16]. Therefore, great efforts have been made to optimize electrode contacts with MoS2 [17–20]. Top contact and edge contact are two fundamental interface geometries between electrodes and 2D materials. Conventional methods usually use 3D metal electrodes to top contact with 2D ⇑ Corresponding authors. E-mail addresses: [email protected] (N. Lin), [email protected] (X. Zhao). http://dx.doi.org/10.1016/j.commatsci.2017.03.004 0927-0256/Ó 2017 Elsevier B.V. All rights reserved.

MoS2. Ohmic contacts can be achieved by choosing metals with suitable work functions for MoS2. Nevertheless, it is difficult to obtain stable and controllable Ohmic contact in experiments. For example, Ohmic and Schottky contacts were found in the MoS2 transistors with Au electrodes by different groups, respectively [5–21]. Theoretical studies have verified that edge contacts between 3D metal electrodes and 2D MoS2 show more advantageous features compared with top contacts, owing to stronger orbital overlaps and better reduction of tunnel barriers [22]. Apart from 3D metal electrodes, 2D metallic TMDs and graphene, which can provide much smaller volumes, are also fabricated as electrode contacted with MoS2. Graphene top contacts perform low contact resistances when they have sufficiently large contact areas [23,24], however, due to the van der Waals gap between graphene and MoS2 [25], the contact resistance increases dramatically as the length of the graphene top contact is reduced below the transfer length to the tens of nm scale [26]. The metallic TMDs electrode forming edge contact with SCTMD can also generate smaller resistance [27,28], however, the metallic TMDs phase is metastable and this is not beneficial for its practical applications. In recent experiments, the new lateral heterostructures, graphene-MoS2 [26], are synthesized in spite of that the two atomic layers show significant crystallographic dissimilarity and

138

J. Sun et al. / Computational Materials Science 133 (2017) 137–144

possess large lattice mismatch (
25%). Moreover, the MoS2-based devices in-plane stitched with graphene electrodes are also fabricated [29]. New graphene edge contact devices show low contact resistance, ohmic behavior at room temperature and the least possible additional volume to the devices. 1D edge contact exhibits significant potential for MoS2 based future device applications. However, there is still a lack of knowledge of the precise edge morphology of the two materials stitched together, which may have

great impact on the electronic and transport properties of the new in-plane heterostructure. Therefore, it is necessary to explore the relationship between the contact geometries and electrical properties of the heterostructure, which may contribute to practical device applications. In this paper, we report a first-principle study on electronic and transport properties of edge contact between graphene and MoS2. We construct four different edge contact configurations of

Fig. 1. Four edge contact geometries of graphene-MoS2 heterostructures, (a) armchair-armchair, (b) zigzag-armchair, (c) armchair-zigzag, and (d) zigzag-zigzag. Shaded parts show the contact regions.

Fig. 2. The relative energy of (a) A-A, (b) Z-A, (c) A-Z, and (d) Z-Z junctions with different lattice constant along the direction vertical to the contact edge. The energy of the most stable configuration has been set to zero.

J. Sun et al. / Computational Materials Science 133 (2017) 137–144

graphene-MoS2 heterostructures and determine the most stable configuration for each heterostructure. Based on the analyses of charge distributions, density of states (DOS), and band alignment, the nature of the graphene-MoS2 heterostructures is revealed. And then, a two-probe system, MoS2 center with graphene electrodes, is applied to clarify the effect of different edge contacts on the electronic transport properties.

2. Computational details The first principles calculations were carried out based on density functional theory (DFT) with the projector-augmented wave (PAW) method, as performed in the Vienna ab initio simulation package (VASP) [30]. The exchange and correlation interactions were treated by generalized gradient approximation (GGA) functional with Perdew, Burke and Ernzerhof (PBE) [31] parameterization. A plane-wave basis set with kinetic energy cutoff of 400 eV

139

was used and the first Brillouin zone was sampled with a 6  4  1 Monkhorst-Pack grid. All structures were relaxed by using the conjugate gradient method until the maximum Hellmann-Feynman forces exerted on each atom less than 0.02 eV/Å. To avoid interaction between neighboring images, a vacuum layer larger than 15 Å was imposed. The magnetism is not considered and spin-unpolarized computations have been performed in our work. The electron transport properties were investigated based on the nonequilibrium Green’s function (NEGF) method implemented in the TRANSIESTA program [32]. The current through the system is calculated using the Landauer–Büttiker formula:

I¼

2e h

Z

TðE;VÞ½f L ðEÞ  f R ðEÞdE

where T(E,V) is the transmission coefficient at energy E and the bias voltage V, fL(E) and fR(E) are the Fermi distribution functions of left and right electrodes. Single-f plus polarization basis set was used

Fig. 3. Isosurface plot of the electron charge density difference for (a) A-A, (b) Z-A, (c) A-Z, and (d) Z-Z. The charge accumulation is represented in yellow and charge depletion is in blue. (e)–(h) Charge transfer value based on Bader charge analysis of A-A, Z-A, A-Z and Z-Z, respectively. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

140

J. Sun et al. / Computational Materials Science 133 (2017) 137–144

and a mesh cutoff was set to be 150 Ry. The k-point samplings for the transmission spectra calculations were 1, 50, and 50 in the x, y, z directions, respectively. 3. Results and discussion At the beginning of our calculations, the unit cell of graphene and MoS2 were optimized with lattice constant of 2.46 Å and 3.18 Å, respectively, which were consistent with previous studies [33,34]. It is known that both graphene and MoS2 can be designed as armchair or zigzag edges, respectively. Therefore, upon formation of a lateral junction between two materials, four different edge contact configurations of graphene-MoS2 heterostructures can be formed, namely, armchair graphene edge contacts with armchair MoS2 edge (A-A) (Fig. 1(a)), zigzag graphene edge contacts with armchair MoS2 edge (Z-A) (Fig. 1(b)), armchair graphene edge contacts with zigzag MoS2 edge (A-Z) (Fig. 1(c)), and zigzag graphene edge contacts with zigzag MoS2 edge (Z-Z) (Fig. 1(d)). We connect 4,9,3,4 graphene rectangular unit cell with 3,4,4,3 MoS2 rectangular unit cell to obtain A-A, Z-A, A-Z, Z-Z configurations with dimensional along contact edge (x direction) 17.08 Å, 22.19 Å, 12.81 Å, 9.86 Å, respectively. The supercell lattice mismatch between graphene and MoS2 are 3.2%, 0.6%, 0.6%, 3.2% for A-A, Z-A, A-Z, Z-Z, respectively. Here, we randomly place graphene and MoS2 ribbons to form the lateral structure. Therefore, the impacts on the electronic properties of this heterostructure induced by the relative position of two ribbons are not considered. No atom defects are introduced in the contact region. As we construct the supercell along y direction, we note that each supercell actually contains two contact edges to satisfy it periodic. In order to ensure there is

no interaction between two edges, we choose supercell of both graphene and MoS2 along y direction with the size larger than 15 Å. The supercells of A-A, Z-A, A-Z and Z-Z along y direction consist of 9, 4, 7, 4 graphene primitive cells connected with 7, 5, 3, 3 MoS2 primitive cells, respectively. Both x and y directions have been satisfied with periodic boundary conditions in our calculations. The most stable configuration of each heterostructure has been determined by calculating the total energy corresponding to the different supercell lattice along the direction vertical to the 1D contact edge (Fig. 2(a)–(d)), considering the occurrences of reconstructions in both graphene and MoS2 edges [35–37]. We individually obtain a local minimum for respectively Z-A, A-Z and Z-Z configurations corresponding to the lattice constant of 35.28 Å, 33.90 Å, and 33.72 Å and consider it as the most stable for each case. However, for A-A structure, only a critical lattice constant (45.07 Å) is obtained and the supercell is destroyed with all atoms out of order as the lattice constant below this value after optimization. This is mainly caused by the enlarged stress arising as decreasing the lattice constant along y direction. Hence, we treat this critical energy as the most stable case. After careful checking the structures of four stable configurations, it is found that graphene and MoS2 are connected with each other through covalent bonds. In addition, it is clearly seen that the hybrid octagon consisting of C, Mo and S atoms exists in A-A structure besides hexagonal rings. A partial ordered pentagons and heptagons arrange in Z-A structure. For A-Z and Z-Z structures, the dominated polygons are mainly not larger than six-membered ring. The discrepancies in the contact region among these heterostructures are expected to lead different electronic structure and transport behavior.

Fig. 4. Total and selected orbital projected densities of states of (a) A-A, (b) Z-A, (c) A-Z, and (d) Z-Z, respectively.

J. Sun et al. / Computational Materials Science 133 (2017) 137–144

141

Fig. 5. The projection of the total DOS on graphene and MoS2 for (a) A-A, (b) Z-A, (c) A-Z, (d) Z-Z. Gra_c (MoS2_c) and Gra_nc (MoS2_nc) represent graphene (MoS2) at and far away from the contact interface, respectively.

In the following part, we examine the charge transfer between MoS2 and graphene by the isosurface plot of electron charge density difference (Fig. 3(a)–(d)). It can be clearly seen that the charge accumulation and depletion mainly occurs at the interface of two materials stitched together for each configuration. And we can see the charge transfer increases in the interface bonding region. Hence the electrons can transfer easily between graphene and MoS2. The charge transfer values are evaluated by Bader charge analysis [38] (Fig. 3(e) and (f)), which can give the charge of each atom and the total charge is calculated by summing over all atoms. It is found that MoS2 and graphene play the role as electron donor and acceptor, respectively. The charge transfer values from MoS2 to graphene are 0.159, 0.153, 0.182 and 0.156 e/Å (total transfer charge value/interface length) for A-A, Z-A, A-Z and A-A, respectively. Specifically, C atoms, which link with Mo atoms, show more obvious charge accumulation than that link with S atoms. The charge depletion value of Mo atoms in contact region is larger than that in far from contact region. Similarly, S atoms in the contact region are populated by fewer electrons compared with those in far from contact region. The electronic structures of these heterostructures are stated by plotting the total and selected orbital projected densities of states (PDOS) as shown in Fig. 4(a)–(d). For all configurations, a remarkable characteristic is that the occurrences of high density of states located around Fermi level. This indicates the semiconductor MoS2 is metalized as contact with graphene due to strong interaction between them in contact region. This phenomenon is similar to MoS2 edge contacting with metals, which is also become metallic due to stronger orbital overlaps and induce high density of states

in the original band gaps [20–22]. As can be seen from the PDOS, these metallic states are mainly constituted by 4d states of Mo atoms coupled with 2p states of S and C atoms. However, the proportion of these states making contribution to the metallic sates around the Fermi level is not uniform for each configuration due to the different contact geometries. From the enlarged PDOS of selected orbital around Fermi level (inset pictures in Fig. 4(a)– (d)), we note that, as for A-A, Z-A and Z-Z configurations, the states at Fermi level are mainly dominated by 4d states of Mo atoms, while for A-Z configuration, the 2p states of S and C atoms nearly have the same contribution as Mo 4d states to the composition of the total states. In addition, the 4d states of Mo atoms toward Fermi levels show more delocalization in A-A, A-Z and Z-Z configurations compared with that in Z-A configuration. This may be a sign of a low resistance Ohmic contact for A-A, A-Z and Z-Z. In order to give further insight into the electronic structures of these heterostructures, we project the total DOS on graphene and MoS2, respectively (Fig. 5(a)–(d)). MoS2 (graphene) at the contact interface and far away from contact interface are represented as MoS2_c (Gra_c) and MoS2_nc (Gra_nc), respectively. We can see that the metallic states at the Fermi level are primarily caused by the MoS2 and graphene in contact interface region due to their direct interaction through covalent bonds and strong states hybridation with each other. MoS2 far away from the contact interface can well preserve its semiconductor nature; meanwhile, the bottom of its conduction band and the top of its valence band shift in various degrees for each configuration probably due to the charge transfer between the contact and noncontact regions.

142

J. Sun et al. / Computational Materials Science 133 (2017) 137–144

In comparison with electronic structure, band offset of the heterostructure is also very important in material and device design [39]. Their precise knowledge is extremely important to engineer electronic and optoelectronic devices. The band alignments of lateral graphene-MoS2 systems are calculated by using the macroscopic averaging method [40]. The macro average electrostatic potential of graphene-MoS2 heterostructure along y direction can be obtained by using the equation

1  VðyÞ ¼ a

Z

yþ2a

y2a

 0 Þdy0 Vðy

 where a is the period length along y direction, VðyÞ is the electroRR  Vðx; y; zÞdxdz, S static potential averaged on xz-plane, VðyÞ ¼ 1S represents the area of the supercell in the xz-plane and Vðx; y; zÞ is calculated by DFT calculations. The macro average electrostatic potential is chosen as a reference, the Fermi level of graphene and the valence band maximum (VBM) and conduction band minimum (CBM) of MoS2 with respect to the electrostatic potential are calculated using their individual supercell. The calculated results are shown in Fig. 6. It can be seen that the macroscopic average potential is discontinuous at the interface of graphene and MoS2, indicating that electrons should overcome a potential barrier to cross from MoS2 section to graphene section. Our PBE calculations indicate that both VBM and CBM of MoS2 locate below the Fermi level of graphene in the lateral heterostructures, that is, CBM of MoS2 is more close to the Fermi level of graphene, which will result in n-type characteristic of

MoS2 FETs. The n-type Schottky contacts also have been predicted in MoS2 edge connecting with Sc and Ti metals, possessing barrier height
0.11 eV and 0.39 eV, respectively [20]. The enlarged conduction band offset (CBO) of graphene-MoS2 has been found in our calculations with the values of 0.46 eV, 0.74 eV, 0.66 eV and 0.64 eV for A-A, Z-A, A-Z, and Z-Z configuration, respectively. The differences of the CBO among four heterostructures are owing to the discrepancies of their contact geometries. We note that each Mo atom in A-A structure bonds with two C atoms, while for other structures, the Mo atoms mainly bond with one C atom. The strong hybrid between Mo and C atom in A-A structure may be the reason that lower band offset is displayed for it. Larger CBO value is easier to hinder the electrons across the graphene-MoS2 interface, which is considered to introduce larger contact resistance. Therefore, we predict A-A and Z-A contact between graphene and MoS2 probably possess the lowest and highest contact resistance among four configurations. Following the studies of electronic structures of graphene-MoS2 complexes, their electron transport properties are investigated by using a two-probe system, as shown in Fig. 7(a), where two semi-infinite electrodes are connected with the scattering region. The calculated transmission spectrum of each configuration at zero bias is shown in Fig. 7(b). It is unexpected that the transmission gap of each configuration (
3 eV) is larger than the intrinsic band gap of pristine MoS2 (
1.6 eV), which is mainly due to the occurrence of band alignments between graphene and MoS2. Previous theoretical studies have mentioned that the transport gap is a sum over electrons and holes Schottky barrier heights [41]. We

Fig. 6. Band alignments of the lateral graphene-MoS2 systems, (a) A-A, (b) Z-A, (c) A-Z, and (d) Z-Z.

J. Sun et al. / Computational Materials Science 133 (2017) 137–144

143

A-Z mildly rises up while that of Z-Z increases more rapidly. This leads the current of Z-Z exceeded that of A-Z as the bias beyond 0.5 V. This different I-V behavior between A-Z and Z-Z is probably due to the different charge transfer across the interface. We note that Bader analysis charge transfer value between graphene and MoS2 in A-Z configuration is larger than that of Z-Z configuration. More charge accumulation or depletion at the contact interface may give rise to electron scattering and is not beneficial for the electrons transport through the contact interface. A clear understanding on the relation between the contact geometry and electronic transport property contributes to practical device applications. In addition, based on these, further efforts can be made to optimize the contact interface and improve the transport behaviors. 4. Conclusions In summary, we have systematically investigated the electronic and transport properties of graphene-MoS2 systems with different contact interfaces. As MoS2 in-plane connecting with graphene, it is found that the electrons transfer from MoS2 to graphene. In addition, MoS2 is metalized due to the occurrence of band gap states of MoS2-graphene heterostructure. These band gap states are mainly localized at the contact interface and dominated by 4d states of Mo atom and 2p states of S and C atoms. Further band alignment investigation shows that the CBM of MoS2 is more close to the Fermi level of graphene with values of 0.46 eV, 0.74 eV, 0.66 eV and 0.64 eV for A-A, Z-A, A-Z, and Z-Z configurations, respectively. Studies on electronic properties clearly show that the carrier transport across a graphene-MoS2 interface is mainly dominated by the detailed structural properties at the interfaces such that the Schottky barrier heights and charge distribution are largely controlled by the contact structure. Hence the contact structure should be carefully considered in order to obtain a low contact resistance. Acknowledgements We acknowledge the National Natural Science Foundation of China (Grant Nos. 21573129 and 21403300), the Natural Science Foundation of Shandong Province (Grant No. ZR2015BQ001), and the General Financial Grant from the China Postdoctoral Science Foundation (Grant No. 2013M531595). The authors also acknowledge a generous grant of computer time from the National Supercomputer Center in Tianjin-TianHe-1(A) and the Norwegian Programme for Supercomputing. References Fig. 7. (a) A Two-probe system where semi-infinite left and right electrode regions are in contact with the central scattering region. (b) Transmission spectrum of each configuration at zero bias. (c) The I-V curve of A-Z and Z-Z configurations.

note that the sum of VB and CB band alignment in our calculations are 2.57 eV, 3.13 eV, 2.97 eV, 2.93 eV for A-A, Z-A, A-Z and Z-Z, respectively, which are consistent with the transmission band gaps. In addition, the locations of VBM (CBM) of four configurations in transmission spectrum show discrepancies due to different contact geometries. To further elucidate the relation between contact geometries and electronic transport properties, we calculate the currentvoltage (I-V) curves for A-Z and Z-Z in low bias (Fig. 7(c)). Unfortunately, A-A and Z-A is not involved here due to their large-size systems (350 atoms for A-A and 384 atoms for Z-A) in order to calculate the electronic transport properties, which are indeed beyond our calculation ability. We can see that the current of

[1] K.S. Novoselov, A.K. Geim, S.V. Morozov, D. Jiang, Y. Zhang, S.V. Dubonos, I.V. Grigorieva, A.A. Firsov, Science 306 (2004) 666–669. [2] A.K. Geim, K.S. Novoselov, Nat. Mater. 6 (2007) 183–191. [3] A.H. Castro Neto, F. Guinea, N.M.R. Peres, K.S. Novoselov, A.K. Geim, Rev. Mod. Phys. 81 (2009) 109–162. [4] S. Das Sarma, S. Adam, E.H. Hwang, E. Rossi, Rev. Mod. Phys. 83 (2011) 407– 470. [5] B. Radisavljevic, A. Radenovic, J. Brivio, V. Giacometti, A. Kis, Nat. Nanotech. 6 (2011) 147–150. [6] K.F. Mak, C. Lee, J. Hone, J. Shan, T.F. Heinz, Phys. Rev. Lett. 105 (2010) 136805. [7] G. Eda, H. Yamaguchi, D. Voiry, T. Fujita, M. Chen, M. Chhowalla, Nano Lett. 11 (2011) 5111–5116. [8] Y. Yoon, K. Ganapathi, S. Salahuddin, Nano Lett. 11 (2011) 3768–3773. [9] Q.H. Wang, K. Kalantar-Zadeh, A. Kis, J.N. Coleman, M.S. Strano, Nat. Nanotech. 7 (2012) 699–712. [10] A.M. Jones, H. Yu, N.J. Ghimire, S. Wu, G. Aivazian, J.S. Ross, B. Zhao, J. Yan, D.G. Mandrus, D. Xiao, W. Yao, X. Xu, Nat. Nanotech. 8 (2013) 634–638. [11] B. Radisavljevic, M.B. Whitwick, A. Kis, ACS Nano 5 (2011) 9934–9938. [12] H. Wang, L. Yu, Y.H. Lee, Y. Shi, A. Hsu, M.L. Chin, L.J. Li, M. Dubey, J. Kong, T. Palacios, Nano Lett. 12 (2012) 4674–4680. [13] S. Das, H.Y. Chen, A.V. Penumatcha, J. Appenzeller, Nano Lett. 13 (2013) 100– 105. [14] I. Popov, G. Seifert, D. Tománek, Phys. Rev. Lett. 108 (2012) 156802.

144

J. Sun et al. / Computational Materials Science 133 (2017) 137–144

[15] J.R. Chen, P.M. Odenthal, A.G. Swartz, G.C. Floyd, H. Wen, K.Y. Luo, R.K. Kawakami, Nano Lett. 13 (2013) 3106–3110. [16] F. Ahmed, M.S. Choi, X. Liu, W.J. Yoo, Nanoscale 7 (2015) 9222–9228. [17] W. Chen, E.J.G. Santos, W. Zhu, E. Kaxiras, Z. Zhang, Nano Lett. 13 (2013) 509– 514. [18] H. Liu, M. Si, Y. Deng, A.T. Neal, Y. Du, S. Najmaei, P.M. Ajayan, J. Lou, P.D. Ye, ACS Nano 8 (2014) 1031–1038. [19] C. Gong, L. Colombo, R.M. Wallace, K. Cho, Nano Lett. 14 (2014) 1714–1720. [20] Z. Li, X. Li, Jinlong. Yang, ACS Appl. Mater. Interfaces 7 (2015) 12981–12987. [21] M. Fontana, T. Deppe, A.K. Boyd, M. Rinzan, A.Y. Liu, M. Paranjape, P. Barbara, Sci. Rep. 3 (2013) 1634. [22] J. Kang, W. Liu, D. Sarkar, D. Jena, K. Banerjee, Phys. Rev. X 4 (2014) 031005. [23] X. Cui, G. Lee, Y.D. Kim, G. Arefe, P.Y. Huang, C. Lee, D.A. Chenet, X. Zhang, L. Wang, F. Ye, F. Pizzocchero, B.S. Jessen, K. Watanabe, T. Taniguchi, D.A. Muller, T. Low, P. Kim, J. Hone, Nat. Nanotech. 10 (2015) 534–540. [24] Y. Liu, H. Wu, H. Cheng, S. Yang, E. Zhu, Q. He, M. Ding, D. Li, N. O. Weiss, Y. Huang, Xiangfeng Duan, Nano Lett. 15 (2015) 3030–3034. [25] Y. Ma, Y. Dai, M. Guo, C. Niu, B. Huang, Nanoscale 3 (2011) 3883–3887. [26] X. Ling, Y. Lin, Q. Ma, Z. Wang, Y. Song, L. Yu, S. Huang, W. Fang, X. Zhang, A.L. Hsu, Y. Bie, Y. Lee, Y. Zhu, L. Wu, J. Li, P. Jarillo-Herrero, M. Dresselhaus, T. Palacios, J. Kong, Adv. Mater. 28 (2016) 2322–2329. [27] R. Kappera, D. Voiry, S.E. Yalcin, B. Branch, G. Gupta, A.D. Mohite, M. Chhowalla, Nat. Mater. 13 (2014) 1128–1134.

[28] H. Yamaguchi, J. Blancon, R. Kappera, S. Lei, S. Najmaei, B.D. Mangum, G. Gupta, P.M. Ajayan, J.L.M. Chhowalla, J.J. Crochet, A.D. Mohite, Nano Lett. 9 (2015) 840–849. [29] M.H.D. Guimarães, H. Gao, Y. Han, K. Kang, S. Xie, C. Kim, D.A. Muller, D.C. Ralph, J. Park, Nano Lett. 10 (2016) 6392–6399. [30] G. Kresse, J. Hafner, Phys. Rev. B 47 (1993) 558–561. [31] J.P. Perdew, K. Burke, M. Ernzerhof, Phys. Rev. Lett. 77 (1996) 3865–3868. [32] M. Brandbyge, J.L. Mozos, P. Ordejón, J. Taylor, K. Stokbro, Phys. Rev. B 65 (2002) 165401. [33] L. Yuan, Z.Y. Li, J.L. Yang, J.G. Hou, Phys. Chem. Chem. Phys 14 (2012) 8179– 8184. [34] J. Kang, S. Tongay, J. Zhou, J. Li, J. Wu, Appl. Phys. Lett. 102 (2013) 012111. [35] J.A. Spirko, M.L. Neiman, A.M. Oelker, K. Klier, Surf. Sci. 542 (2003) 192–204. [36] W. Gao, J. Shu, C. Qiu, Q. Xu, ACS Nano 6 (9) (2012) 7806–7813. [37] G. Levita, P. Restuccia, M.C. Righi, Carbon 107 (2016) 878–884. [38] E. Sanville, S.D. Kenny, R. Smith, G. Henkelman, J. Comput. Chem. 28 (2007) 899–908. [39] A. Franciosi, C.G. Van de Walle, Surf. Sci. Rep. 25 (1996) 1–140. [40] A. Baldereschi, S. Baroni, R. Resta, Phys. Rev. Lett. 61 (1988) 734–737. [41] Y. Pan, S. Li, M. Ye, R. Quhe, Z. Song, Y. Wang, J. Zheng, F. Pan, W. Guo, J. Yang, J. Lu, J. Phys. Chem. C 120 (2016) 13063–13070.