Evolution of the electronic band structure of twisted bilayer graphene upon doping

The electronic band structure of twisted bilayer graphene develops van Hove singularities whose energy depends on the twist angle between the two layers. Using Raman spectroscopy, we monitor the evolution of the electronic band structure upon doping using the G peak area which is enhanced when the laser photon energy is resonant with the energy separation of the van Hove singularities. Upon charge doping, the Raman G peak area initially increases for twist angles larger than a critical angle and decreases for smaller angles. To explain this behavior with twist angle, the energy separation of the van Hove singularities must decrease with increasing charge density demonstrating the ability to modify the electronic and optical properties of twisted bilayer graphene with doping.

layers are typically a few nanometers 38,39 so that the capacitance is significantly larger than that of the Si back gate. Side gate voltages applied to the device were kept to between −1.5 V and 3.0 V to avoid any electrochemical reactions in the polymer electrolyte. The gate voltage was constant during the acquisition of the Raman spectra and then was changed slowly at a rate of 1 mV/s between measurements to ensure the polymer electrolyte remained stable. The leakage current was monitored during the measurements and there was no detectable signal.
The graphene in this study was grown via a low pressure CVD method on 25 μm thick Cu foil placed within a sealed copper pouch 40 (see details in the Method section). Our growth conditions yield large monolayer graphene flakes with smaller multilayer regions originating from the nucleation sites. Figure 1(b) shows an optical microscope image of the graphene flakes after they have been transferred from the Cu foil onto an oxidized Si substrate (285 nm SiO 2 thickness) using a wet transfer method 41 . The number of layers is clearly distinguishable based on the optical contrast. The size of the monolayer region is over 200 μm on a side, while the bilayer is on the order of 25 μm and the multilayer region in the center is smaller than 10 μm. The contrast varies on different regions of the bilayer, which indicates the angle between the first and second layer varies due to different domains in the CVD grown graphene. Specifically, the slightly darker areas (e.g. blue circle in the inset of Fig. 1(b)) correspond to approximately a 12 degree twist angle which absorbs more visible light owing to the enhanced DoS in this energy range 42,43 .

