Giant apparent lattice distortions in STM images of corrugated sp2-hybridised monolayers

We report on the strengths and limitations of scanning tunnelling microscopy (STM) when used for characterising atomic-scale features of quasi two-dimensional materials, such as graphene and single layers of hexagonal boron nitride, which may present strong corrugations when grown epitaxially on a substrate with a lattice mismatch. As a paradigmatic test case, we choose single-layer and bilayer graphene on Ru(0001), because their STM images show both a long-range moiré modulation and complex atomic-scale distortions of the graphene lattice. Through high-resolution STM measurements, we first determine with high accuracy the moiré epitaxial relations of the single layer and the bilayer with respect to the metal substrate. In particular, we also provide direct evidence for the existence of AA-stacked bilayer graphene domains on Ru(0001). We then demonstrate that the local strain distribution, as inferred from the same STM images, can be affected by large errors, so that apparent giant strains arise in some regions of the moiré as an imaging artefact. With the aid of density functional theory simulations, we track down the origin of these fictitious distortions in the high directionality of the graphene π-orbital density combined with the large corrugation of the sample. The proposed theoretical model correctly accounts for the observed dependence of the apparent strain on the STM tip–sample separation and on the different degree of curvature of the second graphene layer with respect to the single layer.


Introduction
After the rise of graphene [1], many other atomically thin materials have been proposed and fabricated [2] in view of their great potential, both as a fundamental physics playground and as a new technological platform. For instance, electronic devices based on the combination of metallic (e.g. graphene) and insulating (e.g. hexagonal boron nitride-h-BN) variants of these quasi two-dimensional materials have already been realised [3] and might provide an alternative to current silicon based integrated circuit technologies. Nanostructure engineering provides a powerful way to tailor the electronic properties of these materials and also to create new functionalities [4]. The intrinsic graphene electronic structure can be modified by introducing a periodic potential through the substrate [5], through a superlattice of adsorbed metal clusters [6], or by engineering a precise strain distribution [7]. The simultaneous presence of strain and local curvature can induce giant pseudomagnetic fields in graphene nanobubbles [8] and could be at the origin of electron-hole puddles in rippled graphene [9]. Novel nanostructures with very large local curvatures and, possibly, large strains are not uncommon in nanostructured sp 2 -bonded atomically thin layers. For instance, bubbles or 'blisters' were created through the controlled delamination of graphene [10], acid treatment of highly oriented pyrolytic graphite [11], rare-gas implantation below h-BN [12], or adsorption of transition metal atoms onto it [13].
An accurate measurement of the strain distribution in such nanostructures is highly desirable, yet there have been only few attempts to date, sometimes with contrasting outcomes. For graphene on Ir(111), no significant stretching or compression of the C-C bonds was reported from low-energy electron diffraction [14]. In the case of graphene on Ru(0001), x-ray [15] and low-energy electron [16] diffraction studies agree on the lateral size of the unit cell, but they derive very different strain distributions. The first study suggests strains up to 7% [15] and chiral lattice deformations [17], while the second concludes that the maximum strain is around 1% [16].
Scanning probe techniques, such as atomic force microscopy (AFM) and scanning tunnelling microscopy (STM), provide a direct visualisation of these systems in real space, so it is tempting to determine the atomic positions from the images and then infer the local strain distribution. Noncontact AFM with functionalised tips yields an unprecedented resolution of individual bonds in aromatic molecules [18,19] and has already been employed to study the surface topography of graphene on Ir(111) with picometer accuracy [14]. However, it was shown for CO-functionalised tips that the interplay of the CO bending and the nonlinear background signal arising from the neighbouring atoms limits the achievable accuracy on the bond length determination [20]. STM images of graphene 'bubbles' fabricated through the controlled oxidation of the Ru(0001) substrate have been used to derive strain maps of the graphene lattice [10]. Such maps reported tensile strains exceeding 10%, which appear to be in contrast with the exceptionally high lateral stiffness of graphene [21,22]. Indeed, it had already been noticed that apparent lattice constant variations of 0.03 nm between regions of opposite curvature show up in STM images of rippled graphene layers deposited on silicon dioxide. The observed effect was attributed to the tilting of the graphene π-orbitals [23], but a complete and quantitative understanding of this phenomenon has remained elusive.
In this work, we select epitaxial graphene on Ru(0001) as an exemplification of a curved sp 2 -bonded network and investigate the moiré structure and the graphene lattice distortions through high-resolution STM images. We first determine with high accuracy the moiré epitaxial relations for both the single layer and the bilayer systems, and provide evidence in favour of AA-stacking of the bilayer. We find that the real unit cells of these structures are much larger than formerly reported, and therefore they are at the limit of incommensurability, very much as in the system graphene/Ir(111) [14,24]. Then we address the atomic-scale features of the graphene lattice and derive the distributions of the C-C bond lengths for both the monolayer and the bilayer. We perform density functional theory (DFT) calculations for these two systems and show that the giant strains evidenced by the STM images are to a large extent fictitious. The high directionality and tilting of the π-orbitals quantitatively account for the observed apparent distortions and explain the detailed features of the apparent bond length distribution, such as the presence of a peak at large tensile strains and its dependence on the STM tip-sample distance. Since imaging artefacts identified in this work generally apply to any curved sp 2 -monolayer, reports on giant atomic distortions from scanning probe imaging of these systems should be judged with great care.
The paper is organised as follows: in section 2, we report high-resolution STM images of graphene on Ru(0001) and of h-BN on Rh(111) exhibiting giant apparent lattice distortions. Then we focus on graphene/Ru(0001) and determine the moiré epitaxial relations of the monolayer and of the AA-stacked bilayer in section 3. Apparent C-C bond length maps and strain distributions derived from the STM images of these two cases are presented in section 4, while the DFT-based modelling of the same systems is given in section 5. Finally, we draw our conclusions in section 6.

