Interlayer interaction in general incommensurate atomic layers

We present a general theoretical formulation to describe the interlayer interaction in incommensurate bilayer systems with arbitrary crystal structures. By starting from the tight- binding model with the distance-dependent transfer integral, we show that the interlayer coupling, which is highly complex in the real space, can be simply written in terms of generalized Umklapp process in the reciprocal space. The formulation is useful to describe the interaction in the two-dimensional interface of different materials with arbitrary lattice structures and relative orientations. We apply the method to the incommensurate bilayer graphene with a large rotation angle, which cannot be treated as a long-range moir\'{e} superlattice, and obtain the quasi band structure and density of states within the first-order approximation.


INTRODUCTION
Recently there has been extensive research efforts on atomically-thin nanomaterials, including graphene, hexagonal boron nitride(hBN), phosphorene and the metal transition dichalcogenides. The hybrid systems composed of different kind of atomic layers have also attracted much attention, and there the interaction at the interface between the different atomic layers often plays an important role in the physical property. In such a composite system, the lattice periods of the individual atomic layers are generally incommensurate due to the difference in the crystal structure and also due to misorientation between the adjacent layers. A well-known example of irregularly stacked multilayer system is the twisted bilayer graphene, in which two graphene layers are rotationally stacked at an arbitrary angle [1]. When the rotation angle is small, in particular, the system exhibits a moiré interference pattern of which period can be much greater than the atomic scale, and such a long-period modulation is known to strongly influence the low-energy electronic motion [2,3,4,5,6,7,8,9,10,11,12,13,14]. Graphene-hBN composite system has also been intensively studied as another example of incommensurate multilayer system, where the two layers share the same hexagonal lattice structure but with slightly-different lattice constants, leading to the long-period modulation even at zero rotation angle [15,16,17,18,19,20,21]. The electronic structure in graphene-hBN system was theoretically studied [22,23,24,25,26,27,28,29,30,31,32,33,34,35], and the recent transport measurements revealed remarkable effects such as the formation of mini-Dirac bands and the fractal subband structure in magnetic fields [36,18,19,20].
The previous theoretical works mainly targeted the honeycomb lattice to model twisted bilayer graphene and graphene-hBN system, and also particularly focus on the long-period moiré modulation which arises when the crystal structures of two layers are slightly different. Then we may ask how to treat general bilayer systems where the lattice vectors of the adjacent layers are not close to each other. In this paper, we develop a theoretical formulation to describe the interlayer interaction effect in general bilayer systems with arbitrary choice of crystal structures and relative orientations. By starting from the tight-binding model with distance-dependent transfer integral, we show that the interlayer coupling can be simply written in terms of a generalized Umklapp process in the reciprocal space. We then apply the method to the incommensurate bilayer graphene with a large rotation angle (θ = 20 • ) which cannot be treated as a long-range moiré superlattice, and obtain the quasi band structure and density of states within the first-order approximation. Finally, we apply the formulation to the moiré superlattice where the two lattice structures are close, and demonstrate that the long-range effective theory is naturally derived.

