Introduction

With the discovery of the half integer quantum Hall effect in graphene1,2 and of topological materials,3,4 Berry phase effects have taken center stage in condensed matter research. In graphene and other crystals with a honeycomb structure, charge carriers have a sublattice and valley degree of freedom, described in a continuum Dirac model by a pseudospin.5,6 As a consequence of this spin-like property, electrons or holes belonging to the two inequivalent valleys (K and K′, see Fig. 1) acquire a Berry phase of π and −π, respectively during a cyclotron orbit. The importance of pseudospin and the Berry phase are most striking when perturbations are smooth on the atomic scale, i.e., sublattice symmetry still holds. In this case, scattering between the two valleys is suppressed and the pseudospin is conserved,5 leading to some important effects that are the hallmark of graphene, such as weak antilocalization,7 the half integer quantum Hall effect1,2 and Klein tunneling.5,8 Importantly, complete backscattering from a state |−q〉 to |q〉 (see Fig. 1) is forbidden, due to pseudospin conservation5 (q is the crystal momentum measured from the K point in the Brillouin zone). This was first shown and explained for metallic carbon nanotubes9,10 and is crucial for the exceptional mobility of graphene.5,11

Fig. 1
figure 1

Pseudospin and intravalley scattering in graphene. a Pseudospin texture around a K point in the first Brillouin zone of graphene. If the perturbation doesn’t distinguish between the two sublattices, it is unable to flip the pseudospin (σ). Therefore, complete backscattering involving small changes in momentum from |−q〉 to |q〉 is forbidden by pseudospin conservation. b Intravalley process responsible for the D′ Raman band of graphene

Full size image

Here we show that scattering on strain fluctuations in graphene can lift the restriction on complete backscattering. This is demonstrated by confocal Raman spectroscopy measurements of crumpled graphene, where we measure a giant increase in the D′ peak intensity. This Raman peak originates from a resonant Raman process which involves intravalley backscattering of charge carriers. The intravalley to intervalley scattering peak intensity ratio is found to be as high as ID′/ID = 30, in contrast to the usual value of ≈0.1.12 Since the strain induced pseudo-magnetic field (Bps) couples to the pseudospin,13,14 the enhancement of the D′ peak at 1620 cm−1 is due to the extra phase acquired by charge carriers undergoing Raman scattering on strain fluctuations. Thus, in contrast to a scalar potential, backscattering of Dirac particles is no longer forbidden. We reproduce our measurement results, using numerical calculation of the double resonant Raman processes.

Results and discussion

Mechanical deformations in 2D materials with a honeycomb atomic structure naturally give rise to a two component pseudogauge field, which is directly proportional to the strain tensor components.15,16,17,18,19 These strain induced fields have a scalar (V(r)) and a vector (A(r)) component,18 being analogous to an electrostatic potential and a magnetic vector potential. The latter having opposite sign in the two valleys and giving rise to a pseudo-magnetic field (Bps).18 For graphene supported on hexagonal BN, Bps is especially strong near bubbles (hundreds of Tesla) and has a major influence on transport properties.20 Furthermore, in the highest mobility heterostructure devices it is very likely that random strain fluctuations are the main factors limiting mobility through intravalley scattering.21 This type of scattering is characterized by small changes in the charge carrier momentum and dominates if the scattering potentials are smooth on the atomic scale, such as charged impurities,22,23,24 or strain fluctuations.21 Such scattering processes are mostly explored in charge transport experiments through weak (anti)localization measurements.7,10,25 In supported graphene, Bps appears due to random strain fluctuations, stemming from non perfect stacking and interaction with the substrate. For SiO2 supported graphene, Bps due to corrugation has values of the order of 1 T. The extra phase aquired by the wave function of the scattered carrier, much like a real magnetic field, suppresses the weak anti-localization effect.25

