This is the submitted version of the article: Sánchez-Ochoa F., Hidalgo F., Pruneda M., Noguez C.. Unfolding method for periodic twisted systems with commensurate Moiré patterns. Journal of Physics Condensed Matter, (2020). 32. 025501: - . 10.1088/1361-648X/ab44f0. Available at: https://dx.doi.org/10.1088/1361-648X/ab44f0 Unfolding method for periodic twisted systems with commensurate Moiré patterns F. Sánchez-Ochoa1 , Francisco Hidalgo2,3 , Miguel Pruneda3 and Cecilia Noguez1 1 Instituto de Fı́sica, Universidad Nacional Autónoma de México, Apartado Postal 20-364, Cd. de México C.P. 01000, Mexico 2 Departamento de Ciencias Básicas, Univerisdad Autónoma Metropolitana-Azcapotzalco, Av. San Pablo 180, Cd. México C.P. 02200, Mexico 3 Catalan Institute of Nanoscience and Nanotechnology (ICN2), CSIC and BIST, Campus UAB, Bellaterra, 08193 Barcelona, Spain E-mail: fsanchez@fisica.unam.mx, cecilia@fisica.unam.mx Abstract. We present a general unfolding method for the electronic bands of systems with double-periodicity. Within density functional theory with atomic orbitals as basis-set, our method takes into account two symmetry operations of the primitive cell: a standard expansion and a single rotation, letting to elucidate the physical effects associated to the mutual interactions between systems with more than one periodicity. As a result, our unfolding method allows studying the electronic properties of vertically stacked homo- or √ heterostructures. √ We apply our method √ (3 ×√ to study √ ( 13 × 13) twisted single-layer graphene, 3) single-layer graphene, √ and (2 3 × 2 3) graphene- ( 7 × 7) tungsten disulfide heterostructure with an interlayer angle of 10.9◦ . Our unfolding method allows observing typical mini gaps reported in heterostructures, as well as other electronic deviations from pristine structures, impossible to distinguish without an unfolding method. We anticipate that this unfolding method can be useful to compare with experiments to elucidate the electronic properties of two-dimensional homo- or heterostructures, where the interlayer angle can be considered as an additional parameter. Keywords: graphene, tungsten disulfide, heterostructures, twisted bilayers, DFT, unfolding Unfolding method for periodic twisted systems 2 1. Introduction folding effects hide the real physical interactions, usually leading to a misinterpretation of the results Vertically stacked two-dimensional (2D) systems with in the Brillouin zone (BZ) of the SC. Indeed, the van der Waals (vdW) interactions have fascinating comparison between the theoretical electronic band physical properties and potential applications in structure of SC and angle-resolved photoemission optoelectronics, electronic transport [1], among other spectroscopy (ARPES) measurements is quite difficult fields. In principle, these 2D-vdW structures should [15]. A successful technique to unveil the change of preserve the physical advantages of their pristine electronic structure, hidden by the use of SC method, monolayers but also allow for achieving unusual and is the band unfolding approach (hereafter, unfolding even superior properties, which cannot be obtained for simplicity), which have been developed as a post- otherwise [2]. Besides the atomic composition of each processing tool for DFT [16, 17, 18, 19, 20, 21] monolayer, the vertical stacking adds a new parameter: and tight-binding calculations as well [22, 23, 24]. the interlayer rotation angle (θ). In general, when two In the unfolding method for DFT, the calculated vertically stacked monolayers of the same composition Kohn-Sham orbitals of SC are unfolded into the (a homostructure) have a θ 6= 0, a Moiré pattern original BZ of the PC. This technique projects the is formed, being its commensurable area inversely electronic states from the BZ of the SC to the BZ of proportional to θ [3]. In contrast, two vertically the reference system by assigning a spectral weight, stacked monolayers of different atomic composition Wk (I), where k is a point on the reciprocal space (heterostructure) also exhibit a Moiré pattern just by of the SC and I denotes I-th electronic band to be the condition of a lattice mismatch of layers [4]. In this unfolded. Consequently, Wk (I) quantifies how much sense, different theoretical and experimental groups of the reference system remains in the perturbed one, have reported new physical properties exhibited by for example, when the translation symmetry is lost. two vertically stacked graphene (G) layers related to θ, However, unfolding methods for DFT calculations have including optical activity[5], topological edge states[6], been applied only for conventional hexagonal SCs and a Mott insulator-to-superconductor transition at a [20, 21, 25]. In contrast, unfolding methods for tight- critical θ [7, 8]. Besides bilayer G, the vertical stacking binding calculations have been improved for hexagonal of transition metal dichalcogenides (TMDCs) with and non-conventional SCs, and applied to twisted G have been considered promising materials for the systems [16, 26]. design of new electronic, optoelectronic, and photonic Therefore, we present an unfolding method to devices, in particular based on tungsten disulfide analyze the electronic structure of 2D-vdW twisted (WS2 ) or molybdenum disulfide (MoS2 ) [5, 9, 10, 11, systems under the scheme of DFT with atomic 12, 13], both layered semiconductors with a direct band orbitals as basis-set. Our method is an extension gap for single layer. of previous work of [27] and [28], allowing the Through first-principles methods, for instance, unfolding of electronic bands of structures constructed based on Density Functional Theory (DFT), it is by conventional and non-conventional SCs. To show possible to study the electronic properties of these the potentiality of the method, we investigate the novel 2D materials, where the use of periodic electronic band structure of√two ideal √ G single-layer boundary conditions is highly desirable to reduce the with periodicities 3 × 3 and 13 × 13. Furthermore, computational effort. However, periodic boundary the electronic band structure of a twisted G-WS2 conditions require matching of the lattice constants of bilayer with θ = 10.9◦ is unfolded using our method, vertically stacked layers to construct a commensurate unveiling the physical consequences of the mutual supercell (SC) of the bilayer heterostructure, such interactions between G and WS2 layers. The content that each layer can be built by a conventional or of the paper is organized as follows. In section 2, we non-conventional SC. By definition, a conventional present a brief description of the unfolding method. SC is constructed just by expansion symmetry Section 3 shows the models studied here to validate operations on the primitive cell (PC), while non- our method. Section 4 summaries the computational conventional SC requires expansion and rotation details employed in the DFT calculations. Then, in operations. In consequence, SC method is a Section 5, we use our method to unfold the electronic common practice in periodic DFT calculations, making band structure of the SCs described above. Finally, in available the study of systems with thousands of section 6 we discuss our main conclusions. atoms, different composition, impurities, vacancies, or combined systems, including vertically stacked 2D- 2. Method vdW heterostructures with θ 6= 0. Nevertheless, the use of SCs leads to an artificial Under the scheme of DFT with a linear combination of folding of the electronic bands as a consequence of atomic orbitals (LCAO) as basis-set and using Bloch’s the reduction of the reciprocal space [14]. Such Unfolding method for periodic twisted systems 3 theorem, the total wavefunction of a system is: carried out in 2D or 3D, respectively. Therefore, the X equivalent u atoms of the SC from PC is easily obtained ψI (k, r) = (exp ik · Rj0 )φj 0 (r)cj 0 i (k); (1) through Eq. (3) and Eq. (4). j0 For 2D hexagonal lattice, any conventional or non- where the sum in j 0 extends to all basis orbitals in conventional SC can be set up using the transformation space, I labels the different bands, and cj 0 i corresponds matrix P, considering c vector as an invariant [29]. to the expansion coefficients of the Kohn-Sham orbitals Thus, Pmn elements in Eq. (4) are P33 = 1 and P13 = obtained by the self-consistent calculation. The base of P23 = P31 = P32 = 0. Finally, the other elements are our unfolding method stands on the calculation of the functions of two integer parameters, m and n to take spectral weight, Wk (I) which quantifies the projection into account expansion and rotation operations. Then, of the electronic state from the BZ of the SC to the the matrix transformation, P, simplifies to: BZ of the PC for the I-th band. Hence, following the   development of Deretzis [27], Wk (I) is obtained as: m n 0     P =  −n m + n 0  . (5) N N 0 0 1 1 X X X Wk (I) =  c̃u+j (k) ×  cu+j (k) , (2) N j j Notice that if m or n = 0, then the 2D hexagonal SC uP C is a conventional one, while a non-conventional SC is where the sum on c coefficients run over all the obtained for m, n 6= 0. For our propose, given a 2D SC, equivalents orbitals j of the N times repeated PC, it is necessary to identify which of the atomic positions and over all atoms u with equivalent positions in SC are equivalent to the PC. This is solved by a process with respect to PC. Here, c is a complex number and c̃ which follows the inverse process described above to set represents its complex conjugate, and N is an integer up a SC from a PC. Additionally, because commonly defined as N = ΩSC /ΩP C , where Ω denotes the volume DFT calculations are carried out including a structural in the real space. relaxation, it is necessary to specify a threshold value Notice that Wk (I) determines the equivalence of to determine whether an atom has deviated from its the electronic states of SC on the corresponding of initial position or not. PC and must obey that 0 ≤ Wk (I) ≤ 1, which is Another important issue for the calculation of true if and only if the basis-set is orthogonal and Wk (I) is the scanning of the BZ along specific paths of normalized [27, 28]. For SC with more than one the SC and their projections on equivalent paths of the periodicity or built up by more than one PC, for BZ of PC [27, 28]. In general, the fractional coordinates instance, any vertically stacked heterostructure, it is (kx , ky , kz ) of any point in the reciprocal space of PC necessary to define on which PC the electronic states transform in similar way as the lattice vectors, or any will be projected. A wrong choice on the PC could point in the real space, do [29]. Therefore, points in induce a bad interpretation on the unfolded bands or the reciprocal space of SC (kx0 , ky0 , kz0 ) are related with even ghost states [26]. Once the PC has been defined, equivalent points in the reciprocal space of PC by an the problem of identify the equivalent u atoms in SC equivalent transformation given by Eq. (3) and the is directly solved because the layer was constructed by transformation matrix defined in Eq. (4): expansion and rotation symmetry operations on the (kx , ky , kz ) = kx0 , ky0 , kz0 P.  PC, which involves an affine transformation between (6) SC and PC. This transformation consists of two parts: For a 2D hexagonal lattice, kz = kz0 = 0 and P is a shift of the origin and a linear part. The former is defined through Eq. (5). given just by a straight translation of the lattice vectors Finally, we can calculate the electronic states of of PC (a, b, c), while the latter implies a change of commensurate 2D Moiré SCs with the new coordinates their orientation, length or both. Through just by a (kx0 , ky0 , 0) based on periodic DFT, and thereafter, we linear transformation, the lattice vectors of the SC (a0 , can apply the unfolding described above to calculate b0 , c0 ) are defined by: Wk (I), obtaining the SC projected band structures a0 , b0 , c0 = (a, b, c) P,  (3) onto the original PC systems. Since we are performing DFT calculations with a non-orthogonal basis-set with P, a transformation matrix, defined as: [30], we employ a Gram-Schmidt process [31] to   P11 P12 P13 orthogonalize the Kohn-Sham orbitals and carry out P =  P21 P22 P23  , (4) the unfolding process through the calculation of Wk (I) P31 P32 P33 via Eq. (2). We want to stress that our unfolding where Pmn are positive integer values and the absolute code reads the whole information of wavefunctions value of the determinant of P corresponds to N . from the corresponding siesta file and, then, performs Notice that Pmn can be parametrized in terms of two the Wk (I) calculation. Beside with the linear scaling or three integer values whether the transformation is [30], siesta code has different advantages: the choice Unfolding method for periodic twisted systems 4 of the basis-set size, different subroutines for writting problem, it is necessary to find SC of G and WS2 the wavefunctions file, and the setting of how many k with similar lattice parameters; then, one of the layers points will be calculated. Putting all this together, our is fixed and the other one is enlarged or reduced to unfolding code is able to work efficiently for twisted assure commensurability. In our case, we fixed WS2 systems with significant number of atoms [27] with as the substrate and G is slightly enlarged or shrunk respect to other unfolding DFT methods based on to match both lattice parameters of SCs. Thus, the plane-waves [20, 21], being able to unfold the electronic in-plane strain (∆) imposed in G and associated to states even from the 1st occupied band up to the I-th commensurability condition is an intrinsic property, desired band and improving the resolution of projected defined as [35]: band structure easily by the number of k points along √ p aG m2 + n2 + mn − aW S2 p2 + q 2 + pq the paths. ∆= p , (9) aW S2 p2 + q 2 + pq 3. Models where m and n indexes correspond to G, while p and q do to WS2 , obeying the following constrains: m ≥ n In order to show the validity and main characteristics and p ≥ q and |∆| ≤ 1.5% [35]. Then, we employ of the present unfolding method, we study three the matrix transformation described in Eq. (5) with structures: a conventional SC of monolayer G, a non- m = 2, n = 2 for G, and p = 2, q = 1 for WS2 , plus the conventional SC of monolayer G, and a twisted bilayer lattice vectors of G and WS2 to generate the periodic G-WS2 heterostructure. First, let us describe the SC of G-WS2 , shown in Figure 2(c). For these indexes, cases of monolayer G. Lattice vectors√ of the PC are ∆ = −1.09%, where the minus sign means compression a = aG (1, 0, 0), b = aG (−1/2, 3/2, 0) and c = of the G layer, the twisted-angle θ = 10.9◦ , the lattice c(0, 0, 1), where aG = 2.44 Å is the lattice parameter parameter of the SC is 8.34 Å, and the total atoms optimized through DFT calculations, which is in very of the SC is 45. Visualization of atomic models was good agreement with other reports [5, 32], and c performed using the vesta program [36]. is the vacuum space imposed for the SC (see the computational details in Section 4). Conventional and 4. Computational details non-conventional √ √ SC of monolayer G have a 3 × 3 and 13 × 13 periodicity, respectively. They are DFT calculations are performed using the siesta constructed using Eq. (3) and Eq. (5) with m = 3 code [30]. Exchange-correlation energy is described and n = 0 for × 3) with N = 9, and m = 3 and √ G(3 √ by local density approximation (LDA) with the n = 1 for G( 13 × 13) with N = 13. Inset of Figure Perdew-Zunger parametrization [37]. Norm-conserving 1(a) and (d) show both SC. Hence, the total number pseudopotentials [38] in their scalar relativistic and of atoms, NT , defined as 2N , are 18 and 26 atoms, non-local form [39] are employed for electron-ion respectively. In general, the rotation angle, θ, between interactions. The expansion of the valence electronic two hexagonal layers is defined as: states is described by a LCAO with an optimized mp + nq + (mq + np) /2 double-ζ polarized (DZP) basis-set plus diffuse orbitals cos θ = √ p , (7) [40, 41, 42]. If the basis-set uses strictly localized m2 + n2 + mn × p2 + q 2 + pq orbitals, it is necessary to include diffuse orbitals for where m and n indexes correspond with one of the a better description of interface and vacuum regions rotated layer, whereas p and q correspond with the in slabs and surfaces [40]. The electronic density other one. Taking one of the layers as a conventional is sampled by a mesh-cutoff energy of 300 Ry and SC, p 6= 0 and q = 0 and Eq. (7) simplifies: converged up to 10−5 . The Monkhorst-Pack scheme mp + (np) /2 [43] is used for sampling the BZ with an optimal cos θ = √ , (8) 25 × 25 × 1 k−grid. Energy mesh-cutoff and k−grid p × m2 + n2 + mn convergence tests were performed leading to the above giving θ = 13.90◦ . mentioned values as the optimal ones. A Fermi-Dirac Now, let us describe the case of twisted G-WS2 smearing of 1 K is adopted to obtain convergence. bilayer. PCs of G and WS2 have hexagonal symmetry; The structural optimization is carried out within the hence, both PCs exhibit the same lattice vectors but conjugate gradient algorithm with a maximum value of different lattice parameter [see Figure 2(a) and (b)]. 0.01 eV/Å for the interatomic forces, as well as 0.5 GPa Letting an optimization of the lattice parameter of in the stress components. A vacuum space of 55 Å is WS2 , we obtain aW S2 = 3.16 Å, which is in very applied to avoid spurious interactions between periodic good agreement with other reports [33, 34]. Since the images along c direction. Finally, all calculations PCs of G and WS2 have different lattice parameters, include a dipole correction, as it was suggested by the construction of a periodic vertically stacked SC of Bengtsson [44], to cancel the artificial electric field G-WS2 is an additional problem. To overcome this Unfolding method for periodic twisted systems 5 3 3 K (a ) 2 (b ) 2 (c ) 1 11 111 1 111 1 1 1 111 11 111 0 0 1 111 M 11 111 1 111 K -1 -1 11 e n e rg y (e V ) e n e rg y (e V ) 111 1 M -2 -2 111 1 111 11 Γ 111 -3 -3 1 010 00 000 0 -4 -4 000 00 000 0 000 -5 -5 0 000 00 000 -6 -6 0 000 00 000 0 -7 -7 000 0 000 00 00 -8 -8 0 -9 -9 Γ K M Γ Γ K M Γ 3 3 (e ) (f) 1 (d ) K 2 2 1 11 111 1 111 1 1 1 111 11 111 M 0 0 1 111 11 111 -1 -1 1 111 e n e rg y (e V ) e n e rg y (e V ) 1 K 111 11 M -2 -2 111 1 111 11 Γ -3 -3 111 1 000 0 000 00 -4 -4 000 0 000 00 -5 -5 000 0 000 0 000 -6 -6 00 000 0 000 -7 -7 00 000 0 000 0 00 -8 -8 0 -9 -9 Γ K M Γ Γ K M Γ Figure 1. Top row: (a) reciprocal space of G(3 × 3) SC. The green (shaded-red) hexagons denotes the first Brillouin zone of G SC (PC). The irreducible Brillouin zone of SC (PC) is enclosed by the filled-green (red) triangle. High symmetry points are coloured in green (red) for SC (PC). The inset shows the stick representation of G(3 × 3) SC layer in real space and the dashed-black lines delimit the size of SC. (b) calculated band structure for G(3 × 3) SC. The plot (c) shows the SC electronic states (gray background) calculated using the red paths in (a) while, after the projection, those SC electronic states with the same symmetry of that reference PC are√plotted√with circles of different size and color, which are proportional to Wk (I). The bottom row depicts the characteristics for G( 13 × 13) SC as the top row. The number of k-points scanned along Γ →K, K→M and M→ Γ were 30, 15 and 30, respectively. The Fermi level is set to 0 eV. Figure 2. Top view of stick representation of (a) G PC, (b) WS2 PC and (c) G/WS2 SC. The PCs in (a) and (b) have the same orientation and indicated by dashed lines. In (c), the G/WS2 SC is highlighted by the green area, while G and WS2 PCs are in red and blue areas, respectively. θ is the angle between the armchair edges of both layers. Color code: gray, yellow, and blue sticks represent C, S and W elements, respectively. Unfolding method for periodic twisted systems 6 induced by the vertical asymmetry of the SC and the SC at the Fermi level around the K point is also periodic boundary conditions. present. Comparing Figure 1(e) and 1(b), there are more electronic bands (gray lines as background) in the 5. Results former at a glance. Nevertheless, when the unfolding process is applied using the new fractional coordinates As a first proof of our unfolding code, we study the K= (2/3, 5/3, 0) and M= (3/2, 1/2, 0) of PC, the fully electronic structure of an ideal G(3 × 3) conventional band structure of PC is again fully recovered, as in layer. The reciprocal space, the first BZ and the Figure 1(f) with Wk (I) = 1. An analysis of unfolded irreducible BZ of the SC and PC are shown in Figure bands of Figure 1(f) reveals the absence of states with 1(a). Also, the inset presents the G(3×3) conventional Wk (I) 6= 1, since the coordinates (kx , ky ) of scanned k- layer. Figure 1(b) shows the electronic band structure points of SC are not the same positions as Bragg points of the SC. The typical linear dispersion of electronic or lines in reciprocal space. Therefore, our unfolding states of the PC at the Fermi level around the K code works for conventional hexagonal SCs, and non- point is also present in SC but around the Γ point. conventional ones as well. Due to the folding process, many electronic states are Now we turn to a more complex system: a merged into the small SC reciprocal space [the green heterostructure built up by of two twisted non- area in Figure 1(a)]. The irreducible BZ of the PC conventional SCs. In Figure 2(c), the √2D commen- √ has the fractional coordinates K= (1/3, 1/3, 0) and surable √ heterostructure √ composed of 2 3 × 2 3 G M= (1/2, 0, 0) in terms of the corresponding reciprocal and 7 × 7 WS2 with an interlayer rotation angle lattice vectors of PC. After the transformation with θ = 10.9◦ is depicted. Hence, the number of times Eqs. (6) and (4), the new fractional coordinates that the PC is repeated on the SC, N in Eq. (2), or of high-symmetry points are K= (1, 1, 0) and M= equivalently the determinant of P is 12 for G, and 7 (3/2, 0, 0) but in terms of the reciprocal lattice vectors for WS2 . The optimized average interlayer distance, d, of SC. With this transformation, the highly symmetric defined as the average distance between carbon atoms points in the reciprocal space of SC can be easily of G and upper sulfur atoms of WS2 is d = 3.28 Å, in identified in the reciprocal space of PC [see Figure good agreement with previous reports [45, 46]. Aver- 1(a)]. When we unfold the electronic states of SC age distances of carbon–carbon (C–C), tungsten–sulfur calculated throughout the red path of G [gray band (W–S) and up and down sulfur atoms (Su –Sd ) are 1.41, structure in Figure 1(c)] onto the PC of G, the fully 2.41 and 3.16 Å, respectively. The irreducible BZ of band structure of the ideal PC is recovered, as it PC of both WS2 and G has the same fractional coordi- is shown by red circles in Figure 1(c), which means nates. Once the transformation is done using Eqs. (6) Wk (I) = 1. Some unfolded wavefunctions in Figure and (5) on the high-symmetry points, the new frac- 1(c) have Wk (I) 6= 1 along M–Γ within -7.5 and - tional positions are K= (0, 2, 0) and M= (1, 1, 0) for G, 4.5 eV range, given by green dots and the yellow dot and K= (1/3, 4/3, 0) and M= (1, 1/2, 0) for WS2 . Mul- along K–M at -7.5 eV, because the wavefunctions are tiplying these new fractional positions by the recipro- degenerated due to scanned k-points match with Bragg cal lattice vectors of SC generate the coordinate points diffractions. However the total summation of Wk (I)s is of the twisted irreducible BZ of PCs in Figure 3(a), ' 1. This characteristic of our unfolding process allows where the red labels identify corresponding BZ of G, seeing two main aspects in comparison with previous while blue labels do with WS2 . unfolding methods with DFT [16, 17, 18, 19, 20, 21] Using the fractional coordinates K= (1/3, 1/3, 0) or tight-binding calculations [22, 23, 24]: first, we and M= (1/2, 0, 0) of the SC, we calculated the can determine a relationship between the reciprocal electronic band structure of the heterostructure space and possible Bragg diffractions in the band G/WS2 scanning the Γ−K−M−Γ path, as it is shown structure; and, second, electronic degeneracy and/or in Figure 3(b). Close to the Fermi level, the typical hybridization can be quantified via Wk (I). These linear energy dispersion associated to the Dirac cones capabilities of our unfolding code are owing to the of G is observed at the Γ point. The preserving of the basis-set used and the normalization process and they linear dispersion in 2D van der Waals heterostructures will be discussed with more details later for the G/WS2 has been reported in other works [32, 47, 35, 48, 49] heterostructure. √ √ by means of ab initio calculations. Furthermore, close We performed the same analysis for 13 × 13 to the Fermi level the degeneracy is broken at the G layer with θ = 16.1◦ . Figure 1(d) sketches the Γ-M path, while it is preserved at the Γ-K path. reciprocal space, the first and the irreducible BZ of For energies lower than -1 eV and greater than +1 SC with the same orientation as those in Figure 1(a). eV there are many electronic states associated to G The inset shows the corresponding SC layer. Again, and WS2 as well, which means that the electronic the typical linear dispersion of electronic states of states of G around the Fermi level fall inside the Unfolding method for periodic twisted systems 7 values are obtained, the effective electronic bands (a ) K (EEBs) for G and WS2 are plotted, as it is shown in K M Figure 4(a) and 4(b), respectively. In general, EEBs show large difference and discontinuities with respect M to the electronic band structure of pristine PC of G K [50, 51] and WS2 [52, 33, 34] as a consequence of Wk (I) M is not exactly equal to 1 in every projected band. Such deviations can be easily observed comparing EEB of Γ G in Figure 5(a) with the band structure of pristine PC of G in Figure 1(c) and 1(f). This means that the symmetry of each layer reconstruction is not equal to the original symmetry of pristine PC. In particular, the translation symmetry is broken in this 2D G/WS2 heterostructure due to interlayer interaction. In the case of G, the unfolding process shows that the Dirac cone around the Fermi level remains 3 unchanged but it is now located at the K point of 2 (b ) PC instead of Γ of SC. Furthermore, it falls inside of the WS2 bandgap. At the Fermi level, there is no 1 bandgap associated to G bands despite the applied 0 symmetrical strain, ∆, which is in agreement with -1 other theoretical results for pristine G undergoes a small symmetrical strain distribution [53]. On the e n e rg y (e V ) -2 other hand, we detected a change in Fermi velocity -3 (v F ) in G by the strain. A slight increase of 1.54% -4 was measured, taking as reference the unfolded band -5 structure of G SC in Figure 1(c). The previous v F -6 value was roughly estimated only for carriers in the valence band. Our results are in good agreement with -7 other experimental and theoretical results [54, 53] with -8 potential applications in strain engineering and flexible -9 electronics. For energies below -3 eV in figure 5(a), Γ K M Γ there are two bands composed of mixed C(2s,2p x ,2p y ) orbitals giving a sp 2 hybridization. These bands are Figure 3. (a) Reciprocal space of the G/WS2 system. The non-perturbed showing a Wk (I) ' 1 along the full purple (joined by the red points), green (green points) hexagons k -path. Thus, since we performed DFT calculations represent the reciprocal space and the first BZ. The irreducible with a localized basis-set, we can extract partial EEBs BZ of SC is represented by a shaded green area. The red and blue hexagons represent the first BZ of PC of G and WS2 , taking into account only specific orbitals or atoms of respectively, and the red (blue) triangle denotes the path of the the SC [27]. This result reflects the usefulness of irreducible BZ for G (WS2 ). K and M labels are in green, red and the present unfolding code based on atomic orbitals blue for the G/WS2 , and both PC of G and WS2 , respectively, and normalized wavefunctions for studying of any while the Γ-point (in black) is a common reference among the three systems. (b) Electronic band structure of G/WS2 with the perturbed periodic system. At the energies plotted in Fermi level set to 0 eV. Figure 5, only bands associated with C(2p z ) orbitals changes drastically. It is observed a splitting bands near to Γ point at energies around -8 eV, which is band gap of WS2 , confirming that the full system is similar to the band structure of graphite [55]. Also, a semimetal. From a simple analysis of Figure 3(b), it different minigaps in the valence bands are observed. is quite difficult to extract more physical information. Such minigaps observed at low energies are in fair Nonetheless, applying our unfolding method to the agreement with experimental results of G adsorbed on electronic band structure of G/WS2 onto each PC, it MoS2 [56, 57], as well as other theoretical calculations is possible to distinguish the physical consequences of [26, 25] with θ = 0◦ . We will discuss more about the mutual interactions among layers. formation of minigaps latter. Another advantage of Let us apply our unfolding code to the G/WS2 the unfolding method for studying these systems in structure to project the electronic states to both PC comparison, for instance, with ARPES [15] technique using the corresponding Γ-K-M-Γ path for G and is that, we can see what happens with the unoccupied WS2 [see Figure 3(a)]. Once the corresponding Wk (I) Unfolding method for periodic twisted systems 8 states as well. In this case, in Figure 4(a) there is path and those located at 0.22 and near to 1/3 and not an equal relation of Wk (I) for each wave function 1/2 starting from M → Γ path are Bragg diffractions. (I ) and k -point, which means there is a particle-hole The rest of minigaps have a different physical origin. asymmetry and, as result of WS2 interaction, the The discontinuity in K-M path looks like an electronic mobility of holes and electrons would behave different resonance showing coupling between monolayers. The in G [24, 58]. Bragg diffraction in WS2 located at 0.6 (1/3) Γ → K Now, let us discuss the unfolding process (M → Γ) path shows minigaps around these points for projecting the electronic states of G/WS2 to the PC energies from -5.5 eV to -7.5 eV. Also, other minigaps of WS2 . The direct-transition electronic band gap are observed near the K (on left-side) and M points (on located at the K point of PC is preserved, which means right-side) with energy -3.5 eV indicated by arrows. A that, under our methodology, the vertical stacking of discontinuity is observed in 0.6 Γ → K path in ∼2.75 G with WS2 does not change the electronic properties eV below the Fermi level. Finally, the number of Bragg of WS2 near the Fermi energy. In the Γ-K path and diffractions in G is greater than WS2 by the size of BZs around -1.5 eV, a small discontinuity is observed in the of PCs. The origin of all different interactions will be last occupied valence band (second unnocupied band), discussed in a future work. green arrows in Figure 4(b). Going down in energy, all We want to stress that, strictly speaking if both electronic bands of WS2 change from -2 eV to -6 eV, layers do not interact between them, Wk (I) must be which means that W(d ) and S(p) orbitals modify their equal to 1 as in Figures 1(c) and 1(f). To corroborate occupation just by G adsorption. Another interesting this, we analyze the G/WS2 heterostructure with the result obtained through the unfolding process is the interlayer electronic interactions almost zero. This energy splitting about ∼1 eV of the band at Γ−K and condition is possible putting both layers separated M−Γ path from -7.5 eV to -8.5 eV [see Figure 4(b)], by an interlayer distance larger than the equilibrium also observed in the PC of G [see Figure 4(a)], which one d = 3.28 Å. Hence, using the optimized G/WS2 means an interlayer interaction repulsion between the structure, we impose an interlayer distance of 5.29 C(2pz ) orbitals of G and S(3p) and W(5d) orbitals Å and, from a single-point DFT calculation, the of WS2 . Furthermore, it is noticeable that there is electronic states of the SC are obtained. Then, we a complementary relationship between the two EEBs applied our unfolding code to project the electronic of G and WS2 around the Γ point for energies from - states to both PC using the corresponding Γ-K-M-Γ 6 eV to -8.5 eV. Wk (I) of G bands [see Figure 4(a)] path for G and WS2 , denoted in Figure 3(a) by red complements the Wk (I) of WS2 [see Figure 4(b)], and blue labels, respectively. Once the corresponding even when the scanned paths are a bit different near Wk (I) values are obtained, the EEBs are plotted. the Γ point. This complementary characteristic of Figure 5(a) shows with red and blue points the EEB Wk (I) between two monolayers is another advantage of G and WS2 projected onto the Γ-K-M-Γ path for of our unfolding process, in comparison with other G, while 5(b) shows the corresponding EEB projected unfolding methods based on Green’s functions for onto the Γ-K-M-Γ path for WS2 . Analyzing the EEBs periodic twisted bilayers [24, 26]. of G and WS2 in Figure 5(a) and 5(b), respectively, Regarding the electronic minigaps observed in they are almost entirely recovered since Wk (I) ' 1 EEBs of G, it is worth pointing out that their number, in the whole energy range and the scanned k path as size and position in this G/WS2 heterostructure with well. This is an outstanding theoretical proof which θ = 10.9◦ is greater than other theoretical reports shows that the strong interlayer interaction, given by when θ = 0◦ [26, 25]. This means that θ has a vdW forces, leads to hybridization and the emerging strong influence on the electronic properties in 2D of minigaps discussed above. However, some Kohn- heterostructures, as in the case of G bilayer [5, 6, 7, 8]. Sham orbitals show Wk (I) < 1, which is associated We have identified two kind of electronic interactions. to interlayer distortions originated by the structural One is related with the electronic gap generation optimization carried out at d = 3.28 Å. Comparing with a direct transition. The other one seems an both Figures 5(a) and (b), the EEBs of G are quite electronic resonance (or discontinuity) between out- different as well as corresponding of WS2 . Such of-plane orbitals of each monolayer. It is possible to differences are originated because in Figure 5(a) the identify both behaviors in the EEBs of G. In Figure states of WS2 are projected onto the BZ of G, while in 4(a), the minigaps are observed along the Γ-K and M- Figure 5(b) the states of G are projected onto the BZ Γ paths and indicated by white arrows, while electronic of WS2 . So, a wrong choice of the k path could induce resonances are detected in the K-M path indicated a misunderstanding or even a bad interpretation of the by green arrows near the M point. Based on the mutual influence between layers. Additionally, in every reciprocal space sketched in Figure 3, we can say that plot in Figure 5 we can observe only the electronic the gaps located in 1/2 and 0.7 starting from Γ → K states from the reference PCs without the presence of Unfolding method for periodic twisted systems 9 0 1 0 1 0 0 . 0 0. 0 0 . 0 0 . 0 0. 0 0 . 1 0 . 1 0. 1 0 . 1 0 . 1 0. 1 0 . 1 0 . 1 0. 1 0 . 1 0 . 2 0. 2 0 . 2 0 . 2 0. 2 0 . 2 0 . 2 0. 2 0 . 2 0 . 2 0. 3 0 . 3 0 . 3 0. 3 0 . 3 0 . 3 0. 3 0 . 3 0 . 3 0. 3 0 . 4 0 . 4 0. 4 0 . 4 0 . 4 0. 4 0 . 4 0 . 4 0. 4 0 . 4 0 . 5 0. 5 0 . 5 0 . 5 0. 5 0 . 5 0 . 5 0 . 5 0. 5 0 . 5 0 . 6 0. 6 0 . 6 0 . 6 0. 6 0 . 6 0 . 6 0. 6 0 . 6 0 . 6 0. 7 0 . 7 0 . 7 0. 7 0 . 7 0 . 7 0. 7 0 . 7 0 . 7 0. 7 0 . 8 0 . 8 0. 8 0 . 8 0 . 8 0. 8 0 . 8 0 . 8 0. 8 0 . 8 0 . 9 0. 9 0 . 9 0 . 9 0. 9 0 . 9 0 . 9 0. 9 0 . 9 1 . 9 1. 0 1 . 0 1 . 0 1. 0 1 . 0 . 0 0 0 . 0 0. 0 0 . 0 0. 0 0 . 0 0. 1 0 . 1 0 . 1 0. 1 0 . 1 0. 1 0 . 1 0. 1 0 . 1 0 . 1 0. 2 0 . 2 0. 2 0 . 2 0. 2 0 . 2 0 . 2 0. 2 0 . 2 0. 2 0 . 3 0. 3 0 . 3 0 . 3 0. 3 0 . 3 0. 3 0 . 3 0. 3 0 . 3 0 . 4 0. 4 0 . 4 0. 4 0 . 4 0. 4 0 . 4 0 . 4 0. 4 0 . 4 0. 5 0 . 5 0. 5 0 . 5 0. 5 0 . 5 0 . 5 0. 5 0 . 5 0. 5 0 . 6 0. 6 0 . 6 0 . 6 0. 6 0 . 6 0. 6 0 . 6 0. 6 0 . 6 0 . 7 0. 7 0 . 7 0. 7 0 . 7 0. 7 0 . 7 0 . 7 0. 7 0 . 7 0. 8 0 . 8 0. 8 0 . 8 0 . 8 0. 8 0 . 8 0. 8 0 . 8 0. 8 0 . 9 0 . 9 0. 9 0 . 9 0. 9 0 . 9 0. 9 0 . 9 0 . 9 1. 9 1 . 0 1. 0 1 . 0 1. 0 1 . 0 . 0 3 3 (a ) (b ) 2 2 1 1 0 0 -1 -1 e n e rg y (e V ) e n e rg y (e V ) -2 -2 -3 -3 -4 -4 -5 -5 -6 -6 -7 -7 -8 -8 -9 -9 Γ 1 /4 1 /2 0 .7 K 1 /2 M 0 .2 2 1 /3 1 /2 2 /3 Γ Γ 1 /3 0 .6 K 2 /3 M 1 /3 0 .6 Γ Figure 4. Calculated EEBs for (a) G and (b) WS2 PCs with an interlayer distance of 3.28 Å. The vertical dashed lines represent Bragg points and lines in the reciprocal space. The fractional values along the paths mean the position at which diffractions occur. Minigaps with an origin different from Bragg reflections are indicated by white arrows, while resonances states are indicated by green arrows. The number of k-points scanned along Γ →K, K→M and M→ Γ were 200, 100 and 200, respectively. The Wk (I) is proportional to symbol size and color. The Fermi level is set to 0 eV. additional bands or the so-called ghost bands [26]. This the interlayer distance is reduced until its optimal is another advantage of the present unfolding method value. The results reveal a strong modulation of the based on atomic orbitals and DFT. Comparing the electronic structure of graphene by tungsten disulfide EEB of strained G in Figure 5(a) with ideal G bands due to the relatively weak van der Waals interactions. in Figure 1(c) around -8 eV and close to Γ point, In graphene, the main contributions are by the orbital a shift energy of the electronic state are observed, normal to the monolayer planes, while the in-plane which is directly related with the strain condition of orbitals form strong bonds without any perturbation. G monolayer. Meanwhile WS2 is a direct bandgap However, around the Fermi level the Dirac cone semiconductor at the K point in its PC, see blue circles remains almost unaffected with small changes in the in Figure 5(b). Fermi velocity of the carriers assisted by symmetric lattice strain. The electronic properties of tungsten 6. Conclusions disulfide are preserved near the Fermi energy, mainly the direct band gap at the same point. The nature of In summary, total energy density functional theory minigaps comes from Bragg conditions in the supercell calculations in combination with an extended unfolding and the interlayer coupling between graphene and band structure approach have been applied to tungsten disulfide monolayers. study conventional and non-conventional graphene supercells, and twisted bilayer of non-conventional 6.1. Acknowledgments graphene and tungsten disulfide supercells. The unfolding method works for periodically twisted Authors acknowledge the partial support from systems, showing that, the full band structure of DGAPA-UNAM project PAPIIT IN109618 and CONA- the primitive cells are recovered if ideal supercells CyT project A1-S-39242. F.S.O. acknowledges to of graphene are studied. Furthermore, we addressed DGAPA- UNAM for postdoctoral fellowship. The nu- the topic of interaction by proximity effects in merical calculations were carried out in clusters Bak- a system of graphene and tungsten disulfide with tum and Flatland at the Institute of Physics (UNAM). double-periodicity. Several advantages of developed method have been presented and the physics of References interlayer interactions is easier to explain. When [1] Geim A K and Grigorieva I V 2013 Nature 499 419 graphene is adsorbed on tungsten disulfide, many [2] Bhimanapati G R, Lin Z, Meunier V, Jung Y, Cha J, Das electronic interactions arise between monolayers when S, Xiao D, Son Y, Strano M S, Cooper V R et al. 2015 ACS nano 9 11509–11539 Unfolding method for periodic twisted systems 10 0 1 0 1 0 0 . 0 0. 0 0 . 0 0 . 0 0 . 0 0. 1 0 . 1 0 . 1 0. 1 0 . 1 0 . 1 0 . 1 0. 1 0 . 1 0 . 1 0 . 2 0. 2 0 . 2 0 . 2 0. 2 0 . 2 0 . 2 0 . 2 0. 2 0 . 2 0 . 3 0. 3 0 . 3 0 . 3 0 . 3 0. 3 0 . 3 0 . 3 0. 3 0 . 3 0 . 4 0 . 4 0. 4 0 . 4 0 . 4 0 . 4 0. 4 0 . 4 0 . 4 0. 4 0 . 5 0 . 5 0 . 5 0. 5 0 . 5 0 . 5 0. 5 0 . 5 0 . 5 0 . 5 0. 6 0 . 6 0 . 6 0. 6 0 . 6 0 . 6 0 . 6 0. 6 0 . 6 0 . 6 0 . 7 0. 7 0 . 7 0 . 7 0. 7 0 . 7 0 . 7 0 . 7 0. 7 0 . 7 0 . 8 0. 8 0 . 8 0 . 8 0 . 8 0. 8 0 . 8 0 . 8 0. 8 0 . 8 0 . 9 0 . 9 0. 9 0 . 9 0 . 9 0 . 9 0. 9 0 . 9 0 . 9 1. 9 1 . 0 1 . 0 1 . 0 1. 0 1 . 0 . 0 0 0 . 0 0. 0 0 . 0 0 . 0 0. 0 0 . 1 0 . 1 0. 1 0 . 1 0 . 1 0. 1 0 . 1 0 . 1 0. 1 0 . 1 0 . 2 0. 2 0 . 2 0 . 2 0. 2 0 . 2 0 . 2 0. 2 0 . 2 0. 2 0 . 3 0 . 3 0. 3 0 . 3 0 . 3 0. 3 0 . 3 0 . 3 0. 3 0 . 3 0 . 4 0. 4 0 . 4 0 . 4 0. 4 0 . 4 0 . 4 0. 4 0 . 4 0 . 4 0. 5 0 . 5 0 . 5 0. 5 0 . 5 0 . 5 0. 5 0 . 5 0 . 5 0. 5 0 . 6 0 . 6 0. 6 0 . 6 0 . 6 0. 6 0 . 6 0 . 6 0. 6 0 . 6 0 . 7 0. 7 0 . 7 0 . 7 0. 7 0 . 7 0 . 7 0. 7 0 . 7 0. 7 0 . 8 0 . 8 0. 8 0 . 8 0 . 8 0. 8 0 . 8 0 . 8 0. 8 0 . 8 0 . 9 0. 9 0 . 9 0 . 9 0. 9 0 . 9 0 . 9 0. 9 0 . 9 1 . 9 1. 0 1 . 0 1 . 0 1. 0 1 . 0 . 0 3 3 (a ) (b ) 2 2 1 1 0 0 -1 -1 e n e rg y (e V ) e n e rg y (e V ) -2 -2 -3 -3 -4 -4 -5 -5 -6 -6 -7 -7 -8 -8 -9 -9 Γ K M Γ Γ K M Γ Figure 5. Calculated EFBSs for (a) G and (b) WS2 PCs. An interlayer distance of 5.29 Å reflects weak electronic interactions between G and WS2 . The projection of G/WS2 SC states onto the WS2 (G) PC using the k-path of G (WS2 ) is plotted on left (right) panel with blue (red) circles. The number of k-points scanned along Γ →K, K→M and M→ Γ were 200, 100 and 200, respectively. The Wk (I) is proportional to symbol size and color. The upper color bar from green to red (yellow to blue) is for Wk (I) of G (WS2 ). The Fermi level is set to 0 eV. [3] Lu C C, Lin Y C, Liu Z, Yeh C H, Suenaga K and Chiu Condens. Matter 25 345501 P W 2013 ACS nano 7 2587–2594 [20] Medeiros P V, Stafström S and Björk J 2014 Phys. Rev. B [4] Dean C R, Wang L, Maher P, Forsythe C, Ghahari F, Gao 89 041407 Y, Katoch J, Ishigami M, Moon P, Koshino M et al. 2013 [21] Medeiros P V, Tsirkin S S, Stafström S and Björk J 2015 Nature 497 598 Phys. Rev. B 91 041116 [5] Kim C J, Sánchez-Castillo A, Ziegler Z, Ogawa Y, Noguez [22] Boykin T B and Klimeck G 2005 Phys. Rev. B 71(11) C and Park J 2016 Nat. Nanotechnol. 11 520 115215 [6] Rickhaus P, Wallbank J, Slizovskiy S, Pisoni R, Overweg [23] Boykin T B, Kharche N, Klimeck G and Korkusinski M H, Lee Y, Eich M, Liu M H, Watanabe K, Taniguchi T, 2007 J. Phys.: Condens. Matter 19 036203 Ihn T and Ensslin K 2018 Nano Lett. 18 6725–6730 [24] Nishi H, Matsushita Y i and Oshiyama A 2017 Phys. Rev. [7] Cao Y, Fatemi V, Fang S, Watanabe K, Taniguchi T, B 95(8) 085420 Kaxiras E and Jarillo-Herrero P 2018 Nature 556 43 [25] Le N B, Huan T D and Woods L M 2016 ACS Appl. Mater. [8] Cao Y, Fatemi V, Demir A, Fang S, Tomarken S L, Luo Interfaces 8 6286–6292 J Y, Sanchez-Yamagishi J D, Watanabe K, Taniguchi T, [26] Matsushita Y i, Nishi H, Iwata J i, Kosugi T and Oshiyama Kaxiras E et al. 2018 Nature 556 80 A 2018 Phys. Rev. Materials 2 010801(R) [9] Roy K, Padmanabhan M, Goswami S, Sai T P, Ramalingam [27] Deretzis I, Calogero G, Angilella G and La Magna A 2014 G, Raghavan S and Ghosh A 2013 Nat. Nanotechnol. 8 EPL 107 27006 826 [28] Haverkort M, Elfimov I and Sawatzky G 2011 arXiv [10] Zhang W, Chuu C P, Huang J K, Chen C H, Tsai M L, preprint arXiv:1109.4036 Chang Y H, Liang C T, Chen Y Z, Chueh Y L, He J H [29] Arnold H 2006 International Tables for Crystallography. et al. 2014 Sci. Rep. 4 3826 Vol. A, Chapter 5.1 (Netherlands: Springer) [11] Georgiou T, Jalil R, Belle B D, Britnell L, Gorbachev [30] Soler J M, Artacho E, Gale J D, Garcı́a A, Junquera R V, Morozov S V, Kim Y J, Gholinia A, Haigh S J, J, Ordejón P and Sánchez-Portal D 2002 J. Phys.: Makarovsky O et al. 2013 Nat. Nanotechnol. 8 100 Condens. Matter 14 2745 [12] Yu W J, Li Z, Zhou H, Chen Y, Wang Y, Huang Y and [31] Press W H, Teukolsky S A, Vetterling W T and Flannery Duan X 2013 Nat. Mater. 12 246 B P 1992 Numerical recipes in Fortran 77, vol. 1 2d ed. [13] Wang Z, Ki D K, Chen H, Berger H, MacDonald A H and (Cambridge: Cambridge Univ. Press) Morpurgo A F 2015 Nat. Commun. 6 8339 [32] Wang Z, Chen Q and Wang J 2015 J. Phys. Chem. C 119 [14] Ashcroft N W and Mermin N D 1976 Solid state physics 4752–4758 (New York: Holt, Rinehart and Winston) [33] Klein A, Tiefenbacher S, Eyert V, Pettenkofer C and [15] Yang H, Liang A, Chen C, Zhang C, Schroeter N B and Jaegermann W 2001 Phys. Rev. B 64 205416 Chen Y 2018 Nat. Rev. Mater. 3 341–343 [34] Yun W S, Han S, Hong S C, Kim I G and Lee J 2012 Phys. [16] Ku W, Berlijn T and Lee C C 2010 Phys. Rev. Lett. 104(21) Rev. B 85 033305 216401 [35] Li L and Zhao M 2014 J. Phys. Chem. C 118 19129–19138 [17] Popescu V and Zunger A 2010 Phys. Rev. Lett. 104(23) [36] Momma K and Izumi F 2011 J. Appl. Crystal. 44 1272– 236403 1276 [18] Popescu V and Zunger A 2012 Phys. Rev. B 85(8) 085201 [37] Perdew J P and Zunger A 1981 Phys. Rev. B 23 5048 [19] Lee C C, Yamada-Takamura Y and Ozaki T 2013 J. Phys.: [38] Hamann D, Schlüter M and Chiang C 1979 Phys. Rev. Lett. Unfolding method for periodic twisted systems 11 43 1494 [39] Kleinman L and Bylander D 1982 Phys. Rev. Lett. 48 1425 [40] Garcı́a-Gil S, Garcı́a A, Lorente N and Ordejón P 2009 Phys. Rev. B 79 075441 [41] Hauschild A, Karki K, Cowie B, Rohlfing M, Tautz F and Sokolowski M 2005 Phys. Rev. Lett. 94 036106 [42] Lee K, Yu J and Morikawa Y 2007 Phys. Rev. B 75 045402 [43] Monkhorst H J and Pack J D 1976 Phys. Rev. B 13 5188 [44] Bengtsson L 1999 Phys. Rev. B 59 12301 [45] Kaloni T, Kou L, Frauenheim T and Schwingenschlögl U 2014 Appl. Phys. Lett. 105 233112 [46] Li S S and Zhang C W 2013 J. Appl. Phys. 114 183709 [47] Di Felice D, Abad E, González C, Smogunov A and Dappe Y 2017 J. Phys. D: Appl. Phys. 50 17LT02 [48] Sa B, Sun Z and Wu B 2016 Nanoscale 8 1169–1178 [49] Wang Y and Ding Y 2015 Phys. Chem. Chem. Phys. 17 27769–27776 [50] Şahin H, Cahangirov S, Topsakal M, Bekaroglu E, Akturk E, Senger R T and Ciraci S 2009 Phys. Rev. B 80 155453 [51] Boukhvalov D, Katsnelson M and Lichtenstein A 2008 Phys. Rev. B 77 035427 [52] Wang Q H, Kalantar-Zadeh K, Kis A, Coleman J N and Strano M S 2012 Nat. Nanotechnol. 7 699 [53] Gui G, Li J and Zhong J 2008 Phys. Rev. B 78 075435 [54] Huang M, Yan H, Heinz T F and Hone J 2010 Nano Lett. 10 4074–4079 [55] Marinopoulos A G, Reining L, Rubio A and Olevano V 2004 Phys. Rev. B 69 245419 [56] Pierucci D, Henck H, Avila J, Balan A, Naylor C H, Patriarche G, Dappe Y J, Silly M G, Sirotti F, Johnson A C et al. 2016 Nano Lett. 16 4054–4061 [57] Coy Diaz H, Avila J, Chen C, Addou R, Asensio M C and Batzill M 2015 Nano Lett. 15 1135–1140 [58] DaSilva A M, Jung J, Adam S and MacDonald A H 2015 Phys. Rev. B 91 245422