key: cord-0901253-doek7so9
authors: Thompson, J. J. P.; Pei, D.; Peng, H.; Wang, H.; Channa, N.; Peng, H. L.; Barinov, A.; Schroter, N. B. M.; Chen, Y.; Mucha-Kruczy'nski, M.
title: Determination of interatomic coupling between two-dimensional crystals using angle-resolved photoemission spectroscopy
date: 2020-07-17
journal: Nature communications
DOI: 10.1038/s41467-020-17412-0
sha: 50b889c37d7c9ff3c4873831f62ff730d91ba6b9
doc_id: 901253
cord_uid: doek7so9

Lack of directional bonding between two-dimensional crystals like graphene or monolayer transition metal dichalcogenides provides unusual freedom in selection of components for vertical van der Waals heterostructures. However, even for identical layers, their stacking, in particular the relative angle between their crystallographic directions, modifies properties of the structure. We demonstrate that the interatomic coupling between two two-dimensional crystals can be determined from angle-resolved photoemission spectra of a trilayer structure with one aligned and one twisted interface. Each of the interfaces provides complementary information and together they enable self-consistent determination of the coupling. We parametrize interatomic coupling for carbon atoms by studying twisted trilayer graphene and show that the result can be applied to structures with different twists and number of layers. Our approach demonstrates how to extract fundamental information about interlayer coupling in a stack of two-dimensional crystals and can be applied to many other van der Waals interfaces.

Following the isolation of graphene (a layer of carbon atoms arranged in regular hexagons) in 2004 [1] , many other atomically thin two-dimensional crystals have been produced and can be stacked in a desired order on top of each other. In contrast to conventional heterostructures, in which chemical bonding at interfaces between two materials modifies their properties and requires lattice matching for stability, stacks of two-dimensional crystals are held together by weak forces without directional bonding. As a result, any two of these materials can be placed on top of each other, providing extraordinary design flexibility [2] [3] [4] . Moreover, subtle changes in atomic stacking, especially the angle between the crystallographic axes of two adjacent layers, can have big impact on the properties of the whole heterostructure, with examples including the observation of Hofstadter's butterfly [5, 6] and interfacial polarons [7] in graphene/hexagonal boron nitride heterostructures, interlayer excitons in transition metal dichalcogenide bilayers [8, 9] , appearance of superconductivity in magicangle twisted bilayer graphene [10, 11] and explicit twistdependence of transport measurements in rotatable heterostructures [12] [13] [14] . Phenomena like these arise because the misalignment of two crystals changes the * M.Mucha-Kruczynski@bath.ac.uk atomic registry at the interface and hence tunes the spatial modulation of interlayer interaction. Consequently, understanding the coupling between two two-dimensional materials at a microscopic level is crucial for efficient design of van der Waals heterostructures.