STM images of graphene/Ru(0001) and h-BN/Rh(111): giant apparent lattice distortions
In this section, we briefly describe two systems for which apparent giant distortions of the sp 2 -bond network can be spotted by eye from high-resolution STM images: graphene on Ru(0001) and h-BN on Rh(111), both shown in figure 1 (see appendix D for sample preparation).
The STM image of monolayer graphene/Ru(0001) in figure 1(a) shows a smoothly corrugated moiré pattern with graphene 'hills' appearing much brighter than the 'valleys', where graphene is much closer to the Ru surface and rather flat. The origin of this arrangement is known [25][26][27] and stems from the change of lattice registry between graphene and substrate, leading to strongly substrate-bonded graphene regions (the 'valleys'), where either of the two C atoms in the graphene unit cell is located atop a Ru atom, alternating with weakly, van der Waals-bonded regions (the 'hills'), where the C 6 rings are centred on the Ru atoms. The first type of graphene/ Ru(0001) stacking appears approximately in 2/3 of the moiré approximate unit cell, while the second in the remaining 1/3, resulting in the observed fraction of dark background to bright hills.
The approximate unit cell, containing only one of such hills, was earlier interpreted as a( ) 11 11 [28][29][30][31] or a( ) 10 10 Ru(0001) surface cell 9 [32,33]. Successive reports have clarified that not all such hills are equivalent [15,16,34]. This fact can also be recognised in figure 1(a), if one concentrates on the hexagonal lattice of depressions corresponding to the C 6 ring centres 10 . From the C 6 rings closest to the hill apexes, marked with red spots in the image, it is evident that adjacent hills are indeed different. The( ) 23 23 structure determined from diffraction techniques [15,16], then confirmed by a combination of STM measurements and ab initio calculations [34], implies four such hills per unit cell. However, our data show no exact repetition even afteŕ ( ) 2 2 such hills, indicating that the real surface unit cell is even larger (see section 3 for more details).
A closer look at the STM image in figure 1(a) (e.g., under grazing angle) reveals that the honeycomb centres of graphene strongly deviate from straight lines. The apparent lattice distortions can be roughly quantified as large as ±10% (see section 4 for a better estimate), which is quite surprising in view of the large stiffness of the C-C bonds in graphene [21,22]. The same kind of distortions can be seen in the STM image of monolayer h-BN on Rh(111) presented in figure 1(b). This system forms a( ) 12 12 moiré pattern and shows alternating regions of stronger and weaker substrate/monolayer bonding [4,35], as in the case of graphene/Ru(0001). However, the relative surface areas of these regions are here inverted 11 with respect to graphene/Ru(0001), and much more abrupt changes from the strongly to the weakly interacting parts of the moiré are observed. In this case, the honeycomb centres of the sp 2 -monolayers deviate even more prominently from straight lines, especially where the high-lying areas bend down to the flat depressions, i.e., where the local curvature is the highest. This can only partially be explained by the lower elastic moduli of h-BN [36] with respect to graphene [22], and points to a correlation between apparent lattice distortions and curvature in the monolayer.
In the following, we will focus on graphene/Ru(0001). After a detailed structural characterisation of the moirés of the monolayer and bilayer systems, we will quantify the atomic-scale distortions of the graphene lattice and identify their origin. pA, T=4.7 K). In both cases, the honeycomb centres are imaged dark. Their positions strongly deviate from straight lines suggesting a large apparent strain. 9 We use here and in the following the surface science convention of labelling the surface unit cell with respect to the substrate primitive vectors. 10 The depressions are attributed to the C 6 ring centres because the contrast between the two C sublattices is much smaller than the one between honeycomb centres and C atoms. 11 This is due to the two disparate atoms in the unit cell [35]. While the B atoms bind weakly to Rh independently of their adsorption site, the N atoms form a strong chemical bond when they are on top of Rh atoms (1/3 of the surface), whereas no such bond is formed in the 2/3 of the surface where N atoms are situated on either of the two threefold Rh hollow sites. Figure 2 shows a single Ru(0001) terrace covered by one monolayer (ML) and two MLs of graphene on the lefthand and on the right-hand sides, respectively. The profile of the 1st ML taken along the red line shows, in addition to the depression of the C 6 -ring centres, a clear contrast of 4±1pm between the two C sublattices. This contrast is inverted as we follow the line scan across the moiré. Thereby it directly reveals the stacking transition from moiré areas where C A atoms are situated on the strongly binding Ru top sites and C B atoms on the more weakly binding threefold substrate hollow sites, to moiré areas where this stacking is reversed.

