A generalized Read–Shockley model and large scale simulations for the energy and structure of graphene grain boundaries

Abstract

The grain boundary (GB) energy is a quantity of fundamental importance for understanding several key properties of graphene. Here we present a comprehensive theoretical and numerical study of the entire space of symmetric and asymmetric graphene GBs. We have simulated over 79 000 graphene GBs to explore the configuration space of GBs in graphene. We use a generalized Read–Shockley theory and the Frank–Bilby relation to develop analytical expressions for the GB energy as a function of the misorientation angle and the line angle, and elucidate the salient structural features of the low energy GB configurations.

---

1 Introduction

Graphene – a two dimensional allotrope of carbon with excellent mechanical and electronic properties – has attracted much attention since it was first produced by direct exfoliation from graphite more than a decade ago.^1–10 Exfoliated graphene is largely monocrystalline; however, exfoliation is not a scalable production technique. Chemical vapor deposition (CVD) is the most widely used scalable production technique, and is being used to produce more than 300 000 m² of graphene annually.^5,6,10 Graphene produced from CVD is polycrystalline,¹¹ and thus it contains intrinsic line defects in the form of grain boundaries (GBs) that have been studied observed experimentally.^12–17 Such GBs have a profound effect on the properties of the polycrystalline materials; for instance, a high GB energy can promote grain growth, or a low GB strength can lead to brittle intergranular fracture. It has also been noted that while some graphene GBs offer minimal resistance to electron transport, other GBs can be highly insulating.¹⁸ Several recent studies have demonstrated that the thermal conductivity of graphene is strongly dependent on the GB structure.^19–22 The mechanical strength of graphene is influenced by GBs,^23–26 which can have a profound influence on sustainability application such as sea water purification by reverse osmosis. Additionally, the GB structure in polycrystalline single-layer graphene strongly modifies electronic transport.^27–31

Thus, it is clear that in order to understand the properties of polycrystalline graphene, it is necessary to characterize graphene GBs. Perhaps the most important property of a GB is its excess energy (per unit length). The excess GB energy, or simply the GB energy, has direct influence on the grain morphology,^32,33 and thus influences all grain morphology dependent properties including strength and transport. It is not feasible to measure the GB energy directly in an experiment. Instead, GB energy in crystalline materials is typically inferred from the equilibrium structure of triple junctions (interface of three GBs).^34–36 Such equilibrium junctions satisfy the Herring equations,³⁴ which can be used to deduce the GB energy if a statistically large amount of experimental data is available.³⁶ No such experimental study has been performed for graphene. Indeed, we (with co-authors) published the first statistically large dataset of high-resolution transmission electron microscopy (HRTEM) observation of graphene GBs only recently.³⁷ Neither us, nor any other group has reported a statistically relevant number of experimental observations of triple junctions in graphene. On the other hand, computer simulations can be used to directly measure the GB energy without resorting to the indirect method of triple junctions.³² Although there have been several numerical studies of graphene GBs, most of these studies have focused on a few special configurations or symmetric GBs, and have not explored the entire configuration space of graphene GBs.^38–48 Here we present a comprehensive numerical and theoretical study of the energy of a very large set of graphene GBs. We have used molecular dynamics (MD) simulations to measure the energy of about 79 000 grain boundary configurations, corresponding to 4122 unique (θ_M,θ_L) points, and spanning the entire space of all possible graphene grain boundaries.

Traditionally, the energy of low angle grain boundaries in bulk materials is understood in terms of Read and Shockley's theory.^49–51 This theory was developed for bulk materials, and is not directly applicable to two dimensional membranes. This is due to the fact that the individual dislocation cores in a two dimensional membrane can buckle out of plane – a mode of relaxation that is not available in their three dimensional counterparts. By buckling out of plane, the dislocation can trade in-plane strain for out of plane bending, thereby lowering its energy considerably.^38,52,53 We use a generalized Read–Shockley theory for two dimensional membranes, and combine it with the Frank–Bilby^50,54–57 equation to elucidate the structure of the energy function of all possible graphene GBs. Further, we develop a theoretical understanding of the salient structural features of graphene GBs. Our results should be applicable to a large class of 2D materials, and will lead to a better understanding of fundamental processes such as grain growth, transport, and strength in these materials.

2 Modeling

2.1 Configuration space of graphene GBs