The impacts of a twisted interface and modulated interlayer coupling on the electronic properties of twodimensional crystals include band hybridization [15] [16] [17] , band replicas and minigaps due to scattering on moiré potential [15, 18, 19] , charge transfer and vertical shifting of bands [17, 20, 21] as well as changes of the effective masses [17, 20] . Variations in the interlayer coupling as a function of the twist angle, θ, were probed for example using photoluminescence, Raman and angleresolved photoemission (ARPES) spectroscopies [20, [22] [23] [24] . Here, we use the last of those methods to image directly the electronic bands in trilayer graphene with one perfect and one twisted interface. From our data, we extract the interatomic coupling, t(r, z), describing coupling between two carbon atoms separated by a vector r 3D = (r, z) = (x, y, z). Such coupling functions, usually based on comparisons to ab initio calculations, can be used to determine electron hoppings in tight-binding [25, 26] and continuum [27, 28] models of corresponding van der Waals interfaces at any twist angle. We show that t(r, z) determined purely by measurements on one of the structures accurately describes electronic dispersions obtained for stacks with different θ and number of layers, providing an experimentally verified set of parameters to model twistronic graphene. Our approach makes use of the fact that a trilayer structure is the thinnest stack that can contain both a perfect and twisted interface. The former, due to translational symmetry, can be straightforwardly described in the real space using t(r, z). At the same time, the impact of the moiré pattern formed at the latter can be captured in the reciprocal space by considering scattering by moiré reciprocal vectors on the momentum-dependent potentialt(q, z) which is a twodimensional Fourier transform F[t(r, z)] of t(r, z) (see the comparison of the two cases in Fig. 1(a) ). As a consequence, this method should enable determination of interatomic couplings for all van der Waals interfaces for which moiré effects were observed.

ARPES of twisted trilayer graphene We grew our graphene trilayers on copper foil using chemical vapour deposition [29, 30] . The inset of Fig. 1(b) shows the intensity map of copper d-band photoelectrons which are attenuated differently by the overlying graphene layers depending on their number. This provides means to identify all of the layers in our stack, shown in the inset with different shades of gray and indicated with the red arrows. As depicted schematically in the main panel of Fig. 1(b) , the bottom two layers form a Bernal bilayer (2L) while the crystallographic axes of the top monolayer (1L) are rotated by an angle θ with respect to those of the layer underneath. As a result, the Brillouin zones corresponding to the bilayer and monolayer are also rotated with respect to each other, Fig. 1(c) . We focus here on the vicinity of one set of the corners of the two Brillouin zones, which we denote K 2 and K 1 , for the bilayer and monolayer, respectively. The separation between these two points, dependent on the twist angle, defines an effective superlattice Brillouin zone, indicated in orange in the inset of Fig. 1(c) .

In Fig. 2 (a), we present ARPES intensity along a cut in the k-space connecting K 2 and K 1 , with the energy reference point set to the linear crossing (Dirac point) at K 1 . Close to each corner, the intensity reflects the low-energy band structures of unperturbed 2L and 1L. Because the bilayer flake is below the monolayer, signal from the former is attenuated due to the electron escape depth effect. In between the two spectra, coupling of the two crystals leads to anticrossings of the bands and opening of minigaps (marked as ε g I and ε g II in the figure). As the size of the superlattice Brillouin zone depends on the twist angle, the energy positions of the minigaps also depend on θ. Moreover, the magnitudes of the minigaps depend on the interlayer coupling between the bilayer and monolayer and also, in principle, vary with θ. However, fundamentally, all of the features in our spectrum originate in interactions between carbon atoms, be it in the same or different layers, at the twisted or aligned interface. This provides us with an opportunity to study the interatomic coupling t(r, z) in carbon materials. (1). Inset shows photoemission intensity from copper substrate which is attenuated by graphene layers above, providing a measure of graphene layer number. The red arrows indicate each of the graphene layers in the trilayer stack and the cyan line corresponds to the distance of 10 µm. (c) Brillouin zones of the Bernal bilayer (black) and rotated monolayer (blue) with bilayer and monolayer graphene low-energy electronic spectra shown in the vicinities of one set of the Brillouin zone corners. The inset depicts in orange the superlattice Brillouin zone and the cyan line indicates the k-space path cuts along which are presented in Fig. 2 (a) and 3.

Parametrizing carbon-carbon interaction potential In order to understand our data, we use a generic Hamiltonian for a van der Waals heterostructure comprised of three layers of the same two-dimensional crystal

In this Hamiltonian, the diagonal block,Ĥ 0 (θ i , ε i ) de-scribes the i-th layer at a twist angle θ i , with on-site energies of atomic sites in this layer, ε i . Here, because only the relative twist between any two adjacent layers is important, we have θ 1 = θ 2 = 0 and θ 3 = θ. Also, our choice of energy reference point is equivalent to ε 3 = 0 and we introduce potential energy difference, 2u = ε 1 − ε 2 , as well as average energy, ∆ = (ε 1 + ε 2 )/2, of layers 1 and 2 (the charge transfer between the copper foil and the graphene layers giving rise to u = ∆ = 0 is discussed in more detail in Ref. 29) . For graphene, the intralayer blocksĤ 0 can be straight-forwardly described using a tight-binding model [31] for a triangular lattice with two inequivalent atomic sites, A and B, per unit cell and nearest neighbour coupling between them γ 0 ≡ −t(r AB , 0), where r AB is a vector connecting neighbouring A and B atoms with the carbon-carbon bond length |r AB | = 1.46Å.

Of more importance for us, however, are the offdiagonal blocksT (θ i − θ i−1 ) which capture the twistdependent interlayer interactions between adjacent layers (we neglect the interaction between the bottom and the top layers which is at least an order of magnitude weaker [32] ). As the bottom two layers are stacked according to the Bernal stacking, a real-space description of the interlayer interaction blockT (0) is possible with the leading coupling t(0, c 0 ) ≡ γ 1 , with interlayer distance c 0 = 3.35 A, due to atoms with neighbours directly above or below them, as shown in Fig. 1 (a) [33] . In contrast, we describe the coupling between the twisted layers, i = 2, 3, in the reciprocal space based on electron tunnelling from a state with wave vector k in layer 2 to a state with wave vector k in layer 3 with the requirement that crystal momentum is conserved [34, 35] , k + G = k + G , where G and G are the reciprocal vectors of layers 2 and 3, respectively. The strength of a given tunnelling process is set by the two-dimensional Fourier transform, F[t(r, z)] =t(q, z), of the real-space coupling t(r, z) so that

where τ = (−|r AB |, 0) andR θ is a matrix of clockwise rotation by angle θ (see Supplementary Note 1 for more details on the construction of the HamiltonianĤ). The uniqueness of a trilayer with one perfect and one twisted interface (as exemplified in Fig. 1 (a) for the case of graphene) lies in the fact that the HamiltonianĤ contains interlayer blocks based on both the real-space (T (0)) and reciprocal-space (T (θ)) descriptions which provide complementary information and at the same time are related to each other because of the Fourier transform connection between t(r, z) andt(q, z). Because of this, comparison of the photoemission data with the spectrum calculated based on Equation (1) provides more information about the interatomic coupling t(r, z) than struc- Table I . Right: Twodimensional Fourier transformt(|q|, c0) of the interatomic coupling, t(|r|, c0), as a function of wave vector |q|.

tures with one type of interface only. For our graphene trilayer, we compute the miniband spectrum ofĤ (see Methods for more details) assuming a Slater-Koster-like two-centre ansatz for t(r, z) [25] ,

where V π and V σ represent the strength of the π and σ bonding [36] , respectively, and α π and α σ their decay with increasing interatomic distance.

In fitting our numerical results to the experimental data in Fig. 2(a) , we first determine the position of 1L Dirac point what sets the ε = 0 reference point. We then use the electronic band gap at K 2 to fix the electrostatic potential 2u and position the bilayer neutrality point halfway in the gap, establishing the potential energy shift ∆. We obtain the in-plane nearest neighbour hopping γ 0 from the slope of the 1L linear dispersion close to the Dirac point at K 1 while the direct interlayer coupling γ 1 is set by the splitting of the 2L lower valence band from the neutrality point at K 2 . Finally, the decay constants α π and α σ are found numerically using the constraints that (i) the magnitudes of the gaps ε g I and ε g II in Fig. 2 (a) match the experimental data and (ii) in the limit of θ = 0,T (θ) from Equation (2) converges to the real-space form ofT (0) as used for coupling between the Bernal stacked layers (see Supplementary Note 2 for further discussion).

The miniband spectrum resulting from our model is shown in red dashed lines in Fig. 2(a) , the functions t(|r|, c 0 ) andt(|q|, c 0 ) are plotted in Fig. 2 (b) and the corresponding values of the parameters γ 0 , γ 1 , α π and α σ are summarized in Table I . The interatomic potential we obtain decays more rapidly in the real space (and hence slower in the reciprocal space) than suggested by computational results [25] . Importantly, parametrization of t(r, z) does not depend on the twist angle and so should be applicable to other graphene stacks with twisted interfaces. It also does not depend on the doping level because, for the relevant range of electric fields, the electrostatic energies ∆ and u do not modify the electron hoppings. At the same time, once these energies are determined for a particular stack, their influence on the band structure (shifting of the positions and magnitudes of anticrossings) is captured through the Hamilto-nianĤ. To confirm applicability of a single parametrization of t(r, z) to different graphene stacks, we compare in Fig. 3 the miniband spectra computed using the parameters from Table I to ARPES intensities measured along a similar K 2 -K 1 k-space cut for, in Fig. 3 (a), a trilayer with θ = 9 • and, in Fig. 3 (b), twisted bilayer with θ = 19.1 • . Our model describes the bands of both of the structures well, despite changes in the twist angle, number of layers, potentials u and ∆ (which vary with growth conditions and thickness of the stack [29] and are determined for each structure individually) and the magnitudes of minigaps.

Probing electron wave function We assess the ac- Table I and shown with red dashed lines) for (a) twisted trilayer with θ = 9 • and (b) twisted bilayer with θ = 19.1 • , both measured along the direction connecting Brillouin zone corners K2 and K1 as shown in Fig. 1 (c) and indicated in the inset. In (a), the grey dashed lines, labelled (I)-(V), indicate energies for which constant-energy ARPES intensity maps are presented in Fig. (4) .

