{'fraction_non_alphanumeric': 0.030945294057271887, 'fraction_numerical': 0.011888877484689862, 'mean_word_length': 4.35682213873362, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Theory of the Raman spectrum of rotated double-layer graphene

We study theoretically the Raman spectrum of the rotated double-layer graphene, consisting of two graphene layers rotated with respect to each other by an arbitrary angle \theta. We find a relatively simple dependence of the Raman G peak intensity on the angle \theta. On the other hand, the Raman 2D peak position, intensity, and width show a much more complicated dependence on the angle \theta. We account for all of these effects, including dependence on the incoming photon energy, in good agreement with the experimental data. We find that it is sufficient to include the interaction between the graphene layers on the electronic degrees of freedom (resulting in the occurrence of van Hove singularities in the density of states). We assume that the phonon degrees of freedom are unaffected by the interaction between the layers. Furthermore, we decompose the Raman 2D peak into two components having very different linewidths; these widths are almost independent of the angle \theta. The change in the intensity and the peak position of one of these two components gives insight into the dependence of the overall Raman 2D peak features as a function of the angle \theta. Additionally, we show regions in the phonon spectrum giving rise to the Raman 2D peak signal. This work provides an insight into the interplay between the mechanical degree of freedom (angle \theta) and the electronic degrees of freedom (singularities in the density of states) in rotated double-layer graphene. Additionally, this work provides a way to establish experimentally the value of the rotation angle \theta using Raman spectroscopy measurement. This procedure becomes even more robust if one repeats the Raman spectroscopy measurement with a different incoming photon energy.


I. INTRODUCTION
The electronic band structure of a single graphene layer near the Fermi level consists of Dirac-cone like structure at the Brillouin zone edge (K point). In this work we study the rotated double-layer graphene (also refereed to as the twisted bilayer graphene) which consists of two single-layers of graphene that are rotated with respect to each other by an arbitrary angle θ. In the special case when θ = 0 • the Dirac cones from the two layers are exactly on top of each other in reciprocal space. However, rotation of one of the graphene layers in real space (θ = 0 • ) is accompanied by a corresponding rotation of its band structure in reciprocal space (around the origin of the reciprocal space). Therefore, when θ = 0 • Dirac cones of the two graphene layers are no longer on top of each other in reciprocal space, but are separated, proportionally to ∼ sin θ/2. Nevertheless the two Dirac cones are still overlapping in a small region of reciprocal space in between the cones. From a perturbation theory argument, one would expect that the interaction between the Dirac cones of the two graphene layers will be particularly strong in region where the Dirac cones are overlapping. Indeed, interaction between the layers in the overlap region opens a hybridization gap and leads to Van Hove singularities in the density of states of the rotated double-layer graphene. Since the position of the overlap region depends on the angle θ, we expect that the rotated double-layer graphene will have an interesting coupling between the mechanical degree of freedom (angle θ) and the electronic degrees of freedom (singularities in the density of states). Many interesting properties of rotated double-layer graphene arise from this tunability, and they have recently been attracting a lot of interest [1][2][3][4][5][6][7][8][9][10][11][12][13][14] .
In this work we study theoretically the influence of the angle θ on the Raman spectrum of rotated doublelayer graphene. Raman spectroscopy is an experimental technique commonly used to characterize carbon based materials as discussed in detail in Ref. 15. Since Raman spectroscopy uses incoming photons with a well defined energy, one can use this spectroscopy to study selectively certain regions of the electronic density of states. Therefore, we can expect an interesting dependence of the Raman signal of the rotated double-layer graphene as the angle θ is varied. Such a dependence of the Raman signal was demonstrated in some recent experimental studies 13,14,[16][17][18][19][20][21] .
The two most prominent Raman signals in graphene based systems are Raman G peak (close to 1600 cm −1 ) and Raman 2D (or G ) peak (close to 2700 cm −1 ). The Raman G peak in the graphene based systems is a simpler process than the 2D peak since it involves creation of just one phonon per each scattered photon. Considering momentum conservation and assuming a negligible momentum of the photon we conclude that the created phonon must occur at the Brillouin zone center. On the other hand, Raman 2D peak involves emission of two phonons per each scattered photon. In this case the momentum conservation implies that the two emitted phonons have arXiv:1304.5002v2 [cond-mat.mes-hall] 1 Nov 2013 arbitrary but opposite momenta. In this work we study both of these Raman peaks (G and 2D) in the rotated double-layer graphene as a function of angle θ. We find a relatively simple dependence of the Raman G peak on the angle θ. Namely, when the incoming photon energy is comparable to the separation between the Van Hove singularities of the rotated double-layer graphene, there is a significant increase (∼70) in the G peak intensity. On the other hand, the Raman 2D peak intensity, position, and width show a much more complicated dependence on the angle θ. All of these features, including the incoming light frequency dependence, are well reproduced in our calculation and agree well with experimental data (detailed comparison is shown in Ref. 13). Furthermore, these results provide a simple way to experimentally determine the angle θ of a rotated double-layer graphene. The angle determination procedure becomes even more robust if one performs Raman spectroscopy with two (or more) different incoming photon energies.
We compute the Raman spectra of the rotated doublelayer graphene using a super-cell tight-binding method. Additionally, for the Raman 2D peak we confirm our findings using a continuum model method. In the super-cell method we choose special values of θ for which there exists a super-periodicity between the two graphene layers. In the super-cell method we treat this enlarged commensurate super-cell as a unit cell of our system. On the other hand, in the continuum model calculation we rely on a simple Dirac equation description of a singlelayer graphene and we add interaction with the other layer in the restricted Hilbert space. These continuum model calculations are less numerically demanding than the super-cell calculations, since they do not rely on the super-periodicity between the two graphene layers. However, we expect that the super-cell tight-binding method is more reliable, and we find that it compares better with experimental data. Unless explicitly mentioned, the results reported here refer to the super-cell tight-binding method.
We provide details of both the super-cell tight-binding method and the continuum model method in Sec. II. In Sec. III we present results of the Raman G and 2D peaks calculations in the rotated double-layer graphene case as a function of the angle θ. We also provide a detailed analysis of these Raman peaks in the rotated double-layer graphene.

II. METHODS
In the Raman process the incoming photon creates a virtual electron-hole pair which then emits (or in some cases, absorbs) a phonon excitation quantum (or a quantum of some other excitation). In the first-order Raman process, the interaction of the single incoming photon results in the emission of a single phonon excitation. In the second-order Raman process two phonons are emitted for each interaction of the incoming photon. For this reason, measurement of the spectrum of the inelastically scattered outgoing photons is a sensitive probe of the electron and phonon degrees of freedom in the sample.
Therefore to describe theoretically the Raman spectrum of the rotated double-layer graphene, we need to know its electron and phonon band structures. Furthermore, we need to evaluate the matrix element of the interaction between the electrons and light, and of the interaction between the electrons and phonons. In the remainder of this section we describe how we computed all of these quantities in the case of the rotated doublelayer graphene. For simplicity, we only compute the Raman intensity in the rotated double-layer graphene relative to the Raman intensity in the single-layer graphene. Therefore, in our calculation we don't include explicitly numerical prefactors common to these two cases.

A. Rotated double-layer graphene unit cell
We start by defining the geometry of the rotated double-layer graphene unit cell.
The single-layer graphene unit cell consists of two carbon atoms (A and B) arranged in a two-dimensional honeycomb lattice. The rotated double-layer graphene consists of a stack of two single graphene layers rotated with respect to each other by some angle θ. We define the angle θ as follows. We start from two identical copies of single-layer graphene, translated along the axis perpendicular to the graphene plane. Next, we perform the rotation of one of the layers by angle θ around the axis passing through one of the carbon atoms. We refer to these two layers either as the top and the bottom layer, or as L = 1 and L = 2 layer.
In the case of our super-cell tight-binding method we work with angles θ for which the resulting double-layer structure is periodic. As shown in Ref. 2 every periodic double-layer structure is characterized by a pair of integers p and q up to a relative translation of layers. Furthermore, angle θ is related to these integers as A few examples of the periodic double-layer structures are shown in Fig. 1, and they all have characteristic moire pattern resulting from the misalignment of two periodic structures. The primitive unit cell of the rotated doublelayer graphene is indicated by black hexagons in Fig. 1.
As can be seen from the Fig. 1, the primitive cell of the rotated double-layer graphene has a much larger area than the single-layer primitive cell containing only two carbon atoms. As shown in Ref. 2, area of the unit cell of the rotated double-layer graphene is N times larger than the single-layer unit cell area, where N is given as Here gcd(a, b) is greatest common divisor of integers a and b. The smallest possible values of N compatible with Eq. 2 are 7, 13, 19, 31, 37, 43, and 49. An N -fold increase of the real space primitive cell is accompanied with a folding of the electron and phonon band structure in the reciprocal space. The folded Brillouin zone area is smaller by a factor of 1/N as compared to the single-layer Brillouin zone. The thick black hexagon in Fig. 2 shows the Brillouin zone of the rotated double-layer graphene for two choices of p and q, while the thick red and blue hexagons indicate the single-layer Brillouin zones of the two individual layers. As one can see from the figure, the rotated double-layer graphene Brillouin zone in these two cases needs to be repeated six [for case from Fig. 2(a)] or twelve [ Fig. 2(b)] additional times in order to cover the same area as the underlying single-layer graphene Brillouin zone. Corresponding realspace super-cells for these two Brillouin zones are shown in panels (g) and (h) of Fig. 1.
Throughout this paper we use a convention in which the wavevector from the Brillouin zone of the rotated double-layer graphene is denoted by k. The reciprocal vector of the rotated double-layer graphene is denoted by G. The wavevector for the Brillouin zone of the two single-layer graphene sheets will be denoted by k and k , for the two sheets respectively. Corresponding reciprocal vectors of the single-layer graphene sheets we will denote as G and G . We always assume that vectors k, k , k , G, G , and G are given in the same coordinate system.

B. Electrons
We describe the electron wavefunctions in the rotated double-layer graphene using a tight-binding model taking into account interaction between the graphene layers. Low-energy electronic excitations in graphene are well described by the carbon π-bonds. For this reason our model includes only one p z orbital per carbon atom which we will denote simply by φ(r), for the carbon atom at the origin. We further assume that the φ(r) orbitals at the two neighboring atomic sites are orthogonal to each other Using orbitals φ(r) we construct the Bloch-like tightbinding basis functions χ kj (r) for each k-vector in the Brillouin zone and for each site j in the rotated doublelayer graphene unit cell (super-cell) as Since rotated double-layer graphene consists of two graphene layers and the primitive unit cell of each graphene layer has two carbon atoms (A and B), the unit cell (super-cell) of the rotated double-layer graphene has 4N carbon atoms. Therefore index j ranges from 1 to 4N . A sum is performed over all lattice vectors R, while the coordinate of j-th orbital in the primitive unit cell is given by the vector t j . The functions χ kj (r) satisfy the periodicity requirement of the Bloch theorem so we can write the m-th electron eigenstate ψ m k (r) simply as a linear combination of the basis functions χ kj (r), The band index m again ranges from 1 to 4N , while only half (2N ) of these bands are assumed to be occupied for undoped systems (as in the single-layer graphene case). Next, by solving the Schrodinger equation for the electrons in the χ kj (r) basis we obtain a set of C m kj coefficients at the each vector k of interest. In order to construct the Schrodinger equation, we use the Slater-Koster parametrization from Ref. 3 fitted to the density functional theory calculation of the rotated double-layer graphene. We also rescale all tight-binding hopping parameters from Ref. 3 by 18% to account for the GW computed self-energy effects 22-25 .

Unfolding of the electron states
We now describe a procedure in which one can rewrite (unfold) the rotated double-layer graphene electron wavefunction in terms of the single-layer graphene basis functions. This procedure will be useful later in the computation of the electron-phonon matrix element of the rotated double-layer graphene.
In a close relation to the Eq. 3 let us now define the following basis functions Here, α = (L, λ) is a composite index where L = 1, 2 indexes the graphene layers and λ = A, B indexes the A and B carbon atoms. The sum over index j in Eq. 5 is performed over all atoms of type α, i.e. over all A or B carbon atoms in either first or second graphene layer. Functions ξ kα (r) respect the periodicity of the singlelayer graphene sheet in the same way that χ kj (r) respect the periodicity of the rotated double-layer graphene. For this reason if we consider indices α from either first or second layer (L = 1 or 2), functions ξ kα (r) become Bloch-like tight-binding basis functions of the first or the second single-layer graphene [in the same way in which χ kj (r) are the basis functions of the rotated double-layer graphene].
Computing the overlap between the rotated doublelayer graphene wavefunction ψ m k (r) and the single-layer graphene basis function ξ κα (r) gives ξ κα |ψ m k = n sc δ κ,k+G P m kGα , where we have defined quantity P m kGα as In Eq. 6, the term δ κ,k+G equals 1 only if κ and k differ by one of the rotated double-layer reciprocal vector G, while n sc is the total number of the super-cells in the entire sample. Since the |ξ κα basis is complete, Eq. 6 implies that the electron wavefunction |ψ m k of the rotated double-layer graphene can be rewritten (unfolded) in terms of the basis functions |ξ k+Gα of the single-layer graphene (here G are the reciprocal vectors of the rotated double-layer graphene). Furthermore, unfolding amplitude of the rotated double-layer graphene electron wavefunction |ψ m k in terms of the single-layer graphene basis function |ξ k+Gα is given by quantity P m kGα defined in Eq. 7.
We perform the unfolding procedure on a fixed set of G vectors, G, such that the following two constraints are satisfied. First, vectors k + G with two different choices of G from the set G never differ from each other either by G or G ( reciprocal vectors of two single-layer graphene sheets), as this would lead to double counting. Second, every unique single-layer wavevector k or k  Fig. 1]. Arrows indicate G vectors from the set G (see main text). Determination of the set of vectors G is relatively easy for these choices of p and q, while it becomes more involved for some other choices since some elements of G must point outside of the single-layer Brillouin zone.
can be written as k + G for some G from G and k from the double-layer graphene Brillouin zone. Black arrows on Fig. 2 indicate two examples of a set of G vectors G satisfying these constraints.