Moiré unit cells of monolayer and bilayer graphene on Ru(0001) from STM
On the 2nd ML, the hexagonally arranged honeycomb centres are imaged as protrusions. As on the 1st ML, both C atoms are individually resolved on this layer as well, but this time as depressions and without a significant apparent height difference between the two sublattices (figure 2, blue line). The equivalent apparent height observed for both C sublattices indicates the existence of AA stacking for graphene bilayers on Ru(0001), and has been reported before [5,37]. In figure 2, the appearance of monolayer and bilayer graphene in a single image also allows us to directly infer the AA-stacking stacking by projecting the honeycomb positions of one layer onto the other (yellow lines). We note that this inference relies on the assumption that the 2nd ML in figure 2 has grown over the 1st one. Large scale images show that 2nd ML areas with AA stacking coexist with those exhibiting AB stacking [5]. However, our sample does not present sufficiently large 2nd ML areas to observe a continuous variation from one type of stacking to the other. We remark that in figure 2(a) the second layer shows a much smaller vertical corrugation than the first uncovered layer. This seems to be at variance with previous STM measurements, showing height profile variations larger than 1Å also for the 2nd ML [37]. pm. The value of the second layer is in good agreement with a recent STM study of the bilayer system [37]. Note that this value also applies to the buried layer since both layers are pseudomorphic and AA-stacked, as inferred above. Although small, the difference between the first and second ML values is significant, showing that the 2nd layer is, as expected, closer to the graphite lattice constant of =  a 246.4 0.2 graphite pm [38], while the 1st layer has an average tensile strain of  ( ) 1.3 0.2 % with respect to the 2nd due to its epitaxy with Ru. Although barely visible in figure 2, the 2nd ML exhibits a moiré pattern as well. It results from the mismatch between Ru and the AA-stacked graphene bilayer. The moiré patterns of both layers are visible in the large scale images in figure 3, emphasised by the superimposed hexagons for clarity. We have determined the moiré unit cells of both layers from the Fourier transforms of these STM images, which are presented in figure 4.
In both cases one can distinguish four of the six first-order spots of graphene, and six first-order spots of the moiré reciprocal lattices. The 1st ML shows many additional higher-order spots originating from the moiré. Labelling the reciprocal lattice vectors of the moiré M and of graphene G, we obtain G/M ratios of 12.57±0.06 and 11.57±0.06 for the 1st and 2nd ML, respectively. The respective angles between M and G are a =  ( ) • 4.0 0.2 and 12  ( ) • 4.0 0.7 . For the 1st layer this rotation can be readily seen in the misalignment between the first-order spot of the graphene atomic lattice (blue square) and the successive-order spots of the moiré lattice (red open circles). A similar rotation has been reported by Borcaet al [39], whereas diffraction studies concluded that the graphene and Ru lattices are aligned [15,16]. For the 2nd ML, the moiré corrugation is too weak to produce multiple-order spots in the Fourier transform, hence the larger error in α.
This collection of observations leads to the reciprocal space structural models sketched in figures 4(c) and (d). The reciprocal lattice vector of the substrate is given by = -S G M [24,39,40]. The commensurate structures agreeing best with experiment are = + n G M M 1 2 with n=12 for the 1st and n=11 for the  12 In order to compensate for the piezo creep along the slow scan direction, the angles between the first order graphene spots were normalised to 60 • .
2ndML. As can be seen from the tables in figures 4(e) and (f), their G/M, α, and a values agree with the measured ones within the error bars. Therefore, our proposed moiré unit cells are´ ( ) 11.57 11.57 R0.3 and  ( ) 10.54 10.54 R0.4 for the 1st and 2nd ML graphene/Ru(0001), respectively (by convention, the angles are those between G and S). Due to the slight rotation, the period of the unit cell defining a perfectly commensurate structure would be rather large, namely, in terms of the substrate primitive vectors,( ) 133 133 for the 1st ML If one ignores the graphene/Ru rotation angle, the 11.57 period of the 1st ML is about half that of thé ( ) 23 23 unit cell derived from x-ray diffraction [15,16]. However, the unrotated structure corresponding to this unit cell has never been observed in real-space. Diffraction techniques average over the two equally abundant rotational domains, leading to the conclusion that the graphene and Ru lattices are aligned [39]. The structure of the AA-stacked 2nd ML is reported here for the first time, and the smallest integer unit cell iś ( ) 21 21 if we disregard the rotation.

