# General theoretical description of ARPES of vdW structures

###### Abstract

We develop a general theory to model the angle resolved photoemission spectroscopy (ARPES) of commensurate and incommensurate van der Waals (vdW) structures, formed by lattice mismatched and/or misaligned stacked layers of two-dimensional materials. The present theory is based on a tight-binding description of the structure and the concept of generalized umklapp processes, going beyond previous descriptions of ARPES in incommensurate vdW structures, which are based on continuous, low energy models, being limited to structures with small lattice mismatch/misalignment. As an example, we study the ARPES bands and constant energy maps for a twisted bilayer graphene structure. The present theory should be useful in correctly interpreting experimental results of ARPES of vdW structures and other system displaying competition between different periodicities, such as density wave phases.

## I Introduction

The development of fabrications techniques in recent years enabled the creation of structures formed by stacked layers of different two-dimensional (2D) materials, referred to as van der Waals (vdW) structures Ponomarenko et al. (2011); Novoselov and Neto (2012); Geim and Grigorieva (2013). By combining layers of materials displaying different properties, it is possible to engineer devices with new functionalities, not displayed by the individual layers. This makes vdW structures very appealing from the applications point of view. As an example, transistors based on graphene and hexagonal boron nitride (h-BN) or semiconducting transition metal dichalcogenides (STMD) Britnell et al. (2012); Georgiou et al. (2012) and photodetectors based on graphene and STMD Britnell et al. (2013); Yu et al. (2013); Massicotte et al. (2015) have already been realized. The properties of a vdW structure depend not only on the properties of the individual layers, but also on how different layers interact with each other. Due to the high crystallographic quality of 2D materials, the interlayer interaction depends crucially on the lattice mismatch and misalignment between different layers. This is clearly exemplified, both experimentally Mishchenko et al. (2014) and theoretically Brey (2014); Amorim et al. (2016a), by the observation of negative differential conductance in vertical graphene/h-BN/graphene tunneling transistors, where the bias voltages at which peak current occurs is controlled by the angle between the misaligned graphene electrodes. A necessary step to fully understand and take advantage of vdW structures is to characterize their electronic properties, and how these are dependent on lattice mismatch/misalignment.

Angle resolved photoemission spectroscopy (ARPES) is an extensively
used tool to characterize the electronic degrees of freedom of materials
Damascelli et al. (2003); Kordyuk (2014); Moser (2017). In crystals, ARPES
is generally understood as a direct probe of the electronic band structure
of occupied states over the Brillouin zone. Nevertheless, even in
a perfect crystal, where the notions of reciprocal space and Brillouin
zone are well defined, this picture might breakdown, as the ARPES
response is weighted by matrix elements which describe the light induced
electronic transition between a crystal bound state to a photoemitted
electron state. For bands that are well decoupled from the remaining
band structure, the ARPES matrix elements are featureless, and indeed
ARPES can be seen as a direct probe of the band structure. However,
exceptions to this can occur and the matrix elements can impose selections
rules on the transitions. Two well known examples where this occurs
are graphite Shirley et al. (1995) and graphene Bostwick et al. (2007a, b); Mucha-Kruczyński et al. (2008); Gierz et al. (2010); Jung et al. (2010); Puschnig and Lüftner (2015); Lee et al. (2017).
In these materials part of the Fermi surface is not observed in constant
energy ARPES maps.^{1}^{1}1As shown in Ref. Gierz et al. (2010) it is actually possible to observe
the full Fermi surface of graphene in ARPES by using polarized
incident light. This phenomena can be explained theoretically if scattering
of the photoemitted state by the graphene lattice is taken into account.
Such effects will not be considered in the present paper. This effect is due to the ARPES matrix elements, which suppress the
signal from some parts of the band structure.

Another case where the ARPES matrix elements should play an important role is in systems with competing periodicities, such as materials with charge density wave (CDW) phases Voit et al. (2000); Che (1995). In this case it is easy to understand how the interpretation of ARPES as a direct probe of the band structure can break down. Let us suppose that for a given material in the normal, undistorted phase, ARPES accurately maps the electronic band structure. Suppose that the system undergoes a transition into a commensurate CDW, with a larger unit cell. If the distortion is small, the bands in Brillouin zone will be weakly perturbed apart from back folding into the new, smaller, Brillouin zone due to the enlarged unit cell. If the CDW perturbation is weak, by an adiabatic argument, the ARPES mapped bands observed in both phases must be essentially unchanged. This means that the signal of the back folded bands must be very weak, and the observed ARPES bands will mostly follow the bands of the undistorted phase in an extended zone scheme. The suppression of the back folded bands is due to the ARPES matrix elements. This was exemplified in Ref. Voit et al. (2000) in a simple one-dimensional tight-binding model. These effects may also be relevant when interpreting ARPES experiments in cuprates, for which a hidden density wave state has been proposed Chakravarty et al. (2003).