Interlayer Hamiltonian for general incommensurate atomic layers
We consider a bilayer system composed of a pair of two-dimensional atomic layers having generally different crystal structures as shown in Fig. 1. We write the primitive lattice vectors as a 1 and a 2 for layer 1 andã 1 andã 2 for layer 2, which are all along in-plane (x-y) direction. The reciprocal lattice vectors are defined by G i andG i for layer 1 and 2, respectively, so as to satisfy a i · G j =ã i ·G j = 2πδ ij . The area of the unit cell is given by S = |a 1 × a 2 | and S = |ã 1 ×ã 2 | for layer 1 and 2, respectively. Without specifying any details of the model, we can easily show that an electronic state with a Bloch wave vector k on layer 1 and one withk on layer 2 are coupled only when where G = m 1 G 1 + m 2 G 2 andG =m 1G1 +m 2G2 are reciprocal lattice vectors of layer 1 and 2, respectively. This is regarded as a generalized Umklapp process between arbitrary misoriented crystals, and it can be easily understood by considering the wave decomposition in the plain wave basis as follows. A Bloch state on layer 1 (say φ (1) k ) is expressed as a summation of e i(k+G) over the reciprocal vectors G, and one on layer 2 (φ (2) k ) is expressed as a summation of e i(k+G) overG. Also, Hamiltonian of the total system consists of Fourier components of G andG. As a result, the matrix element φ In the following, we actually calculate the matrix elements for generalized Umklapp processes using the tight-binding model. We assume that a unit cell in each layer contains several atomic orbitals, which are specified by X = A, B, · · · for layer 1 andX =Ã,B, · · · for layer 2. The lattice positions are given by R X = n 1 a 1 + n 2 a 2 + τ X (layer 1), where n i andñ i are integers, and τ X and τX are the sublattice position inside the unit cell, which can have in-plane and out-of-plane components. When the atomic layers are completely planar and stacked with interlayer spacing d, for instance, we have τ X · e z = 0 for layer 1 and τX · e z = d for layer 2, where e z is the unit vector in z-direction. Let us define |R X ≡ φ X (r − R X ) as the atomic state of the sublattice X localized at R X . The atomic orbital φ X may be different depending on X. We assume the transfer integral from the site R X to RX is written as −TX X (RX − R X ), i.e., depending on the relative position RX − R X and also on the sort of atomic orbitals of X andX. The interlayer Hamiltonian to couple the layer 1 and 2 is then written as When the superlattice period is huge, constructing Hamiltonian in the real space bases becomes hard because it requires the relative inter-atom position for every single combination of the atomic sites of layer 1 and layer 2. On the other hand, the interlayer coupling is described in a simpler manner in the reciprocal space. We define the Bloch basis of each layer as where k andk are the two-dimensional Bloch wave vectors parallel to the layer, N (Ñ ) is the number of unit cell of layer 1(2) in the total system area S tot . Although the layer 1 and layer 2 are generally incommensurate, we assume the system has a large but finite area S tot = N S =ÑS to normalize the wave function. Then we can show that the matrix elements of U between Bloch bases can be written as, which is non-zero only when the generalized Umklapp condition Eq. (1) is satisfied. Here t(q) is the in-plane Fourier transform of the transfer integral defined by where zX X = (τX −τ X )·e z , and the integral in r is taken over the two-dimensional space of S tot . G andG run over all the reciprocal lattice vectors of layer 1 and 2, respectively. Since tX X (k) decays in large k, we only have a limited number of relevant components in the summation of Eq. (5). Eq. (5) is derived in a straightforward manner as follows. By inserting U in Eq. (3) to the definition of UX X , we have By applying the inverse Fourier transform of Eq. (6), to Eq. (7), the second summation is transformed as where we used in the last equation, Using Eqs. (7) and (9), we have which is Eq. (5). In the summation in R X in the first line, we used the transformation