curacy of our parametrization of the interatomic potential, t(r, z), further by modelling directly the ARPES intensity data (we use approach developed in Ref. 37 and applied to the graphene/hexagonal boron nitride heterostructure in Ref. 38 ; see Methods and Supplementary Note 3 for further details). In graphene materials, interference of electrons emitted from different atomic sites within the unit cell provides additional information about the electronic wave function [37] . This is best visualized by ARPES intensity patterns at constant electron energy, which we present, both as obtained experimentally (top row) and simulated theoretically (bottom row), in Fig. 4 for the trilayer sample with θ = 9 • and energies indicated with grey dashed lines in Fig. 3 . For the map at the energy ε = 0, the two spots of high intensity indicate the positions of the valleys K 1 and K 2 . For energies 0 < ε < −0.6 eV, the bilayer and monolayer dispersions are effectively uncoupled. The crescent-like intensity pattern in the vicinity of K 1 reflects the pseudospin of n = 1 (evidence of Berry phase of π [39] ) of electrons in monolayer graphene. In contrast, in bilayer graphene, the low-energy band hosts massive chiral fermions [40] with pseudospin n = 2 so that the outer ring pattern in the vicinity of K 2 displays two intensity maxima, feature best visible in panel (II) . Because in our model all electron hoppings are generated naturally by t(r, z), agreement of our ARPES simulation with experimental data provides confirmation that our model and parametrization of the interatomic coupling t(r, z) leads to the correct band structure. Finally, panels (III)-(V) in Fig. 4 show the constant-energy maps in the vicinity of the minigaps which open due to hybridization of the bilayer and monolayer bands. The merging of 1L and 2L contours in panel (III) leads to a van Hove singularity and an associated peak in the electronic density of states, similarly to the case of twisted bilayer graphene [15] and discussed also for twisted trilayer graphene [29] (in the latter, the position of the van Hove singularity is established by tracking the minigap; the former is caused by saddle points in the electronic dispersion as the bands flatten at the anticrossings and so every minigap is accompanied by a van Hove singularity). Overall, our simulated patterns correctly reflect the evolution of the minigap as a function of energy and wave vector as well as the measured photocurrent intensity.

