This is the submitted version of the article: Guster B., Robles R., Pruneda M., Canadell E., Ordejón P.. 2 × 2 charge density wave in single-layer TiTe2. 2D Materials, (2019). 6. 015027: - . 10.1088/2053-1583/aaf20b. Available at: https://dx.doi.org/10.1088/2053-1583/aaf20b 2D Mater. 2 ×2 Charge Density Wave in single-layer TiTe 2 1 1 1 Bogdan Guster , Roberto Robles , Miguel Pruneda , Enric 2 1 Canadell , Pablo Ordejón 1 Catalan Institute of Nanoscience and Nanotechnology (ICN2), CSIC and BIST, Campus UAB, Bellaterra, 08193 Barcelona, Spain 2 Institut de Ciència de Materials de Barcelona (ICMAB-CSIC), Campus Bellaterra, 08193 Barcelona, Spain E-mail: pablo.ordejon@icn2.cat Abstract. A density functional theory study concerning the origin of the recently reported 2 ×2 charge density wave (CDW) instability in single-layer TiTe 2 is reported. It is shown that, whereas calculations employing the semi-local functional PBE favor the undistorted structure, the hybrid functional HSE06 correctly predicts the 2 ×2 distortion, . The study suggests that the magnitude of the semi-metallic overlap between the valence band top at Γ and the conduction band bottom at M is a key factor controlling the tendency towards the distortion. It is also shown that tensile strain stabilizes the 2 ×2 CDW, and we suggest that this fact could be further used to induce the instability in double-layers of TiTe 2 which in the absence of strain remain undistorted. The driving force for the CDW instability seems to be the same phonon mediated mechanism acting for single- layer TiSe2 , although in single-layer TiTe 2 the driving force is smaller, and the semimetallic character is kept below the transition temperature. Keywords : Single layer dichalcogenides, charge density waves, density functional theory, exchange-correlation funtionals, titanium ditelluride, 2D materials. 2×2 Charge Density Wave in single-layer TiTe 2 2 1. Introduction successfully predict that the 2 ×2 CDW structure is more stable than the non-distorted structure for single- Transition metal dichalcogenides of the groups IV layer TiSe2 [18, 19, 13]. Overall, these observations and V rank among the most controversial materials suggested the hypothesis that something really new exhibiting charge density wave (CDW) instabilities [1, and challenging is at work in single-layer TiTe 2 [8]. 2]. The possibilities of strong or weak electron- However, one should note that the CDW phonon coupling scenarios in group V 2 H -MX2 (M transition occurs at 100 K in single-layer TiTe 2 [8] = Nb, Ta; X = S, Se) and either phonon mediated but at a considerably higher temperature, 232 K, in or excitonic mechanisms in group IV 1 T -TiSe2 have single-layer TiSe2 [11]. Hence, the driving force for the been discussed for decades [1]. Many of these distortion must be considerably weaker in single-layer systems also exhibit superconductivity (SC) under TiTe2 . Before concluding that a new scenario is needed certain conditions and the competition between the to grasp the origin of the unexpected 2 ×2 CDW in this two instabilities remains an important question still material, one should wonder about the appropriateness unanswered [3, 4]. These materials are built from MX 2 of the so far successful GGA-type of DFT approaches layers interacting through weak van der Waals forces to the CDW instabilities in single-layer group IV and thus are easily exfoliated [5]. Consequently, they and V dichalcogenides. Such an appraisal is needed oer the possibility to examine the above mentioned because it impinges on very fundamental questions issues at the two-dimensional (2D) limit as well as by concerning CDW instabilities at the 2D limit. Note smoothly varying the density of carriers through gate that it has been recently shown [20] that GGA-type doping. This is at the origin of the huge revival of functionals like PBE [21] overestimate the overlap interest recently raised by these materials [6, 7, 8, 9, between the Ti 3d and Se 4p levels in bulk 1T -TiSe2 . 10, 11]. This can be corrected by using hybrid functionals Indeed, intriguing dierences of these few-ake or like HSE06 [22, 23] leading to an improvement of even single-layer materials with their bulk counterparts the electronic description of bulk TiSe 2 [20]. In the have been discovered. Recent reports on the existence following we report a DFT study of the likeliness of a of a very weak pseudo-gap at the Fermi energy in 2×2 CDW in single-layer TiTe 2 employing both PBE single-layer NbSe2 [6, 12] or the possible occurrence and HSE06 functionals which provides useful insight on of incommensurate modulations for slightly electron the origin of the CDW instability in single-layer TiTe 2 . doped TiSe2 crystals of thicknesses less than 10 nm [10, 13] make clear that we are still far from a full 2. Results and Discussion understanding of the physics of CDW materials and more particularly when the screening is reduced. An isolated layer of 1T -TiTe2 is made of an hexagonal In this context, the recent report of a 2 ×2 CDW lattice of Ti atoms in an octahedral environment of in single-layer TiTe2 by Chen et al. [8] came as a Te atoms (Fig. 1a). The repeat unit of the hexagonal very intriguing surprise. Since long [14, 15, 16] it has bulk crystal structure contains just one of these layers. been known that bulk 1T -TiTe2 does not exhibit the A detailed description of our calculation method is 2×2 CDW that occurs in isostructural 1 T -TiSe2 [17] presented in Appendix A. Let us start our analysis In addition, the 2×2 CDW is not observed anymore by briey considering the PBE description of the in double-layer TiTe2 [8]. In contrast, the 2×2 CDW electronic structure. The optimized a cell parameter instability is observed in TiSe 2 from the single- is 3.804 Å for the single layer and 3.815 Å for the layer, for ultrathin lms with up to six layers [18], bulk, in good agreement with the bulk experimental and for the bulk crystal [17]. The experimental value of 3.777 Å [24]. Within a single-layer there indication of the occurrence of the CDW in single- are several Te...Te contacts shorter than the sum of layer TiTe2 is even more surprising when considering the van der Waals radii so that the valence bands, that rst-principles density functional theory (DFT) which have their maximum at Γ and are mostly built calculations found that single-layer TiTe 2 shows no from Te 5p orbitals, are considerably wide and overlap tendency to distort towards the 2 ×2 CDW structure with the bottom part of the Ti 3 d bands, which have at the generalized gradient approximation (GGA) their minima at the M point (Fig. 1b). As expected level [8]. Yet, calculations of the same quality from the fact that in the bulk there are short Te...Te 2×2 Charge Density Wave in single-layer TiTe 2 3 Figure 1. Single-layer TiTe2 . (a) Side and top views of a layer. (b) Calculated GGA band structure with (black lines) and without (red dots) spin-orbit coupling. Γ = (0, 0, 0), M= (1/2, 0, 0) and K= (1/3, 1/3, 0) in units of the reciprocal hexagonal lattice vectors. (c) GGA phonon dispersion in the Γ - M segment of the Brillouin Zone (BZ) for the undistorted structure. (d) Frozen-phonon total energy calculation as a function of the soft phonon mode amplitude at M, calculated with the GGA functional. Energies are given in meV per formula unit relative to the undistorted phase. All results have been obtained using the PBE functional. contacts in the direction perpendicular to the layer, the Table 1. Evolution with tensile strain of: (i) the energy semimetallic overlap is 19% larger in the bulk. Note dierence between the undistorted structure and the relaxed that, as it is clear from Fig. 1b, inclusion of spin-orbit 2×2 CDW structure, and (ii) displacement of the Ti atoms coupling eects does not have any noticeable eect in the plane parallel to the layer from their position in the and therefore it will not be considered anymore in the undistorted phase. All values obtained using the PBE functional. following. Thus, according to the PBE calculations single-layer TiTe2 is a semimetal exactly as the bulk. Tensile Strain (%) ∆E (meV/f.u.) ∆ Ti(Å) In contrast with the case of single-layer TiSe 2 , for 0 0.00 0.000 which the same type of calculations led to a phonon 1 -0.001 0.007 with imaginary frequency at the M point [13], our 2 -0.09 0.057 calculations with the GGA functional for single-layer 3 -0.29 0.061 TiTe2 (Fig. 1c) show no phonons with imaginary 4 -0.47 0.067 frequency, in agreement with those of Chen et al. [8]. 5 -0.74 0.080 Although there is some remnant of the instability at M 6 -1.27 0.095 (notice the optical branch that disperses downwards 7 -2.04 0.118 and shows the lowest frequency at the M point, reminiscent of the mode that becomes unstable for strain (allowing the atomic positions to relax for each TiSe2 ), there is no denite indication of a phonon strain applied). The strain is dened as s = δm/m0 instability that may lead to the 2 ×2 distortion of where m0 is the unstrained cell parameter and δm+m0 the structure. Detailed structural optimizations of the strained cell parameter. Thus, positive values 2×2 supercells conrmed this result. However, as correspond to tensile strains. As shown in Table 1, indicated by the frozen-phonon total energy calculation a tensile stress as small as 1% is sucient to make as a function of the soft phonon mode amplitude at the 2×2 CDW structure slightly more stable than M (Fig. 1d), the GGA potential energy surface is the undistorted one. For strains larger than 2%, the extremely at. Under such circumstances, even if 2×2 CDW is denitely favored. When this happens, strictly speaking the PBE calculations disagree with energy gaps open at the crossings of the folded the experimental results in that no tendency toward valence and conduction bands although, as shown in the 2×2 CDW distortion is found, the results are Fig. 2, the distorted structure is still semimetallic somewhat inconclusive and a closer look is needed. and only for relatively large tensile strains (between The very at frozen-phonon energy curve of 5 and 6%) a clear band gap occurs around the Fermi Fig. 1d) suggests that small external perturbations level. Thus, it appears that GGA-PBE description could be able to change the relative stability of the of single-layer TiTe2 is only consistent with the the undistorted and 2×2 CDW structures. This could real situation when a slight tensile strain in imposed occur, for instance, by the eect of strain. Thus, in the calculation. The main eect of the strain we studied the evolution of the band structure and is a decrease of the intralayer Te...Te short contacts the relative stability of the undistorted 1 ×1 and which leads to a decrease of the Te 5 p bandwidth 2×2 CDW structures as a function of biaxial tensile 2×2 Charge Density Wave in single-layer TiTe 2 4 Figure 2. Calculated GGA band structures for the stable 2 ×2 CDW structure of single-layer TiTe 2 under tensile biaxial strain calculated with the PBE functional: (a) 0%, (b) 3% and (c) 6% biaxial strain. Γ = (0, 0, 0), M = (1/2, 0, 0) and K = (1/3, 1/3, 0) in units of the reciprocal hexagonal lattice vectors. and, consequently, of the semimetallic overlap. Only HSE06 level, unstrained single-layer TiTe 2 is indeed when this overlap decreases with respect to the PBE unstable towards the 2×2 CDW distortion. By description of the system the 2 ×2 CDW becomes more relaxing the structure around the minimum of the stable. This observation is consistent with the fact that frozen-phonon curve we have obtained an energy gain when Te...Te interlayer interactions come into play in of 3.8 meV/formula unit. This stabilization energy the bulk or even in the double-layer, the semimetallic is much lower than the value obtained for a TiSe 2 overlap increases and the 2 ×2 CDW is not observed monolayer using a PBE functional, 6 meV/f.u [13]. anymore. Hellgren et al. [20] showed that for bulk TiSe 2 hybrid functionals predict a much higher stabilization energy than the PBE functional. Thus, it is understandable that according to our PBE type studies single-layer TiTe2 does not tend to experience a 2 ×2 CDW instability. Overall, our data suggest a weak driving force for the distortion. Figure 3. Energy dierence (in meV per formula unit) between the undistorted and 2 ×2 CDW structures of an unstrained TiTe 2 single-layer according to DFT calculations using the hybrid HSE06 functional These results are reminiscent of the above mentioned work concerning the functional type dependence of the semimetallic overlap of bulk TiSe2 [20] and prompted us to reconsider the stability Figure 4. Calculated band structure for TiTe 2 single-layer using the PBE (a) and hybrid HSE06 (b) functionals. The size of the unistorted vs. 2×2 CDW structures using the of the green and red circles is proportional to the Te and Ti hybrid type functional HSE06 [22, 23]. Shown in character, respectively. Fig. 3 is a frozen-phonon calculation of the energy dierence between the undistorted structure and the Since this is a signicant result let us consider 2×2 CDW structure following the soft phonon mode in more detail the region of the semimetallic overlap. distortion. The curve clearly shows that, at the A fatband analysis of the PBE and HSE06 band 2×2 Charge Density Wave in single-layer TiTe 2 5 structures is shown in Figs. 4a and b, respectively. the nature of the band structure which entirely applies A close inspection at the region around Γ points out to the HSE06 one reported in Fig. 4b. some clear dierences between the results of both functionals. For PBE, above the (partially empty) top of the valence band at Γ there is a non-degenerate band and a pair of degenerate bands at a higher energy. The order is the opposite for the HSE06 functional (Fig. 4b), for which the two degenerate bands are lower in energy. Around the Fermi level a heavy mixing of the t2g levels of the Ti atom and the 5 p levels of the Te atoms occurs. Assuming a local coordinate system in which the three-fold symmetry axis of the octahedron occurs along the z direction, the non- degenerate Ti-based band is essentially built from the Ti dz2 orbital, whereas the top of the valence band, which is doubly degenerate at Γ, is mainly built from the Ti dx2 −y2 and dxy orbitals which are somewhat tilted because the plane of the Ti atoms is not a symmetry plane. The single band at Γ is slightly higher in energy than the doubly-degenerate pair because of the slight rhombohedral distortion. The upper, doubly degenerate, Te-based bands are built from the Te px and py orbitals. Having identied the nature of these levels it is immediately clear that the Te-based doubly degenerate set is higher in energy for the PBE-type calculation. This means that, in the absence of Ti-Te hybridization, the Te px /py valence bands would raise up to higher energies at Γ than in the HSE06-type calculations. Consequently, when the hybridization is switched on, a larger number of electron and hole carriers is induced. Indeed, a careful look at the band structures in Fig. 4 shows that there is a signicantly larger overlap between the top of the valence band Figure 5. HSE06 band structure (a) and Density of States around Γ and the bottom of the conduction band (DOS) (b) for single-layer TiTe 2 . Broken red lines correspond to the 2×2 CDW structure and black continuous lines to the at M for the PBE functional than for HSE06. We undistorted structure. The bands are represented in the BZ of note that the PBE-type band structure for single-layer the 2×2 structure and therefore are folded with respect to that TiSe2 , which provides a satisfactory description of the in Fig. 1b. The origin of the energy scale is the Fermi level of relative stability of the undistorted and 2 ×2 CDW in the undistorted structure. this system [13, 19], shows exactly the same topology and particularly, the same band ordering as the HSE06 The HSE06 band structure and density of states band structure of Fig. 4b. for both the undistorted and the 2 ×2 CDW structures As soon as one moves along the Γ-M line (i.e. are reported in Figs. 5. As shown in both gures the along the a* direction), the only symmetry element stabilizing 2×2 CDW distortion opens a gap at the preserved is the symmetry plane perpendicular to the Fermi level. The analysis of the nature and origin layer and going along the a* direction. One of the of the distortion goes along the same lines presented two Ti-based doubly degenerate levels at Γ and the Ti in detail in our recent work for the TiSe 2 single-layer dz2 mix and interact with one of the Te p bands near Γ and will not be repeated here [13]. However, we note leading to the slowly descending band from Γ-M which that the non-activated conductivity is kept in single- is associated with the electron pockets near M . With layer TiTe2 below the CDW transition temperature [8] the local system of axes mentioned above, the crystal so that no gap should occur at the Fermi level. orbitals around M are almost exclusively made of tilted Consequently, the HSE06 calculations exaggerate the Ti dx2 −y2 orbitals (i.e. a mixing of Ti dx2 −y2 and Ti dz2 tendency towards the transition. As mentioned above, which leads to the tilting of the orbital) exactly as we the PBE calculations suggest that under reasonable recently reported for the TiSe 2 single-layers [13]. We tensile strain a 2×2 CDW is favored without the refer the reader to this work for a detailed analysis of development of a band gap. Only when the strain is relatively large (i.e. around 5-6%, see Fig. 2) a 2×2 Charge Density Wave in single-layer TiTe 2 6 band gap really opens. The HSE06 band structure of the distorted phase and, as a consequence, a of Fig. 5a is similar to that in Fig. 2c corresponding noticeable band gap of more than 0.1 eV is induced to a 6% tensile strength. We believe that the HSE06 at the Fermi level. This is in contrast with the functional exaggerates the stability of the 2 ×2 CDW metallic character of the TiTe 2 single-layers below and that the real situation concerning the weaker the transition temperature. Interestingly, PBE type stabilization of the CDW phase and the lack of a full calculations for the case where the single-layer is bandgap opening would rather correspond to that of subject to a slight tensile strain also favor the the PBE functional with a moderate tensile strain. 2×2 distortion while keeping the semimetallic overlap. This observation, together with previous results on The study suggests that the magnitude of the semi- bulk TiSe2 [20] suggest that a predictive description metallic overlap is a key factor controlling the of these TiX2 systems is attainable by using the tendency towards the distortion and consequently only HSE06 functional and tuning the actual contribution functionals describing such overlap very accurately can of the exact exchange. However, from the viewpoint provide a truly predictive description of the electronic of the physical origin of the CDW we do not nd structure. any noticeable dierence with the previously reported According to the present study the mechanism of analysis of the 2×2 CDW instability of TiSe2 single- the CDW instability in single-layer TiTe 2 seems to layers. be the same phonon mediated mechanism acting for To further assess our conclusion we carried out single-layer TiSe2 [13] although now the driving force HSE06 calculations for a TiTe 2 double-layer. Two is smaller and the semimetallic character is kept be- single-layers with the optimized structure were placed low the transition temperature. As mentioned above, with the interlayer bulk distance. The stabilization the magnitude of the semimetallic overlap seems to be energy of the 2×2 CDW is reduced to practically one- one of the key factors in controlling the likeliness of the half the value in the single-layer. Taking into account 2×2 CDW. Taking into account that the overlap should the small (and exaggerated) value for the single-layer, increase in these TiX2 systems when the number of the driving force for the 2×2 CDW distortion in the short Te...Te contacts increase or when they become TiTe2 double-layer must be extremely small or most stronger, the overlap should increase from single-layers likely nil, as experimentally found. to bulk and from X=S to X=Te. Since the instabil- Our study points out an interesting possibility. ity is not observed in bulk TiS 2 nor in double-layer Since tensile strain has been found to be a useful TiTe2 , it seems that only a relatively narrow range of technique to induce modications in single-layer or semimetallic overlaps is associated with the instabil- few-ake materials [25, 26, 27, 28], it is possible that ity. In this respect, a signicant result of the study is the 2×2 CDW can be induced in double-layers or that tensile strain stabilizes the 2 ×2 CDW distortion triple-layers of TiTe2 by using a small tensile strain. ‡ in single-layer TiTe2 . This could be used to induce the Another useful hint provided by our study is that the instability in double- or triple-layers of TiTe 2 which in 2×2 CDW in single-layer TiTe 2 may be more stable the absence of strain remain undistorted, to induce a and result with a band gap under tensile strain or may stronger distortion leading to the creation of a band be suppressed under a slight compresive strain. gap in single-layer TiTe2 or most likely to suppress the 2×2 CDW under a small compression. Such studies 3. Conclusions could provide useful insight on the CDW mechanism of group IV 1T -TiX2 phases. A density functional theory study concerning the origin of the 2×2 CDW distortion recently reported This work was supported by Spanish MINECO experimentally for single-layer TiTe 2 has been carried (Grants FIS2015-64886-C5-3-P and FIS2015-64886- out. This report is surprising because neither double- C5-4-P, and the Severo Ochoa Centers of Excellence layer nor bulk TiTe2 exhibit the 2×2 distortion Program under Grants SEV-2013-0295 and SEV-2015- and a PBE-based DFT study predicts that the 0496), Generalitat de Catalunya (Grant 2017SGR1506 undistorted structure is also more stable for the single- and the CERCA Programme) and by the European layer. Our study shows that, whereas calculations Union H2020-EINFRA-5-2015 MaX Center of Excel- employing the semi-local functional PBE favor the lence (Grant No. 676598). undistorted structure, the hybrid functional HSE06 correctly predict the 2×2 distortion. However, the Appendix A. Computational details HSE06 calculations seem to exaggerate the stability ‡ Note that tensile strain has been shown to decrease the The geometrical optimizations, electronic and phononic semimetallic overlap and eventually lead to a band gap without band structures were carried out using a numerical changing the size of the unit cell both in single-layer [29] and atomic orbitals density functional theory (DFT) [31, bulk 1T -TiS2 ) [30]. 2×2 Charge Density Wave in single-layer TiTe 2 7 32] approach implemented in the Siesta code [33, [10] Li L J, O'Farrell E C T O, Loh K P, Eda G, Özyilmaz B 34]. The Perdew-Burke-Ernzerhof (PBE) functional and Castro-Neto A H 2016 Nature 529 185 [11] Chen P, Chan Y H, Fang X Y, Zhang Y, Chou M, Mo was used to account for the exchange-correlation S K, Hussain Z, Fedorov A V and T-C C 2015 Nature energy [21]. The core electrons have been re- Communications 6 8943 placed by norm-conserving scalar relativistic pseu- [12] Silva-Guillén J A, Ordejón P, Guinea F and Canadell E dopotentials [35] factorized in the Kleinman-Bylander 2016 2D Materials 3 035028 [13] Guster B, Canadell E, Pruneda M and Ordejón P 2018 2D form [36]. We include the 3p shell of Ti explicity Materials 5 025024 in the valence, as semicore states. We have used a [14] de Boer D K G, van Bruggeb C F, Bus G W, Coehoorn R, split-valence double-ζ basis set including polarization Haas C and Sawatzky G A 1984 Phys. Rev. B 29 6797 functions [37]. The non-linear core-valence exchange- [15] Koike Y, Okamura O, Nakanomyo T and Fukase T 1983 J. Phys. Soc. Jpn. 52 597 correlation scheme [38] was used for all elements. In [16] Allen P B and Chetty N 1994 Phys. Rev. B 50 14855 the direction normal to the single-layer we chose a vac- [17] Di Salvo F J, Moncton D E and Waszczak J V 1976 Phys. uum space of 50 Å in order to avoid possible interac- Rev. B 14 43214328 tions between the layer and its images. In the case of [18] Chen P, Chan Y H, Wong M H, Fang X Y, Chou M Y, Mo S K, Hussain Z, Fedorov A V and Chiang T C 2016 Nano geometrical optimization calculations, the atomic co- Lett. 16 63316336 ordinates were relaxed until the forces on them were [19] Singh B, Hsu C H, Tsai W F, Pereira V M and Lin H 2017 below 10−5 eV/Å. In all calculations, we use a cuto Phys. Rev. B 95 245136 [20] Hellgren M, Baima J, Bianco R, Calandra M Mauri F and of 2500 Ry for the real space integrals, and a toler- Wirtz L 2017 Phys. Rev. Lett. 119 176401 ance of 10−7 and 10−6 eV on the density matrix and [21] Perdew J P, Burke K and Ernzerhof M 1996 Physical Review the total energy, respectively, for the convergence of Letters 77 38653868 the SCF cycle. To sample the Brillouin cell for the [22] Heyd J, Scuseria G E and Ernzerhof M 2003 J. Chem. Phys. 118 82078215 electronic states, a Monkhorst-Pack [39] k-point grid [23] Krukau A V, Vydrov O A, Izmaylov A F and Scuseria G E of 300×300×1 was used for the undistorted minimum 2006 The Journal of Chemical Physics 224106 cell and it was scaled accordingly where supercell cal- [24] Arnaud Y and Chevreton M 1981 J. Sol. State Chem. 39 culations were performed. 2230 [25] Choi S M, Jhi S H and Son Y W 2010 Nano Lett. 10 3486 All HSE06 hybrid functional [22, 23] calculations 3489 were performed using the VASP code [40]. Core [26] Johary P and Shenoy U B 2012 ACS Nano 6 54495456 electrons were treated by means of the projector [27] Biele R, Flores E, Ares J R, Sánchez C, Ferrer I J, Rubio- augmented wave method [41, 42] including semicore Bollinger G, Castellanos-Gomez A and D'Agosta R 2018 Nano Research 11 225232 states for Ti. We used a planewave basis set with an [28] Roldán R, Chirolli L, Prada E, Silva-Guillén J A, San-Jose energy cuto of 330 eV. A Monkhorst-Pack [39] k-point P and Guinea F 2017 Chem. Soc. Rev. 46 43874399 grid of 24×24×1 was used to sample the Brillouin zone. [29] Xu C, Brown P A and Shuford K L 2015 RSC Adv. 5 83876 [30] Samanta A, Pandey T and Dingh A K 2014 Phys. Rev. B. 90 174301 References [31] Hohenberg P and Kohn W 1964 Physical Review 136 B864 B871 [1] Rossnagel K 2011 J. Phys.:Condens. Matter 23 213001 [32] Kohn W and Sham L J 1965 Physical Review 140 A1133 [2] Wilson J, Di Salvo F and Mahajan S 1975 Adv. Phys. 24 A1138 117201 [33] Soler J M, Artacho E, Gale J D, García A, Junquera [3] Castro Neto A H 2001 Physical Review Letters 86 4382 J, Ordejón P and Sánchez-Portal D 2002 Journal of [4] Noat Y, Silva-Guillén J A, Cren T, Cherkez V, Brun C, Physics: Condensed Matter 14 27452779 Pons S, Debontridder F, Roditchev D, Sacks W, Cario [34] Artacho E, Anglada E, Diéguez O, Gale J D, García A, L, Ordejón P, García A and Canadell E 2015 Phys. Rev. Junquera J, Martin R M, Ordejón P, Pruneda J M, B 92(13) 134510 Sánchez-Portal D and Soler J M 2008 Journal of Physics: [5] Bhimanapati G R, Lin Z, Meunier V, Jung Y, Cha J, Das S, Condensed Matter 20 064208 Xiao D, Son Y, Strano M S, Cooper V R, Liang L, Louie [35] Troullier N and Martins J L 1991 Physical Review B 43 S G, Ringe E, Zhou W, Kim S S, Naik R R, Sumpter 19932006 B G, Terrones H, Xia F, Wang Y, Zhu J, Akinwande [36] Kleinman L and Bylander D M 1982 Physical Review D, Alem N, Schuller J A, Schaak R E, Terrones M and Letters 48 14251428 Robinson J A 2015 ACS Nano 9 1150911539 [37] Artacho E, Sánchez-Portal D, Ordejón P, García A and [6] Ugeda M M, Bradley A J, Zhang Y, Onishi S, Chen Y, Soler J M 1999 Physica Status Solidi (b) 215 809817 Ruan W, Ojeda-Aristizabal C, Ryu H, Edmonds M T, [38] Louie S G, Froyen S and Cohen M L 1982 Physical Review Tsai H Z, Riss A, Mo S K, Lee D, Zettl A, Hussain Z, B 26 17381742 Shen Z X and Crommie M F 2016 Nature Physics 12 92 [39] Monkhorst H J and Pack J D 1976 Physical Review B 13 [7] Sugawara K, Nakata Y, Shimizu R, Han P, Hitosugi T, Sato 51885192 T and Takahashi T 2016 ACS Nano 10 13411345 [40] Kresse G and Furthmüller J 1996 Physical Review B 54 [8] Chen P, Pai W W, Chan Y H, Takayama A, Xu C Z, Karn 1116911186 A, Hasegawa S, Chou M Y, Mo S K, Fedorov A V and [41] Blöchl P E 1994 Physical Review B 50 17953 Chiang T C 2017 Nature Communications 8 516 [42] Kresse G and Joubert D 1999 Physical Review B 59 1758 [9] Xi X, Zhao L, Wang Z, Berger H, Forró L, Shan J and Mak K F 2015 Nature Nanotechnology 10 765769