Confocal Raman spectroscopy is another powerful tool to investigate strain fluctuations in graphene. Mechanical deformation induced softening or hardening of the phonon mode energy is detectable through the shift, splitting and broadening of the G and 2D peaks.21,26,27 However, until now direct detection of the scattering on strain fluctuations has been lacking. Both small (intravalley) and large momentum (intervalley) scattering is measurable separately via the double resonant D′ and D peaks at ≈1620 and ≈1350 cm−1. Lattice defects produce both intravalley as well as intervalley scattering, giving contribution to both D′ and D peaks, smooth defects which are less efficient at producing large momentum change, mostly contribute to the D′ peak. Numerical calculations by Venezuela et al.28 show that for closely packed alkali metals on graphene, acting as smooth Coulomb scattering centers, the D′ peak should be more intense than the D peak. However, the peak itself should be undetectably small. Furthermore, for lattice defects, able to induce both types of scattering, the measured ratio of intensities is ID′/ID = 0.1.12 This is much smaller than the expected value of ~0.5 from analytical theory of the double resonant processes, suggesting that pseudospin effects could play an important role, particularly in intravalley processes.29

The double resonant D′ process is sketched in Fig. 1b. After the creation of the electron/hole pair the largest contribution to the D′ intensity is given by processes involving scattering by both an electron and a hole.28,30 Backscattering of the electron and/or hole involve an elastic defect scattering with small momentum change and a scattering involving an LO phonon along the \(\overline {{\mathrm{\Gamma M}}}\) direction.28,31 It has been shown previously that the process involves a small portion of phonon phase space, as well as relatively small regions of the Dirac cone.28,30 Since the biggest contribution to the intensity involves backscattering within a single valley along \(\overline {{\mathrm{\Gamma M}}}\), this process necessarily involves a flip in pseudospin. Indeed it is suspected by Rodriguez-Nieva et al.29 that the pseudospin related phases of the excited electron and hole play a dominant role in suppressing the backscattering necessary for the D′ peak. Since Coulomb scatterers do not change the phase of the charge carrier wave functions, the resulting suppression of backscattering straightforwardly explains the immeasurably small calculated intensity of the D′ peak.28 The situation changes drastically if the scattering potential has a vector component, i.e., there is a sizable pseudo-magnetic field involved in the defect scattering. If the charge carriers stay within the same valley, this changes the phase of the electron or hole, similarly to a real magnetic field, enabling backscattering.

To detect Raman scattering from strain fluctuations with a sizable Bps, we have measured confocal Raman maps on crumpled graphene flakes exfoliated onto a SiO2 surface. We have used to our advantage that during the exfoliation some flakes tend to crumple (see Fig. 2a). Crumpling enhances Bps caused by the strain fluctuations by at least a factor of 100 compared to the ~1 T25 resulting from surface roughness of SiO2. Additionally, we found the same results on samples that were crumpled using mild annealing and crumpling by a tungsten tip (see Methods). We have measured a total of ten crumpled graphene samples, among which six show an anomalously high and dispersive D′ peak (see Supplementary information S1). Here, we illustrate the effect through the example of a characteristic sample (Fig. 2), while data for other samples can be found in the Supplementary information.

Fig. 2
figure 2

Intravalley scattering in crumpled graphene. a Optical microscope image of an exfoliated graphene layer on SiO2. b AFM image of crumpled area. Position within flake shown by red rectangle in a. c Raman spectroscopy map of the D′ peak intensity in a crumpled graphene area in the position shown by the red rectangle. Scalebars in b, c 1 μm. d Raman spectrum showing a large D′ peak measured at the crumpled area shown in b, c by green circle, Insets: larger magnification AFM image of the crumpled area and zoom of the D peak region of the spectrum. Scalebar: 500 nm. e D′ peak of the area shown in d, measured with three different excitation wavelengths. Colors correspond to the respective laser color (red: 633 nm, green: 532 nm, blue: 488 nm). Black spectrum is measured using 532 nm excitation on the uncrumpled part of the flake, shown by the black circle in c. f Change in the peak position (Δω) of the 2D, G, and D′ peaks, as a function excitation energy (EL), for the area marked by the green circle in b–d. Measured dispersions (ΔωD′/ΔEL) of the experimental peaks are shown. Calculated D′ peak position as a function of laser energy is shown by orange triangles

Full size image

