Real-time observation of epitaxial graphene domain reorientation

Graphene films grown by vapour deposition tend to be polycrystalline due to the nucleation and growth of islands with different in-plane orientations. Here, using low-energy electron microscopy, we find that micron-sized graphene islands on Ir(111) rotate to a preferred orientation during thermal annealing. We observe three alignment mechanisms: the simultaneous growth of aligned domains and dissolution of rotated domains, that is, ‘ripening'; domain boundary motion within islands; and continuous lattice rotation of entire domains. By measuring the relative growth velocity of domains during ripening, we estimate that the driving force for alignment is on the order of 0.1 meV per C atom and increases with rotation angle. A simple model of the orientation-dependent energy associated with the moiré corrugation of the graphene sheet due to local variations in the graphene–substrate interaction reproduces the results. This work suggests new strategies for improving the van der Waals epitaxy of 2D materials.

A ttempts to grow large, structurally pristine twodimensional (2D) materials such as graphene and transition metal dichalcogenides are often frustrated by the fact that different rotational variants nucleate and grow, leading to a polycrystalline film [1][2][3][4][5] . This is not surprising given the weak film-substrate interactions inherent in the van der Waals epitaxy of these materials, which disables typical epitaxial control such as lattice matching. It is surprising, however, that rotational disorder in graphene films greatly varies with the substrate. For example, a single orientation prevails on Ru(0001) 6 and, under optimized growth conditions, on Ge(110) 7 . Graphene on Cu spans a small range of orientations 2 , whereas graphene on Pd(111) 8 and Ir(111) 4,9,10 can adopt many orientations. What then controls island orientation? The relevant factors are not understood. Consequently, current efforts to improve film quality focus on reducing the nucleation density [11][12][13] or using substratespecific recipes such as selective etching steps 14 . A more general approach requires understanding what dictates island orientation during van der Waals epitaxy.
Here we use low-energy electron microscopy (LEEM) to observe in real time the evolution of graphene domains during thermal annealing. We show that large, micron-scale graphene domains on Ir(111) can reorient themselves to a preferred orientation by three mechanisms after nucleation. We estimate that domains aligned with the Ir(111) lattice are energetically favourable by as little as B0.1 meV per C atom. A simple model of the orientation-dependent energy associated with the moiré corrugation of the graphene sheet due to local variations in the graphene-substrate interaction reproduces the results. These results provide valuable guidance for the synthesis of 2D materials on current and future substrates: annealing immediately after nucleation will improve quality and this effect will be greater for substrates that induce larger corrugations in the 2D sheets.