Quantifying the STM apparent graphene strain
Following the procedure outlined in appendix A, we constructed a bond-length map of the graphene honeycomb lattice (see figure 5(a)) from the apparent positions of the C 6 -ring centres of figure 1(a). The C-C bond lengths are colour-coded with blue compressive and red tensile apparent strain, and the map is overlaid onto the original moiré pattern. The approximate( ) 11 11 unit cell is shown in yellow. The spatial strain distribution clearly follows the periodicity of the underlying moiré pattern. The tensile strain is located prevalently on top of the hills, while the strongly interacting regions present more variability, from slightly elongated to strongly compressed bonds. The normalised bond-length histogram in figure 5(b) exhibits a shoulder at around 7% apparent tensile strain and the extrema exceed ±10% apparent strain. The C-C bonds that give rise to the shoulder are in fact those lying in the hill regions of the moiré indicated by the red circle in figure 5(a). To elucidate this point we have collected the lengths of those high-lying bonds into a separate histogram (see figure 5(b)), which indeed explains the shoulder in the original histogram. The observation of such a high level of strain is puzzling given the large in-plane stiffness of graphene. It can only be obtained under very large external stress [41]. This makes highly strained geometries energetically unfavourable. In addition, having highly strained hill regions, where the graphene-metal interaction is weak, is also counter-intuitive. Similar observations were erroneously attributed to real strain in the past [10]. The first indication that these bond lengths are to a large extent the result of an imaging artefact is inferred from the evolution of the bond-length histogram in figure 5(b) with tip-sample distance. Retracting the tip by 1Åfrom the surface leads to higher apparent strain and a shift of the shoulder in the distribution towards larger bond lengths. Errors arising from the in-plane projection are more than an order of magnitude smaller than the observed apparent strain (details in appendix C). It is evident that a direct identification of the atomic positions as local extrema in the constant current contours is not possible in strongly rippled monolayers of graphene [23] and other sp 2 -hybridised materials, with serious consequences in the interpretation of data from scanning probe microscopy techniques. Nevertheless, the interplay between atomic geometry and electronic structure properties that underlies such imaging artefact still needs to be unraveled.