Twist angle determination by STM and Raman measurements
While previous studies have employed contrast measurements to identify the twist angle 44,45 , we find Raman spectroscopy to be a more precise and convenient method. Raman spectra taken at the locations marked in the inset of Fig. 1(b) are plotted in Fig. 1(c). There is no D peak visible in the spectra either before or after the polymer Scientific REPORTS | 7: 7611 | DOI:10.1038/s41598-017-07580-3 electrolyte measurements which indicates that the measured area is free of defects 18 . In tBLG, both the Raman G and 2D peaks have been found to vary with the twist angle as measured with transmission electron microscopy 23,24 . Remarkably, on the slightly darker region (blue curve), which has a twist angle of 12 degrees, the G peak intensity is 25 times larger than on the lighter region that has a twist angle less than 3 degrees (green curve). Both bilayer spectra are reliably distinguishable from that of monolayer graphene marked by the black diamond (black curve). For the 532 nm laser in our study, the G peak is enhanced for twist angles from 10 to 15 degrees.
The twist angle between the two layers of tBLG leads to a moiré pattern, whose wavelength depends inversely on the twist angle θ. Therefore, measuring the wavelength of the moiré pattern using STM gives a direct measurement of the twist angle. All topography measurements were performed in a low-temperature ultrahigh vacuum STM at a temperature of 4.5 K. These measurements were performed before applying the polymer electrolyte to the device. Figure 1(d) shows a typical topography image of tBLG where a 1.19 ± 0.03 nm moiré pattern is visible. The upper inset of Fig. 1(d) is a zoomed in topography image showing the graphene atomic lattice. A Fourier transform of the topography image gives the wavevectors of the atomic lattice and moiré pattern. The lower inset of Fig. 1(d) is a FFT of the topography with the atomic lattice and moiré pattern wavevectors marked by red and blue circles respectively. The ratio between the atomic lattice length a and the moiré pattern wavelength D gives the twist angle θ between the two lattices using a/D = 2 sin(θ/2).
Using the above procedure, more than 15 different tBLG flakes were imaged using STM to determine their twist angles. Figure 2(a) plots the measured Raman G peak area of these tBLG flakes normalized to monolayer graphene. There is a clear enhancement of the G peak area near a twist angle of 12 degrees. Our data agrees with previous results 23,24 , where the twist angle was determined by transmission electron microscopy. The enhancement at a specific angle range arises due to the enhanced joint density of states at the photon energy of the laser excitation. The twist angle for the black points was measured using STM topography. The uncertainty in angle of the black points represents the variation in the atomic lattice length a and moiré pattern wavelength D measured by STM. The blue curve is the theoretical calculation within the angle range where the G peak is enhanced showing that the maximum enhancement occurs at 12.6°. The additional red points were only measured by Raman spectroscopy and their twist angles were determined as discussed in the text. The error bars in G peak area represent the standard deviation of the area of a Lorentzian fit to the G peak. Inset shows Raman pathways that contribute to the enhanced G peak. (b) Raman R peak position as a function of twist angle. The error bars in R peak position represent the uncertainty of the peak position when fit with a Lorentzian. (c),(d) Raman intensity near the G peak as a function of wavenumber and gate voltage for 11.1° (c) and 14.1° (d) respectively.
Scientific REPORTS | 7: 7611 | DOI:10.1038/s41598-017-07580-3 To quantify the enhancement, the Raman G peak area as a function of twist angle is calculated using 46 eV is the photon energy of the laser excitation, E ab (k) is the energy separation of electronic states a and b, E ph = 0.196 eV is the G peak phonon energy, γ is the inelastic scattering rate for which we use 0.15 eV, M el-ph is the electron-phonon coupling matrix element, and M op is the optical transition matrix element. We follow ref. 1 to calculate the electronic band structure of tBLG, which retains the linear dispersion near the Dirac points and exhibits vHs and splitting where the bands from the two layers meet. We include only the four lowest energy bands to perform the calculation as depicted in the inset of Fig. 2(a). The largest enhancement occurs when the photon energy of the laser is on resonance with the energy separation of a large number of electronic states as denoted by the green arrows in the inset of Fig. 2(a). M el-ph is taken as a constant in the calculation as we are only interested in changes in the band structure due to twist angle variations. Furthermore, the twist angle corresponding to the maximum enhancement depends on laser photon energy and therefore M el-ph can not be responsible for the observed enhancement 23,24 . M op along the line passing through the two Dirac points is modeled by a phenomenological Lorentzian as in ref. 24 with no k-space dependence in the orthogonal direction. The blue curve in Fig. 2(a) plots the result of the calculation which matches the experimentally measured data.
The center position of the peak in the calculated curve is at 12.6 degrees, which is defined as the critical angle. Using a different photon energy shifts the critical angle. The calculation allows a Raman spectroscopy measurement of the G peak area to determine two possible twist angles for the tBLG flake. To distinguish between the two possible twist angles, we use the Raman R peak associated with the transverse optical (TO) phonon, which emerges from a double-resonance intervalley process mediated by the periodic potential 26 . In the angle range where the G peak is enhanced, the R peak is clearly visible, as shown in Fig. 1(c). The R peak position decreases linearly with increasing twist angle, −11.7 ± 0.8 cm −1 /degree, as shown in Fig. 2(b) which allows the twist angle to be uniquely identified. Again, the black circles correspond to tBLG flakes whose twist angles were measured by STM. The red circles correspond to tBLG flakes whose twist angles were inferred from Raman spectroscopy measurements of the G peak area. Using a combination of the R peak position and the G peak area, the twist angle can be determined solely from Raman measurements.

