This is the submitted version of the article: Vincent T., Panchal V., Booth T., Power S.R., Jauho A.-P., Antonov V., Kazakova O.. Probing the nanoscale origin of strain and doping in graphene-hBN heterostructures. 2D Materials, (2019). 6. 015022: - . 10.1088/2053-1583/aaf1dc. Available at: https://dx.doi.org/10.1088/2053-1583/aaf1dc Probing the nanoscale origin of strain and doping in graphene-hBN heterostructures Tom Vincent1, 2, Vishal Panchal1i, Tim Booth3, Stephen R. Power4, 5, 6, Antti-Pekka Jauho3, Vladimir Antonov2, 7 & Olga Kazakova1 1. National Physical Laboratory, Teddington, TW11 0LW, United Kingdom; 2. Royal Holloway, University of London, Egham TW20 0EX, United Kingdom; 3. Center for Nanostructured Graphene (CNG), DTU Nanotech, Department of Micro- and Nanotechnology, Technical University of Denmark, DK-2800 Kongens Lyngby, Denmark; 4. Catalan Institute of Nanoscience and Nanotechnology (ICN2), CSIC and The Barcelona; Institute of Science and Technology, Campus UAB, 08193 Bellaterra, Catalunya, Spain; 5. Universitat Autònoma de Barcelona, 08193 Bellaterra, Catalonia, Spain; 6. School of Physics, Trinity College Dublin, Dublin 2, Ireland. 7. Skolkovo Institute of Science and Technology, Nobel str. 3, Moscow, 143026, Russia Email: olga.kazakova@npl.co.uk Abstract We use confocal Raman microscopy and modified vector analysis methods to investigate the nanoscale origin of strain and carrier concentration in exfoliated graphene-hexagonal boron nitride (hBN) heterostructures on silicon dioxide (SiO2). Two types of heterostructures are studied: graphene on SiO2 partially coved by hBN, and graphene fully encapsulated between two hBN flakes. We extend the vector analysis methods to produce spatial maps of the strain and doping variation across the heterostructures. This allows us to visualise and directly quantify the much-speculated effect of the environment on carrier concentration as well as strain in graphene. Moreover, we demonstrate that variations in strain and carrier concentration in graphene arise from nanoscale features of the heterostructures such as fractures, folds and bubbles trapped between layers. For bubbles in hBN-encapsulated graphene, hydrostatic strain is shown to be greatest at bubble centres, whereas the maximum of carrier concentration is localised at bubble edges. Raman spectroscopy is shown to be a non-invasive tool for probing strain and doping in graphene, which could prove useful for engineering of two-dimensional devices. Introduction Two-dimensional (2D) materials, such as graphene, demonstrate a great potential for device fabrication due to their high mobility, extremely low thickness, high strength and flexibility.[1] The ability to stack different 2D materials into van der Waals heterostructures with novel properties creates further opportunities for engineering devices with custom-tailored properties.[2] Graphene sheets encapsulated in hexagonal boron nitride (hBN) have been shown to have a strongly enhanced carrier mobility compared to bare graphene and are also protected from atmospheric adsorbants.[2] Precise knowledge of strain and doping in graphene and its heterostructures is crucial for tailored i Formerly. Now at: Bruker Nano Surfaces UK, Coventry, CV4 9GH, United Kingdom. device performance. This is often complicated at the nanoscale, where graphene strain and doping variations alter carrier mobility[3], Fermi level[4] and optoelectronic properties[5]. It has also been shown that non-uniform strain can induce strong pseudomagnetic fields greater than 300 T[6], which may in turn provide a platform to manipulate the sublattice[7–9] and valley[10,11] degrees of freedom. The sensitivity of graphene’s physical properties to strain opens up the possibility of using deliberately induced strain as a method of controlling various parameters in graphene devices, this has been referred to as “straintronics”.[12] However, strain can also occur as an uncontrolled artefact of many of the processes involved in fabricating graphene devices, including deposition on a substrate[3,13], assembly into van der Waals heterostructures[14] and thermal annealing.[15] For this reason, quantitative methods for determining strain variation are vital. Raman spectroscopy has proved a useful tool for studying strain and doping in graphene[13,16–18]. Its characteristic Raman peaks (G and 2D) are affected by both strain and doping, which means that the basic Raman analysis does not allow a straightforward separation of intertwined strain and charge effects. Lee et al. first demonstrated that the effects of strain (𝜀) and hole doping (𝑛) can be optically separated from each other by correlation analysis, enabling their quantification.[16] Here, we adapt and further develop the aforementioned approach to the study of two graphene-hBN heterostructures on SiO2: A. A simple pristine heterostructure consisting of a sheet of single layer graphene (SLG) partially covered by a flake of multilayer hBN; B. A heterostructure formed of a sheet of SLG encapsulated between two multilayer hBN flakes, in which bubble-like structures formed during fabrication. Moreover, we extend the strain-doping analysis to produce spatial maps of the strain and doping variation. Using this method, we unambiguously link nanoscale variations in strain and doping to local features and defects in the heterostructures, such as fractures, folds, bubbles and edges. Theory Quantifying strain and doping from graphene’s Raman spectra is complicated. The Raman shifts of the G and 2D peaks, 𝜔𝐺 and 𝜔2𝐷 , depend on both 𝜀 and 𝑛.[16] If only one peak is considered, it is therefore not possible to determine 𝜀 or 𝑛, unless the other quantity is known beforehand[16–18]. Using a number of prior assumptions and experimental data sets, Lee et al. proposed a method[16] that uses correlation analysis of both 𝜔𝐺 and 𝜔2𝐷 to separately determine 𝜀 and 𝑛, without prior knowledge of either, which we discuss below. It has been shown both experimentally[19–25], and theoretically[26–29], that for SLG with a constant 𝜀 and varying 𝑛, or vice versa, the mutual position of the G and 2D peaks for a given Raman spectrum in 𝜔𝐺 -𝜔2𝐷 space will approximate a straight line. Figure 1a shows a schematic illustration of this. In these coordinates, the relative distance between an experimental point and the point corresponding to pristine (i.e. unstrained and undoped) graphene can be decomposed into two vector components, 𝒗𝜺 and 𝒗𝒏 , each parallel to one of the two straight lines. Here, 𝒗𝜺 is the shift due solely to strain and 𝒗𝒏 is the shift due solely to hole doping. The blue circle in Figure 1a shows an example of a pair of 𝜔𝐺 and 𝜔2𝐷 values, representing a typical Raman spectrum. The distance from the point corresponding to pristine graphene (charge-neutral and unstrained), shown as a red circle, is decomposed into vectors 𝒗𝜺 and 𝒗𝒏 , shown as black arrows. Once these vectors have been obtained, the values of 𝜀 and 𝑛 can be found by comparing the shifts with known reference values.[16,18] In general, the method is only valid for the case of undoped or p-type graphene. Although both n- and p-type doping cause an increase in the value of 𝜔𝐺 , for n-doped graphene the point (𝜔𝐺 , 𝜔2𝐷 ) does not move in a straight line[16], making the vector analyses overcomplicated and unreliable. The curved line means that it is not possible to use vector decomposition to unambiguously determine strain and electron doping. The straight and curved trajectories shown by p- and n-doped graphene, respectively, are also shown in Figure 1a. Figure 1. (a) A schematic representation of the vector decomposition method[16] for separating the effects of strain and doping based on Raman peak shift in graphene. The red circle represents pristine graphene (i.e. unsupported and charge neutral), the blue circle is an example of an experimental point. The red and blue lines show experimental trajectories for p-type and n-type graphene, respectively[16]. (b, c) The resulting values of hydrostatic strain and doping, respectively, for allowed combinations of G and 2D peak positions according to the Mueller method[18] of strain-doping decomposition. Values above the line of 𝑛=0 are not treated as valid by the model, and so are left blank. In both figures, the dashed and solid black lines represent 𝜀=0 and 𝑛=0, respectively. The point where they cross corresponds to pristine, undoped and unstrained graphene.[18] Other prerequisites for a successful analysis include graphene with a low defect density[16], because high defect densities can also cause peak shifts in graphene which would interfere with the separation of strain and doping contributions.[30] For nominally undoped graphene, it is usual to observe a scatter in the G peak frequency of ±1 cm-1. This means that the model may return negative values of 𝑛 for particular Raman measurements, while the mean doping level of many individual spectra remains at zero. However, by using confocal Raman and keeping track of the location where each spectrum was taken to produce spatial maps of 𝜀 and 𝑛, we show that this scatter is not random, but instead correlated to nanoscale features of the heterostructures under investigation. Lee’s analysis is further complicated by the fact that different types of strain produce different slopes for 𝒗𝜺 [16,18,21–23,26,27], originating from biaxial and uniaxial strain as well as the coexistence of both strain components. Additionally in the case of uniaxial strain, the slope also varies with the strain orientation[16,18]. If the type of strain is unknown, uncontrolled, or if there is a mixture of different types of strain, as in the heterostructures studied in this work, Lee’s method is no longer valid. This is because an unknown error in the slope of 𝒗𝜺 will lead to an unknown error in the resulting value of 𝜀. Mueller et al. proposed a modification to Lee’s method that allows for arbitrary strain configurations to be determined along with doping.[18] Arbitrary strains can be written in terms of the components of a biaxial strain tensor: 𝜀𝑥𝑥 (𝒓) 𝜀𝑥𝑦 (𝒓) 𝜀(𝒓) = ( ), (1) 𝜀𝑦𝑥 (𝒓) 𝜀𝑦𝑦 (𝒓) then further decomposed into a hydrostatic component, defined as: 𝜀ℎ = 𝜀𝑥𝑥 + 𝜀𝑦𝑦 , (2) and a shear component, defined as: 2 𝜀𝑠 = √(𝜀𝑥𝑥 − 𝜀𝑦𝑦 ) + 4𝜀𝑥𝑦 2 , (assuming 𝜀 𝑥𝑦 = 𝜀𝑦𝑥 ). (3) Using Mueller’s terminology, hydrostatic strain refers to an isotropic expansion of the lattice, which causes a shift in 𝜔𝐺 and 𝜔2𝐷 , and shear strain refers to a change in the shape of the lattice, which leaves the area of a unit cell unchanged and results in a splitting of the G and 2D peaks.[18] . It is worth noting that this definition of hydrostatic strain yields values twice as high as the corresponding biaxial strain (𝜀𝑏𝑖𝑎𝑥𝑖𝑎𝑙 = 𝜀𝑥𝑥 = 𝜀𝑦𝑦 ). By choosing the slope of 𝒗𝜺 corresponding to changing 𝜀ℎ with a constant 𝑛, the hydrostatic strain component can be found from the Raman shift of the G and 2D peaks. To obtain a single value of Raman shift from peaks that may be split into two by shear strain, the mean centre points of the split peaks are used, 𝜔 ̅̅̅̅ 𝐺 and ̅̅̅̅̅. 𝜔2𝐷 Then the peak shift contributions due to 𝑛 and 𝜀ℎ can be separated using vector decomposition as described in the Lee model. The value of 𝑛 can be found by comparing the relevant peak shift components to known reference values. The value of 𝜀ℎ can be found from the G peak shift due to strain according to: Δ𝜔𝐺ℎ = −𝜔𝐺0 𝛾𝐺 𝜀ℎ , (4) where Δ𝜔𝐺ℎ is the G shift due solely to 𝜀ℎ , 𝜔𝐺0 is the unstrained G frequency and 𝛾𝐺 is the Grueneisen parameter for the G peak, which quantifies how much the peak shifts with strain.[18,26] If the splitting due to 𝜀𝑠 is large enough to be measured, 𝜀𝑠 can be determined according to: Δ𝜔𝐺 = 𝜔 ̅̅̅̅𝛽 𝐺 𝐺 𝜀𝑠 , (5) where Δ𝜔𝐺 is the splitting of the G peak and 𝛽𝐺 is the shear deformation potential for the G mode[18,31]. In practice, Δ𝜔𝐺 may be too small compared to the peak width to allow the split peaks to be resolved. To calculate 𝜀 and 𝑛, we follow the approach used by Mueller[18] and assume 𝜔𝐺 =1583 cm-1 and 𝜔2𝐷 =2678 cm-1 for the point corresponding to pristine graphene (for a laser excitation wavelength of 532 nm), 2.21 and 0.55, for the gradients of 𝒗𝜺 and 𝒗𝒏 , respectively, and a Grueneisen parameter for the G mode of 𝛾𝐺 =1.8. The specific reference values used to calculate 𝑛 were taken from Froehlicher and Berciaud’s work, backgating a monolayer graphene device to control its carrier density.[25] Heat-maps showing how the resulting values of 𝜀 and 𝑛 vary across 𝜔𝐺 -𝜔2𝐷 space are given in Figures 1b and 1c. Results Heterostructure A consists of a sheet of single layer graphene, ~70 μm in size, partially covered by exfoliated hBN of around 20 nm in thickness and with lateral dimensions of ~50×30 μm. Atomic force microscopy (AFM) scans (Figures 2a and b) reveal one dimensional (1D) features on both the graphene and hBN covered graphene. These appear to be a mixture of small fractures and folds, which most likely occurred during the exfoliation or stacking of the heterostructure materials. Figure 2b shows a close-up of a fracture and a fold. Whereas the fracture has homogeneous width along the whole length, the fold starts tightly compacted where it contacts the fracture but spreads out as it gets further away. Detailed line profiles taken from the features in Figure 2b are plotted in Figures 2c and 2d. The bottom of the fracture in Figure 2c is ~0.6 nm lower than the average level of the graphene. This suggests a monolayer nature of the graphene. Figure 2. Pristine graphene-hBN heterostructure A. (a) An AFM scan of the overall heterostructure. The image shows 1D features on the surface of the hBN. (b) A higher resolution AFM scan, taken from the region indicated by the white square in 2a and indicating that 1D features include both fractures and folds in the graphene. The red rectangles show where line the profiles shown in 2c and 2d were taken. (c) A line profile taken across the fracture as shown in 2b. The dashed line indicates the level of the surrounding graphene layer. The inset shows a close up of the fracture itself. (d) A line profile taken across the fold as shown in 2b. (e) A map of the Raman 𝐼𝐷 /𝐼𝐺 ratio shows higher defect density along the fractures in graphene. Some fractures exist beneath the hBN and correspond to surface features visible in the AFM. Black regions correspond to areas with no graphene. (f, g) Histograms of the Raman 𝐼2𝐷 /𝐼𝐺 ratio and the width of the 2D peak, 𝛤2𝐷 . The blue and red areas in the histograms correspond to exposed and hBN-covered graphene, respectively. A Raman map of the heterostructure was taken, and Lorentzians were fitted to the graphene 2D, G and D peaks, as well as the characteristic hBN peak at around 1367 cm-1 for every spectrum. The fitted intensities of the 2D and hBN peaks were then used to determine points corresponding to graphene and hBN areas. Figure 2e shows a map of the D to G peak intensity ratio, 𝐼𝐷 /𝐼𝐺 , which is used as an indicator of lattice defect density in graphene. The near-zero values on most of the flake indicate a very low defect density.[30] The areas previously identified as fractures show a high defect density, which is to be expected if the lattice is broken at these points. Some features, which are seen on top of the hBN in the AFM, correlate to fractures in the graphene, which propagate from the hBN covered area to the area of bare graphene. To better verify the thickness and quality of the graphene, we plotted histograms of the 2D to G peak intensity ratio, 𝐼2𝐷 /𝐼𝐺 , and the width of the 2D peak, 𝛤2𝐷 , in Figures 2f and 2g. Both parameters are altered by the presence of hBN on top of graphene, so bare graphene on SiO2 and hBN covered graphene were plotted separately, in blue and red, respectively. For the bare graphene, the values of 𝐼2𝐷 /𝐼𝐺 cluster around 2.2, and the values of 𝛤2𝐷 around 27 cm-1, which are typical values for SLG on SiO2.[16] We can assume that the part of the same graphene flake that is covered by hBN is of the same intrinsic quality. A scatter plot demonstrating the distribution of 𝜔𝐺 and 𝜔2𝐷 across heterostructure A is shown in Figure 3a. The axes corresponding to 𝑛=0 and 𝜀ℎ =0 are added for clarity (see also Figure 1a). By colouring the points corresponding to bare graphene and hBN/graphene differently, it is clear that the presence of the hBN divides the scatter plot into two main populations: the hBN covered points are shifted in the direction corresponding to more compressive strain and decreased doping. Within both separate populations, there are well defined clusters spread out parallel to the hydrostatic strain direction and less well defined clusters aligned with the doping direction. This indicates that there is a spread of values of both strain and doping in the heterostructure, but that most of the variation in Raman shift is caused by the strain. Figure 3. Vector decomposition analysis applied to heterostructure A. (a) A scatter plot showing the distribution of G and 2D positions for the heterostructure. The blue and red points correspond to Raman spectra taken from bare and hBN- covered graphene on SiO2, respectively. The solid and dashed lines indicate 𝑛=0 and 𝜀ℎ =0, respectively. The point where they cross corresponds to pristine, undoped and unstrained graphene. (b) A map of the hydrostatic strain distribution across the heterostructure. (c) A map of the doping variation across the heterostructure. In both Raman maps, black regions correspond to areas with no graphene. The arrows indicate the region of folded graphene detailed in Figure 2b. A map of the hydrostatic strain calculated for each pixel in the Raman map is shown in Figure 3b. The graphene beneath the hBN is significantly more compressively strained, at ~-0.1% strain, than the uncovered graphene, which shows tensile strain, at ~0.5%. This is to be expected, given that the heterostructure was stacked at a high temperature and then allowed to cool to ambient temperature. The different thermal expansion coefficients of hBN and graphene cause them to contract differently, effectively putting a strain on both materials.[15] At the most prominent fractures, the strain is relaxed to nearly zero, which, again, is to be expected, as fracturing by its nature is a mechanism that takes a system from a state of high strain to a relaxed state. At the upper right edge of the graphene, there is a region with less tensile strain than the rest of the bare graphene. The close proximity to the large fracture is likely to cause relaxation of a significant area of the graphene flake nearby. On the other hand, at the graphene fold (shown in Figure 2b) and the area of graphene spreading out from the fold, there is a higher degree of tensile strain than in the surrounding area of graphene. It should be noted that while the area of high tensile strain is directly localised at the fold, its shape does not coincide with the fold, i.e. it protrudes significantly wider than the topological dimensions of the fold. Thus, it is noteworthy that these important variations in the physical properties (i.e. strain) of the material are clearly observed using the Raman vector decomposition method whilst they are only partly visible in the AFM topography. Figure 3c shows the corresponding map for the hole concentration of heterostructure A. The averaged doping of graphene protected by the hBN layer is ~1×1012 cm-2, which is significantly lower than the value of ~3×1012 for the exposed graphene. This clearly demonstrates the protective properties of the hBN layer, which shields the underlying graphene from environmental adsorbates, the primary source of charge carriers (primarily p-type) in exfoliated graphene. Interestingly, the same decrease in the carrier concentration is observed both at folds (and the associated spread areas) and fractures in graphene. The same reduction of doping along 1D features is seen both in hBN covered graphene and bare graphene, indicating that in this case the effect is not due to the environmental doping but has an intrinsic nature (e.g. graphene substrate interaction). Doping is also reduced at those areas where a fold spreads out into the surrounding graphene, this is most likely due to limited adhesion to the substrate, preventing charge transfer from SiO2 to graphene. We further discuss the results obtained using heterostructure B, which consists of a single layer of graphene encapsulated between two flakes of hBN, both of ~20 nm in thickness. Figure 4a shows an optical image of the heterostructure with annotations explaining its structure. During the transfer process, bubbles were introduced between the layers, most likely due to stacking at an insufficiently high temperature.[32] Many of these bubbles are large enough to be clearly observed both optically and using AFM techniques, as shown in Figure 4e, even through the 20-nm-thick hBN on the top of the heterostructure. Again, a Raman map of the heterostructure was taken, and Lorentzians were fitted to the graphene 2D, G and D peaks, and the characteristic hBN peak. The fitted intensities of the 2D and hBN peaks were used to determine the points corresponding to graphene and hBN. In the measured part of this heterostructure, the whole graphene flake is covered or encapsulated by hBN. The map of the 𝐼𝐷 /𝐼𝐺 ratio in Figure 4d shows a low average defect density. There are high 𝐼𝐷 /𝐼𝐺 values at the edge of the graphene sheet inside the hBN, which is to be expected. There are also 1D features, which spread from the upper right of the flake (invisible in optical and AFM images), where graphene is in contact with the SiO2 substrate (Figure 4a). These can be attributed to fractures introduced during the stacking of the heterostructure. Overall, the bubbles are not visible in the defect map, which indicates that the graphene structure remains intact within the bubbles. It is noteworthy, however, that the defective area in the lower right part of the heterostructure corresponding to several bubbles was introduced by excessive laser heating in previous extensive experiments. Figure 4. Heterostructure B (hBN encapsulated graphene). (a) An optical image of the whole heterostructure. The lower hBN flake is outlined in blue, the upper one in red, and the encapsulated graphene in yellow. The dark spots in the centre of the heterostructure are contamination on the surface. The solid white rectangle indicates the region from which the Raman maps in this work were taken. (b, c) Histograms showing the distributions of 𝐼2𝐷 /𝐼𝐺 , and 𝛤2𝐷 . (d) A map of the Raman 𝐼𝐷 /𝐼𝐺 ratio shows higher defect density at the edge of the encapsulated flake. Two bubbles in the lower right part of the heterostructure show a high defect density caused by excessive laser heating in previous Raman experiments. (e) An AFM scan of the lower part of the heterostructure. The red bar indicates a bubble from which line profiles were extracted (shown in Figure 5d). The histograms in Figures 4b and 4c show values of I2D/IG and Γ2D clustered at ~3.2 and ~22.5 cm-1, respectively. These compare favourably to the values from hBN covered graphene from heterostructure A, which confirms that the graphene is single layer, and that the prerequisites for the model to be valid are met. A scatter plot showing the distribution of 𝜔𝐺 and 𝜔2𝐷 across heterostructure B is presented in Figure 5c. The axes corresponding to 𝑛=0 and 𝜀ℎ =0 are added again to help visualise changes due to doping and strain. The points on the scatter plot are much more densely packed than they were for heterostructure A, implying a smaller overall variation in both doping and strain. This is a reflection of the fact that the majority of the graphene flake is fully encapsulated, as well as the smaller size of the area under inspection. The average carrier concentration is close to zero, as can be seen from the fact that the points are clustered around the line of 𝑛=0, indicating that graphene is close to charge neutrality, as would be theoretically expected for fully encapsulated graphene, i.e. in the absence of charge transfer both from the environment and substrate. As such, we observe the expected scatter of ±1 cm-1. The maps of the separated doping and hydrostatic strain across the heterostructure are presented in Figures 5a and 5b. The native strain in the encapsulated graphene is slightly more compressive compared to freestanding graphene, with the strain varying from ~-0.06% to -0.03%. Within the bubbles, the graphene becomes nearly strain-free or slightly tensile strained. The damaged bubbles in the lower right part of the heterostructure show a tensile strain of ~0.06%. In all cases, the strain value peaks at the bubbles’ centres, i.e. the strain distribution has a simple dome shape. This is to be expected for bubbles in 2D materials[23,33]. At the 1D features shown in the defect density map in Figure 4d, there may be a small relaxation of the compressive strain, however the strain variations that correlate to bubbles make this difficult to judge. At the edges of the graphene flake, there is a higher tensile strain of ~0.06%, except for the small area of graphene in the top right corner, which is located directly on SiO2 rather than on the lower hBN flake. It is interesting to note that there are a relatively large number of small bubbles clear from the strain map that are not visible in the optical image or AFM (Figures 5a, 4a, and 4e, respectively). Figure 5. Vector decomposition analysis applied to heterostructure B. (a) A map of the hydrostatic strain distribution across the heterostructure. (b) A map of the doping variation across the heterostructure. In both Raman maps, black regions correspond to areas with no graphene. (c) A scatter plot showing the distribution of G and 2D peak positions for the heterostructure. The solid and dashed lines indicate 𝑛=0 and 𝜀ℎ =0, respectively. The point where they cross corresponds to pristine, undoped and unstrained graphene. (d) AFM, hydrostatic strain and doping profiles taken across the example bubble indicated by white bars in 5a and 5b, as well as by the red bar in Figure 4e. The encapsulated graphene is nearly charge-neutral, which is to be expected as hBN shields graphene not just from environmental adsorbates, but also from doping from the substrate. Due to the aforementioned scatter about the point of zero doping, some points are returned by the model with negative values of 𝑛. These negative values are difficult to interpret. They might be an indication that the point corresponding to pristine graphene (i.e. the point of 𝑛=0), may need to be defined more precisely in order to quantify such small variations in doping. The origin point was taken from literature, and was calculated from an experiment in which graphene was back-gated to control its carrier density while Raman measurements were performed[18,25]. It is also possible that a small native strain in the graphene in Ref. 25 may have caused a slight offset to the position of the origin point, which becomes significant at such low doping levels. However, Figure 5b shows that rather than being random, the observed scatter is correlated to nanoscale features of the heterostructure. The ability to resolve such features in the scatter clearly demonstrates the capability of the confocal Raman method to probe carrier concentration variation over nanoscale domains. The bubbles are a feature that have a clear effect on the carrier concentration, however the profile of the charge distribution across a bubble looks very different to that for strain. The charge distribution shows a crater like shape with rings of positive charge forming along the bubble edge. Line profiles, taken across a bubble in the strain and doping maps together with an AFM height profile, are shown in Figure 5d. Unlike the strain, which peaks at the centre of the bubble, where the AFM also shows the greatest height, the charge concentration is greatest on the sides of the bubble, where the steepest slope in the bubble wall is observed. It is worth noting that bubbles in graphene are buried under the ~20 nm-thick top layer of hBN, which masks the exact topography of the bubbles as measured by AFM. It is worth considering the possible mechanisms that can lead to such charge density fluctuations in the bubble structures. It has been predicted that nanobubbles can induce pseudomagnetic fields in graphene[6], which can lead to the emergence of pseudo-Landau levels and cause charge density variations within the bubbles.[34–37] However, fitting the height and strain profiles of the bubbles in this work to simple membrane models[34,38] suggests that the resulting pseudomagnetic fields are too weak (~10 mT) to induce the strong spatial fluctuations seen in the experiment. In the absence of pseudomagnetic effects, a sharp interface at the edge of the bubble can induce electronic standing waves within the bubble[34,35], similar to a quantum corral.[39] Certain modes in such systems exhibit ring-like features[40], but we note that a wider range of modes with different features, including central peaks, should occur for bubbles of different sizes, whereas the experimental system studied here only shows ring-type features. Furthermore, the Fermi wavelengths required for such long-range oscillations are inconsistent with density fluctuations on the order ~1×1012 cm-2 in the experimental system. Sharp interfaces can also give rise to significant localised states in their own right[34,35] due to discontinuities in the strain profile, but the bubbles in our system do not present sharp enough interfaces when approximated using the membrane model. Furthermore, no discontinuities are noticed in the strain profiles near the bubble edges in Fig 5. A recent work by Huang et al. reported qualitatively similar features in Raman on graphene bubbles[41], which they attribute to optical standing waves forming inside the bubbles at certain heights. However, the bubbles studied by Huang et al. are on the order of hundreds of nanometres high and several microns in radius, much larger than those studied in this work, which we measure from AFM to be 18±3 nm in height and 600±100 nm in radius. It is therefore unlikely that the doping variations we see are due to the same mechanism. Increased reactivity is predicted when graphene is strained[42], and a consequent charge transfer between graphene and various adsorbants could lead to charge density fluctuations. However, the increase in charge density in our case does not coincide with the maximum strain, and there are no defect signatures in the Raman data. Ruling out these possibilities, the charge density features could be caused by modulation of the graphene-hBN interaction near the cavity edge, possibly due to varying interlayer separation[43] or edge potentials at hBN edges.[44] Conclusions We have shown that nanoscale variations in strain and doping in graphene-hBN heterostructures are not random, but instead arise from local features in the heterostructures. Fractures, folds, edges and bubbles are all examples of nanostructures affecting these variations. If strain engineering of 2D materials is to become viable, we will need a greater understanding of what part these features play in determining strain and doping at a more global level. Confocal Raman equipped with the vector analysis method can provide a purely optical non-invasive tool for probing strain and doping in graphene devices at the nanoscale, and can reveal variations in the physical parameters, which are not accessible by methods such as AFM and optical microscopy alone. Methods Raman measurements. The Raman measurements were performed using a Renishaw inVia confocal Raman microscope, using a 532 nm excitation laser, with ~10 mW power incident on the sample, and a 1800 line/mm diffraction grating. To determine whether the peak splitting is negligible, the polarisation configuration of the laser and detector also needs to be taken into account. With linear polarisation, the relative intensity of the two components of the split peaks varies with the angle of polarisation. If one component is very low in intensity, it might be difficult to observe splitting even if it is non-negligible.[18] To reduce this dependency, we use a quarter-wave plate to produce a circular polarisation in the incident laser beam and we leave the scattered beam free of polarisation optics. This results in a measurement setup that has only a small polarisation angle preference introduced by the spectrometer grating, typically on the order of 10%. In the samples studied here we did not observe any peak splitting, so to find the mean Raman shifts a single Lorentzian was fitted to each peak. Atomic force microscopy measurements. AFM measurements were performed with a Bruker Dimension Icon scanning probe microscope, using Peak Force tapping mode with Bruker PFQNE-AL probes. Heterostructure fabrication. The van der Waals heterostructures studied in this work were assembled using a method described by Pizzocchero et al.[32] According to the method, the hBN and graphene are stacked on SiO2 at a temperature of 110°C, so it is to be expected that some strain will be introduced as the structure cools to ambient temperature, due to the different thermal expansion coefficients of the materials.[15] The bubbles seen in heterostructure B are consistent with bubbles described by Pizzocchero et al. and are most likely the result of stacking the hBN and graphene at a temperature lower than is recommended.[32] Acknowledgements This project has received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement GrapheneCore2 785219 number. The work has been also financially supported by the Department for Business, Energy and Industrial Strategy though NMS funding. The work has been partly realized within the Joint Research Project 16NRM01 GRACE: Developing electrical characterisation methods for future graphene electronics. This project has received funding from the EMPIR programme co-financed by the Participating States and from the European Union’s Horizon 2020 research and innovation programme. S.R.P. acknowledges funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skodowska-Curie grant agreement No 665919 and from the Irish Research Council under the laureate awards programme. ICN2 is funded by the CERCA Programme/Generalitat de Catalunya and supported by the Severo Ochoa programme (MINECO, Grant. No. SEV-2013-0295). The Center for Nanostructured Graphene (CNG) is sponsored by the Danish Research Foundation, Project DNRF103. The authors are grateful to Eli Castañon for fabrication of a heterostructure and Irina Grigorieva for useful discussions. References [1] Butler S Z et al 2013 Progress, challenges, and opportunities in two-dimensional materials beyond graphene ACS Nano 7 2898–926 [2] Novoselov K S, Mishchenko A, Carvalho A and Castro Neto A H 2016 2D materials and van der Waals heterostructures Science (80-. ). 353 aac9439 [3] Couto N J G, Costanzo D, Engels S, Ki D K, Watanabe K, Taniguchi T, Stampfer C, Guinea F and Morpurgo A F 2014 Random strain fluctuations as dominant disorder source for high-quality on-substrate graphene devices Phys. Rev. X 4 1–13 [4] Huang M, Yan H, Heinz T F and Hone J 2010 Probing strain-induced electronic structure change in graphene by Raman spectroscopy Nano Lett. 10 4074–9 [5] Fei Z et al 2016 Ultraconfined Plasmonic Hotspots Inside Graphene Nanobubbles Nano Lett. 16 7842–8 [6] N. Levy, S. A. Burke, K. L. Meaker, M. Panlasigui, A. Zettl, F. Guinea, A. H. Castro Neto M F C 2010 Strain-Induced Pseudo–Magnetic Fields Greater Than 300 Tesla in Graphene Nanobubbles Science (80-. ). 544–7 [7] Schneider M, Faria D, Viola Kusminskiy S and Sandler N 2015 Local sublattice symmetry breaking for graphene with a centrosymmetric deformation Phys. Rev. B - Condens. Matter Mater. Phys. 91 1–5 [8] Settnes M, Power S R and Jauho A P 2016 Pseudomagnetic fields and triaxial strain in graphene Phys. Rev. B 93 1–8 [9] Georgi A et al 2017 Tuning the Pseudospin Polarization of Graphene by a Pseudomagnetic Field Nano Lett. 17 2240–5 [10] Qi Z, Bahamon D A, Pereira V M, Park H S, Campbell D K and Neto A H C 2013 Resonant Tunneling in Graphene Pseudomagnetic Quantum Dots Nano Lett. 13 2692–7 [11] Settnes M, Power S R, Brandbyge M and Jauho A-P 2016 Graphene Nanobubbles as Valley Filters and Beam Splitters Phys. Rev. Lett. 117 276801 [12] Si C, Sun Z and Liu F 2016 Strain engineering of graphene: a review Nanoscale 8 3207–17 [13] Neumann C et al 2015 Raman spectroscopy as probe of nanometre-scale strain variations in graphene Nat. Commun. 6 1–7 [14] Ahn G, Kim H R, Ko T Y, Choi K, Watanabe K, Taniguchi T, Hong B H and Ryu S 2013 Optical probing of the electronic interaction between graphene and hexagonal boron nitride ACS Nano 7 1533–41 [15] Bao W, Miao F, Chen Z, Zhang H, Jang W, Dames C and Lau C N 2009 Controlled ripple texturing of suspended graphene and ultrathin graphite membranes Nat. Nanotechnol. 4 562–6 [16] Lee J E, Ahn G, Shim J, Lee Y S and Ryu S 2012 Optical separation of mechanical strain from charge doping in graphene Nat. Commun. 3 1024 [17] Banszerus L, Janssen H, Otto M, Epping A, Taniguchi T, Watanabe K, Beschoten B, Neumaier D and Stampfer C 2017 Identifying suitable substrates for high-quality graphene-based heterostructures 2D Mater. 4 025030 [18] Mueller N S et al 2017 Evaluating arbitrary strain configurations and doping in graphene with Raman spectroscopy 2D Mater. 5 015016 [19] Das A et al 2008 Monitoring dopants by Raman scattering in an electrochemically top-gated graphene transistor Nat. Nanotechnol. 3 210–5 [20] Das A, Chakraborty B, Piscanec S, Pisana S, Sood A K and Ferrari A C 2009 Phonon renormalization in doped bilayer graphene Phys. Rev. B - Condens. Matter Mater. Phys. 79 1– 7 [21] Ding F, Ji H, Chen Y, Herklotz A, Dörr K, Mei Y, Rastelli A and Schmidt O G 2010 Stretchable graphene: A close look at fundamental parameters through biaxial straining Nano Lett. 10 3453–8 [22] Metzger C, Rémi S, Liu M, Kusminskiy S V., Castro Neto A H, Swan A K and Goldberg B B 2010 Biaxial strain in graphene adhered to shallow depressions Nano Lett. 10 6–10 [23] Zabel J, Nair R R, Ott A, Georgiou T, Geim A K, Novoselov K S and Casiraghi C 2012 Raman spectroscopy of graphene and bilayer under biaxial strain: Bubbles and balloons Nano Lett. 12 617–21 [24] Bruna M, Ott A K, Ijäs M, Yoon D, Sassi U and Ferrari A C 2014 Doping dependence of the Raman spectrum of defected graphene ACS Nano 8 7432–41 [25] Froehlicher G and Berciaud S 2015 Raman spectroscopy of electrochemically gated graphene transistors: Geometrical capacitance, electron-phonon, electron-electron, and electron- defect scattering Phys. Rev. B - Condens. Matter Mater. Phys. 91 1–17 [26] Mohiuddin T M G et al 2009 Uniaxial strain in graphene by Raman spectroscopy: G peak splitting, Grüneisen parameters, and sample orientation Phys. Rev. B - Condens. Matter Mater. Phys. 79 1–8 [27] Mohr M, Maultzsch J and Thomsen C 2010 Splitting of the Raman 2D band of graphene subjected to strain Phys. Rev. B - Condens. Matter Mater. Phys. 82 1–4 [28] Pisana S, Lazzeri M, Casiraghi C, Novoselov K S, Geim A K, Ferrari A C and Mauri F 2007 Breakdown of the adiabatic Born-Oppenheimer approximation in graphene Nat. Mater. 6 198–201 [29] Lazzeri M and Mauri F 2006 Nonadiabatic Kohn anomaly in a doped graphene monolayer Phys. Rev. Lett. 97 29–32 [30] Ferrari A C and Basko D M 2013 Raman spectroscopy as a versatile tool for studying the properties of graphene Nat. Nanotechnol. 8 235–46 [31] Thomsen C, Reich S and Ordejón P 2002 Ab initio determination of the phonon deformation potentials of graphene Phys. Rev. B - Condens. Matter Mater. Phys. 65 1–3 [32] Pizzocchero F, Gammelgaard L, Jessen B S, Caridad J M, Wang L, Hone J, Bøggild P and Booth T J 2016 The hot pick-up technique for batch assembly of van der Waals heterostructures Nat. Commun. 7 11894 [33] Khestanova E, Guinea F, Fumagalli L, Geim A K and Grigorieva I V. 2016 Universal shape and pressure inside bubbles appearing in van der Waals heterostructures Nat. Commun. 7 1–10 [34] Settnes M, Power S R, Lin J, Petersen D H and Jauho A-P 2015 Bubbles in graphene - a computational study J. Phys. Conf. Ser. 647 012022 [35] Settnes M, Power S R, Lin J, Petersen D H and Jauho A-P 2015 Patched Green’s function techniques for two-dimensional systems: Electronic behavior of bubbles and perforations in graphene Phys. Rev. B 91 125408 [36] Carrillo-Bastos R, Faria D, Latgé A, Mireles F and Sandler N 2014 Gaussian deformations in graphene ribbons: Flowers and confinement Phys. Rev. B 90 041411 [37] Moldovan D, Ramezani Masir M and Peeters F M 2013 Electronic states in a graphene flake strained by a Gaussian bump Phys. Rev. B 88 035446 [38] Kim E-A and Castro Neto A H 2008 Graphene as an electronic membrane EPL (Europhysics Lett. 84 57007 [39] Crommie M F, Lutz C P and Eigler D M 1993 Confinement of Electrons to Quantum Corrals on a Metal Surface Science (80-. ). 262 218–20 [40] Jolie W, Craes F, Petrović M, Atodiresei N, Caciuc V, Blügel S, Kralj M, Michely T and Busse C 2014 Confinement of Dirac electrons in graphene quantum dots Phys. Rev. B 89 155435 [41] Huang Y et al 2018 Raman Spectral Band Oscillations in Large Graphene Bubbles Phys. Rev. Lett. 120 186104 [42] Wu Q et al 2013 Selective surface functionalization at regions of high local curvature in graphene Chem. Commun. 49 677–9 [43] Yankowitz M, Watanabe K, Taniguchi T, San-Jose P and LeRoy B J 2016 Pressure-induced commensurate stacking of graphene on boron nitride Nat. Commun. 7 13168 [44] Zheng F, Sasaki K, Saito R, Duan W and Gu B-L 2009 Edge States of Zigzag Boron Nitride Nanoribbons J. Phys. Soc. Japan 78 074713