C. Phonons
In this work we assume that the interaction between the two layers of graphene does not affect the phonon band structure of the rotated double-layer graphene. Nevertheless, working in the Brillouin zone of the rotated double-layer graphene we need to take into account that the set of the single-layer phonons at wavevectors q + G are folded to the wavevector q for all G from G. Unlike for the electron states, which do get affected by the interaction between the two graphene layers, here the unfolding procedure corresponds simply to the relabeling of states. In the folded double-layer Brillouin zone we denote phonon states with the wavevector label q and the branch label ν. On the other hand, in the unfolded singlelayer Brillouin zone, this same phonon would be labeled with wavevector q + G with vector G chosen from the set G. Therefore labels ν and G are interchangeable for a unique phonon branch (corresponding for example either to the G or the 2D peak) of the single-layer graphene.
In our calculations of the Raman G peak we use the phonon frequency of the G phonon to equal 1561 cm −1 , as found in Ref. 26. The G phonon atomic displacement pattern can uniquely be determined from the representation theory analysis of the graphene space group.
Calculation of the Raman 2D peak is more complex than that of the G peak, since it involves emission of two phonons with arbitrary (and opposite) momentum. For this reason we need information about the 2D phonon frequencies in a relatively large region of the phonon Brillouin zone close to the K-point (phonons far away from the K-point give negligible contribution to the Raman 2D peak). It was shown in Ref. 27 that the Raman 2D peak profile relies strongly on the compensation between the phonon and the electron trigonal warpings, which are shown to be of a different sign. Therefore, the Raman 2D spectrum in graphene is very sensitive to the details of the phonon band structure. For this reason we use as an input to our calculations a high-order polynomial fit to the computed 2D phonon frequencies from Ref. 27, and we include the computed trigonal warping effect. Atomic displacements of 2D phonons are inferred only at the K (and K') points of the Brillouin zone by the representation theory analysis. The atomic displacement pattern of the phonons near the K (or K') point are approximated by the displacement pattern at K (or K') point.

D. Electron-light interaction
In the weak field approximation 28 interaction of electrons with light is given by the following operator where A is the vector potential of the electromagnetic field of light in the Coulomb gauge ( = 1). Therefore the matrix element of this operator between two electron states f | and |i up to a constant equals P · f |∇|i , where P is the polarization direction of the incoming or outgoing light. Expressing electron states f | and |i in terms of the basis functions φ(r) we are left with computing φ|∇|φ where φ| and |φ are the tight-binding basis orbitals φ(r) at the different atomic sites. We work under approximation 29 that the matrix element φ|∇|φ is exactly zero between φ| and |φ that are not on the first-neighbor sites in the same graphene layer. Additionally, assuming p z -like symmetry of the φ(r) orbitals one can easily show that under this approximation all matrix elements φ|∇|φ can be determined up to a single constant prefactor.

E. Electron-phonon interaction
Interaction between electrons and phonons H ph qν is described in terms of the deformation potential δV qν . The deformation potential is defined as a change in the effective potential experienced by electrons as a result of the phonon excitation from the ν-th phonon branch with a wavevector q. Similarly as in the case of an electron-light interaction, using symmetry and taking into account the interaction between the nearest neighbors, the electronphonon matrix elements can be computed up to a constant prefactor for both G and 2D phonon modes.
Since we assumed no changes to the phonon band structure coming from the interaction between the graphene layers, we compute the electron-phonon matrix element using the same electron-phonon interaction operator as in the case of a single-layer graphene. We start from the electron-phonon matrix element between any two rotated double-layer graphene states k| dl and |k + q dl (dropping electron band indices), Next we express the rotated double-layer states in the basis of the single-layer states using Eqs. 6 and 7 (we absorb coefficients P m kGα into . . . | sl and | . . . sl for simplicity), Furthermore, here we also drop sum over α for simplicity. Remembering that the branch index ν for the folded phonon band structure is just a relabeling of vectors G 3 from the set G we can write the electron-phonon matrix element as Next, as we mentioned earlier, we assume that the electron-phonon interaction operator in the rotated double-layer graphene is the same as in the single-layer graphene. For this reason, the electron-phonon matrix operator conserves the crystal momentum of the singlelayer graphene. Therefore, only one of the G 3 vectors will give a non-zero contribution to the electron-phonon matrix element,

F. Raman intensity
Using standard perturbation technique methods 27,30 , one can show that the intensity of the outgoing photon at frequency ω out for the first-order Raman process can be computed as while for the second-order Raman process it is given by Here frequency of the ν-th (µ-th) phonon branch with the momentum q is denoted as ω ν q (ω µ q ), while the incoming light frequency is denoted as ω in . Furthermore, here for simplicity we always assume that the phonon-dependent terms (phonon frequencies, electron-phonon matrix elements) appearing in the first-order Raman process are due to the G mode, while those appearing in the secondorder Raman process are due to the 2D mode. Scattering amplitudes K ν AB and K qνµ ABC are summed over all virtually excited states A, B, and C. Sum in Eqs. 13 and 14 is performed coherently over the electron states and incoherently over the phonon states. Delta functions ensure the conservation of energy.
In both Eqs. 13 and 14 we focus only on the processes involving emission, not absorption, of phonons, and we work at zero temperature. Furthermore, we are neglecting the momentum of the photon. Therefore to conserve total momentum, the emitted phonon in the firstorder Raman process must have zero momentum. In the second-order process momentum of one phonon (q) must be compensated by that of the other phonon (−q). For this reason, the first sum in Eq. 13 is performed over zero-momentum phonons, from arbitrary phonon branch ν. Similarly, first sum in the Eq. 14 is performed over all pairs of phonons with momenta q and −q, from possibly different phonon branches ν and µ.
Scattering amplitudes K ν AB and K qνµ ABC are most easily represented graphically using Feynman diagrams as in Figs. 3 and 4. Explicit expressions for these diagrams can be found in Refs. 30 and 27, here we provide as an example contribution from Fig. 3(a) for the first order Raman process, Similarly, we also provide an explicit expression for the contribution of the second order Raman process from Fig. 4(a),  We find that it is important to include all of the firstorder diagrams for the Raman G peak (as shown in Fig. 3). For the 2D peak we include only diagrams shown in Fig. 4, since other permutations give much smaller contributions in the case of a single-layer graphene, see Ref. 27 for more details.

G. Remaining parameters
Here we discuss remaining parameters and calculation details used in this work. For the electron and hole linewidth γ appearing in Eqs. 15 and 16 for the Raman G and 2D peak intensity we use γ 2 = 190 meV and γ 2 = 201 meV for the 1.96 eV and 2.41 eV photon energy calculation respectively, independent of electron wavevector k. We choose this value of electron and hole linewidth in order to reproduce the Raman G peak enhancement factor (discussed later, in Sec. III B) consistent with experiment done with 1.96 eV incoming photon energy 13 . Nevertheless, we find this value to be somewhat consistent with the sum of linewidths coming from the electronphonon 27 (32 and 43 meV for 1.96 eV and 2.41 eV photon energy calculation respectively) and electron-electron interaction (∼ 100 meV 31 ). Using the electron linewidth coming just from the electron-phonon interaction (as done in Ref. 27) would have resulted in a much larger Raman G peak enhancement. The linewidth γ 2 used for the 2.41 eV incoming photon energy calculation was computed from the value for the 1.96 eV incoming photon energy by including the difference in the estimated electronphonon linewidths (43 meV−32 meV=11 meV).
To speed up the convergence of the Raman calculation in the case of a rotated double-layer graphene we interpolate various electron and phonon quantities from a coarser reciprocal space grid onto a finer grid. Furthermore, we neglect Raman amplitudes for which unfolding intensity, electron-light matrix element, or electronphonon matrix element fall below a certain threshold value. We check that our results are fully converged with respect to this threshold.

H. Continuum model method
In this work, we also make use of the continuum model developed in Ref. 5, in order to confirm our super-cell tight-binding based calculation of the Raman 2D peak. The continuum model, as compared to the super-cell tight-binding method, uses an electron wavefunction in the restricted Hilbert space. The Brillouin zone folding in the super-cell tight-binding calculation implies that the interaction between the graphene layers introduces hybridization between states at N different wavevectors (state with wavevectors k is hybridized with states k+G, G ∈ G). On the other hand, in our continuum model, we select only a subset of vectors G from G for which the interlayer hybridization is the strongest. In particular, a state with wavelength k in layer L = 1 hybridizes only with the electron states in layer L = 2 with wavelength k + G . In the case of our continuum model calculation, we consider three reciprocal vectors G of layer L = 1 for which the norm |k + G | is minimized. This approximation can be justified with a perturbative calculation, as in Ref. 5.
Furthermore, as compared to the super-cell tightbinding calculation, in our continuum model calculation we are using a simple parametrization 5 of the interaction strength between the graphene layers that depends on only one parameter, c. Additionally, in our continuum model we neglect trigonal warping of the single-layer graphene band structure, and assume perfectly linear Dirac cone band structure parametrized with the band velocity v F . Following Ref. 5, we take v F = 10 6 m/s and c = 0.11 eV. Since we neglect the trigonal warping effect, Raman G peak intensity vanishes in the continuum model. This is because the electron-phonon matrix elements change sign under the operation K + κ → K − κ in the continuum model, 32 where K is the Dirac point of single-layer graphene. The product of the electronlight matrix elements and the energy denominators in Eq. 15 do not change under this operation in the continuum model. The sum over scattering amplitudes in Eq. 13 therefore vanishes in the continuum model without trigonal warping. For this reason we use continuum model only to compute the Raman 2D peak.
The electron-light matrix element, the electron-phonon matrix element, and the Raman intensity in the continuum model are computed as in the super-cell based method. We use the same value of the electron linewidth in the two calculations.

III. RESULTS AND DISCUSSION
In this section we present results of our calculations of the Raman spectra in the rotated double-layer graphene.

A. Electronic structure
We start with a discussion of the electronic structure of the rotated double-layer graphene. Density of states for varying angles θ are given in Fig. 5 with thin gray lines, while that for the single-layer graphene is given with a thick black line. The density of states of the singlelayer graphene in this range of energies linearly increases with the energy as one moves away from the Fermi level (Fermi level is at the zero energy in Fig. 5). This linear dependence of the density of states originates from the well known Dirac cones at the Brillouin zone corners of the single-layer graphene band structure.
The two graphene layers in the rotated double-layer graphene are rotated with respect to each other by an angle θ. For this reason, the Dirac cones from each layer are not exactly on top of each other (in the reciprocal Step-like singularity arises from the energy maximum or minimum (as a function of momentum) while the logarithmic divergence arises from the energy saddle point. As angle θ is increased, these singularities move further away from the Fermi level (compare lighter and darker gray lines in the figure). Thick black line is showing the density of states of a single-layer graphene, multiplied by two, so that it can be compared more easily to the rotated double-layer graphene case. space) but are instead rotated with respect of each other by the angle θ. Therefore, the two Dirac cones are overlapping only in a small region of the reciprocal space, and position of this overlap in the reciprocal space depends on the angle θ. In this overlap region interaction between the two layers opens a hybridization gap, which in turn leads to the occurrence of prominent Van Hove singularities both in the occupied and empty states, whose position again depends on θ. For example in θ = 6.01 • case Van Hove singularities occur near ±0.5 eV, while for the θ = 13.17 • case Van Hove singularities occur near ±1.0 eV.

B. Raman G peak
As shown in the Fig. 5, the energy at which the Van Hove singularities occur in the rotated double-layer graphene depends on the angle θ. For larger values of θ, the Van Hove singularities occur further away from the Fermi level. In particular, for the larger value of angle θ the Van Hove singularity of the occupied states are moved to lower energies, while those of the empty states are moved to the larger energies. When separation between the Van Hove singularities of the empty and occupied states matches the incoming photon energy, we expect to see changes of the rotated double-layer graphene Raman spectrum. Angle θ for which the incoming photon For range 30 • < θ < 60 • we find almost the same Raman G peak intensity for θ = 30 • + ∆ as for the θ = 30 • −∆ case. Intensity is measured relative to a single-layer graphene. We find ∼ 70 fold enhancement in the Raman G peak intensity for 1.96 eV incoming photon energy near the critical angle, 10 • . This enhancement shifts to the higher angles θ for higher incoming photon energy (2.41 eV), consistent with shift in the Van Hove singularity. At angles away from both sides of the critical angle, we find Raman G peak enhancement close to 2 (as would be expected in the limit of no interaction between the layers). Comparison with experimental data (in good agreement with our calculation) is shown in Ref. 13. energy is close to the separation between the Van Hove singularities we will refer to as the critical angle.
The computed Raman G peak intensity in the rotated double-layer graphene is given in Fig. 6 as a function of angle θ for two different incoming photon energies (black and red line). The Raman G peak intensity in Fig. 6 is given in terms of the intensity of a single-layer graphene. We find a ∼ 70 fold enhancement of the Raman G peak intensity at angles θ close to the critical angle, 10 • (1.96 eV incoming photon energy, black line in Fig. 6). At angles below and above this critical angle we find that the Raman G peak enhancement factor is close to 2. Therefore, in that region of angles θ Raman signal of the rotated double-layer graphene is almost the same as that of two independent graphene sheets. Furthermore, we also find that the G peak enhancement shifts to the higher angles θ with higher incoming photon energy (red line in Fig. 6). This behavior we attribute to the shift in the energy of the Van Hove singularity as a function of angle θ, as observed already in Fig. 5.
Unlike the Raman 2D peak, the Raman G peak in graphene is a single phonon process and therefore its width and peak position depend solely on the phonon lifetime and frequency. We assumed in our calculation that the phonon lifetime and frequency are not affected by the interaction between the two graphene layers. For this reason, Raman G peak width and position are inde- pendent of the angle θ, in agreement with experimental observations in Ref. 13.
We find a very strong dependence of the Raman G peak enhancement at the critical angle on the effective electron and hole linewidth γ appearing in Eq. 15. Dependence of the Raman G peak enhancement at the critical angle on the value of parameter γ is shown in Fig. 7. Dotted line in Fig. 7 is a fit to the function ∼ γ −2 . As already mentioned in Sec. II G due to this strong dependence of the Raman G peak enhancement on γ, we have chosen value of γ which gives Raman G peak enhancement in agreement with experiment at 1.96 eV incoming photon energy. Nevertheless, the value of γ we obtained is consistent with that obtained from the electron-phonon and electron-electron interaction estimates.
In the Raman calculations of phonon excitations it is a common 15 practice to neglect the A 2 term in the electron-light interaction Hamiltonian (as in Eq. 8). However, Ref. 33 claims that under certain conditions A 2 terms are important for the Raman G peak process. Since these conditions are not satisfied in a typical experimental situation (ω in ∼ 2 eV, as in Ref. 13, and assuming γ = 0 would lead to a divergent G peak enhancement, see Fig. 7) we refer inclusion of A 2 term to the future work, as it would go beyond the scope of this manuscript.

Influence of coherence
We find a large influence of coherence (interference) in the calculation of the Raman G peak. (Similar observation was found in Ref. 33.) This is true both for the coherence between different Feynman diagrams (shown in Fig. 3) and for the coherence between different electronic states appearing in Eq. 13. Influence of both of these coherences is illustrated in Fig. 8. Figure 8 shows four different ways the sum given in Eq. 13 is performed. Horizontal axis of Fig. 8 shows the difference in electronic energies (∆E) appearing in the energy denominator as in Eq. 15. The vertical axis of Fig. 8 shows value of the Raman G peak intensity, if the sum in Eq. 13 is performed over all pairs of electronic states with energy separation up to ∆E. Dotted lines in the Fig. 8 show the Raman intensity for the Raman G peak if the coherent sum appearing in Eq. 13, | AB K ν AB | 2 , is replaced with an incoherent sum ( AB |K ν AB |) 2 . Solid lines show results for when the sum is performed coherently, as in Eq. 13. Additionally, dotted lines are downscaled 300 times in Fig. 8 so that they can be compared more easily to the coherent result. Blue (gray) lines in Fig. 8 represent Raman G peak intensity when the sum in Eq. 13 is performed only over two Feynman diagrams shown in Fig. 3 (a) and (b), while black lines shows results for all twelve first-order diagrams in Fig. 3.
From Fig. 8 we can reach several conclusions about the Raman G peak in graphene. First, we find that the coherence in Eq. 13 between different electronic states leads to the suppression of the Raman G peak intensity by more than 300 times. Second, in order to achieve the fully converged result, we find that the sum in Eq. 13 has to be performed over electron-hole pairs separated in energy more than the incoming photon energy ω in (1.96 eV in the case of Fig. 8). This is especially true for the coherent calculation (solid lines). Third, we find that if the sum in Eq. 13 is performed only up to the energies close to the incoming photon energy ω in , that the sum is dominated by two diagrams shown in Fig. 3 (a) and (b). However, if the sum in Eq. 13 is continued to the energies larger than ω in , Feynman diagrams from Fig. 3 (c) to (l) start to dominate (compare solid blue and black lines in Fig. 8).

C. Raman 2D peak
Similarly as in the case of the Raman G peak, we expect to see changes in the Raman 2D peak when the angle θ is close to the critical angle. In fact we find an even more complicated dependence of the 2D Raman peak on angle θ than that of the Raman G peak.
Comparing the first-order Raman calculation (as for the Raman G peak) given in Eq. 13 to the second-order Raman calculation (as for the 2D peak) given in Eq. 14 we see that in the latter case the sum is performed over all phonon momenta q in the entire phonon Brillouin zone. Phonons at different momenta q have different frequency, ω ν q , which in general would lead to the Gaussianlike spread in the Raman intensity I 2 (ω out ), even if the phonon lifetime is infinite. This observation is not true for the Raman G peak since it involves only a single phonon frequency ω ν 0 , and therefore its Lorentzian-type width comes solely from the finite phonon lifetime, and its peak position is determined by ω ν 0 . The super-cell tight-binding method computed Raman 2D peak position, intensity, and width in the rotated double-layer graphene are given in Fig. 9. For all three features of the 2D peak we find a complex variation as a function of the angle θ, especially so near the critical angle (∼ 10 • for 1.96 eV incoming photon energy). Similarly as in the case of the Raman G peak, we find that these features shift to the larger angle θ if the incoming photon energy is increased. Again, as in the case of the Raman G peak, this behavior is consistent with the angle θ dependent position of the Van Hove singularities shown in Fig. 5.
The position of the Raman 2D peak in Fig. 9 is indicated relative to the single-layer graphene case at the same incoming photon energy (since even for the singlelayer case the Raman 2D peak position depends on the incoming photon energy). We find that the position of the Raman 2D peak of the rotated double-layer graphene is shifted to the larger energies with respect to the single-layer graphene case. The observed shift is nonmonotonic, starting out small (∼ 4 cm −1 ) at large angles (> 20 • ). Close to the critical angle ∼ 10 • (for 1.96 eV incoming photon energy) the shift in the peak position increases to ∼ 14 cm −1 and is followed by a steep drop to ∼ 4 cm −1 at about 7 • . At even lower angles (< 7 • ) there is a steep rise in the 2D peak position.
The intensity of the Raman 2D peak is somewhat less complicated than the peak position and the peak width.
The Raman 2D peak intensity shows almost a step-like change close to the critical angle ∼ 10 • (for 1.96 eV incoming photon energy), having an intensity comparable to two independent single-layers at higher angles, and ∼ 4 times smaller intensity at the smaller angles.
The width of the Raman 2D peak at angles above 15 • (for 1.96 eV incoming photon energy) is comparable to that of a single-layer graphene, ∼ 31 cm −1 . At the smaller angles (< 15 • ) there is a sharp increase in the Raman 2D peak width. Additionally, close to 8 • Raman 2D peak width suddenly jumps to 52 cm −1 . Below 8 • there is again a non-monotonic behavior of the width, starting with a decrease followed by a sharp increase below 6 • .
The results of the continuum model calculation of the Raman 2D peak in Fig. 10 show similar overall features as the super-cell tight-binding calculations. The angle dependence of the peak position, intensity and width follow the same trends in both calculations, but the numerical values are somewhat different. In addition, there are some spurious features present in Fig. 10 that are not present in the super-cell tight-binding calculation. For example, the Raman 2D peak intensity and width in the region from θ = 5 • to θ = 15 • show some small features not present in the super-cell tight-binding calculation. We expect that these differences are occurring due to the approximations introduced into the continuum model calculation (see Sec. II H). In particular, the lack of the trigonal warping in the continuum model becomes especially important at large energy of the incoming photons and for large angle θ. Additionally, the reduction of the Hilbert space becomes more important at low angles θ.

Procedure to experimentally determine angle θ
The dependence of the position, the intensity, and the width of the Raman 2D peak on the angle θ provide a simple route to experimentally determine angle θ. However, for some range of values of θ the position and the width of the Raman 2D peak depend non-monotonically on the angle θ. Naively, one would expect that this would make it impossible to uniquely determine θ in that range of angles. Nevertheless, combining all three properties of the Raman 2D peak (position, intensity, and width) make it easier to uniquely assign angle θ. Furthermore, combining Raman measurements at two different incoming photon energies gives additional way to uniquely determine the angle θ even in the region where the position and the width of the Raman 2D peak depend non-monotonically on θ. For example, if one measures for the incoming photon energy of 1.96 eV change in the Raman 2D peak position of 8 cm −1 , according to the black line in Fig. 9 this measurement can correspond to angle θ of ∼ 5 • , ∼ 10 • ,  Fig. 14. Intensity is defined as the area under the peak (not peak height). Width is defined as the full width at half of the peak maximum (FWHM). Peak intensity and peak position are defined relative to a single-layer graphene. See main text for more details. Comparison with experimental data (in good agreement with our calculation) is shown in Ref. 13. or ∼ 15 • . However if one repeats the measurement on the same sample with a larger incoming photon energy (for example 2.41 eV as shown by the red line in Fig. 9) and the change in the Raman 2D peak position becomes smaller than 8 cm −1 , angle θ can be assigned uniquely to ∼ 10 • . On the other hand, if the Raman 2D shift becomes larger than 8 cm −1 then θ is either ∼ 5 • or ∼ 15 • . Since these two angles are quite far apart (by construc- FIG. 10. (Color online.) Position, intensity, and width of the Raman 2D peak using a more approximate method (continuum model). Conventions are the same as in Fig. 9, but the range of vertical scales is not the same as in Fig. 9. tion), other Raman 2D features like intensity or width can be used to determine which of the two angles should be assigned. Similar procedure can also be used for the non-monotonic dependence of the Raman 2D peak width.

Decomposition into contributing phonons
According to the Eq. 14 second order Raman process (as for the Raman 2D peak) can be decomposed into a decoherent sum of contributions coming for the pair of phonons (q, µ) and (−q, ν) with opposite momenta q and possibly different phonon branches µ and ν. Since the phonon branches of the 2D mode arise from the Brillouin zone folding, branch indices µ and ν can be relabeled with the rotated double-layer graphene reciprocal vectors G (as discussed in Sec. II C). Figure 11 shows regions of the phonon Brillouin zone which contribute the most to the Raman 2D peak, for varying angle θ. Contributions from the phonon pair (q, µ) and (−q, ν) is equally distributed among the un- Color scale for each panel is scaled individually to the largest intensity for that panel, since otherwise the overall intensity of the panels for small angles would be too small (see Fig. 9 showing decrease of the Raman 2D peak intensity at small angles θ).  folded vectors q+G with reciprocal vector G corresponding to both µ and ν. Brillouin zone of both bottom (solid line) and top (dashed line) single-layer graphene are indicated with black lines. For simplicity, only contributions from phonons in one graphene layer are shown in Fig. 11, and only the region close to the Brillouin zone corner (K point) is shown.
For large values of angle θ [for example θ = 27.80 • in Fig. 11(h)] we find a characteristic triangular region (red) in the phonon Brillouin zone around the K-point with the largest contribution to the Raman 2D peak. Similar behavior we find in the calculation of a single-layer graphene, as consistent with the decomposition found in Ref. 27. At angles smaller or equal to the critical angle, this triangular region is significantly modified. Largest modification we find when the K-point of the Brillouin zone of the top graphene layer is overlapping with the triangular region in the bottom layer [see for example Fig. 11(e)]. As shown in Fig. 12 this modification occurs precisely at the critical angle, at which the Dirac cones in the electron Brillouin zone are overlapping.

Peak substructure, two Gaussian components of the 2D peak
Our calculations show that the profile of the Raman 2D peak [I 2 (ω out ) in Eq. 14] can be well fitted with two Gaussians with varying position, intensity, and width of each Gaussian function (compare black and yellow line in Fig. 13). We find this to be true both for the singlelayer graphene and for the rotated double-layer graphene. For the single-layer graphene importance of using two Gaussians as opposed to only one is more subtle. However, for the rotated double-layer graphene just below the critical angle, positions of these two Gaussians are somewhat apart from each other, leading to the more pronounced two-peak feature. Similar feature has been found in the experimental measurements, near the critical angle 34 . Furthermore, these two Gaussian components of the Raman 2D peak behave differently as a function of angle θ which will be of interest in analyzing angle θ dependent data for the rotated double-layer graphene.
First, let us analyze these two Gaussian components in the case of a single-layer graphene. We find that these Gaussian components in this case are centered around nearly the same frequency (difference is only 3.5 cm −1 at 1.96 eV incoming photon energy) and have nearly the same intensity. Additionally, we find that the width of one Gaussian component (narrow component) is 30 cm −1 while the width of the other Gaussian component (broad component) is almost two times larger, 59 cm −1 . Figure 14 shows the position, the width, and the intensity of these two Gaussian components in the case of a rotated double-layer graphene (broad and narrow Gaussian components are shown with different color in Fig. 14). Data in Fig. 14 is shown for the super-cell tightbinding calculation, but similar results are obtained with Quite surprisingly, we find that the broad Gaussian component of the Raman 2D peak in the rotated double layer graphene is nearly independent of the angle θ. There is an overall decrease in the intensity of the broad component below the critical angle (∼ 10 • ) but the changes in the position and the width are almost negligible.
For the narrow Gaussian component in the rotated double layer graphene we again find that its width almost does not depend on the angle θ. On the other hand, the peak intensity and the peak position of the narrow component show a drastic change below the critical angle (∼ 10 • ). In particular, exactly at the critical angle the narrow component nearly vanishes. Below the critical angle (5 • < θ < 10 • ) the narrow component reappears but with significantly lower peak position (−3 cm −1 below the critical angle as compared to 18 cm −1 above the critical angle). At the even lower angle (θ < 5 • ) the narrow component nearly disappears once again.
This appearance and disappearance of the narrow component gives an insight into the complex behavior of the overall position, intensity, and width of the Raman 2D peak (black line in Fig. 14). For example, the overall increase in the width of the Raman 2D peak near the critical angle (∼ 10 • ) can be explained by the disappearance of the narrow Gaussian component at the same angle. Similarly, reappearance of the narrow component with lower frequency below the critical angle (5 • < θ < 10 • ) explains the overall change in the peak position of the Raman 2D peak. Additionally, reappearance of the narrow component at the lower frequency than the broad component is consistent with the experimentally observed two peak structure of the Raman 2D peak in the same range of angles θ.
It is tempting to interpret the broad and narrow Gaus- sian components of the 2D peaks as coming from the corners of the triangular region (inner phonons, Ref. 27) in Fig. 11 and from the triangular faces (outer phonons) respectively. Indeed, similar two-peak feature of the Raman 2D peak has been found in Ref. 27, but for significantly larger incoming photon energies (3.8 eV). These two features of the Raman 2D peak were denoted as 2D + (inner) and 2D − (outer) in Fig. 26 of Ref. 27. However, origin of the two peak features we find here is decidedly different. We demonstrate this by taking our single-layer graphene calculation and considering only small slices (in certain region of angles around the K-point) of the triangular regions in the phonon Brillouin zone either near the triangular corners or faces. We find in both cases that the two-Gaussian peak feature persists, with similar fitting parameters. Instead, we find that this two-peak structure of the Raman 2D peak originates from the sum over electronhole pair states in Eq. 14 (not different phonon states as for the feature found in Ref. 27). In particular, we find that the electron-hole pairs which are separated by the energy close to the incoming light energy give rise to the narrow component of the 2D peak from Fig. 14, while the higher energy electron-hole pairs give rise to the broad component from Fig. 14. More specifically, for the incoming photon energy of 1.96 eV, we find that the narrow component of the 2D peak originates from the electron-hole pairs separated up to ∼ 2.1 eV. Electronhole pairs between ∼ 2.1 and ∼ 2.6 eV give rise to the broad component.

Influence of interlayer interaction on electron wavefunctions and eigenenergies
The tight-binding model of the rotated double-layer graphene used in our study is based on a Slater-Koster parametrization from Ref. 3. This parametrization assigns a hopping term to any pair of p z orbitals on two carbon atoms. These carbon atoms can either be in the same, or two different graphene layers. Therefore, if we set to zero all hopping terms between pair of carbon atoms in the different graphene layers (interlayer hopping), we can effectively turn off the interaction between the two graphene layers.
The effect of allowing the electron interlayer hopping in our calculation is twofold. Firstly, it affects electron wavefunctions. The change in the electron wavefunctions modifies electron-light and electron-phonon matrix elements, which in turn changes Raman intensity of both G and 2D peak, as given for example in the numerators of Eqs. 15 and 16. Secondly, interlayer hopping affects electron eigenenergies. Electron eigenenergies in turn affect Raman G and 2D intensities through the denominators in for example Eqs. 15 and 16. Figure 15 shows which features of the Raman 2D peak can be explained solely by the influence of the interlayer hopping on the electron wavefunctions, and which by the influence on the electron eigenenergies. Dotted red (blue) line in Fig. 15 shows the Raman 2D peak position, intensity, and width for the calculation in which the interlayer hopping is given only for the electron eigenenergies (electron wavefunctions). Solid black line in these graphs are the same as in Fig. 9, showing the results of the full Raman 2D peak calculation (with interlayer hopping considered both for electron eigenenergies and wavefunctions).
From Fig. 15 we conclude that the position of the Raman 2D peak is almost completely determined by the influence of the interlayer hopping on the electron eigenergies. On the other hand, intensity of the Raman 2D peak is determined by the interlayer hopping influence on the electron wavefunctions. Finally, increase in the width of the Raman 2D peak at low angles θ is well described by the influence of the interlayer hopping on the electron wavefunctions. However, influence of the interlayer hopping on the electron wavefunctions does not reproduce feature in the Raman 2D peak width near the critical angle (∼ 10 • ).

D. Limit of small and limit of large angles
Here we discuss properties of the Raman 2D and G peaks of the rotated double-layer graphene in the limit of small (close to 0 • ) and large (close to 30 • ) angles θ. For the Raman G peak we find that in both limits (0 • and 30 • ) intensity of the G peak is similar to that of a single-layer graphene (multiplied with number of layers in the rotated double-layer graphene, two). In fact, for the entire range of angles θ, except close to the critical angle, we find that the Raman G peak intensity is similar to that of a single-layer graphene (times two).
The situation with the Raman 2D peak is again more complicated. Figure 16 shows calculated Raman 2D profiles for the rotated double-layer graphene (black) shifted for clarity in the vertical direction proportionally to the value of the angle θ. The Raman 2D profile of the single-layer graphene (multiplied by two) is indicated with thicker red line in Fig. 16. From Fig. 16 one can see that the Raman 2D spectrum of the rotated double-layer graphene above θ ≈ 15 • is already converging towards that of a single-layer graphene (red).
On the other hand, in the limit of a small angle θ (close to 0 • ) Raman 2D peak intensity of the rotated doublelayer graphene is significantly smaller than that at the larger angles, or that of the single-layer graphene. We find similar reduction in intensity in the case of the AB (blue in Fig. 16) and the AA (green in Fig. 16) stacked double-layer graphene. Additionally, peak position and width for small angles θ are qualitatively similar to that of the AB and AA stacked double-layer graphene. Similarity with the AB and AA stacked double-layer graphene is not unexpected since the rotated double-layer graphene in the limit of very small angles θ is composed of a hexagonal super-periodic arrangements of AB and AA stacked regions. This pattern is already visible to some degree on Fig. 1(b) for the case of θ = 9.43 • and is even more prominent at smaller angles θ.
However, in the sharp contrast to the AB and AA stacked double-layer graphene, we find no prominent multi-peak structure in the case of the rotated doublelayer graphene in the limit of a very small angle θ. Furthermore, double peak structure discussed earlier in Sec. III C 3 is of a different origin, and separation in frequency between the two Gaussian components is much smaller.

IV. SUMMARY AND OUTLOOK
In this work we provided a theoretical description of the two most prominent Raman signals in rotated doublelayer graphene (G peak and 2D peak). We find a relatively simple dependence of the Raman G peak intensity on the angle θ. On the other hand, position, intensity, and width of the Raman 2D peak as a function of angle θ is much more complex. All of our findings are in good agreement with available experimental data 13 . We trace the origin of the complex dependence of the Raman 2D peak signal on the angle θ by decomposing the Raman 2D peak into two Gaussian components with quite different widths that are nearly independent on the angle θ. In fact, strong dependence of the intensity and position of one of the components is responsible for the overall changes to the Raman 2D peak.
Additionally, we discuss importance of coherence in the Raman G peak calculation. We analyze both coherence over the various electron-hole pairs, and coherence over the various Feynman diagrams contributing to the [I2(ωout = ω) from Eq. 14] for the rotated double-layer graphene (thin black lines), the single-layer graphene multiplied by two (thick red line), the AB stacked double-layer graphene (blue), and the AA stacked double-layer graphene (green). The rotated double-layer graphene spectra are shifted in the vertical direction, proportionally to the angle θ, for clarity. Raman 2D profile of the single-layer graphene (red) is shifted vertically proportional to θ = 30 • .
Raman G peak. In the case of the Raman 2D peak we analyze regions of the phonon Brillouin zone contributing to the Raman signal, and explore the influence of the interlayer interaction on the electron wavefunctions and eigenenergies.
Our study provides a way to experimentally determine angle θ of the rotated double-layer graphene using only the Raman spectroscopy measurement. Angle determination becomes even more robust if one repeats Raman spectroscopy measurement with a different incoming photon energy, as discussed in Sec. III C 1. Finally, this work provides an insight into the coupling between the mechanical degree of freedom (angle θ) and the electronic degrees of freedom (singularities in the density of states) in the rotated double-layer graphene. We expect similar effects to occur if even more layers of graphene are stacked on top of each other, or if different graphenelike two-dimensional systems are stacked on top of each other.