Charge doping induced by a polymer top gate
To study the charge density dependence of the Raman G peak, we use a voltage on the polymer gate, V g to induce charge carriers in the tBLG. As the charge density changes, the G peak behaves very differently depending on the twist angle as shown in Fig. 2(c) and (d). These data are from two different tBLG flakes which both have approximately 12 times enhancement at zero gate voltage but two different twist angles. In both cases, the position of the G peak shifts to higher wavenumbers as seen in monolayer graphene 47,48 and previous Raman studies of tBLG 20,37 . Similarly, the 2D peak position for all tBLG flakes shows the same behavior as monolayer graphene where it is initially constant then shifts to higher wavenumbers for hole doping and lower wavenumbers for electron doping 20,37 . While these aspects of the measurements are the same for all twist angles, the intensity of the G peak shows a strong twist angle dependence. As the gate voltage increases the G peak intensity initially increases but then quickly decreases in Fig. 2(c) while it continually increases in Fig. 2(d). To understand why these two different flakes behave so differently with gate voltage, we must examine the effect of the voltage on the band structure of the tBLG.
In monolayer graphene, the voltage on the polymer gate acts to both shift the Fermi level and produce a potential between the graphene and electrode. The Fermi level is related to the charge density n by , where ħ is the reduced Planck constant and v F = 1.05 × 10 6 m/s is the Fermi velocity of monolayer graphene 39 . The potential between the gate electrode and graphene is given by where ε PE is the dielectric constant of the polymer electrolyte, E is the electric field and C g is the geometric capacitance per unit area of the Debye layer in the polymer electrolyte. Combining these two effects gives, Here V g is measured relative to the voltage when the Fermi level is at the Dirac point by finding at which gate voltage the G peak position is a minimum. To determine the geometric capacitance, we measured the G peak position for monolayer graphene as a function of the applied gate voltage. The shift of the G peak ΔΩ G is proportional to the Fermi energy, namely, ΔΩ G = αE F , where α is 42 cm −1 eV −1 taken from ref. 49. For each device, ΔΩ G versus gate voltage was measured to determine C g . For the different devices measured, the capacitance per unit area varied within the range 0.7-2.2 μF cm −2 , which is about two orders of magnitude larger than that of the Si back gate. The geometric capacitance per unit area for bilayer graphene areas was taken to be the same as that of nearby monolayers since the polymer electrolyte was uniformly distributed over these length scales. The doping at the position of the laser spot was directly measured and therefore global spatial variations do not affect the results.
In tBLG, there is an additional effect because the top layer is more heavily doped than the bottom layer due to its proximity to the polymer electrolyte 37 . This asymmetric doping causes an interlayer potential Δφ to develop. A simple model is adopted from ref. 37 to calculate the charge densities n T and n B in the top and bottom layers respectively. The tBLG is modeled as a parallel plate capacitor with two weakly coupled monolayers separated by an interlayer distance d 0 = 0.34 nm. This model is valid as long as the doping level is small compared to the separation of the vHs as is the case in our experiment. The top and bottom layers experience different electric fields, given as B tBLG , where ε tBLG is the dielectric constant of the tBLG. Similar to equation (2) in the monolayer case, for bilayer we use the charge density and electric field of the top layer to obtain, The interlayer potential energy eΔφ is equal to the energy difference between the Dirac points of the top and bottom layers ϕ In the parallel plate capacitor model, the interlayer potential is simply the electric field between the layers E B times their separation d 0 . Namely, φ ∆ = en C / B tBLG , where ε = C d / tBLG tBLG 0 is the effective interlayer capacitance per unit area. Combining the above two equations we obtain,   n T and n B can be solved for any value of V g using equations (3) and (4) and the measured C g . We take = × − C 8 10 tBLG 6 F cm −2 in our calculation 50 . In the inset of Fig. 3(a), we use these equations to convert the applied gate voltage, V g to the charge density n T (black) in the top layer, n B (red) in the bottom layer and n Total (blue) for the tBLG corresponding to Fig. 2(c). Once the conversion from gate voltage to charge density is known, the area of the G peak is fit for each gate voltage with a single Lorentzian to obtain the blue curve in Fig. 3(a). A similar procedure is used for all other tBLG flakes to obtain the remaining curves in Fig. 3(a) and (b).