Our parametrization of t(r, z) is applicable to a wide range of twist angles, including the magic-angle regime [10, 34] as well as the 30 • -twisted bilayer graphene quasicrystal [41, 42] . To mention, it yields the k-space interlayer coupling at the graphene Brillouin zone corner K,t(|K|, c 0 ) = 0.11 eV. This agrees with the values used in effective models of the low-twist limit of twisted bilayer graphene [27, 34, 35, 43] which requiret(|K|, c 0 ) as the only parameter. Overall, our form of t(r, z) decays more rapidly in the real space (and hence slower in the reciprocal space) than usually assumed. This might explain the discrepancy between theory and experimental ARPES intensities of Dirac cone replicas observed for the case of 30 • -twisted bilayer graphene in Ref. 41 .

As we have shown, the same interatomic coupling t(r, z) can be used in graphene structures with different number of layers as, similarly to the case of perfect graphite and other layered materials, coupling to the nearest layer dominates the interlayer couplings. The continuum approach has been applied extensively to model the graphene/graphene interface, including to predict the existence of the magic angle [34] . Hence, in Supplementary Figure 1 , we use our results to simulate ARPES spectra for twist angles in the vicinity of the magic angle, θ ≈ 1.1 • , and show qualitative agreement with the recent experimental data [44, 45] . The continuum model was also used successfully to interpret experimental observations in graphene on hexagonal boron nitride [5] as well as homo-and heterobilayers of transition metal dichalcogenides [46, 47] . Our approach allows for experimental parametrization of the interatomic coupling t(r, z) for each of these interfaces as well as for others for which influence of neighbouring crystals can be approximated by considering the harmonics of the moiré potential [43, [48] [49] [50] [51] [52] . To comment, previous studies suggest that adapting our model to stacks of transition metal dichalcogenides requires taking into account changes in the interlayer distance as a function of the twist angle [20] . Moreover, in contrast to graphene, for which the part oft(q, z) most relevant to modelling twisted interfaces is that for q pointing to the Brillouin zone corner, q ≈ K, for transition metal dichalcogenides more significant changes due to interlayer coupling occur in the vicinity of the Γ point. In multilayers of 2H semiconducting dichalcogenides MX 2 (M = Mo, W, and X = S, Se), coupling of the degenerate states at the Γ point built of transition metal d z 2 and chalcogen p z orbitals leads to their hybridization and splitting which drives the direct-to-indirect band gap transition [53, 54] . Using the form of t(r, z) suggested in Ref. [26] for chalcogen p z -to-p z hopping (which dominates the interlayer coupling) in transition metal disulfides and diselenides, we computed the correspondingt(q, z) and obtained an estimate of t(Γ, c X−X ) ∼ 1.2 eV for interlayer nearest neighbour distance between chalcogen sites, c X−X ≈ 3Å. Taking into account the fractional contribution of the p z orbitals to the top valence band states at Γ in a monolayer [26] , we obtain coupling between two such states in bilayer ∼ 0.4 eV. This, in turn, suggests band splitting of ∼ 0.8 eV, in qualitative agreement with observations [53] [54] [55] . This supports the idea that our model can accurately describe and parametrize interatomic coupling between materials other than graphene.

Experimentally, our approach requires fabrication of trilayer (or thicker) stacks with one twisted and one perfect interface in order to benefit from the complementarity of the information obtained from self-consistent realand momentum-space description of the interfaces. However, to note, building on the observations of superconductivity in magic-angle twisted bilayer graphene [10, 11] , structures containing both a twisted and a perfect interface like twisted trilayer graphene [56, 57] , double bilayer graphene [58] [59] [60] [61] [62] [63] or double bilayer WSe 2 [64] recently attracted attention on its own due to observation of correlated electronic behaviour. Our approach provides one of the avenues to build an experimentally validated singleparticle base to study such effects. It could be, in principle, also applied to stacks of different materials, as long as one of the interfaces is commensurate and can be described in the real space in a tight-binding-like fashion. Finally, apart from continuum models, the interatomic coupling t(r, z) can also be used directly in large scale tight-binding calculations for commensurate twist angles [25, 26, [65] [66] [67] .

The ARPES measurements were performed at the Spectromicroscopy beamline at the Elettra synchrotron (Trieste, Italy). Before measurements, the samples were annealed at 350 • for 30 minutes. The experiment was then performed at a base pressure of 10 −10 mbar in ultrahigh vacuum and at the temperature of 110 K. We used photons with energy of 74 eV and estimate our energy and angular resolution as 50 meV and 0.5 • , respectively. For each sample, we determined the twist angle θ by measuring the distance between the Brillouin zone corners K 2 and K 1 which depends on the twist angle,

Further comments on experimental analysis of ARPES intensity are provided in Supplementary Note 4.

Theoretical calculations We write the Hamiltonian H in Equation (1) in the basis of sublattice Bloch states constructed of carbon p z orbitals φ(r 3D ) [31] ,

where k is electron wave vector, X = A, B is the sublattice, R l are the lattice vectors of layer l and τ X,l points to the site X in layer l within the unit cell selected by R l . We include in the basis all states coupled to k throughT (θ) which are less than a distance 28π 3 √ 3r AB sin θ 2 away from it, compute the matrix elements ofĤ in this truncated basis and diagonalize the resulting matrix numerically. In order to simulate the ARPES intensity, we project the eigenstates of the moiré Hamiltonian,Ĥ, on a plane-wave-like final state (see Supplementary Note 3 for more details and Ref. 38 for a detailed discussion of this approach for the case of graphene on hexagonal boron nitride). We determine the broadening of the ARPES signal as well as the decay constant for the intensity of Bernal bilayer signal by fitting to the experimental data. 

The authors declare no competing interests.

The data used in this study are available from the University of Bath data archive at https://doi.org/10.15125/BATH-00864 [68] .

As mentioned in the main text, we aim to construct our HamiltonianĤ based on blocks corresponding to intra-and interlayer couplings,

Above, the diagonal blocksĤ 0 (θ i , ε i ) describe consecutive graphene layers (starting from the top left for layer 1) at a twist angle θ i and with on-site energies of atomic sites in this layer, ε i . Without loss of generality, we choose θ 1 = 0 (we measure twists with respect to the crystallographic directions of layer 1) and, from our structure, θ 2 = 0 immediately follows.

In contrast, we set the reference point for energy as the on-site energies of layer 3 equal to zero, ε 3 = 0. From the definitions of the energies ∆ and u in the main text, we get

Consider now the sublattice Bloch states constructed of carbon p z orbitals φ(r 3D ) ≡ φ(r, z),

where k = (k x , k y ) is electron wave vector, X = A, B is the sublattice, R l are the lattice vectors of layer l, τ X,l points to the site X in layer l within the unit cell selected by R l and z l defines the position of layer l along the z-axis. In the basis {|k, A l , |k, B l }, the intralayer blocks take a familiar tight-binding form [1] ,

where a = √ 3|r AB | is the lattice constant of graphene. Similarly, the interlayer coupling block between layers 1 and 2 which are perfectly aligned and stacked according to the Bernal stacking [2] ,T

Let us now consider the interface between layers 2 and 3. The positions of the j-th and m-th sublattices in these layers can be expressed as R j = n 1 a 1 + n 2 a 2 + τ j,2 ,

where a 1 and a 2 are the primitive lattice vectors of layer 2 and n 1 , n 2 , n 1 , n 2 ∈ Z. For no twist, the two layers would be Bernal-stacked and that prescribes the relation between τ j,2 and τ m,3 . Here, we choose

Following earlier work [3, 4] , the interlayer coupling matrix element 3 k , m, |T |k, j 2 between sublattice Bloch states on layers 2 and 3 can be written as

The matrix element above can be understood in the following way: (1) By writing all four matrix elements in the form of a 2 × 2 matrix, we obtain Eq. (2) of the main text,T

where we only use τ ≡ τ A,2 = −(0, a/ √ 3).

In order to obtain the energy dispersion for a given k, we include into our Hamiltonian all states coupled to k throughT (θ) which are less than a distance 28π (8) can be grouped according to the same coupling magnitudes,t(k + G, c 0 ) ≡t(|k + G|, c 0 ) ≈t(|K + G|, c 0 ).

Because of the rapid decay oft(k, c 0 ) as a function of |k|, the sum converges. Taking note of the partial cancellation of the phase factors, in the limit θ → 0 equivalence of the descriptions in the real and reciprocal spaces requires that

to reproduce the interlayer coupling in Bernal bilayers. As γ 1 ≡ t(0, c 0 ), this condition provides us with additional constraints on the behaviour of t(r, z) and its Fourier transform t(k, z).

whereφ p e / , p ⊥ e / = dr dz e − i p e ·r e − i p ⊥ e z φ(r, z),

is the Fourier transform of the p z orbital φ(r, z).

Due to the rotational symmetry of the p z orbital,φ p e / , p ⊥ e / =φ |p e / |, p ⊥ e / . Moreover, for the given photon energy, ω, and work function, W , we have p ⊥ e |p e | so that φ R θ l (k + g) − G l , p ⊥ e / can be approximated by a constant and dropped. Finally, in this work we only study points for which G l = 0. As a result,

where we used the fact that z l = lc 0 . We combine both phases exp(iϕ X,l ) and exp(− i p ⊥ e lc 0 ) into a single factor,

which we fit to experiment. To note, the experimental data in this work suggests that the phase difference between the aligned graphene layers is approximately e iπ . Finally, we model the Dirac delta in Supplementary Equation (13) with a Lorentzian

with half-width-half-maximum γ (the value of which, γ ≈ 0.17 eV, we obtain by comparison to the experimental data).

In most cases, our graphene samples are electron-doped so that the Dirac and neutrality this distance is directly governed by the twist angle).

As part of our data analysis, we extract the differences between on-site potentials on different graphene layers. For Bernal bilayer, we investigate the intensity at the corner of its Brillouin zone as a function of energy as the gap in the electronic dispersion is equal to 2|u| [2] . The neutrality point is located halfway in the gap. In turn, the energy shift between the neutrality point and the position of the Dirac point of layer 3 defines ∆.

Finally, in our theoretical analysis, we assume such an orientation of the graphene layers that one of the Brillouin zone corners of the first layer is located at K = 4π 3a (1, 0). In comparison, depending on the experimental geometry, the Brillouin zones of the graphene layers can as introduced in the main text and to be compared with experimental data shown in Ref. [11] and [12] , respectively. Above, we used our model with the parametrization of t(r, z) as presented in the main text. We choose the spectral broadening as well as wave vector and energy ranges to match the experimental data (unfortunately, neither of the experimental works zoom in on the flat band features). The results are in qualitative agreement with the presented experimental photoemission intensity. We do not make comparisons to the constant energy maps presented in [11] as it is unclear to us exactly what the symmetrization procedure is that is used to generate these maps.

Electric field effect in atomically thin carbon films

Van der Waals heterostructures

2D materials and van der Waals heterostructures

Van der Waals heterostructures and devices

Cloning of dirac fermions in graphene superlattices

Hofstadter's butterfly and the fractal quantum hall effect in moiré superlattices

Emergence of interfacial polarons from electron-phonon coupling in graphene/h-BN van der Waals heterostructures

Strong interlayer coupling in van der waals heterostructures built from single-layer chalcogenides

Interlayer valley excitons in heterobilayers of transition metal dichalcogenides

Correlated insulator behaviour at halffilling in magic-angle graphene superlattices

Unconventional superconductivity in magic-angle graphene superlattices

Resistivity of Rotated Graphite-Graphene Contacts

Twistable electronics with dynamically rotatable heterostructures

Tunable crystal symmetry in graphene-boron nitride heterostructures with coexisting moiré superlattices

Evidence for Interlayer Coupling and Moiré Periodic Potentials in Twisted Bilayer Graphene

Direct Observation of Interlayer Hybridization and Dirac Relativistic Carriers in

Graphene/MoS2 van der Waals Heterostructures

Determination of band offsets, hybridization, and exciton binding in 2D semiconductor heterostructures

Band Alignment and Minigaps in Monolayer MoS2-Graphene van der Waals Heterostructures

Direct observation of minibands in a twisted graphene/WS2 bilayer

Direct Measurement of the Tunable Electronic Structure of Bilayer MoS2 by Interlayer Twist

Strong interlayer hybridization in the aligned SnS2/WSe2 hetero-bilayer structure

Tailoring the Electronic Structure in Bilayer Molybdenum Disulfide via Interlayer Twist

Probing the Interlayer Coupling of Twisted Bilayer MoS2 Using Photoluminescence Spectroscopy

Evolution of interlayer coupling in twisted molybdenum disulfide bilayers

Localization of Dirac Electrons in Rotated Graphene Bilayers

Ab initio tight-binding Hamiltonian for transition metal dichalcogenides

Graphene Bilayer with a Twist: Electronic Structure

Generic miniband structure of graphene on a hexagonal substrate

Substrate Doping Effect and Unusually Large Angle van Hove Singularity Evolution in Twisted Bi-and Multilayer Graphene

A review of chemical vapour deposition of graphene on copper

The electronic properties of graphene

Gate-induced interlayer asymmetry in ABA-stacked trilayer graphene

Landau-level degeneracy and quantum hall effect in a graphite bilayer

Moiré bands in twisted double-layer graphene

Interlayer interaction in general incommensurate atomic layers

Simplified LCAO Method for the Periodic Potential Problem

Characterization of graphene through anisotropy of constant-energy maps in angle-resolved photoemission

Moiré miniband features in the angle-resolved photoemission spectra of graphene/hBN heterostructures

Experimental observation of the quantum Hall effect and Berry's phase in graphene

Unconventional quantum Hall effect and Berry's phase of 2π in bilayer graphene

Dirac electrons in a dodecagonal graphene quasicrystal

Quasicrystalline 30 • twisted bilayer graphene as an incommensurate superlattice with strong interlayer coupling

Ab initio theory of moiré superlattice bands in layered twodimensional materials

Visualization of the flat electronic band in twisted bilayer graphene near the magic angle twist

Direct evidence for flat bands in twisted bilayer graphene from nano-ARPES

Resonantly hybridized excitons in moiré superlattices in van der Waals heterostructures

Evidence for moiré excitons in van der Waals heterostructures

Moiré minibands in graphene heterostructures with almost commensurate √ 3× √ 3 hexagonal crystals

Anomalous Light Cones and Valley Optical Selection Rules of In-terlayerExcitons in Twisted Heterobilayers

Topological mosaics in moiré superlattices of van der Waals heterobilayers

Interlayer hybridization and moiré superlattice minibands for electrons and excitons in heterobilayers of transition-metal dichalcogenides

Topological Insulators in Twisted Transition Metal Dichalcogenide Homobilayers

Direct Measurement of the Thickness-Dependent Electronic Band Structure of MoS2 Using Angle-Resolved Photoemission Spectroscopy

Direct observation of the transition from indirect to direct bandgap in atomically thin epitaxial MoSe2

Electronic properties of single-layer and multilayer transition metal dichalcogenides MX2

Electrically tunable correlated and topological states in twisted monolayer-bilayer graphene

Tunable van Hove Singularities and Correlated States in Twisted Trilayer Graphene

Spin-polarized Correlated Insulator and Superconductor in Twisted Double Bilayer Graphene

Correlated Insulating States in Twisted Double Bilayer Graphene

Correlated states in twisted double bilayer graphene

Tunable correlated states and spinpolarized phases in twisted bilayer-bilayer graphene

Tunable correlation-driven symmetry breaking in twisted double bilayer graphene

Density-Wave States in Twisted Double-Bilayer Graphene

Interaction effects and superconductivity signatures in twisted double-bilayer WSe2

Lattice relaxation and energy band modulation in twisted bilayer graphene

Minimum model for the electronic structure of twisted bilayer graphene and related structures

Formation of Bloch Flat Bands in Polar Twisted Bilayers without Magic Angles

Dataset for article "Determination of interatomic coupling between two-dimensional crystals using angle-resolved photoemission spectroscopy

Supplementary References

The electronic properties of graphene

Landau-Level Degeneracy and Quantum Hall Effect in a Graphite Bilayer

Moiré bands in twisted double-layer graphene. Proceedings of the National Academy of

Interlayer interaction in general incommensurate atomic layers

Characterization of graphene through anisotropy of constantenergy maps in angle-resolved photoemission

Moiré miniband features in the angle-resolved photoemission spectra of graphene/hBN heterostructures

Graphene Sublattice Symmetry and Isospin Determined by Circular Dichroism in Angle-Resolved Photoemission Spectroscopy

Visualizing Electronic Chirality and Berry Phases in Graphene Systems Using Photoemission with Circularly Polarized Light

Illuminating the dark corridor in graphene: Polarization dependence of angle-resolved photoemission spectroscopy on graphene

Direct measurement of quantum phases in graphene via photoemission spectroscopy

Direct evidence for flat bands in twisted bilayer graphene from nano-ARPES

Visualization of the flat electronic band in twisted bilayer graphene near the magic angle twist

Determination of interatomic coupling in twisted graphene using angle-resolved photoemission spectroscopy Thompson et al.

Using Fermi's golden rule, we write ARPES intensity as [5, 6] Iwhere M f,i is the matrix element describing transition of the electron from the initial state in the crystal in band i to the final state f , ω is the energy of the incident photon, ε i,k is the energy of an electron in the crystal in band i and with wave vector k, ε p e is the energy of the photoelectron with momentum p e and W is the work function of graphene. Within the dipole approximation,where is the vector potential of the incident photon,p is the momentum operator, |final stands for the final state of the photoelectron and |k, i denotes the wave function of the electron in the crystal. The latter is a linear combination of the sublattice Bloch states, Supplementary Equation (2), corresponding to wave vectors connected by a superlattice reciprocal vector g = G − G,with the coefficients c g,i X,l (k) provided by diagonalization of the HamiltonianĤ, Supplementary Equation (1) . Here, g is the moiré reciprocal superlattice vector. We approximate the final state with a plane wave (justified for incident photon energies above 50 eV [7] ) with momentum p e = (p e , p ⊥ e ), so thatFollowing from the chiral properties of the graphene wave function, the light-matter interaction ·p leads to an angle-dependent phase difference between the distinct atomic orbitals in the matrix element, e iϕ X,l [8] [9] [10] (note that this includes the effect of changing the sign of γ 1 [10] ). Hence,e iϕ X,l c g,i X,l (k) e i p e ·r e i p ⊥ e z |R θ l (k + g), X l = g l,X,G l e iϕ X,l c g,i X,l (R −θ l (p e / + G l ) − g)e iG l ·τ X,l e − i p ⊥ e z lφ R θ l (k + g) − G l , p ⊥ e / ,