Graphene on Ru(0001): A corrugated and chiral structure

We present a structural analysis of the graphene/Ru(0001) system obtained by surface x-ray diffraction. The data were fit using Fourier-series expanded displacement fields from an ideal bulk structure, plus the application of symmetry constraints. The shape of the observed superstructure rods proves a reconstruction of the substrate, induced by strong bonding of graphene to ruthenium. Both the graphene layer and the underlying substrate are corrugated, with peak-to-peak heights of (0.82 +/- 0.15) A and (0.19 +/- 0.02) A for the graphene and topmost Ru-atomic layer, respectively. The Ru-corrugation decays slowly over several monolayers into the bulk. The system also exhibits chirality, whereby in-plane rotations of up to 2.0 degrees in those regions of the superstructure where the graphene is weakly bound are driven by elastic energy minimization.


3
argue that this symmetry breaking is driven by energy minimization based on elastic energy considerations.

Experimental
Sample preparation and the SXRD measurement setup at the Surface Diffraction Station of the Materials Science Beamline, Swiss Light Source, have already been detailed in [18]. It was demonstrated from simple simulations of the 25/23 superstructure rod (SSR) that the substrate must also be corrugated, since oscillations with the appropriate periodicity of approximately 1.0 out-of-plane substrate reciprocal lattice units (r.l.u., 2π/c) on the SSRs only start to appear if one includes a corrugation of the substrate. Here, we present further SXRD data from the same sample, which in addition to the SSRs now includes in-plane data.
Because of the very large number of atoms involved in the superstructure, it is impossible to fit each atomic position individually. Instead, we have parametrized the structural model using a small set of physically reasonable parameters. The in-plane and out-of-plane deviations of the atomic positions from an ideal flat structure of the graphene and of the uppermost layers of the Ru substrate are described by a 2D Fourier-series expansion. We truncate this series after the fourth Fourier component, since higher orders could not be resolved in the diffraction data. The displacement field of the system is allowed to adopt the lower p3 symmetry, since this is the lowest symmetry still compatible with the apparent measured sixfold diffraction symmetry, which only arises because of the superposition of the two possible terminations of the hexagonal close-packed (hcp) substrate [23]. Because the p3 symmetry allows chiral structures, we have to sum over the signals from domains of each enantiomer and assume a 50% distribution.
Details of the implementation of the Fourier expansion and of symmetry constraints are given in the appendix. Here, we discuss only those aspects that are needed to understand the results. First, it is important to note that because the 25-on-23 structure contains 2 × 2 corrugation periods, only the even Fourier components, that is, the second and fourth, must be considered. This is also demonstrated by the absence of signal at the 22/23, 24/23, . . . SSRs. For each atom within the supercell, the in-plane and out-of-plane deviations x, y and z are described by the two Fourier components. In total, both graphene and ruthenium require nine fitting parameters each in order to describe their corrugations.
In addition to the 18 corrugation parameters we introduce a factor, λ, which describes an exponential decay of the substrate corrugation amplitude with substrate depth z This decay applies to all the three amplitudes used for the description of the substrate displacement function. We fix the minimum distance from the substrate to graphene layer, d C−Ru , as 2.0 Å [20,24], since our model is relatively insensitive to this parameter within physically sensible limits (±0.1 Å). The parameter d Ru 1 −Ru 2 is the distance between the first and second Ruatomic layers. Lastly, a global scaling factor S is required, resulting in a total of 21 free-fitting parameters.
We begin by defining regions of the supercell, where we consider a flat graphene layer lying commensurably 25-on-23 on top of a flat Ru substrate (see figure 1). The gray-shaded region in figure 1(b) indicates where the first of the two C atoms within a 'normal' graphene unit cell sits on top of an Ru atom of the topmost substrate atomic layer (red atoms), and the second atom sits on top of an Ru atom from the second substrate atomic layer (green atoms, the hcp position). Henceforth, we refer to this as the (top, hcp) region. Using the same arguments, the red area is the (hcp, fcc) region and the green one is the (fcc, top) region [25]. Fitting 6 was performed using GenX [26], an optimization program using the differential evolution algorithm, which helps avoid getting trapped in local minima [27]. The errors of the fitted parameters are estimated by an increase in the GOF of 5%.
We fit d Ru 1 −Ru 2 to the CTR data alone (figure 2(a)), as this is sensitive to small differences in the interplanar spacing of the topmost two Ru-atomic layers but is largely insensitive to the form of the weakly scattering superstructure. The best fit had an R-factor of 5.2%, for d Ru 1 −Ru 2 = (2.080 ± 0.003) Å, which should be compared to a bulk value of 2.141 Å. This equates to a contraction of 2.8%, in agreement with the papers [28][29][30].