Raman G peak area dependence on charge doping
The Raman G peak area of tBLG with twist angles smaller than 3 degrees or bigger than 16 degrees do not depend on the charge density as shown by the orange markers in Fig. 3(a) and (b) respectively. So, we focus on the angle range from 10 to 15 degrees where the G peak area depends on charge density. The G peak area of tBLG with a twist angle equal to the critical angle decreases monotonically for both the electron and hole doped regions as shown by the purple markers in Fig. 3(a). For all angles smaller than the critical angle, similar behavior is observed as shown by the other markers in Fig. 3(a). In contrast, for angles larger than the critical angle, the area initially increases with doping before decreasing at higher doping levels as shown in Fig. 3(b). The black curve in Fig. 3(b) corresponds to a twist angle of 13.4 ± 0.1 degrees. Its maximum G peak area occurs at a total charge density of 5.3 ± 0.9 × 10 12 cm −2 . This behavior is observed for either electron or hole doping. The total charge density required to maximize the G peak area for a slightly larger twist angle is much larger. The maximum G peak area occurs at 13 ± 2 × 10 12 cm −2 for the green curve which corresponds to a twist angle of 13.7 ± 0.1 degrees. The maxima for the blue and red curves are not yet reached within the charge density range plotted. The charge density when the maximum G peak area occurs as a function of twist angle is shown by the black markers in Fig. 4(a).
To understand this asymmetric angle dependence, we examine the modification of the tBLG band structure upon doping. In undoped tBLG as shown in the inset of Fig. 2(a), the Dirac points of the two layers are at the Fermi level since there is no charge doping or interlayer potential. The two vHs are therefore aligned in momentum exactly between the two Dirac points. However, when the sample is doped, the two Dirac points shift away from the Fermi level asymmetrically 1,[33][34][35][36][37] . When this occurs, the two vHs are displaced horizontally in momentum as seen in Fig. 4(b) and (c) for electron and hole doping respectively. Therefore, as the tBLG is doped the energy separation of the bands splits into two different contributions. The joint density of states can be defined as . Figure 4(d) shows the JDOS as a function of energy and charge density consider- ing only the interlayer potential imposed by the asymmetric doping in the two layers. Half of the JDOS becomes lower in energy with increasing doping (E red branch denoted by red arrows in Fig. 4(b) and (c)) while the other half becomes higher in energy (E blue branch denoted by blue arrows in Fig. 4(b) and (c)). Variation of twist angle shifts the E red and E blue branches uniformly up or down in energy but does not change their separation. So, for all angles there should be a charge density when the laser photon energy is resonant with a peak in the JDOS and an enhanced Raman G peak will occur. We use equation (1) to calculate the G peak area where the part ΔE(k) depends on charge density since the band structure is changing. The red curve in Fig. 4(a) is the theoretical calculation for the charge density when the Raman G peak area is a maximum considering only the interlayer potential imposed by the asymmetric doping in the two layers as in Fig. 4(d). Indeed, we find the maximum of the G peak area occurs at finite charge densities for both smaller and larger twist angles due to the resonant scattering pathways discussed above. Near the critical angle, the maximum G peak area occurs at zero charge density because of the large initial drop in the JDOS with doping which overwhelms the resonant effect. The charge densities calculated are symmetric about the critical angle and are much bigger than experimentally observed. Therefore, the interlayer potential alone can not explain the asymmetric angle dependence that is observed in the experiment.
To explain the results, there must be some mechanism that breaks the symmetry between the smaller and larger twist angles.
To explain the data shown in Fig. 3(a) and (b) for different twist angles, we must have a change in the band structure that depends on charge doping. In monolayer graphene, the Fermi velocity reduces upon charge doping due to electron-electron interactions [50][51][52][53] . Theoretical work for tBLG has also predicted a reduction 33 but it has not been experimentally explored. A reducing Fermi velocity causes the slope of the bands in Fig. 4(b) and (c) to decrease with doping, which results in a shortening of all scattering pathways. This causes both peaks in the JDOS to move toward lower energy as shown in Fig. 4(e). It prevents the E blue branch for smaller twist angles from increasing to match E ex . Likewise, it causes the E red branch for bigger twist angles to decrease faster and match E ex at lower charge densities. Therefore, the Raman G peak area for twist angles smaller than the critical angle drops monotonically within the charge density range in our experiment since no scattering pathway will be resonant with the excitation photon energy or the effect will be overwhelmed by the large initial drop in the JDOS with doping. On the other hand, the Raman G peak area for bigger twist angles increases to a maximum when E red matches E ex , overcoming the effect of the large initial drop in the JDOS with doping, then it drops at higher doping levels as E red becomes smaller than E ex . Here we employ changes in the Fermi velocity as a fitting parameter to explain the qualitative difference between the Raman results as a function of charge density for different twist angles. The blue curve in Fig. 4(a) is the theoretical calculation for the charge density when the maximum G peak area occurs considering both the interlayer potential and a reduction of Fermi velocity. We use a simple linear reduction of the Fermi velocity with total charge doping in the calculation as 6 m/s is the Fermi velocity of tBLG at zero charge density which matches the critical twist angle observed experimentally. The slope of the reduction, α, is adjusted to match the turning charge densities for the larger twist angles. A reduction of around 2 percent of the Fermi velocity at the largest charge density obtained in the experiment is needed to match our data. While the exact functional dependence of the Fermi velocity could be something other than a linear decrease with increasing charge density, we have chosen to use this as it is the simplest form and qualitatively explains our data. In particular, as shown by the red line in Fig. 4(a) a non-changing band structure gives qualitatively different results from our experiments and is symmetric about the critical angle. A change in the band structure which reduces the energy of the scattering pathways is necessary to explain our experimental results which are asymmetric about the critical angle.
To confirm our model for the reduction in Fermi velocity with increasing doping, we calculate the variation of the G peak area with charge density. This is done using equation (1) with a changing separation of the energy levels due to the charge density and a reduction of the Fermi velocity. Figure 3(c) and (d) plot numerical calculations of the Raman G peak area variation with doping for the smaller and larger twist angles corresponding to Fig. 3(a) and (b). The calculations agree very well with the experimentally measured curves. They exhibit the distinct difference in density dependence of the smaller and larger twist angles. Furthermore, they demonstrate that the reduction of Fermi velocity plays an important role in the Raman signal of tBLG with doping. Without the decrease of Fermi velocity with doping, the curves for smaller and larger twist angles would be identical.
In summary, modification of the electronic band structure of tBLG upon doping is realized using a polymer electrolyte top gate. Dirac points from the two layers are shifted apart by an interlayer potential Δφ because of the asymmetric doping causing the two vHs to be misaligned horizontally in momentum. Distinct behavior with charge density is observed for twist angles smaller and larger than the critical angle due to changes in the electronic band structure. From our modeling, a small reduction in the Fermi velocity with doping causes significant changes in the Raman G peak area. We hope that our results will stimulate further theoretical work into the modification of the band structure of twisted bilayer graphene with charge doping. Our results demonstrate the rich physical properties of tBLG and how they can be modified with doping paving the way for tailoring the electronic band structure of tBLG by external potentials.