In vdW structures, such competing periodicities naturally occur due to the lattice mismatch between different layers. Therefore, a general theory capable of correctly taking into account ARPES matrix elements is essential to interpret ARPES data from vdW structures. We point out that the modelling of ARPES in twisted bilayer graphene Pal and Mele (2013) and graphene/hBN Mucha-Kruczyński et al. (2008) structures has been considered previously in the literature. However the models employed relied on effective low energy, continuous descriptions of the system, which are only valid for small misalignment angles. As the field of vdW structures develops, a more general and flexible approach is required. The goal of this work is to develop a general framework to theoretically model ARPES, which is valid for both commensurate and incommensurate structures formed by arbitrary materials and with arbitrary lattice mismatch/misalignment.

The structure of this paper is as follows. In Section II, we review the description of vdW structures based on tight-binding models and generalized umklapp processes developed in Refs. Bistritzer and MacDonald (2010); Koshino (2015). We use this description of vdW structures to compute the ARPES matrix elements for an arbitrary structure in Section III. We apply the general framework to model ARPES in twisted bilayer graphene in Section IV. Conclusions are drawn in Section V. For completeness and the convenience of the reader, in the Appendix A we briefly review the derivation of the ARPES intensity within the non-equilibrium Green’s function approach.

## Ii Tight-binding description of vdW structures

A first step in modelling ARPES of vdW structures is to describe the electronic states bound to the structure. Following Refs. Bistritzer and MacDonald (2010); Koshino (2015) we employ a tight-binding model to describe the bound states. We will focus on bilayer structures, where each layer as a periodic structure with Bravais lattice sites given by , with for the top and bottom layers, respectively. The single particle Hamiltonian of the structure is given by

(1) |

where and are the tight-binding Hamiltonians of the isolated top and bottom layers and () describes the hopping of electrons from the bottom (top) to the top (bottom) layer. More concretely, the intralayer terms are written as

(2) |

and the interlayer terms as

(3) |

with . In the previous equations, the indices run over lattice sites and the indices run over sublattice, orbital and spin degrees of freedom. creates an electron in state , a localized Wannier state in layer , lattice site and sublattice site . The Wannier wavefunction in real space reads where is the Wannier wavefunction or type centered on the origin. and are intralayer hopping terms, which we assume to be invariant under translations by lattice vector of the respective layer, while and are interlayer hopping terms, which describe the coupling between the two layers. It is convenient to express the electronic operators in term of Fourier components

(4) |

where belongs to the unit cell of the reciprocal lattice of layer and is the number of unit cells in layer . Notice that if is a reciprocal lattice vector of layer , i.e. , then . These states bring the Hamiltonians of the isolated layers to a block diagonal form,

(5) |

where . For the interlayer term, we assume a two-center approximation for the hopping elements and write them in terms of its Fourier transform components Koshino (2015) as (focusing on the term)

(6) |

where is the area of the unit cell of layer . With this we can write as

(7) |

where are reciprocal lattice vectors of layer . is written in a similar way. Equation (7) tells us that states of the two layers with crystal momentum and are only coupled provided and exist, such that the generalized umklapp condition is satisfied Koshino (2015).

In an extended Brillouin zone scheme, the generalized umklapp condition can be satisfied if for each and we write and for any . This fact motivates us to look for eigenstates of the complete Hamiltonian Eq. (1) of the form

(8) |

with and coefficients that are to be determined. We now introduce a convenient compact notation. We define as a vector with entries given by ; in the same way we introduce the matrices with entries , with entries and defined in a similarly way. This allows us to write the eigenvalue problem which determines the eigenstates and energies of the vdW structure as

(9) |

where are the eigenvalues,

(10) | ||||

(11) |

and the Hamiltonian matrix is written as

(12) |

where is a block diagonal matrix