Irregularly stacked honeycomb lattices
We apply the general formulation obtained above to the irregularly stacked bilayer graphene with an arbitrary rotation angle. Here we consider a pair of hexagonal lattices as shown in Fig.  2(a), which are stacked with a relative rotation angle θ and interlayer spacing d. We define the primitive lattice vectors of layer 1 as a 1 = a(1, 0) and a 2 = a(1/2, √ 3/2) with the lattice constant a, and those of the layer 2 asã i = R a i where R is the rotation matrix by θ. Accordingly, the reciprocal lattice vectors of layer 1 (G 1 , G 2 ) and layer 2 are related byG i = R G i . The atomic positions are given by for X = A, B (layer 1) andX =Ã,B (layer 2), where Here we take the origin at an A site, and define τ 0 as the relative in-plane translation vector of the layer 2 to layer 1.
To describe the electron's motion, we adopt the single-orbital tight-binding model for p z atomic orbitals. Then TX X (R) does not depend on indexesX and X, and it is approximately written in terms of the Slater-Koster parametrization as, [37] The typical parameters for graphene are a ≈ 0.246nm, d ≈ 0.335 nm, V 0 ppπ ≈ −2.7 eV, V 0 ppσ ≈ 0.48 eV, and r 0 ≈ 0.184a [11]. Once the transfer integral T (R) between the atomic sites is given, one can calculate the in-plane Fourier transform t(q), and then obtain the interlayer Hamiltonian by Eq. (5). Since the transfer integral is isotropic along the in-plane direction, we can write t(q) = t(q) with q = |q|. Fig. 2(a) illustrates the lattice structure at rotation angle θ = 20 • . Fig. 2(b) shows the Brillouin zones of the two layers in the extended zone scheme, where the blue symbols (filled circles, triangles and open circles) represent the positions of k + G for some particular k (here chosen as the zone corner K + ) with several G's. Figs. 2(c) plots the corresponding positions of these symbols inside the first Brillouin zone of layer 2. This indicates the wave numbers of layer 2,k = k+G−G, which are coupled to k of layer 1 under the condition of Eq. (1). The amplitude of the coupling is given by t(k + G), and it solely depends on the distance to each symbol from the k-space origin in Fig. 2(b). In the present parameter choice, the amplitudes for filled circles, triangles and open circles are t(K) ≈ 110 meV, t(2K) ≈ 1.6 meV, and t( √ 7K) ≈ 0.062 meV respectively, where K = |K + | = 4π/(3a). The couplings for other k-points are exponentially small and negligible.
When the lattice vectors of the two layers are incommensurate (i.e., do not have a common multiple) as in this example, we cannot define the common Brillouin zone nor calculate the exact band structure, since the interlayer interaction connects infinite number of k-points in the Brillouin zones of layer 1 and layer 2. But still, we can obtain an approximate band structure, considering only the first-order interlayer processes while neglecting multiple processes. Let us consider a particular wave vector k of layer 1, and take allk's in layer 2 which are directly coupled to k as in Fig. 2(c). By neglecting exponentially small matrix elements, we can construct a finite Hamiltonian matrix including only the bases of a single wave vector k of layer 1 and severalk's in layer 2. By diagonalizing the matrix, we obtain the energy eigenvalues ε nk labeled by the index n. We then define the spectral function contributed from layer 1 as where g (1) nk is the total wave amplitudes on layer 1 in the state |nk . The density of states contributed from layer 1 is expressed as where the summation is taken over the first Brillouin zone of layer 1. We perform the exactly same procedure for the layer 2 by considering a singlek of layer 2 and all k's in layer 1 coupled tok, and obtain the spectral function and the density of states for the layer 2. The total density of states of the system is given by D 1 + D 2 . Fig. 3 plots the total density of state D(ε) calculated for the incommensurate bilayer graphene with θ = 20 • . The red dashed curve is the density of states of decoupled bilayer graphene (i.e., twice of monolayer's). We actually see additional peak structures due to the interlayer interaction, and these features are consistent with the exact tight-binding calculation for a commensurate bilayer graphene at a similar rotation angle [11]. Fig. 4 illustrates the Fermi surface reconstruction at several different energies. The right panel in each row shows the layer 1's spectral function A 1 (k, ε) in presence of the interlayer coupling. The left panel shows the Fermi surface in absence of the interlayer coupling, where the black curves represent the equienergy lines for the layer 1's energy dispersion ε 1 (k), and the pink curves are for the layer 2's dispersion with k-space shift, ε 2 (k + G). Thickness of the pink curves indicate the absolute value of the interlayer coupling t(k + G).  Fig. 4(a) shows a typical case (ε = 2.5 eV) where we observe small band anticrossing at the intersection of the Fermi surfaces of layer 1 and layer 2. Fig. 4(b) (ε = −2.85 eV) is for the energy at which the density of states exhibits a dip [ Fig. 3] . There the interlayer coupling is relatively strong, and indeed we see that the original triangular Fermi pockets of the individual layers are strongly mixed and reconstructed into a large single Fermi surface. Note that the coupling strength depends not only on t(k + G), but also the relative phase factors between different sublattices. At even lower energy ε = −5 eV [ Fig. 4(b)] beyond the van-Hove singularity, the layer 1 and the layer 2 have almost identical Fermi surface surrounding Γ point in absence of the interlayer coupling, and they are hybridized into a pair of circles with different radii corresponding to the bonding and anti-bonding states. This feature is also observed in the density of states [ Fig. 3] as a large split of the band bottom. Such a splitting does not occur in the positive energy region, because there A and B sublattices in the same layer have the opposite signs in the wave amplitude, so that the interlayer mixing vanishes due to the phase cancellation.
In the above approximation, we neglect the second order processes such thatk points of layer 2 (linked from initial k of layer 1) are further coupled to other k ′ of layer 1. We can include such higher-order processes up to any desired order as follows. To include the second order process in calculating A 1 (k, ε), for example, we take allk of layer 2 which are coupled to k of layer 1 and also take all k ′ of layer 1 which are coupled tok, and then construct a Hamiltonian matrix with a larger dimension. From the obtained eigenstates, we calculate the spectral function using Eq. (16), but then g (1) nk should be the total wave amplitudes from layer 1 at the original k (0-th order), without including k ′ .