Methods
Raman measurement setup. Raman spectroscopy was performed on a home built Raman system using a 532 nm Nd:YAG laser under ambient environment and room temperature. The laser was focused onto the sample by a 50× objective with NA = 0.5, which yielded a spot size smaller than 1.5 μm. The power of the laser was less than 1 mW in order to avoid damaging the sample. The reflected light passed through a 532 nm high-pass filter before being dispersed by a 600 lines/mm grating. The resulting spectrum was imaged on a thermoelectrically cooled CCD giving a spectral resolution of about 1 cm −1 .
Scientific REPORTS | 7: 7611 | DOI:10.1038/s41598-017-07580-3 Polymer electrolyte preparation. The polymer electrolyte was composed of lithium perchlorate and polyethylene oxide (PEO) with a weight ratio of 1 to 8 dispersed in methanol and stirred for several hours at 50 °C until thoroughly dissolved. A drop of the ion gel was put onto the graphene device and baked for 5 mins at 100 °C to evaporate the methanol.
Graphene sample synthesis. The graphene in this study was grown via a low pressure CVD method on 25 μm thick Cu foil placed within a sealed copper pouch 40 . The pouch was heated to 1040 °C and annealed for 1 hour under a H 2 gas flow rate of 4 sccm at a pressure of 60 mTorr. Then a CH 4 flow rate of 1.3 sccm was added to the system for 30 minutes to grow the graphene. The system was quickly cooled down to 350 °C with the same gas flows as during growth, and then the gas was stopped.