(13) |

with similarly defined, and is a dense matrix

(14) |

with . For a commensurate structure, there exists and such that and the sums in Eq. (8) become finite, and consequently the matrix is finite. In this case, the energies and eigenstates of for restricted to the Brillouin zone of the commensurate structure provide us the full spectrum and eigenstates of the bilayer structure. For an incommensurate structure, the sums in Eq. (8) involve an infinite number of terms and the corresponding Hamiltonian matrix is infinite. However, an approximation of the eigenstates and energies of the incommensurate structure can still be obtained by suitably truncating , by considering a finite number of reciprocal lattice vectors and . Even in the case of a commensurate structure, the supercell might be very large, , leading to a very large matrix , and in that situation it might still be beneficial to compute the eigenstates and energies of the system approximately by truncating . In the following, we will generally refer to the approximate eigenvalues as band structure, even in the case where we have an incommensurate structure and the concept of Brillouin zone no longer applies.

We will now prove a formal relation satisfied by the solutions of Eq. (9) that will be useful in the next section: Assume that

(15) |

is an eigenstate of with eigenvalue . Then

(16) |

is an eigenstate of with the same eigenvalue. This allows us to identify

(17) | ||||

(18) |

This statement can be proved by looking at the structure of . First, we notice that

(19) | ||||

(20) |

With these we can write

(21) |

For a commensurate structure, the set is finite and periodic modulo reciprocal lattice vectors of the commensurate lattice. Therefore coincides with appart from a reordering of the vectors. Assume that this reordering is implemented by a permutation matrix such that . Since both the Hamiltonian and the vector are reordered in the same way, we can write

(22) |

For an incommensurate structure, the set
is infinite and therefore, appart from a reordering we have that
and coincide^{2}^{2}2Notice that this is only true formally for the incommensurate case,
not being true for any finite truncation of the matrix and we also obtain Eq. (22). This proves
our statement for both the commensurate and the incommensurate case.
Following the same argumentation it can also be formally shown that

(23) | ||||

(24) |

In this section, the coupling between the two layers was assumed to only give origin to hopping terms between the two layers, not affecting the intralayer Hamiltonians and , which preserve the translational symmetry of the isolated layers. Besides this effect, there is also the possibility of one of the layers inducing a potential to which the electrons in the other layer will be subjected to Ortix et al. (2012); Yankowitz et al. (2012); Wallbank et al. (2013). The coupling between the layers can also lead to structural relaxation, which gives origin to a modulation of the intralayer hoppings due to the displacement of the atomic positions San-Jose et al. (2014); Slotman et al. (2015); Jung et al. (2015, 2017). Although we will not explore those effects in the present work, we note that these corrections will have a spatial modulation given by and can therefore be incorporated in the present formalism by including off-diagonal terms in the matrices and of the form

(25) |

where describes the potential and intralayer hopping modulation on the top (bottom) layer.

## Iii ARPES in vdW structures

We wish to model an experimental situation where the incident electromagnetic field is monochromatic with frequency , and the electron detector is placed at position , far away from the crystal sample, collecting electrons emitted with energy along the direction . In this situation, the energy resolved ARPES intensity can be evaluated from Adawi (1964); Mahan (1970); Schaich and Ashcroft (1970, 1971); Caroli et al. (1973); Feibelman and Eastman (1974); Pendry (1976) (see also the Appendix A for a brief derivation)

(26) |

where is the chemical potential, the index runs over crystal bound states, with corresponding wavefunction and energy ; is the spectral function, which in the non-interacting limit reduces to ; and are the ARPES matrix elements. These are given by

(27) |

where is the screenedFeibelman and Eastman (1974) vector potential and is the complex conjugate of a photoemitted state with energy , which is the solution of the Lippmann-Schwinger equation

(28) |

where , with , is the free space electronic Green’s function (which is explicitly given by Eq. (81)) and is the crystal potential. Using integration by parts, we can write the ARPES matrix element as

(29) |