AFM measurements show that the graphene layer consists of a network of folds on it’s surface. Measuring the Raman spectra with a spatial resolution of ~500 nm and excitation wavelength of 532 nm, we map the intensity of the D′ peak (Fig. 2c) within the crumpled area marked in Fig. 2a. The most striking feature of the map is that the crumpled regions show increased D′ intensity, the region with the strongest enhancement marked by a green circle. To shed more light on the Raman scattering, we plot the complete Raman spectrum of this region in Fig. 2d. The most surprising feature of the spectrum is the extremely high D′ peak, having 20% the intensity of the G peak. This is unprecedented because the D′ peak is only observed when the sample contains strong lattice defects, such as vacancies, sp3 type defects of grain boundaries,12,28,32 and it is always accompanied by the D peak, with its intensity being mostly ~10% of the D intensity (ID′/ID = 0.1).12,33 This result has to be considered in the light that in the area where the spectrum was measured there are no graphene edges. As expected for defect free graphene, the D peak intensity is just barely larger than the background (see inset in Fig. 2d). Furthermore, the graphene flake shows the pristine Raman spectrum of graphene in the uncrumpled areas. Ratios ID′/ID = 10 have been also found in samples, where the laser spot contains both crumpled areas and edges (see Supplementary information S1). The ID′/ID maps of two additional samples can be seen in Fig. 3.

Fig. 3
figure 3

ID′/ID intensity maps of crumpled graphene. a, b Optical microscope images of two different crumpled graphene layers on SiO2. c Map of the ID′/ID ratio of the sample shown in a. d Map of the ID′/ID ratio of the sample shown in b. Maps are measured using 532 nm excitation. Additional data regarding the samples and further samples can be found in the Supplementary information S1. Black scalebars: 10 μm

Full size image

It is known that overlapping graphene layers can produce a non dispersive peak ranging from 1540 up to 1630 cm−1.34,35,36 The non dispersive nature of this peak is due to the fact that the LO phonon from around the \({\mathrm{\Gamma }}\) point taking part in the process is selected by the mismatch angle of the two graphene layers and not by the laser energy.35 Therefore, if the mismatch angle between the overlapping graphene layers is in the 4° to 6° range, this so called R′ peak could be mistaken for the D′ peak. To make sure that the peak around 1620 cm−1 is indeed the D′ peak, we measured the change in the peak position with changing excitation wavelength. The data for the presented sample can be seen in Fig. 2e, with the graph color corresponding to the excitation laser color. The measured dispersion is 7.8 cm−1/eV and is similar to the values measured on other samples (see Supplementary information S2). It is slightly lower than the value of ~10 cm−1/eV expected for the D′ peak,32,33 see also our calculation of the strain induced D′ peak dispersion in the same plot (see Fig. 2f). The reason for the lower dispersion value might be the fact that the D′ peak is mixed with the non dispersive R′ peak in certain areas. An example where this effect is strong is shown in Fig. S5 of the Supplementary information. To reduce the possible interference from the R′ peak in our measurements, we have prepared an additional set of samples. These graphene samples had their top and bottom covered in a 5 nm thick poly-vinyl alcohol (PVA) film. Crumpling was achieved by using a sharp tungsten needle and a micro-manipulator stage. We found similar D′ peak intensities in the crumpled areas as in the case of other samples (see Supplementary Figs. S6 and S7). These experiments rule out any major effect from the R′ peak, since the PVA layer on the graphene prevents any atomically clean interface to form between the overlapping graphene.

Within the literature there have been some observations of a barely measurable D′ peak on graphene supported on nanosized pillars37 and nanoparticles,38 exhibiting wrinkling. However, it is unclear what the D to D′ intensity ratio is in these experiments. In experiments of graphene wrinkling on a polymer, the D′ peak is obscured by the presence of the polymer substrate.39

Next, we reproduce our measured D′ peak intensities by numerical calculations of the double resonant processes.28,40 To explore the pseudo-magnetic field within crumpled graphene and to construct a theoretical model for the Raman scattering potential, we study the wrinkling of graphene, using the LAMMPS molecular dynamics code41 (for details see Supplementary information S4). We consider two model geometries that make up a crumpled graphene sheet.42 One involves a single fold (Fig. 4b), seen all over the crumpled sample (Fig. 2b), the other involves a double fold (Fig. 4c), constructed by creating a second fold in a singly folded graphene. The double fold shows a complex strain pattern, containing both compression and tension and an increase in the strain by a factor of 10 compared to the single fold, as evidenced by the color scale on the atoms in Fig. 4b, c.