Fig. 1a shows a general GB at the interface of two grain [capital G, script] and [capital G, script]′, with lattice vectors (v₁,v₂) and (v′₁,v′₂), respectively (the inset in Fig. 1b shows the lattice vectors and the primitive unit cell of graphene). GBs in graphene are characterized by two angles, namely, the misorientation angle θ_M, and the line angle θ_L. The misorientation angle characterizes the relative rotation of the two grains, i.e., v′_i = R_θMv_i, where R_θM represents a positive rotation by θ_M, whereas the line angle characterizes the deviation of the GB from the line of symmetry between the two grains. For any given (θ_M,θ_L) pair there is a third degree of freedom given by the relative sliding of the grains along the GB, however, we choose the sliding that gives the lowest GB energy, thereby effectively eliminating this degree of freedom. Due to the symmetries of the graphene lattice, the space of unique (θ_M,θ_L) pairs is reduced to a triangular area,³⁷ as shown in Fig. 1b. Commensurate GBs (CSL) exist at certain special values of (θ_M,θ_L), while an approximately commensurate GB can be constructed at any (θ_M,θ_L).^37,50,58 We simulate all commensurate GBs with a repeat length less than 2000 Å.³⁷ Further, we grid the (θ_M,θ_L) space in steps of 0.5° with approximately commensurate GBs. For each unique (θ_M,θ_L) pair, several simulations have to be performed to explore the relative sliding between the two grains; in all we have simulated and evaluated the energy of over 79 000 GB configurations corresponding to 4122 unique (θ_M,θ_L) pairs. The details of the GB configurations and structures used in this study can be found in ref. 37. The excess energy per-unit-length of a GB is calculated as γ(θ_M,θ_L) = (E_total − n_atomsE_bulk)/l_GB, where E_total is the net energy of the configuration, n_atoms is the number of atoms in the configuration, E_bulk is the energy per-atom in the reference crystal (=−7.81 eV for the AIREBO potential), and l_GB is the length of the GB. We use the AIREBO potential^59,60 as implemented in the LAMMPS code,⁶¹ and all our GB configurations are thoroughly relaxed, and allow for out of plane deformations (see the Methods section for details). The GB structures used in this study are available online.⁶²

Fig. 1 (a) Two grains [capital G, script], [capital G, script]′ and the GB interface (black line) between them. The half line angle, θ_L/2, measures the deviation of the interface from the symmetric GB (dashed black line). The Burger's vector needed to complete the interfacial Burgers circuit is shown by the red arrow. (b) The space of unique GBs in graphene. The colored (yellow) triangular region contains all unique (θ_M,θ_L) pairs (up to symmetry). The red circles and the black dots show all commensurate, and approximately commensurate GBs, respectively, with a repeat distance smaller than 2000 Å that were simulated in this study. The inset shows the lattice vectors and the primitive unit cell of graphene.

2.2 Dislocation model for GBs

Before discussing the GB structures and energy in detail, we present a dislocation based model for the GBs. This model will be used to elucidate the structure of the GBs, and to derive functional forms for the GB energy. The Frank–Bilby equation can be used to calculate the interfacial Burger's vector (per-unit-length) of the geometrically necessary dislocations for a GB with misorientation θ⁰_M and line angle θ⁰_L. This is the Burger's vector required to close the Burger's circuit shown in Fig. 1a (the red arrow), and is given by n⁰ =2 sin(θ^0′_M/2)(cos(θ⁰_L/2)e₁ − sin(θ⁰_L/2)e₂), where, θ^0′_M = θ⁰_M for 0 < θ⁰_M ≤ 30°, while θ^0′_M = 60° − θ⁰_M for 30° < θ⁰_M ≤ 60°, and, e_1,2 are the unit vector parallel and perpendicular to v₁, respectively (ESI Section S1†). Let the energy of this GB be γ⁰ = γ(θ⁰_M,θ⁰_L). Now consider a GB near this configuration, with the perturbed misorientation and line angle given by (θ⁰_M + δθ_M,θ⁰_L + δθ_L). The perturbation in the density of interfacial Burger's vector is given by δn = (∂n⁰/∂θ⁰_M)δθ_M + (∂n⁰/∂θ⁰_L)δθ_L. We assume that this change in Burger's vector density is accommodated by well separated (1,0) dislocations introduced along the boundary. There are three independent (1,0) dislocations in the reference crystal, with Burgers vectors in the directions R_{0°,60°,120°}e₁; thus δn = δn₁e₁ + δn₂R_60°e₁ + δn₃R_120°e₁, giving

δn₁e₁ + δn₂R_60°e₁ + δn₃R_120°e₁ = (∂n⁰/∂θ⁰_M)δθ_M + (∂n⁰/∂θ⁰_L)δθ_L, (1)

where δn₂ is the perturbation in the density of Burger's vector due to dislocations in the R_60°e₁ direction, etc. The above vector equation provides two constraints on the perturbation of the density of three independent (1,0) dislocations, thus leaving the system indeterminate. We obtain one more condition by writing a perturbed GB energy and minimizing it with the above constraints. Since the perturbation is small, the new dislocations introduced into the GB are well separated and do not interact. In the traditional Read–Shockley theory, the energy of an isolated dislocation has a divergent logarithmic term.^49,56 This term is due to the fact that, in a bulk material, the long range strain field of an isolated dislocation decays with distance as 1/r. However, it is known that in a two dimensional membrane the bending stiffness is small, and it is energetically favorable to trade long range strain for out of plane deformation, thereby removing the logarithmic term from the energy of the isolated dislocation core.^52,53 It can be shown that for two dimensional membranes, each isolated dislocation costs a finite (constant for a given GB) amount of energy.^52,53 Thus, we can write the following minimization problem for the perturbation δnwhere G is the shear modulus, μ is the Poisson's ratio, b is the Burger's vector, and c_i's are dimensionless constants for a given configuration (θ⁰_M,θ⁰_L) representing the energy required to embed the (1,0) dislocations into the GB. The validity of such a perturbational form for the GB energy has been tested numerically and experimentally.^51,63 If c_i's are known, then the minimization of the perturbed energy in eqn (2) with respect to the dislocation density δn_i can be performed analytically. However, the coefficients c_i's are unknown functions of (θ⁰_M,θ⁰_L), and thus a close form solution to the minimization problem is infeasible. Yet, we will show that considerable insight can be gained from this formulation.