where we have neglected the contribution due to ,
where is the total charge in the crystal, which is a common
approximation Feibelman and Eastman (1974); Pendry (1976); Lüders et al. (2001). In the
same spirit, we neglect effects of screening in
and assume it to be described by a plane wave ,
where is the amplitude,
is the wavevector, satisfying
with the speed of light, and
is the polarization vector for the polarizations. For
simplicity we will also approximate
by a plane wave Shirley et al. (1995),
which will greatly simplify the form of ,
while providing a non-trivial description of the ARPES matrix elements.^{3}^{3}3This approximation is sometimes insufficient to explain the observed
data as shown in Ref. Gierz et al. (2011). In this case, multiple
scatterings of the photoemitted electron by the lattice become important
and a full solution of the Lippmann-Schwinger equation (28)
is required. With these approximations we obtain

(30) |

and the ARPES matrix element becomes proportional to the Fourier transform of the crystal bound states.

In order to obtain the ARPES intensity for a bilayer vdW structure we need to evaluate Eq. (30) with the crystal bound state given by Eq. (8). In real space we have that

(31) |

and the ARPES matrix elements Eq. (30) become

(32) |

and we is the transferred momentum, with and indicating the components parallel and perpendicular to the axis, and is the Fourier transform of the Wannier wavefunctions. Using the relations given by Eqs. (17), (18), (23) and (24) together with the in-plane momentum conserving Kronecker symbols, we can write

(33) | ||||

(34) |

and the ARPES matrix elements can be rewritten as

(35) |

where we have defined

(36) |

and assumed that the total area of both layers is the same, that is Koshino (2015). It is still necessary to evaluate . Assuming that the Wannier functions are well localized they can be written in a separable form as Shirley et al. (1995)

(37) |

where is the radial wave function and are real spherical harmonics, where

(38) | ||||

(39) |

and being associated Legendre polynomials. Using the plane wave expansion

(40) |

where is a spherical Bessel function, we can write as

(41) |

where and we have used the orthogonality of real spherical harmonics. For the case in which are given by hydrogen-like wavefunctions

(42) |

where , , is the effective Bohr radius, with the Bohr radius and the effective nuclear charge, and is a generalized Laguerre polynomial; then can be evaluated analytic as is given by Bransden and Joachain (1982); Moser (2017)

(43) |

where and is a Gegenbauer polynomial. Analytic expressions for are also available if are approximated by Slater type Belkić and Taylor (1989) or Gaussian type Yükçü (2016) orbitals.

Summarizing the results of this section, we can write the ARPES intensity corresponding to photoemitted electrons with energy emitted along direction , due to an incident electromagnetic field with frequency , wavenumber and polarization vector as

(44) |

where, , , ,

(45) |

and in order to include broadening effects we can approximate the spectral function by a Lorentzian

(46) |

with the broadening factor.

We make two remarks regarding Eq. (45). First, we point out that its form depends on the chosen convention to define the Fourier components of the electronic operators in Eq. (4). It is also possible to work with an alternative convention, where

(47) |

The operators and are related via . Using the convention of Eq. (47), Eq. (31) would be written as

(48) |

where and . With these definitions, would be obtained from Eq. (45) by replacing . Secondly, we would like to emphasize that the factors of in Eq. (45) ensure that ARPES intensity has the same symmetries under rotations along the direction as the bilayer vdW structure .

We would like to point out that the present formalism can also be applied to model ARPES in any system subject to competing periodicities, such as systems displaying CDW Flicker and van Wezel (2016).

## Iv An application: ARPES in twisted bilayer graphene

We now apply the general formalism developed in the previous section to the case of twisted bilayer graphene. The Bravais basis vectors for the bottom layer are written as

(49) | ||||

(50) |

where is the lattice parameter, and the positions of the A and B atoms are given by

(51) | ||||

(52) |

The top layer is rotated with respect to the bottom one by an angle of such that , for , and , for , where is the rotation matrix

(53) |

and is the separation between the two layers. The corresponding reciprocal space basis vectors are given by

(54) | ||||

(55) |

for the bottom layer, and for the top layer we have . In Fig. 1 we show the first Brillouin zone of both layers. It is also shown the Brillouin zone of the Moiré superlattice, whose associated reciprocal lattice basis vectors are given by Amorim et al. (2016b). We will describe each individual graphene layer within the nearest-neighbour tight-binding model for orbitals Castro Neto et al. (2009), which reads

(56) |

where is the nearest-neighbout hopping and we have written . For the interlayer coupling, we have , where the function can be written in terms of Slater-Koster parameters as

(57) |

where is the distance between the two atoms with the distance projected in the plane, and we use the parametrization of Refs. Moon and Koshino (2013); Koshino (2015):