Results and discussion
The starting model for the search of all the other parameters was a strained 25-on-23 flat graphene layer lying commensurably on a flat ruthenium bulk structure. The best fit for the SSR and in-plane data has an R-factor of 13.4% (figures 2(b) and (c)). The peak-to-peak corrugation height of graphene is (0.82 ± 0.15) Å, in agreement with the papers [15,17,31], whereas that of the uppermost Ru-atomic layer is (0.19 ± 0.02) Å and is out-of-phase with respect to the graphene corrugation (figure 3). The exponential decay length of λ = (7.0 ± 0.4) Å means that there is still approximately a tenth of the distortion of the first Ru-atomic layer at a depth of four Ru-atomic layers. This strongly supports the idea of a chemisorbed graphene layer with significant interaction with the substrate [7,24], [32][33][34].  Details of the final structure are summarized in figure 4. Figure 4(a) shows a clear corrugation of the graphene with the hills lying in the weakly bound (hcp, fcc)-region. The hills have a triangular shape, in remarkable agreement with earlier STM data [15,16]. Although  in-plane movements of up to (0.25 ± 0.03) Å of the graphene are observed (Figure 4(b)), the bond lengths are distorted by less than 0.1 Å. This requires a twisting motion and indeed the in-plane movements exhibit a chiral signature, in which the largest movements occur at the steepest flanks of the hills, as one might expect, based on simple elastic strain considerations. Note that this feature emerged naturally from the fitting and was not implemented a priori into the model. The biggest rotation angle of the hexagons is 2.0 • found on the flanks as well as on top of the hills.
The elastic energy was calculated to test the physical validity of the presented parametrization approach and the resulting model. It takes into account the in-plane and outof-plane displacements of surface atoms from their 'ideal' positions due to the 25/23 surface reconstruction. From our model, we calculate an elastic energy [35][36][37] due to strain of 9.3 eV per supercell, assuming zero strain for a flat 25-on-23 graphene layer 7 . Fitting the data to the higher p3m1-symmetry results in an increase in elastic energy by 83%, while the R-factor of 14.7% is significantly higher than that for the p3-symmetry. Even if we were to assume zero strain for a flat graphene layer having the bulk graphite in-plane lattice constant (figure 4(c)), this has no significant influence on the energy difference between the two different symmetry models. A histogram of all the bond lengths in the graphene superstructure demonstrates that the implementation of the lower p3-symmetry allows the bond lengths to be more preserved relative to bulk graphite.
Very recently, an independent study using low-energy electron-diffraction (LEED) [38] has been published on the same Ru single crystal using the same graphene preparation and characterization, where the authors claim a corrugation of the graphene layer of 1.5 Å and a corrugation of the topmost ruthenium layer of 0.23 Å. In that study, the system is described by a p3m1-symmetry and the unit cell is cut down to one of the four inequivalent sub-unit cells. Both of these measures were taken in order to reduce computational time. An SXRD simulation of the coordinates extracted from the LEED study performed led us to a similarly high R-factor of 34.0% to that of the fit results of the LEED analysis. The reason for the discrepancies, which are far outside the error bars, are not yet resolved, although possible explanations are the already-mentioned restriction to p3m1-symmetry and a 12-on-11 superstructure-a full dynamical scattering LEED calculation of the system with p3-symmetry is at present beyond computational capabilities. In addition, the fact that LEED only probes the topmost layers, while SXRD demonstrates that significant vertical displacements occur down to at least four atomic layers of the Ru substrate, might also play an important role.

Summary and conclusion
In summary, we have determined the graphene/Ru(0001) structure in unsurpassed detail. This was only possible by adopting a parametric Fourier description of the superstructure using only a small number of physically reasonable parameters. Up to the mirror-symmetry breaking, the final model agrees excellently with previous STM studies. We find a graphene and ruthenium corrugation peak-to-peak height of (0.82 ± 0.15) Å and (0.19 ± 0.02) Å, respectively. The ruthenium corrugation is out-of-phase with that of graphene and decays exponentially down to a depth of several ruthenium layers. Importantly, we have also discovered the new and 8 potentially highly significant property of areal chirality in the in-plane movements, which are most evident on the flanks of the hills of the corrugation. We propose that this symmetrybreaking phenomenon is induced by elastic energy minimization of the graphene layer. To test the validity of this, we calculated the elastic energy of the graphene superstructure to be 9.3 eV, less than two-thirds of that for the p3m1 case.
In the following, the implementation of the symmetry constraints and the Fourier expansion to the graphene-on-ruthenium model will be briefly described. The displacement dr of an atom sitting at point r is expressed by its two-dimensional Fourier series where s and t ∈ {0, 2, 4} are the orders, A i s,t and B i s,t are the Fourier coefficients, φ i s,t are the phases of the corrugation and i ∈ {x, y, z}. Note that the phase of the out-of-plane displacements influences the valley and hill shapes and positions of the corrugation allowed by the p3symmetry ( figure A.1). Since by equating A i and B i to The rotation operators used for the description of the p3-symmetry are R 1 and R 2 , which in a hexagonal coordinate system describe a 120 • rotation counterclockwise and clockwise around the origin, respectively (figure A.2); they are defined by It can be easily shown that R 1 = R −1 2 ≡ R. Note that dR has to fulfill the p3-symmetry constraint, which results in  The p3-symmetry constraint operators. R is defined as a rotation by 120 • counterclockwise around the origin, while R −1 is the rotation clockwise by 120 • around the origin  Regarding the considered Fourier components in the analysis, the zeroth order is the 23/23 reflection, the first-and third-order components correspond to the 24/23 and 26/23 systematic absences, respectively, and the 25/23 to the second-order component. Hence the fourth order refers to the 27/23 reflection, and along the h-direction (equivalent to the k-direction) one can limit (s, t) = (2, 0) and (s, t) = (4, 0). For the sake of simplicity, we describe here only the implementation to the second order.
From equation (A.5) and table A.1, one can derive the following expressions for the single components of dr that describe the displacement field. We do not include the orders (s, t) for the sake of simplicity. Hence, we obtain six fitting parameters for the displacement field dr, namely A x , A y , A z , B x , B y and B z . Since we lock the phase for the fourth order to be the same as that for the second order, there will be nine fitting parameters. The effect the fourth-order harmonic with a locked phase has on the structure is shown in figure A.3. For an amplitude up to 0.25 of that of the second harmonic, the low regions in the structure will be flattened out.