Fig. 4
figure 4

Strain within graphene folds, deformations. a 10 × 10 nm graphene piece used in the molecular dynamics calculations.41 b Single graphene fold along the zigzag direction. c Doubly folded graphene. Color scale on atoms shows the local strain, calculated by taking the average of the relative nearest neighbor distance changes. Positive values (yellow) correspond to tensile strain. d–f Pseudo-magnetic field magnitude in the structures shown in b, c. d Pseudo-magnetic field generated by the σ-π overlap17 within the fold seen in b. e Bps magnitude due to the hopping modulation within the fold seen in (b). f Bps magnitude due to hopping modulation within the double fold seen in c

Full size image

Having determined the geometry of the crumpled graphene, it becomes possible to calculate the the pseudo-magnetic field Bps, from the atomic positions. It can be understood, as a result of the local modification of inter-atomic hopping, either through changes in bond length or bond angles.18 Additionally, Bps can be interpreted as the consequence of Dirac quasiparticles existing in a curved space-time.18,43,44,45,46 Here, we calculate Bps from the modulation of the hopping parameter, as described by Guinea et al.47 This model reproduces the magnitude of Bps measured via scanning tunneling microscopy.13,48,49 Computing Bps, directly from the modulation of the atomic positions50 we find a maximum Bps of 20 T in the single fold, and as expected from the difference in strain, the double fold shows a Bps in the 200 T range (Fig. 4e, f). This Bps originates from the modulation of the nearest neighbor hopping energy due to strain.15,16,47 However, in cases where the curvature of the graphene sheet is large, there is another sizable component to the pseudo-magnetic field, due to the hybridization of the σ and π bonds.17 The total vector potential is the sum of the hopping induced A(r) and the curvature induced Aσπ(r). The experimentally relevant part of the vector potential,13 the pseudo-magnetic field, is then formed by the rotor of the sum of these two contributions: \(B_{{\mathrm{ps}}} = \nabla \times ({\boldsymbol{A}} + {\boldsymbol{A}}^{\sigma \pi })\). Using the formula for Aσπ from Rainis et al.,51 we calculate the Bps for our fold. For a radius of curvature R of a fold parallel to the zigzag direction, we have \(A_x^{\sigma \pi }({\boldsymbol{r}}) = 3\varepsilon _{\pi \pi }a^2{\mathrm{/}}8R({\boldsymbol{r}})^2\) and \(A_y^{\sigma \pi }({\boldsymbol{r}}) = 0\), where \(\varepsilon _{\pi \pi } \approx 3{\kern 1pt} {\mathrm{eV}}\)17,51 and a = 1.42 Å. In our calculations R is around 3 Å, corresponding to the 2–3 Å measured by Annett et al.,52 while Rainis et al.51 calculate 7 Å. As an example, in the case when R = 4 Å the maximum curvature induced Bps is around 200 T, an order of magnitude larger than the hopping induced one (see Fig. 4d, e), in accordance with the finding of Rainis et al.51 Both contributions to Bps are also dependent on the orientation of the fold within the graphene lattice, with the maximum Bps present in folds parallel to the zigzag direction and zero Bps for armchair.51,53 For the double fold, the hopping induced Bps is in itself in the 200 T range, due to the larger strains (Fig. 4f). A large collection of double or multiple folds may be necessary to create the large Bps values needed to observe the enhanced D′ peak, such as in the region shown by a green circle in Fig. 2b, c. Using the field theoretical approach in curved space43,44,45,46 to calculate Bps, we find values of similar magnitude. However, in our experiment we lack the spatial resolution to be able to distinguish between the two possible interpretations of Bps.54 For further discussion see Section S6 of the Supplementary information.