DFT modelling of the apparent graphene strain
In this section we show through DFT simulations of the graphene/Ru(0001) system how the electronic structure properties of the corrugated graphene affects the STM images and we propose a simple theoretical model that accounts for the main features of the imaging artefact. The calculations were carried out within a generalisedgradient approximation [42] plus semi-empirical corrections to van der Waals interactions [43] using the planewave pseudopotential code QUANTUM ESPRESSO [44]. The graphene/Ru(0001) moiré pattern was described in an approximate way using a( ) 12 12 graphene/( ) 11 11 Ru periodic simulation cell which contains a single moiré hill, as shown in figure 6(c). Since the differences among the four protrusions in a more realistić ( ) g 25 25 /Ru( ) 23 23 simulation cell were shown to be rather small [26,34], our choice constitutes a good approximation for the present purpose. Further details about the simulation cell and the computational parameters are given in appendix B.
In figure 5(c), we report the distributions of the actual (three-dimensional) C-C bond lengths obtained from the atomic positions in the optimised graphene/Ru(0001) structure, and of their projections on the Ru(0001) plane. The actual bond lengths from DFT indicate strains well below 3%, in stark contrast with the 10% maximal strain obtained from the analysis of the STM images. The geometrical projection of the C-C bonds on the Ru(0001) plane exclusively leads to an apparent reduction of bond length, and this applies only to a small fraction of the bonds (see also appendix C). Moreover, the symmetry assumptions in the bond-length extraction procedure from STM images introduce only minor changes in the histograms, as we show in appendix A. Therefore, the origin of apparent giant strains in these STM images should lie in a peculiar electronic structure effect, which invalidates the customary assumption that local extrema in the constant-current STM profile can be directly identified with atomic positions.
In this perspective, we considered a theoretical model based on the high directionality of the π-orbitals in graphene and other sp 2 -hybridised monolayers. In conjugated carbon systems, the orthogonality of the πorbital with respect to the σ-bonds can be expressed through the π-axial vector p v sketched in the inset of figure 6(b). It can be defined as the vector forming equal angles with the three C-C bonds, so it can be easily determined from the coordinates of the conjugated C atom and its three nearest neighbours [45]. The angle between p v and the σ-bonds can be written as p q + 2 p , where q p is the so-called pyramidalisation angle and q = 0 p for pristine graphene [45]. Hence, q p relates with the amount of sp 2 -to-sp 3 re-hybridisation in graphene [46], while p v gives the approximate orientation of the π-orbital and provides a direct link between the π-orbital density and the atomic geometry. As we will show later in this section (see figure 6(c) and its discussion, as well as appendix B), the experimental images can be reproduced quite well by only considering the DFT charge density associated to the graphene π-orbitals. We now demonstrate that the π-vector correctly describes the orientation of the π-orbital density also in our system, and that it can therefore be used to establish a connection between the apparent atomic positions in STM images and the actual ones.
In figure 6(a), we show an array of p v vectors computed for a row of C atoms in the optimised graphene/Ru (0001) structure, superimposed to the electronic charge density associated to the π-orbitals. It can be seen from the figure that the arrows representing the p v closely follow the maxima of the π-orbital density of graphene. Since the variation of the graphene work-function is small across the graphene/Ru(0001) unit cell [26], we here take all p v vectors with a fixed length ℓ, comparable to the experimental tip-sample distance, and identify their extremities with the bright spots associated to the corresponding C atoms in constant-current STM. It is already obvious from figure 6(a) that the curvature of graphene on Ru(0001) introduces strong deviations of p v with respect to the substrate surface normal and a misalignment between p v vectors of neighbouring C atoms, resulting in electron density maxima that are laterally shifted with respect to the real C position and an apparent shortening or stretching of C-C bonds.
In order to quantify this effect, we built histograms of the bond lengths obtained from these apparent C positions for two different lengths of p v , representing tip-sample distances of , respectively (see also further below). The histograms, shown in figure 6(b), are much broader than that built from the actual DFT bond lengths and shown previously in figure 5(c), and are now comparable to the experimental histograms in figure 5(b). The apparent strain exceeding 10% agrees well with the experimental STM observation, as does the appearance of a shoulder at large strains, which is resolved as a peak in the theoretical histograms. The origin of this shoulder is the same as in the experiments, namely the subset of C-C bonds lying in the hill regions of the moiré (shaded areas in the histograms). In spite of its simplicity, this model captures other important features of the experimental histogram, such as the shift of the shoulder toward larger bond lengths, as well as the general broadening of the distribution with increasing tip-sample separation, showing that the 'tilting' of the π-orbital density due to the corrugation of graphene is the essential ingredient of the apparent strain in the STM images.
Further support for this interpretation is given by a similar analysis applied to the STM data for the 2nd ML graphene/Ru(0001) system and the corresponding atomistic DFT simulation. We used the same procedure outlined in appendix A to extract the apparent C-C bond-length map from an STM image of the 2nd ML graphene/Ru(0001), where now the ring centres are identified with the local maxima of the STM contrast (after subtraction of the long-wavelength moiré). The resulting bond-length histogram, shown in figure 7(a), displays a much narrower distribution compared to the 1st ML case in figure 5(b). The significant reduction of the The arrows are the π-axial vectors of the carbon π-orbitals, and the colours represent a contour plot of the LDOS integrated from −2.5eV to E F . (b) Histograms of apparent bond lengths obtained from the π-axial vector model applied to the DFT graphene geometry at 2.5 and 3.5Åtip-sample distance. The inset shows a schematic picture of a π-axial vector constructed from the atomic positions of a conjugated C atom and its three C-nearest-neighbours. (c) Simulated STM image (top) and apparent height profile (bottom) from LDOS integrated from −2.5eV to E F . The profile is taken along the direction shown, while the level of electron density has been chosen to give the same corrugation as in experiment. apparent lattice distortions in the 2nd ML with respect to the 1st ML can be readily rationalised within our model as a reduced tilting of the π-orbitals resulting from the smaller moiré corrugation of the 2nd ML. Figure 7(b) shows an apparent bond length histogram obtained from the π-axial vector model applied to the topmost layer of a bilayer graphene/Ru(0001) system as described in the DFT simulations (see appendix B). The simulated histogram reproduces the narrowing of the bond-length distribution, as well as the drop in the abundance of bonds with tensile strain seen in the experimental histograms in figure 7(a). The different methods used to compute C-C bond lengths can account for both the slight downward shift of the main peak in the DFT histogram for the 2nd ML with respect to the 1st ML, not seen in the experimental histograms, and the lack of a tail at compressive strains in the 2nd ML, which is instead present in the experimental histogram (a more detailed discussion can be found in appendix A).
We already pinpointed in figure 6(a) the correspondence between the orientation of the π-axial vector and that of charge density lobes of the π-orbitals. This suggests that an alternative simulation approach of the STM image based on the charge density, such as, for instance, the widely-adopted Tersoff-Hamann method [47], would lead to equivalent conclusions. Unfortunately, it turns out that this method fails in reproducing the experimental STM images presented in this work, e.g., figures 1(a) and 2(a), hinting at more complex tip-sample interactions [48] which are beyond the scope of our analysis. Nevertheless, we here present a supplemental analysis based on the DFT charge density of the π-orbitals for the simulated 1st ML case. In fact, as explained in more detail in appendix B, we find good agreement with the experimental STM images presented here when we integrate the local density of states (LDOS) down to about −2.5eV below the Fermi level, corresponding to the range of energies of the π-orbitals. In figure 6(a), the integrated LDOS is plotted on a plane normal to the Ru(0001) surface, together with the atomic geometry and the π-axial vectors for a subset of C atoms. A simulated constant-current STM image is then presented in figure 6(c), together with a representative linescan extracted from the same. Both the image and the linescan reproduce the main experimental features of figure 2(b), such as the sharp contrast between atomic positions and ring centres, the slight asymmetry of the two graphene sublattices in the covalently-bound region, and the contrast inversion between them upon their stacking transition on the Ru surface. The isosurface of electron density that gives an apparent surface corrugation on the graphene valleys comparable to experiment lies approximatively 2.5Å above the graphene layer, justifying our choice of d tip .
In conclusion, the approximate method used here to compute the direction of the π-orbitals is indeed accurate and the π-axial vector model captures the main features of the experimental STM images and bondlength histograms of graphene on Ru(0001). This model provides a clear interpretation of the apparent strain observed in the STM-imaged graphene network as a fictitious distortion related to the rippling of graphene/Ru (0001) and to the high directionality of its π-orbitals. Thanks to its simplicity, it can easily be applied to quantify the apparent distortions of other sp 2 -hybridised materials.