3 Results and discussions

3.1 Structure and energy of symmetric GBs

Graphene GBs are known to be composed of rings of five and seven carbon atoms (apart from the usual hexagonal rings).^12,37,38 These pentagon–heptagon pairs form the cores of the dislocations with the shortest Burger's vector. Fig. 2a and b show the dislocation core with Burger's vector (1,0) and (1,0) + (0,1).³⁸ Fig. 2e–i show segments of symmetric GBs with θ_M = 10°, 20°, 30°, 40°, and 50°. It can be seen that the GBs are composed of (1,0) or (1,0) + (0,1) dislocation cores. The GBs with low misorientation (θ_M = 10°, 20°) are composed of (1,0) dislocations whose Burger's vectors are aligned in a single direction. Analogously, the GBs with large misorientation (θ_M = 50°) are composed of (1,0) dislocations whose Burger's vectors alternate by 60° in their orientation. As mentioned previously, the graphene lattice has three independent (1,0) dislocations whose Burger's vectors are rotated by 60° with respect to each other (the system b_1,2,3 in Fig. 2c and d, for example). Thus, while in principle a general symmetric GB can have (1,0) dislocations with three different orientations, the energy minimizing configurations of low angle symmetric GBs have their (1,0) dislocations aligned along just one direction; furthermore, the high angle symmetric GBs have dislocations aligned with two of the three possible directions. These structural features can be explained by considering the symmetric low and high angle GBs as perturbations about the pristine crystal (obtained at θ_M = 0°, 60°), and solving the minimization problem given by eqn (2) (details in ESI Section S2†). The essential insight is that for low angle GBs the perturbed interfacial Burger's vector is almost parallel to the lattice vector v₁ (or equivalently to the lattice Burger's vector b₁, see Fig. 2c). Thus, the energy minimizing configuration results when all the (1,0) cores are aligned with v₁ (equivalently, b₁). Aligning the core with R_60°e₁ instead, for example, would need twice the number of dislocations, and hence would cost twice the amount of energy, and thus would be suboptimal. Similarly as it can be seen that for high angle GBs the perturbed interfacial Burger's vector is almost perpendicular to v₁ (equivalently b₁, Fig. 2d), it is energetically not beneficial to have a dislocation Burger's vector aligned with v₁ (or b₁). Hence, the high angle GBs have dislocations with Burger's vectors aligned with b₂ and b₃ only.

Fig. 2 The structure and crystallography of graphene GBs. (a and b) show the dislocation core of isolated (1,0) and (1,0) + (0,1) dislocations, respectively. The color indicates the excess energy per atom in units of eV. (c and d) show isolated dislocations at a low angle (θ_M = 10°) and high angle (θ_M = 50°) symmetric GB. The shortest interfacial Burger's vector is indicated by n. (e)–(i) show sections of the lowest energy GBs at θ_M = 10°, 20°, 30°, 40°, 50°; color scheme same as (a and b).

Fig. 3a shows the numerically measured energy function for symmetric GBs, i.e., γ(θ_M,0). It is well known that the GB energy has cusps at special high symmetry (low Σ CSL, where CSL stands for the ‘Coincident Site Lattice’, and Σ denotes the ratio of the volume of the unit cell of the CSL to that of the regular lattice, see ref. 50 for a detailed discussion of CSL) boundaries,^50,64–67 and these can be seen clearly in Fig. 3a and 4a. There are two prominent cusps for graphene:³⁹ first at the Σ₇(θ_M = 21.78°, θ_L = 0°) GB, and the second at Σ₁₃(θ_M = 32.2°, θ_L = 0°) GB. Apart from these, there are the obvious families of cusp singularities at θ_M = 0°, 60°. The Σ_7,13 GBs are strong local energy minima of the GB energy. These minima arise due to favorable interactions between the dislocation cores. For intermediate values of misorientation (15° ≲ θ_M ≲ 45°), the density of required dislocations is high, and the individual cores cannot be well separated. We note that for isolated (1,0) as well as (1,0) + (0,1) cores, there is compression at the tip of the leading pentagonal ring and dilation at the tail of the trailing heptagonal ring (seen by the relative shortening and stretching of the bonds, most clearly visible in Fig. 2a and b).²³ This local straining leads to significant out of plane buckling near the dislocation, as seen in Fig. 3d.³⁸ However, as the dislocation density increases with increasing θ_M, and two (1,0) or (1,0) + (0,1) approach each other, their strain fields cancel, and there is a reduction in the elastic energy of the system. This cancellation of strain fields can be inferred from Fig. 3b and c, where it can be seen that the Σ_7,13 GBs have almost no out of plane buckling, because the strain fields cancel out very effectively in these GBs with tightly arranged dislocations. On the other hand, at higher θ_M (note that the dislocation density peaks at θ_M = 30°), the increased density of the dislocations leads to higher energy per-unit-length. Thus, there is a competition between the energy increase due to higher dislocation density, and energy decrease due to dislocation interaction. It can be seen that initially the GB energy increases with θ_M, thus the energy increase dominates over the energy reduction. However, the reduction becomes significant, and the net energy starts to decrease at about θ_M = 18°. This reduction in energy reaches a first optimum for the Σ₇ = (θ_M = 21.78°) CSL GB (Fig. 3b and e) where all (1,0) dislocation pairs are aligned, and there is a separation of exactly 1 carbon–carbon bond between them. At this optimal configuration there is significant reduction in the elastic energy, resulting in the first cusp in the GB energy (Fig. 3a). Increasing the misorientation θ_M further initially leads to an increase in the GB energy. This is due to the fact that a higher dislocation density pushes dislocations closer; however, geometrically it is still favorable to have all dislocations aligned in the same direction. Thus, creating a (1,0) + (0,1) pair incurs an energy penalty. However, with further increase in θ_M, it becomes progressively more favorable for the individual dislocations to stagger and merge to form (1,0) + (0,1) cores. This process leads to a reduction in energy starting at about θ_M = 24°. An optimal configuration is reached at the Σ₁₃ = (θ_M = 32.2°) CSL GB where all (1,0) + (0,1) line up perfectly (Fig. 3c and f), and leads to a large reduction in the elastic energy, resulting in the second, deeper cusp in the GB energy. On increasing θ_M further, the net dislocation density decreases and the (1,0) + (0,1) cores separate, ultimately resulting in the behavior for large misorientations discussed previously.

Fig. 3 (a) The energy of all simulated symmetric GBs in units of eV Å⁻¹. The filled black circles show the simulation data, while the solid line is a fit to eqn (3). The filled red circles show the magnitude of the fitting error, which is smaller than 0.06 eV Å⁻¹ everywhere. (b and c) show the Σ_7,13 GBs, and (d) shows a symmetric GB with θ_M = 10°, colored by the net out of plane displacement in unit of Å. (e)–(g) show the same GBs colored by excess energy per-atom in units of eV.

Fig. 4 (a) Energy of all simulated symmetric and asymmetric GBs in units of eV Å⁻¹. The surface shows data from simulation, while the grid is a fit to eqn (4). (b) The error in fit shown in units of eV Å⁻¹. (c) and (d) The basis functions, |h₁(θ_M,θ_L,θ_M^c) + h₂(θ_M,θ_L,θ_M^c)|, |h₁(θ_M,θ_L,θ_M^c)| + |h₂(θ_M,θ_L,θ_M^c)| used to fit the cusp singularity at θ_M^c = 21.78° in eqn (4).

Having understood the most salient features of the GB structure, we now turn our attention to the GB energy. A simple analysis of eqn (2) shows that for 2D materials the GB energy function has a absolute value (| |) type singularity at the cusps (ESI Section S2†). Thus, we propose the following functional form for the energy of the symmetric GBs:

where p_i, a_i are dimensionless fitting parameters. The overall factor of |sin 3θ_M| gives the correct asymptotic form at the cusps at θ_M = 0°, 60°. The first term (p_i's) fits the smooth variation in the energy function. The second term (a_i's) fits the cusps at angles . Although any desired number of cusps can be included, we include the two prominent cusps at the Σ_7,13 CSL GBs, thus n_c = 2, and . Note that this form satisfies all symmetry requirements, namely a period of 120°, and even reflection symmetries about θ_M = 0° (i.e., γ_sym(θ_M) = γ_sym(−θ_M) = γ_sym(120° + θ_M)), and has the correct asymptotic form near the cusp singularities. We do not include the n = 0, 1 terms because the corresponding harmonics are included in the expression for the cusps. The solid line in Fig. 3a shows a fit of eqn (3) to the simulation data with n = 4, giving a total of just 5 fitting parameters. The values of these parameters at the best fit are p₂ = −3.70 × 10⁻², p₃ = 6.18 × 10⁻³, p₄ = 1.99 × 10⁻², a₁ = 8.91 × 10⁻², a₂ = 2.00 × 10⁻¹, while G, μ are measured to be 325.68 GPa, and 0.318, respectively, from separate MD simulations (Methods section). The maximum absolute error for the fit is 0.026 eV Å⁻¹. It is clear that the theory provides an excellent fit to the data with a minimal number of fitting parameters. Including the higher harmonics (bigger n) does not result in a significant improvement in the results (ESI Fig. S1†). We find that fitting a Fourier series without including the correct asymptotic form of the cusps results in very poor performance; fits with as many as 50 free Fourier components are needed for an accuracy similar to our fit with 5 parameters (ESI Fig. S2†). Finally, we note that our functional form is reminiscent of the form used by Sethna and Coffman,⁵⁸ however, their form, while more pedagogical, had several redundant parameters (a total of 16 parameters, as opposed to our 5), and thus did not provide a minimal description of the GB energy.

3.2 Energy of asymmetric GBs

Fig. 4a shows the numerically measured energy function γ(θ_M,θ_L). It can be seen that the variation of γ(θ_M,θ_L) in the θ_L direction is significantly smaller than that in the θ_M direction. This is a direct consequence of the fact that the magnitude of the interfacial Burger's vector, n, is given by 2 sin(θ_M/2) and is independent of θ_L. Since the energy is largely a function of the magnitude of the interfacial Burger's vector, it follows that the energy variation in the θ_L direction is smaller. However, the cusps in the symmetric energy function γ_sym(θ_M,0°) become ridges in the full energy function γ(θ_M,θ_L). Thus, even though the energy function varies slowly in the θ_L direction, its structure is made interesting by the presence of these ridges, particularly due to the symmetry requirement γ(θ_M,θ_L) = γ(60° − θ_M,60° − θ_L) which makes the ridges turn (or vanish). The numerical data suggest that there are two kinds of ridges: one that join cusps at (θ⁰_M,0) to its periodic counterpart at (60° − θ⁰_M,60°), and another that continues almost straight in the θ_L direction without bending. Based on an analysis of eqn (2), we propose the following form for the general GB energy to captures these features (see ESI Section S2† for details)where p_ij, a^t_i, a^s_i are fitting parameters, I_ij ≡ (1 + (−1)^i+j)/2 is an indicator function that is 1 if both i, j are even or odd, and zero otherwise, , and . The indicator function is needed to make sure that the symmetry requirement γ(θ_M,θ_L) = γ(θ_M + 60°,θ_L + 60°) is satisfied. Note that the functional form satisfies all other symmetry requirements as well (overall period of 120° in θ_M,θ_L, even mirrors at θ_M,θ_L = 0° and 60°, i.e., γ(θ_M,θ_L) = γ(−θ_M,θ_L) = γ(θ_M,−θ_L) = γ(120° + θ_M,θ_L) = γ(θ_M,120° + θ_L) = γ(60° − θ_M,60° − θ_L)). We also set p₀₀ = p₁₀ = p₀₁ = p₁₁ = 0 because the constant term is not needed, and the other harmonics are contained in h′_is. The a^t_i terms model ridges that turn, while a^s_i terms model the ridges that remain straight. Taking n, m = 4, 4 and k = 2 as before, gives a fit with 15 free parameters, which is presented in Fig. 4a. The values of the parameters for the best fit can be found in the ESI Section S3.† The maximum absolute error of fitting is 0.07 eV Å⁻¹, indicating a good quality fit. Fig. 4b shows the fitting error. As before, adding further harmonics does not improve the quality of the fit significantly.

4 Conclusion

We find that the Read–Shockley type dislocation model provides an accurate description of the structure and energy of graphene GBs. The functional forms for energy derived on the basis of this formulation are numerically efficient, containing just 5 fitting parameters for the symmetric GB energy, and 15 fitting parameters for the entire GB space. The absolute error in our fits is smaller than 0.07 eV Å⁻¹ everywhere. We find that main source of this error is the limited size of our simulation cells (due to computational limitations). It can be seen in Fig. 4b that the largest error is concentrated in narrow bands around θ_M = 0°, 60°, 32.2°. These are high symmetry configurations, with θ_M = 0°, 60° being perfect crystals, and θ_M = 32.2° being the Σ₁₃ GB. The GBs vicinal to these high symmetry configurations have structures that are nominally the same as the high symmetry configuration, plus additional (or missing) dislocations separated by large distances (∼|b|/δθ_M). These well separated ‘perturbations’ produce out of plane distortions that need a very large cell to relax completely. Our simulations use a 1000 Å wide cell (in the direction perpendicular to the GB), and while we see significant decrease in energy over small cells (we have studied cells with widths of 50–1000 Å) due to relaxation, even the 1000 Å wide cell is not sufficiently large enough to fully relax the GB energy and obtain the infinite cell size limit. Note that this problem exists mostly for GBs that are vicinal to high symmetry configurations.

Our model for GB energy is based on a ‘small perturbation’ approach; however, all the evidence presented so far provides only indirect validation of the perturbation idea. The perturbation model can be supported by inspecting GBs in the vicinity of the high symmetry Σ_7,13 boundaries. We consider both symmetric and asymmetric perturbations, as shown in Fig. 5. This figure (plots a–d) shows that GBs in the vicinity of the Σ₇ GB have basically the same structure as Σ₇, plus an occasional extra (or missing) dislocation, as the case might be. The asymmetric perturbation ((δθ_L ≠ 0)) sometimes results in a faceted boundary (Fig. 5d and h), with the facet locally following the high symmetry GB. The kinks joining the facets are composed of extra dislocations that are not present in the high symmetry GB. The same observations are true for the GBs vicinal to the Σ₁₃ GB. It is remarkable that the GB generation algorithm is able to capture all the features expected from well annealed graphene GBs.

Fig. 5 Perturbations about the symmetric Σ_7,13 GBs. The top panel shows a zoom-in of the atoms in the white box in the panel directly below. Color shows the excess energy per-atom in units of eV. (a) Several repeats of the coincident site lattice (CSL) unit cell of the Σ₇ GB. (b)–(d) Perturbations about the Σ₇ GB with (δθ_M,δθ_L) = (−0.78°,0°), =(0.72°,0°), and =(0°,4.84°), respectively. (e) Several repeats of the CSL unit cell of the Σ₁₃ GB. (f)–(h) Perturbations about the Σ₇ GB with (δθ_M,δθ_L) = (−2.2°,0°), =(1.8°,0°), and =(0°,6.4°), respectively.

We have found that all of the approximately 79 000 lowest-energy configuration GBs that we have simulated consist of only pentagon–heptagon pairs, and the usual hexagonal rings. No other geometric configurations were observed for the lowest energy boundaries. For example, the 5-8-5 configurations that have been previously observed experimentally at grain boundaries^28,68,69 were not found in our structures. From an energy point of view, the 5-8-5 defects are vacancy defects and thus should be precluded from the ground state structures. However, non-equilibrium structures, such as the 5-8-5 defects, can indeed be captured by our algorithm if the Hamiltonian (eqn (2) in ref. 37) is not driven to its minima (the convergence criteria could be suitably relaxed, or Metropolis sampling could be performed at a suitably defined “temperature”). Further, the absence of such defects from our GBs is consistent with the fact that, to the best of our knowledge, such defects have not been observed in free-standing graphene films. Rather, they have been observed either in films on substrate or in free-standing graphene films after electron beam irradiation has modified their structure.⁷⁰ The ubiquity of pentagon–heptagon pairs in graphene grain boundaries is consistent with our HRTEM study of 176 boundaries,³⁷ the majority of which did not contain any rings of more than 7 or less than 5 carbon atoms. Further, the GB generation algorithm is able to capture faceting where appropriate.

To conclude, the main contribution of this work is to develop a fundamental understanding of the structure and energy of the entire space of graphene GBs. We have developed analytical expressions for GB energy as functions of the misorientation and line angle that can be readily used in future calculations of grain growth or other GB related phenomena.³² We hope that our analysis will pave the way for a deeper understanding of GB interfaces in graphene and other 2D materials.

5 Methods

The GB structures used in this study were generated by using the algorithm introduced in ref. 37, and are available online.⁶² In order to minimize boundary effects, we used GB structures that were 1000 Å wide, thus the GB was 500 Å from the boundaries. All GB structures used in this study are periodic in the direction parallel to the GB, and had a maximum length of 2000 Å. The atoms were allowed to relax in the out of plane direction without any constraint; a plot of the maximum out of plane displacement as a function of the misorientation and the line angle can be found in ESI Section S4.† The simulations were done in the LAMMPS⁶¹ code with the AIREBO interatomic potential.^59,60 Atoms in a strip of width 10 Å on the edges of the sample were held fixed at their ideal lattice positions during the simulations in order to reduce the boundary effects. Each sample was prepared by first relaxing the atoms with the conjugate gradient (CG) algorithm so that the force on each atom was less than 0.01 eV Å⁻¹. The sample was then held at 300 K for 10 picoseconds by simulating the NVT ensemble (fixed Number of atoms, Volume, and Temperature). This step was used to introduce any out of plane deformation that might have been missed by the CG algorithm. Finally, the atoms were again relaxed to within a residual force of 0.01 eV Å⁻¹ with the CG algorithm. The atoms at the strips on the edges were held fixed throughout these steps. As mentioned previously, the energy of the GB was measured as γ(θ_M,θ_L) = (E_total − n_atomsE_bulk)/l_GB, where E_total is the net energy of the configuration, n_atoms is the number of atoms in the configuration, E_bulk is the energy per-atom in the reference crystal (=−7.81 eV for the AIREBO potential), and l_GB is the length of the GB. The atoms that were held fixed and not allowed to relax were not included in the energy calculations. The linear elastic constants G and μ were calculated with the widely used technique of imposing small deformations on a relaxed bi-periodic graphene crystal, measuring the energy, and fitting the measured energy to the energy expression from linear elastic theory. This method yields G = 325.68 GPa, and μ = 0.318. The fits of the measured GB energy to eqn (3) and (4) were done with a standard least-squares algorithm.

Acknowledgements

A. S. thanks R. O. Ritchie for hosting him at LBNL, and the Miller Institute for Basic Research in Science at University of California Berkeley for providing financial support via the Miller Research Fellowship. The work was supported at LBNL by the Mechanical Properties of Materials Program (KC-13) funded by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences, Materials Sciences and Engineering Division under contract no. DE-AC02-05CH11231. C. O. acknowledges the Molecular Foundry, supported by the Office of Science, Office of Basic Energy Sciences, of the U.S. Department of Energy under Contract no. DE-AC02-05CH11231.

References

1. A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov and A. K. Geim, Rev. Mod. Phys., 2009, 81, 109–162.
2. A. K. Geim and K. S. Novoselov, Nat. Mater., 2007, 6, 183–191.
3. K. Novoselov, A. K. Geim, S. Morozov, D. Jiang, M. Katsnelson, I. Grigorieva, S. Dubonos and A. Firsov, Nature, 2005, 438, 197–200.
4. A. A. Balandin, S. Ghosh, W. Bao, I. Calizo, D. Teweldebrhan, F. Miao and C. N. Lau, Nano Lett., 2008, 8, 902–907.
5. X. Li, W. Cai, J. An, S. Kim, J. Nah, D. Yang, R. Piner, A. Velamakanni, I. Jung, E. Tutuc, S. K. Banerjee, L. Colombo and R. S. Ruoff, Science, 2009, 324, 1312–1314.
6. S. Bae, H. Kim, Y. Lee, X. Xu, J.-S. Park, Y. Zheng, J. Balakrishnan, T. Lei, H. R. Kim, Y. I. Song, Y.-J. Kim, K. S. Kim, B. ÃŰzyilmaz, J.-H. Ahn, B. H. Hong and S. Iijima, Nat. Nanotechnol., 2010, 5, 574–578.
7. Y. Liu, X. Dong and P. Chen, Chem. Soc. Rev., 2012, 41, 2283–2307.
8. F. Bonaccorso, Z. Sun, T. Hasan and A. Ferrari, Nat. Photonics, 2010, 4, 611–622.
9. X. Sun, Z. Liu, K. Welsher, J. Robinson, A. Goodwin, S. Zaric and H. Dai, Nano Res., 2008, 1, 203–212.
10. W. Ren and H.-M. Cheng, Nat. Nanotechnol., 2014, 9, 726–730.
11. O. V. Yazyev and Y. P. Chen, Nat. Nanotechnol., 2014, 9, 755–767.
12. P. Y. Huang, C. S. Ruiz-Vargas, A. M. van der Zande, W. S. Whitney, M. P. Levendorf, J. W. Kevek, S. Garg, J. S. Alden, C. J. Hustedt, Y. Zhu, J. P. Park, P. L. McEuen and D. A. Muller, Nature, 2011, 469, 389–392.
13. J. An, E. Voelkl, J. W. Suk, X. Li, C. W. Magnuson, L. Fu, P. Tiemeijer, M. Bischoff, B. Freitag, E. Popova and R. S. Ruoff, ACS Nano, 2011, 5, 2433–2439.
14. K. Kim, Z. Lee, W. Regan, C. Kisielowski, M. Crommie and A. Zettl, ACS Nano, 2011, 5, 2142–2146.
15. S. Kurasch, J. Kotakoski, O. Lehtinen, V. SkaÌĄkalovaÌĄ, J. Smet, C. E. Krill III, A. V. Krasheninnikov and U. Kaiser, Nano Lett., 2012, 12, 3168–3173.
16. H. I. Rasool, C. Ophus, Z. Zhang, M. F. Crommie, B. I. Yakobson and A. Zettl, Nano Lett., 2014, 14, 7057–7063.
17. Y. Tison, J. Lagoute, V. Repain, C. Chacon, Y. Girard, F. Joucken, R. Sporken, F. Gargiulo, O. V. Yazyev and S. Rousset, Nano Lett., 2014, 14, 6382–6386.
18. O. V. Yazyev and S. G. Louie, Nat. Mater., 2010, 9, 806–809.
19. A. Bagri, S.-P. Kim, R. S. Ruoff and V. B. Shenoy, Nano Lett., 2011, 11, 3917–3921.
20. A. Y. Serov, Z.-Y. Ong and E. Pop, Appl. Phys. Lett., 2013, 102, 033104.
21. K. L. Grosse, V. E. Dorgan, D. Estrada, J. D. Wood, I. Vlassiouk, G. Eres, J. W. Lyding, W. P. King and E. Pop, Appl. Phys. Lett., 2014, 105, 143109.
22. P. Yasaei, A. Fathizadeh, R. Hantehzadeh, A. K. Majee, A. El-Ghandour, D. Estrada, C. Foster, Z. Aksamija, F. Khalili-Araghi and A. Salehi-Khojin, Nano Lett., 2015, 5, 4532–4540.
23. Y. Wei, J. Wu, H. Yin, X. Shi, R. Yang and M. Dresselhaus, Nat. Mater., 2012, 11, 759–763.
24. L. Yi, Z. Yin, Y. Zhang and T. Chang, Carbon, 2013, 51, 373–380.
25. J. Zhang, J. Zhao and J. Lu, ACS Nano, 2012, 6, 2704–2711.
26. R. Grantab, V. B. Shenoy and R. S. Ruoff, Science, 2010, 330, 946–948.
27. J. C. Koepke, J. D. Wood, D. Estrada, Z.-Y. Ong, K. T. He, E. Pop and J. W. Lyding, ACS Nano, 2013, 7, 75–86.
28. C. Ma, H. Sun, Y. Zhao, B. Li, Q. Li, A. Zhao, X. Wang, Y. Luo, J. Yang, B. Wang and J. Hou, Phys. Rev. Lett., 2014, 112, 226802.
29. C. Ma, H. Sun, H. Du, J. Wang, A. Zhao, Q. Li, B. Wang and J. Hou, Nanoscale, 2015, 7, 3055–3059.
30. A. W. Cummings, D. L. Duong, V. L. Nguyen, D. Van Tuan, J. Kotakoski, J. E. Barrios Vargas, Y. H. Lee and S. Roche, Adv. Mater., 2014, 26, 5079–5094.
31. K. W. Clark, X.-G. Zhang, I. V. Vlassiouk, G. He, R. M. Feenstra and A.-P. Li, ACS Nano, 2013, 7, 7956–7966.
32. Y. Mishin, M. Asta and J. Li, Acta Mater., 2010, 58, 1117–1151.
33. A. Rollett, D. Srolovitz and M. Anderson, Acta Metall., 1989, 37, 1227–1240.
34. C. Herring, Phys. Rev., 1951, 82, 87–93.
35. G. Hasson, J.-Y. Boos, I. Herbeuval, M. Biscondi and C. Goux, Surf. Sci., 1972, 31, 115–137.
36. B. L. Adams, S. Ta'Asan, D. Kinderlehrer, I. Livshits, D. Mason, C.-T. Wu, W. Mullins, G. Rohrer, A. Rollett and D. Saylor, Interface Sci., 1999, 7, 321–337.
37. C. Ophus, A. Shekhawat, H. Rasool and A. Zettl, Phys. Rev. B: Condens. Matter Mater. Phys., 2015, 92, 205402.
38. O. V. Yazyev and S. G. Louie, Phys. Rev. B: Condens. Matter Mater. Phys., 2010, 81, 195420.
39. T.-H. Liu, G. Gajewski, C.-W. Pao and C.-C. Chang, Carbon, 2011, 49, 2306–2317.
40. J. M. Carlsson, L. M. Ghiringhelli and A. Fasolino, Phys. Rev. B: Condens. Matter Mater. Phys., 2011, 84, 165423.
41. Y. Liu and B. I. Yakobson, Nano Lett., 2010, 10, 2178–2183.
42. S. Malola, H. Häkkinen and P. Koskinen, Phys. Rev. B: Condens. Matter Mater. Phys., 2010, 81, 165447.
43. J. Zhang and J. Zhao, Carbon, 2013, 55, 151–159.
44. A. Romanov, A. Kolesnikova, T. Orlova, I. Hussainova, V. Bougrov and R. Valiev, Carbon, 2015, 81, 223–231.
45. Z. Zhang, Y. Yang, F. Xu, L. Wang and B. I. Yakobson, Adv. Funct. Mater., 2015, 25, 367–373.
46. X. Zhang, Z. Xu, Q. Yuan, J. Xin and F. Ding, Nanoscale, 2015, 7, 20082–20088.
47. J. Han, S. Ryu, D. Sohn and S. Im, Carbon, 2014, 68, 250–257.
48. Z.-L. Li, Z.-M. Li, H.-Y. Cao, J.-H. Yang, Q. Shu, Y.-Y. Zhang, H. J. Xiang and X. G. Gong, Nanoscale, 2014, 6, 4309–4315.
49. W. Read and W. Shockley, Phys. Rev., 1950, 78, 275.
50. A. P. Sutton and R. W. Balluffi, Interfaces in crystalline materials, Clarendon Press, 1995.
51. W. Gui-Jin and V. Vitek, Acta Metall., 1986, 34, 951–960.
52. H. S. Seung and D. R. Nelson, Phys. Rev. A: At., Mol., Opt. Phys., 1988, 38, 1005–1018.
53. S. Chen and D. Chrzan, Phys. Rev. B: Condens. Matter Mater. Phys., 2011, 84, 214103.
54. F. Frank, Acta Metall., 1953, 1, 15–21.
55. B. Bilby, R. Bullough and E. Smith, Proc. R. Soc. London, Ser. A, 1955, 263–273.
56. D. Hull and D. J. Bacon, Introduction to dislocations, Pergamon Press Oxford, 1984.
57. J. Hirth, R. Pond, R. Hoagland, X.-Y. Liu and J. Wang, Prog. Mater. Sci., 2013, 58, 749–823.
58. V. R. Coffman and J. P. Sethna, Phys. Rev. B: Condens. Matter Mater. Phys., 2008, 77, 144111.
59. S. J. Stuart, A. B. Tutein and J. A. Harrison, J. Chem. Phys., 2000, 112, 6472–6486.
60. D. W. Brenner, O. A. Shenderova, J. A. Harrison, S. J. Stuart, B. Ni and S. B. Sinnott, J. Phys.: Condens. Matter, 2002, 14, 783.
61. S. Plimpton, J. Comput. Phys., 1995, 117, 1–19.
62. A. Shekhawat and C. Ophus, https://sites.google.com/a/lbl.gov/atomic-structure-repository/.
63. N. A. Gjostein and F. Rhines, Acta Metall., 1959, 7, 319–330.
64. D. Wolf, Acta Metall., 1990, 38, 781–790.
65. D. M. Duffy and P. W. Tasker, Philos. Mag. A, 1983, 47, 817–825.
66. S. Tsurekawa, T. Tanaka and H. Yoshinaga, Mater. Sci. Eng., A, 1994, 176, 341–348.
67. V. V. Bulatov, B. W. Reed and M. Kumar, Acta Mater., 2014, 65, 161–175.
68. J. Lahiri, Y. Lin, P. Bozkurt, I. Oleynik and M. Batzill, Nat. Nanotechnol., 2010, 5, 326–329.
69. B. Yang, H. Xu, J. Lu and K. Loh, J. Am. Chem. Soc., 2014, 136, 12041–12046.
70. J. Kotakoski, J. Meyer, S. Kurasch, D. Santos-Cottin, U. Kaiser and A. Krasheninnikov, Phys. Rev. B: Condens. Matter Mater. Phys., 2011, 83, 245420.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/c6ra07584c

---