Long-period moiré superlattice
When the primitive lattice vectors of layer 1 and those of layer 2 are slightly different, the interference of two lattice structures give rise to a long-period moiré pattern, and then we can use the long-range effective theory to describe the interlayer interaction. [2,3,4,5,6,7,8,9,10,11,12,13,14] Here we show that the long-range effective theory can be naturally derived from the present general formulation, just by assuming that the two lattice structures are close  1 (k, ε) in presence of the interlayer coupling. The left panel shows the Fermi surface in absence of the interlayer coupling, where the black curve is for the layer 1's energy dispersion ε 1 (k), and the pink curve is for the layer 2's dispersion with k-space shift, ε 2 (k + G). Thickness of the pink curve indicates the absolute value of the interlayer hopping t(k + G). to each other. We define a linear transformation with a matrix A that relates the primitive lattice vectors of layer 1 and 2 asã Correspondingly, the reciprocal lattice vectors becomẽ to satisfy a i · G j =ã i ·G j = 2πδ ij . When the lattice structures of the two layers are similar, the matrix A is close to the unit matrix. The reciprocal lattice vectors of the moiré superlattice is given by small difference between G i andG i as When the two layers are identical and rotationally stacked with a small angle, for example, the matrix A is given by a rotation R. Noting (R † ) −1 = R, we have When layer 1 and layer 2 have different lattice constant as in the graphene-hBN bilayer, the matrix A is given by the combination of the isotropic expansion and the rotation as αR, where α is the lattice constant ratio. Since [(αR) † ] −1 = α −1 R, we have When the specific form of the matrix A is given, we can immediately derive the interlayer matrix elements for the long wavelength components using Eq. (5) as, where m 1 and m 2 are integers.
For example, let us derive the interlayer Hamiltonian of the irregularly stacked bilayer graphene with a small rotation angle. Since the low-energy spectrum of graphene is dominated by the electronic states around the Brillouin zone corners K and K ′ , we consider the matrix elements for initial and final k-vectors near those points. The K and K ′ points are located at K ξ = −ξ(2G 1 + G 2 )/3 for layer 1 andK ξ = −ξ(2G 1 +G 2 )/3 for layer 2, where ξ = ±1 are the valley indexes for K and K ′ , respectively. When we start from k = K + , for example, a typical scattering process is illustrated in Figs. 5(b) and (c) in a similar manner to Fig. 2. There the electron at K + in layer 1 is coupled to K + + m 1 G 1 + m 2 G 2 in layer 2 with absolute amplitude t(K + +m 1 G 1 +m 2 +G 2 ). As already argued in the previous section, the coupling amplitudes for filled circles, triangles and open circles in Figs. 5(b) and (c) are t(K) ≈ 110 meV, t(2K) ≈ 1.6 meV, and t( √ 7K) ≈ 0.062 meV respectively, and the couplings to other k-points are negligibly small. The matrix element changes when the initial vector k is shifted from K + , but we neglect such a dependence assuming k is close to K + . As a result, we obtain the interlayer Hamiltonian of near K ξ from Eq. (23) as where r is the in-plane position, ω = exp(2πi/3) and τ 0 is set to 0. The total Hamiltonian is written in the basis of {A, B,Ã,B} as where H 1 and H 2 are the intralayer Hamiltonian of layer 1 and 2, respectively, defined by with Pauli matrices σ x and σ y , and the graphene's band velocity v. If we only take the terms with t(K) in the matrix U , the expression becomes consistent with the previous formulation for graphene-graphene bilayer. [2,7,8,11] Note that the Hamiltonian matrix depends on the actual choice of K points out of the equivalent Brillouin zone corners; where r 0 = (A − 1)τ 0 . This reflects the fact that relative sliding between two layers leads to a shift of the moiré interference pattern in the real space. The only exception is when the two layers shares the same lattice vectors G i =G i , where G M i vanishes so that τ 0 cannot be eliminated by shifting the origin. In this case, the energy band actually becomes different depending on the sliding vector τ 0 . For the graphene bilayer case, in particular, the expression of U with G i =G i becomes equivalent to that of regularly-stacked graphene bilayer with a interlayer sliding. [38,39,40]

Conclusion
We theoretically studied the interlayer interaction in general incommensurate bilayer systems with arbitrary crystal structures. Using the generic tight-binding formulation, we demonstrate that the interlayer coupling in the reciprocal space is simply expressed in terms of a generalized Umklapp process. We applied the formula to the incommensurate honeycomb lattice bilayer with a large rotation angle, which cannot be treated as a long-range moiré superlattice, and actually obtain the quasi band structure and density of states within the first-order approximation. Finally, we apply the formulation to the moiré superlattice where the two lattice structures are close, and derive the long-range effective theory with a straightforward calculation.