Conclusion
We demonstrated the appeal and the limitations of STM when studying sp 2 -hybridised monolayers. The epitaxial relations between these layers and their substrates (unit cell, rotation) or between subsequent layers (stacking) can be determined with high accuracy, and even the sublattice symmetry-breaking can clearly be resolved. However, when such layers present strong corrugations or bubbles, giant apparent lattice distortions can arise due to the directionality of their π-orbitals that dominate the STM imaging process. In the case of graphene on Ru(0001), we showed that this apparent local strain was mainly virtual, and under certain imaging conditions up to one order of magnitude larger than the actual strain. These effects can quantitatively be reproduced and understood in DFT simulations, therefore enabling a precise deconvolution into real and fictitious features of the STM images. Analogous apparent distortions are expected also in AFM images, where the atomic contrast depends even more strongly on the directions of dangling bonds. Figure A1 illustrates how we derive C-C bond-length maps (figure 5) from STM images with atomic resolution ( figure 1(a)). We first remove the long-wavelength corrugation of the moiré pattern by Fourier filtering the raw data in figure A1(a) leading to (b). On this image, the local minima, representing the graphene honeycomb centres, are recognised by our home written image processing program and fitted by gaussians to obtain their precise ( ) x y , coordinates. The resulting positions are overlaid onto the STM image as green crosses in figure A1(c). The C-C bond lengths are obtained from these positions assuming a perfect honeycomb lattice, where the C-C bonds lie on the lines joining each ring centre with its second-nearest neighbours and have a length of one third of the ring centre distance. The resulting bond-length map with colour code as in figure 5 of the main text is shown in figure A1(d) for a portion of figure A1(a). In the bottom-right corner, one bond is highlighted in yellow, together with the two ring centres used for the determination of its length.
Since the real graphene moiré lattice on Ru(0001) is not a perfect honeycomb lattice, we address the error induced by the above symmetry assumption by analysing the DFT data for the 1st ML. In figure A2(a) we compare C-C distances derived directly from the targeted C positions in the DFT π-orbital model with those derived indirectly from apparent ring-centre positions. The latter are obtained as centre of mass of the targeted positions of the 6 C atoms forming the ring. The mapping between direct and indirect bond reconstruction shows correspondence to a very good approximation. In particular, no systematic bias is introduced by the indirect method (i.e., through the ring-centre positions). The large majority of the bond lengths lie within a 5% error interval (grey), while only about 30 out of the 432 C-C bond lengths in the supercell show a larger error. In particular, the indirect method does not lead to an overestimation of the maximal values of apparent strain (both compressive and tensile), but rather slightly mitigates the effect. This is confirmed by comparing the histograms of bond lengths obtained from the two different methods, as shown in figure A2(b) for the 1st ML. The indirect method (green) induces a shift toward lower strains of the peak at high tensile strain and a slight decrease of the number of bonds at compressive strains, but the overall appearance and the width of the histogram is very close to that obtained through the direct method (blue). In figure A2(c), we compare bond-length histograms of the two methods for the 2nd ML. Also in this case, the direct and indirect methods lead to a similar shape and width of the histogram, with the appearance of a small spurious tail at compressive strains (also observed in the experimental histogram in figure 7(a)). We conclude that these differences are very small compared to the large apparent lattice distortions inferred from experiment.
DFT calculations were performed with the PBE gradient-corrected functional [42] using the QUANTUM ESPRESSO electronic structure code [44]. The interactions of nuclei and core electrons with the valence electrons were described through ultrasoft pseudopotentials [49,50], while the electronic wave functions and charge densities of the valence electrons were expanded in plane wave basis sets with kinetic energy cutoffs of 25 and 150 Ry, respectively. We also checked that increasing those cutoffs (up to 35 and 250 Ry, respectively) does not influence significantly the atomic geometry of the relaxed structures. The semi-empirical corrections to the van der Waals interactions proposed by Grimme [43] were included in order to ensure a more realistic description of the geometry [51].
The moiré structure of the 1st ML and the 2nd ML were simulated with a periodic supercell containing a four-layer-thick( ) 11 11 Ru(0001) slab overlaid with one or two( ) 12 12 graphene layers, respectively. Modeling the moiré as a single protrusion periodically repeated is sufficient for our purposes, since Iannuzzi and  [34]. Here, we used the experimental Ru lattice constant of 2.706Å, which is close to the theoretical one (2.746Å). Given the approximate periodicity, the graphene lattice constant (2.463 Å within PBE, very close to its experimental value of 2.461 Å) was stretched by less than 1% in order to match the lateral size of the( ) 11 11 Ru slab. We used a Methfessel-Paxton smearing of the occupations of 0.27eV and a Γ-point sampling of the Brillouin zone. This setup was used in previous works [26,51] and was shown to provide a good accuracy for this large supercell [51]. The Ru atoms in the two bottom layers of the slab were kept fixed in their bulk positions. The remaining atomic positions were relaxed until the corresponding forces dropped below 0.026eV Å −1 . The resulting corrugation of the 1st ML is about 1.33Å. Previous calculations employing the same functional, but only three Ru layers [51], found a slightly smaller value (1.195Å). Our test calculations with three and five Ru layers resulted in corrugations of 1.19Å and 1.33Å, respectively, indicating that at least four layers are needed to reach convergence.
Finally, we discuss the simulation of the STM images from the DFT electronic densities. In the standard Tersoff-Hamann method for STM image simulations [47], one generally assumes that the tunnelling current is proportional to the integrated LDOS of the sample at the location of the STM tip. Only electronic states with binding energies between the Fermi level, E F , and the tip-sample bias voltage (about −100 meV in our case) are usually included in the LDOS integral. We found that the experimental STM images presented in this work (e.g. Figures 1(a) and 2(a)) cannot be reproduced by following this recipe. STM images obtained by integrating the LDOS between −0.1 eV and E F (not shown) present a marked transition from a honeycomb pattern in the hill regions to a hexagonal pattern in the covalently bound regions of the moiré where the interaction with the Ru surface breaks the equivalence of the two graphene sublattices. This result is in agreement with previous theoretical works employing the same simulation cell and the same bias as in the present measurements [25,26], or a different bias [51]. This strong sublattice asymmetry in the covalently-bound graphene regions was also observed in previous STM experiments [52,53], but is only barely perceptible in figures 1(a) and 2(a). The reason for the asymmetry is the much lower DOS at the Fermi level for the C atom atop the surface Ru atom with respect to the C atom on the hollow site, which makes only the latter appear in usual STM images. To reproduce the peculiar features of our experimental setup, we had to extend the range of the DOS integral to energies down to about −2.5eV with respect to the Fermi level. In this way, the peak in the DOS of the atop C atom is captured and sublattice unbalance mitigated. Further extending the integration limit to lower energies would completely remove the sublattice asymmetry, which can yet be distinguished in the experimental line scans of figure 2(b). The inadequacy of the Tersoff-Hamann procedure in this specific case is even more obvious for the 2nd ML system, where the contrast inversion seen in the experiments could not be reproduced by any choice of the integration range, pointing to a higher-order tip-sample interaction mechanism [48].
Appendix C. Shortening factor in the 2D-projected C-C distances The 2D projection of the C-C bond network of the corrugated graphene layer onto a surface plane normal to the Ru[0001] direction can result in an apparent shortening of C-C distances, especially for the strongly tilted bonds. Here, we estimate the maximum shortening by taking into account the geometrical shape of the Rusupported graphene layer as inferred from the STM images and DFT simulations. Figure A2. Deriving C-C distances from ring centres versus apparent C positions. (a) Apparent C-C distances as derived from ring centre positions versus corresponding targeted C-C distances, as obtained for the 1st ML in the DFT model. The targeted C-C distances are built using the π-orbital model with a π-vector length of 2.5Å, while the ring centres are obtained as centres of mass of the targeted C positions. (b) Histograms of C-C distances in (a), as extracted from C positions (blue) or from ring centres (green). (c) As in (b), but for the 2nd graphene ML on Ru(0001).