We mention that here we assume Bps is not homogeneous on the scale of the magnetic length, thus Landau quantization48 due to strain is not expected. As an example, the magnetic length at 200 T is lB = 1.7 nm. This length scale is more than a factor of 3 larger than the size of areas with high magnetic field, see Fig. 4f, making the formation of Landau levels impossible. Even though Landau levels are not expected to form, it has been shown that folded graphene areas with rapidly changing Bps can host bound states of Dirac fermions, as shown experimentally55 and theoretically.56,57

Having quantified the Bps magnitude in folded, crumpled graphene, we turn our attention to numerically calculating the double resonant Raman processes in the presence of strain. Calculations are performed similarly to Venezuela et al.28 and Kürti et al.40 For this we need to find the electron/hole—defect Hamiltonian: Hdef. We model the Bps seen in crumpled graphene by a simple model, using a Gaussian deformation of the form: \(h \cdot {\mathrm{exp}}\left( { - \frac{{x^2 \,+\, y^2}}{{b^2}}} \right)\), where h is the height and b is the width of the Gaussian. For a Gaussian, the analytical form of the vector potential \(|{\boldsymbol{A}}| \propto h^2/b^2\) is well known.58 By choosing the height h and width b of the bump to be 5 and 20 Å respectively, the Bps within the bump area reaches a maximum of 200 T, switching sign on the nanometer scale (see Fig. 5a). Both the magnitude and the length scale is similar to the values found within the LAMMPS calculations. In atomic units, the Hamiltonian describing the defect is \(H_{{\mathrm{def}}} = {\boldsymbol{A}}({\boldsymbol{r}}) \cdot \nabla + V({\boldsymbol{r}})\), where A is the vector potential and we also include V, the scalar potential generated by the mechanical deformation.18 In order to calculate the scattering matrix element \(\langle {\boldsymbol{k}}|H_{{\mathrm{def}}}|{\boldsymbol{k}} \pm {\boldsymbol{q}}\rangle\) between two states with wave vectors k and k ± q we need to calculate A( ± q) and V( ± q):28

$$\begin{array}{*{20}{l}} {V( \pm q) = \mathop {\sum}\limits_l {\kern 1pt} e^{ \pm i{\boldsymbol{qr}}_l}V({\boldsymbol{r}}_l)} \hfill \\ {{\boldsymbol{A}}( \pm q) = \mathop {\sum}\limits_l {\kern 1pt} e^{ \pm i{\boldsymbol{qr}}_l}{\boldsymbol{A}}({\boldsymbol{r}}_l)} \hfill \end{array}$$
(1)
Fig. 5
figure 5

Calculated resonant Raman spectra. a Model used in Raman calculations. Gaussian of shape: \(h \cdot {\mathrm{exp}}\left( { - \frac{{x^2 \,+\, y^2}}{{b^2}}} \right)\), where h is the height and b is the halfwidth of the deformation (inset). Color plot shows the calculated pseudo-magnetic field pattern for a deformation having b = 2 nm and h = 0.5 nm. b Plot of |A(k)|2 in the first Brillouin zone (white hexagon). c Calculated D and D′ peak intensity for a Gaussian deformation seen in a. d Ratio of D′ and D intensities (ID′/ID) as a function of excitation energy and Gaussian width b. The Gaussian amplitude is h = b/4, in order to keep the magnitude of the mechanical strain within the Gaussian constant, only changing the spatial extent of the deformation. Color scale is logarithmic

Full size image

Figure 5b shows |A(k)|2 in the first Brillouin zone of graphene for varying sizes of the deformation. While varying the width, the height to width ratio (h/b) was kept constant at 0.25, in order to keep the maximum value of |A(k)| constant.58