Results
Three mechanisms of graphene domain reorientation. The first mechanism identified, shown in Fig. 1a,b, is that some orientations grew, while others shrank simultaneously. Figure 1c compares the two frames and highlights the changes after 155 min. R0 and R11 domains grew, while the R12 domain simultaneously dissolved. The second process appears on further annealing. After 200 min, the reflectivity of the R12 domain circled in Fig. 1b changes in time. This transformation occurs over B80 min as shown by the recorded reflectivity in Fig. 1e. Selected-area low-energy electron diffraction (LEED; Fig. 1f) reveals that this entire region has undergone a continuous change in orientation from R12 to R4. Interestingly, the initial R12 domain dissolves as the newly created R4 domain grows. The third process is shown in Fig. 1g. As the first two processes occur, the boundary between the R11 and R12 domain, denoted by the red arrows, move into the R12 domain (see also Supplementary Movie 1). We performed 20 such annealing experiments at constant temperature in which these domain reorientation mechanisms were observed. More highly rotated domains always dissolved, while less-rotated domains grew. Within those 20 experiments, the continuous lattice rotation was observed five times and domain boundary motion was observed only once. Typically, an hour-long annealing step eliminated half of the rotational variants.
We now turn to the cause of the simultaneous growth and dissolution. In general, growth or dissolution is dictated by the surface adatom concentration adjacent to the island, c, relative to the equilibrium adatom concentration for the graphene domain, c eq . Growth occurs when c4c eq and dissolution occurs when (a) LEEM image of a sub-monolayer graphene film on Ir(111) at 1,050°C (75 mm field of view (FOV)). The spatial extent of individual domains of a given orientation can be deduced from differences in the reflectivity within the graphene islands. Rotation angle is denoted by the number following the letter R, as determined by selected-area low-energy electron diffraction (LEED coc eq . For an ideal adatom lattice gas, the chemical potential difference associated with growth or dissolution is DmEkT ln c/c eq , where k is the Boltzmann constant and T is the absolute temperature 15 . LEEM reflectivity measurements show that the adatom concentration around each island is the same throughout the field of observation (tens of microns), which is consistent with attachment-limited growth of graphene on Ir (refs 10,16). Therefore, the simultaneous growth of a domain of one orientation and dissolution of another with a different orientation indicates that the equilibrium adatom concentration of a given domain depends on its orientation relative to the substrate and, correspondingly, that the chemical potential is orientation dependent as illustrated in Fig. 2a. The relative chemical potential of different rotations can be immediately determined by combining the experimental results, as shown in Fig. 2b. All observations indicate that domains aligned with the Ir substrate, R0, have the lowest chemical potential and that the chemical potential increases as the angle of rotation increases from R0. This picture is consistent with the observed spontaneous change in domain orientation shown in Fig. 1d. A domain can lower its chemical potential by reducing its rotation angle. With its new orientation, the R4 domain has a lower c eq compared with the R12 domain, which results in the growth of the R4 domain after the transformation, while the remaining R12 domain continues to dissolve.
Temperature-dependent ripening. The simultaneous growth and dissolution, that is, 'ripening', is sensitive to temperature. As one sees in Fig. 3, graphene islands were imaged as the temperature was varied in steps of B5 K. At lower temperatures, all domains grew, while at high temperatures domains shrank. By adjusting the temperature, we could change the fraction of domains which grew. This temperature dependence is due to the segregation of the dissolved carbon in the Ir bulk to the surface when the amount of carbon dissolved in the bulk becomes greater than the solubility limit. At lower temperatures, bulk segregation causes the C adatom chemical potential to be above the equilibrium value of all the islands. As the temperature increases, the segregation flux to the surface decreases and the adatom concentration decreases below the equilibrium value.
What is the magnitude of the driving force for a rotated domain RX to realign itself to R0? To estimate it, we use the temperature-controlled ripening experiments described above to determine the differences in equilibrium adatom concentrations between different rotations. From these adatom concentrations, we then extract the difference in chemical potentials of the islands. In previous work, the velocity of R0 domains was precisely determined as a function of the adatom concentration and temperature 10     ARTICLE where c eq,R0 is the equilibrium adatom concentration for the R0 island, n is the cluster size for attachment to a kink and is equal to four, and B(T) ¼ exp(31-(33,800±1,300)/T) nms À 1 is the kinetic coefficient 10 . Converting the A(t) plot of Fig. 3c into a growth velocity allows us to determine the adatom concentration when the R14 and R8 domains are stationary and in equilibrium with the adatoms. If a rotated domain, RX, is stationary while a standalone R0 island continues to grow, the difference in chemical potential can be estimated from the adatom concentrations by When the R14 domain was stationary, the R0 velocity was 540 ± 36 nm h À 1 (1,000°C), where the error is from the s.d. of the dA/dt line fit. When the R8 domain was stationary from 2,800 to 3,200 s (1,030°C), the R0 velocity decreased to 144±36 nm h À 1 . Using equation (2), we find that the R14 domain has a chemical potential of 0.07 ± 0.01 meV per C atom above R0, while R8 is 0.01±0.01 meV per C atom above R0. Five temperature-controlled experiments were conducted in the same manner. The estimated differences in chemical potential relative to R0 are shown in Fig. 3d. The results were reproducible across the five experiments and did not depend on the shape or size of the islands (see Supplementary Note 1). Thus, we attribute the differences in chemical potential to differences in binding strength of the various rotations. The chemical potential increases with rotation from R0, as expected. The trend is not linear, with small differences in chemical potential for small rotations and larger increases in chemical potential as the rotation angle increases. Due to the uncertainty in B(T), the absolute scale of the chemical potential may be a factor of two larger or a quarter smaller.
These measurements allow us to estimate the difference in chemical potential that cause the observed grain boundary motion in Fig. 1g. It is extremely small, 0.007 meV per C atom. An alternative, albeit less direct, argument based on the boundary curvature induced by edge pinning gives a similar estimate (see Supplementary Note 2). We note that we only observed the boundary move in the case of one degree difference in rotation: evidently domain boundary mobility must decrease with relative misorientation and limit this mechanism of coarsening. On the other hand, the continuous lattice rotation was seen five times and was active when grain boundary motion was not.
Moiré corrugation model. These small energy differences seem beyond the capabilities of ab initio techniques such as density functional theory (DFT) to compute directly. To understand the origin of the orientation-dependent chemical potential differences and the extent to which they can be generalized, we make a model of the energy differences due to the corrugation of the graphene sheet, which is known to significantly change with orientation. In incommensurate, undistorted moirés, all equivalent positions of C atoms with respect to the substrate are occupied with equal probability, independent of the layer's rotation 18 , suggesting an equal energy for all rotational variants. Changing lattice misfit with the substrate does not alter this conclusion; so it is doubtful that lattice misfit determines the preferred orientation. However, it is experimentally observed that graphene sheets are not undistorted. Instead, experimental results 4,19-21 and DFT 22 have shown that graphene sheets on metals are corrugated. Scanning tunnelling microscopy images of R14, R19 (ref. 9) and R30 (ref. 4) graphene on Ir(111) reveal that, in addition to the large-scale modulation with moiré periodicity, the graphene topography also exhibits a pronounced fine-scale corrugation (in the case of R30 with a periodicity of B5 Å), which is absent in the R0 phase. We have used DFT to compute the structure as a function of orientation (see Supplementary Fig. 1 for the optimized structures). As seen in Fig. 4a, we reproduce these experimental trends.
The reason for this corrugation is that the preferred distance of the graphene sheet from the substrate varies locally, depending on the lateral position of the C atoms with respect to the substrate (for example, whether the C atom is at an atop position, a threefold hollow site and so on). To maximize its binding to the substrate, the graphene sheet flexes slightly to follow a contour of minimum C-Ir potential. Due to the bending rigidity of the graphene, the C atoms are further away from their optimal positions as the wavelength of the corrugation shortens. So the short-wavelength corrugations observed experimentally for rotated graphene would lead to an energy cost.
To estimate the magnitude of this effect, we construct a simple model for the observed corrugations. For simplicity, we assume that the graphene is incommensurate with respect to the Ir substrate, that is, that the possible energy reduction due to locking into a (higher order) commensurate lattice is not sufficient to overcome the energy cost of straining the graphene from its preferred lattice constant. We assume that the preferred distance z 0 of C atoms from the Ir(111) substrate is represented by a 2D sinusoid, which has opposite extrema at the atop sites and hollow sites. We take the force on each C atom to be proportional to how far the atom is displaced from the preferred distance and assume that the spring constant k is uniform. In the relaxed configuration, this force is balanced by the forces caused by bending the graphene sheet. In the continuum limit, the energy cost E b of periodic corrugations is 23,24 where l is the mean bending rigidity, B1.4 eV 24 , and K is the mean curvature. Our model approximates the local curvature K i at a given atom i as being proportional to the height difference Dz i between atom i and the average of its surrounding three nearest neighbours. This gives a force on each atom equal to 3 ffiffi ffi 3 p lDz i =d 2 , where d is the C-C distance in graphene (for a derivation, see Supplementary Note 3). The equilibrium sheet corrugation is given when this bending force balances the displacement force.
This model has only two parameters, the amplitude of the sinusoidal preferred C atom distance and the spring constant k. To apply this model to graphene on Ir(111), we fit these two parameters to our DFT calculations for R0 shown in Fig. 4a. DFT is reliable for R0 corrugations because it reproduces those from experiments 22 . From the DFT calculations, we obtain the equilibrium corrugation and the forces on each C atom for a sheet constrained to be flat. We then adjust the amplitude of the sinusoidal preferred separation distance and the spring constant k to exactly reproduce the amplitude of the corrugation of the relaxed sheet and the root mean square of the forces on the flat sheet. Given these parameters, we predicted the corrugation for domains rotated by 7.5, 15 and 30 degrees. Height profiles extracted from our moiré corrugation model along the closepacked zigzag direction are shown in Fig. 4b.
The agreement between this prediction and DFT calculations is striking. It reproduces the position of all C atoms in the profiles in Fig. 4b to within 0.05 Å, with most atoms within 0.01 Å. The good agreement makes it plausible that it correctly captures the energy due to C atom displacement and bending of the graphene sheet with changing orientation. And indeed, as shown in Fig. 4c, despite its extreme simplicity, the model correctly reproduces the observed stability of R0 and it easily accounts for the magnitude of the measured energy differences. That the model tends to overestimate the energy differences is not surprising given that the model does not allow lateral relaxation of atomic positions and the factor of two in the uncertainty of the scale of the experimental data (Fig. 3d). (Such small energy differences are generally difficult to compute directly with DFT because of the large moiré unit cells; however, we have used DFT to estimate the binding energy difference between R30 and R0 and find it to be less than B1 meV per C atom, consistent with our model and experiment. See Methods and Supplementary Table 1.) Thus, the driving forces for the post-nucleation reorientations that we observe can be explained by the orientation dependence of the moiré corrugation. This driving force applies to all three observed reorientation mechanisms.

Discussion
Our results have significant implications on the nature of rotational disorder of the graphene islands on various substrates. Scanning tunnelling microscopy shows that graphene on Au(111) 3 and Cu(111) 25 exhibits the same trend of increasing nearest-neighbour corrugation with rotation angle as seen on Ir(111). One would expect the R0 orientation to be energetically preferred. As expected, like on Ir(111) 16 , the majority orientation on Cu(111) 2 and Au(111) 26 is indeed R0. This is especially noteworthy for Au(111), where R0 has a lattice mismatch of 17% 26 and clearly supports the conclusion that the moiré corrugation dictates the resulting island orientation. This raises the possibility that in the systems where good alignment is observed in as-grown films, for example, Cu(111) 2 , Ru(0001) 6 , Au(111) 26 and Ge(110) 7 , such alignment occurs shortly after nucleation by the continuous lattice rotation mechanism when the islands are still small. Furthermore, these results imply that the rotational disorder of other 2D materials (hexagonal-BN, transition metal dichalcogenides) can be controlled only when grown on substrates where corrugations exist, for example, on metals [27][28][29][30] .
In summary, we have identified the energetically preferred graphene orientation on Ir(111). Domains can evolve to this orientation during thermal annealing by three distinct mechanisms. The driving force for the alignment can be as small as 0.1 meV per C atom, much less than the total binding strength of graphene to metal surfaces. We propose that the origin of the preferential alignment is caused by the varying degree to which C atoms can attain the preferred distance from the substrate: the graphene-bending rigidity prevents the sheet from following the short-wavelength corrugations inherent in highly rotated domains. Since graphene corrugations of similar size are ubiquitous, the driving force for the coarsening on Ir reported here is expected to be similar in other systems. Despite this small driving force for rotational alignment, our results indicate that post-nucleation annealing of 2D materials can still improve rotational order. This strategy means that high quality 2D materials can be grown on more substrates than is now recognized. Our understanding of the driving force for reorientation can guide the search for new substrates.

Methods
Experimental. Graphene islands were nucleated by heating the Ir substrate between 800 and 900°C and introducing ethylene gas into the LEEM chamber. After nucleation, the ethylene flow was turned off. The temperature was measured by a W-Rh thermocouple welded directly to the substrate.
Determining domain orientation. The domain orientation relative to the Ir(111) lattice was identified by selected-area LEED. The moiré LEED spots are sensitive to small changes in domain orientation and were used to determine the orientation to within ± 0.5°. The distance between each pair of moiré superstructure LEED spots was measured and normalized by the distance between the corresponding pair of first-order graphene LEED spots. This ratio was computed for the three pairs of spots per LEED pattern and then averaged. Over 47 LEED patterns were analyzed and orientations were assigned based on this ratio. This ratio is easily determined using a simple vector model and there was excellent agreement between the experimental and expected ratios.
Velocity measurement. The average R0 growth velocity was computed by measuring island area, A, and perimeter, p, with time: u ave ¼ (1/p ave )dA/dt, where dA/dt is the slope of a line fit to the area versus time data for the time interval in question, and p ave was the average perimeter over the time interval. Density functional theory. All DFT calculations were performed with the optB86b-vdW functional with spin-polarization and using the projectoraugmented-wave GW pseudopotentials supplied by VASP. A cutoff energy of 400 eV for the plane-wave basis and the Brillouin zone was sampled using 1 Â 1 Â 1 and 2 Â 2 Â 1 k-point grids. The optB86b-vdW exchange-correlation functional was employed to account for van der Waals interactions that are required for a qualitative description of the intermolecular interactions between the graphene sheet and the Ir(111) surface. The functional is of the form: x is the exchange energy, E LDA c is the local density approximation to the correlation energy and E vdW À DF c is the non-local correlation energy term.
The binding energies per C atom of the systems were determined by the following equation: where E G/Ir111 is the total energy of the system (R0, R7.5, R15, R22.5 or R30 on the Ir(111) substrate), E G and E Ir(111) are the total energies of the optimized isolated graphene layer and the Ir(111) surface, and N is the total number of C atoms. The graphene deformation energies (bending energies) were determined from the difference between the single-point energy of the optimized graphene sheet on the Ir(111) surface and the fully optimized graphene sheet within the same unit cell. Two k-point grids were used to compute graphene deformation energies and binding energies per C atom. The larger 2 Â 2 Â 1 k-point grid was used to test whether the predicted energies from the 1 Â 1 Â 1 k-point grid were converged with respect to k-space. The results indicate that the computed graphene deformation energies are converged; however, the binding energies are not. This conclusion was reached by comparing the difference in energies with respect to k-space (DE/Dk). For the deformation energies, DE/Dk is much smaller than the magnitude of the values. On the other hand, the relative binding energies between the different orientations (R0, R7.5, R15, R22.5 and R30) are smaller than DE/Dk. On the basis of this analysis, it is concluded that the predicted ordering of the deformation energies is statistically significant, but the relative binding energies between orientations is not significant. As a result, the orientation with the greatest binding energy per C atom cannot reliably be determined from these results. Increasing k-space much further is computationally intractable; moreover, each system has different lattice dimensions that affect the relative k-point spacing. Computationally determining differences of one-tenth of a meV is difficult given the various numerical approximations used within DFT (that is, k-points and energy cutoffs).