The calculated D, D′, and 2D peaks for a Gaussian with h = 5 Å and b = 20 Å is shown in Fig. 5c. We reproduce ID′/ID = 32 in accordance with the experimental spectrum of Fig. 2d. As a crosscheck to our calculations we also compute the D and D′ intensity for a “hopping” defect, as implemented by Venezuela et al.,28 modeling a lattice defect in graphene. For the hopping defect, the D′ peak intensity is 11% of the D intensity, as expected for lattice defects, such as vacancies.12 Comparing the relative intensity with respect to the 2D peak, it is clear that the D′ peak is showing a measurable intensity, as opposed to the D′ peak generated by Coulomb scatterers.28 The D′ intensity is generated almost completely by the vector potential A, with the scalar potential giving only a slight (less then 10−4) contribution (see Supplementary information S5). This can be easily explained if we consider the scattering matrix element between states \(|\langle {\boldsymbol{k}}|H_{{\mathrm{def}}}|{\boldsymbol{k}} \pm {\boldsymbol{q}}\rangle |^2 = |H_{{\mathrm{def}}}({\boldsymbol{q}})|^2{\mathrm{cos}}^2(\theta _{{\boldsymbol{k}},{\boldsymbol{k}} \pm {\boldsymbol{q}}}{\mathrm{/}}2)\), where θk,k ± q is the angle between the the initial and the scattered state.9 For backscattering (θk,k ± q = π), if Hdef is of purely scalar character, the matrix element is zero. However, if Hdef has a vector potential component it can be nonzero. This is because A acts on the pseudospin of graphene, as evidenced by the pseudo-Zeeman effect.13 Thus, it is able to change the phase of the charge carriers.

The vector potential A has a slight contribution to the D peak as well, as can be seen in Fig. 5c. This is due to the non zero scattering potential at large k values, around the K points. Due to the Fourier transformation in Eq. (1), the contribution at large k increases as we make the Gaussian narrower. By decreasing b, but keeping the aspect ratio the same, we can see a marked increase in the scattering potential |A(k)|2 at the K points, increasing the D peak intensity. This can be observed in Fig. 5d, where we plot the ID′/ID ratio on a log scale, as a function of b and excitation energy. Although it is clear that the parameter b is a tuning knob for the ID′/ID ratio, it is not possible to attribute a certain b value to the measured ratio, because the crumpling patterns can have a more complex Bps pattern than what is assumed in the Gaussian model. Additionally, the calculated model does not include lattice defects, which could be present in the experiment, such as edges. For the map seen in Fig. 2b, the D peak intensity is too low to prepare a meaningful ID′/ID map. For other samples, the map of this ratio can be found in the Supplementary information.

Our experiments show a first example of resonant Raman scattering on potentials that are smooth on the atomic scale. In the lack of intervalley scattering, the pseudo-magnetic field created by the strain induced vector potential shifts the phase of the electron or hole wave function, in addition to the Berry phase. This lifts the restriction on backscattering, enabling an enhancement in the intravalley backscattering rate by orders of magnitude, evidenced by the enhancement of the D′ peak in the Raman spectrum. The drastic increase in intravalley scattering, as a consequence of strain fluctuations, may be an important factor in lowering the mobility of certain CVD grown graphene samples, especially if wrinkles and folds are introduced during the transfer process. Furthermore, it is well known that defects such as grain boundaries, vacancies, etc. can have a specific ID′/ID fingerprint.12 Such defects also distort the graphene lattice around them59 and this local strain field can be specific to the type of defect. Our results show that if we want to understand the origin of this Raman fingerprint, scattering on the defect induced local strain fields has to be taken into account. If experimentally, a more controlled folding of graphene can be achieved, the Bps present in these folds could be used to steer and guide55 valley polarized60,61 Dirac fermions.

Methods

Graphene samples were prepared by micromechanical exfoliation from natural graphite, purchased from NGS Trading & Consulting GmbH. Raman measurements were carried out using a Witec 300rsa+ confocal Raman system, using 488, 532, and 633 nm excitation lasers. Laser power was kept at 0.5 mW for all lasers.

Crumpling of the graphene layers has been induced either in the exfoliation stage, as with the sample in the main text and the sample shown in Fig. S4 of the Supplementary information. However, this method has a low yield. It is known that annealing can induce folding of graphene.52 Therefore, to increase the yield of the crumpled graphene within our samples, we have used annealing to induce wrinkling and crumpling. Samples 1–4 have been prepared by placing them onto a hotplate, preheated to 160 °C for 15 min, after which the samples are removed and placed onto a metal surface to allow quick cooling. Additionally, the crumpling can be also induced by the tip of a tungsten needle.