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

Ashivni Shekhawat*abc, Colin Ophusd and Robert O. Ritchieac
aDepartment of Materials Science and Engineering, University of California Berkeley, USA. E-mail: shekhawat.ashivni@gmail.com
bMiller Institute for Basic Research in Science, University of California Berkeley, USA
cMaterials Sciences Division, Lawrence Berkeley National Laboratory, Berkeley, USA
dNational Center for Electron Microscopy, Molecular Foundry, Lawrence Berkeley National Laboratory, Berkeley, USA

Received 23rd March 2016 , Accepted 27th April 2016

First published on 29th April 2016

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[thin space (1/6-em)]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[thin space (1/6-em)]000 m2 of graphene annually.5,6,10 Graphene produced from CVD is polycrystalline,11 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.18 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,34 which can be used to deduce the GB energy if a statistically large amount of experimental data is available.36 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.37 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.32 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[thin space (1/6-em)]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–Bilby50,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 (v1,v2) and (v1,v2), 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., vi = RθMvi, 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,37 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 Å.37 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[thin space (1/6-em)]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) = (EtotalnatomsEbulk)/lGB, where Etotal is the net energy of the configuration, natoms is the number of atoms in the configuration, Ebulk is the energy per-atom in the reference crystal (=−7.81 eV for the AIREBO potential), and lGB is the length of the GB. We use the AIREBO potential59,60 as implemented in the LAMMPS code,61 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.62
image file: c6ra07584c-f1.tif
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 θ0M and line angle θ0L. This is the Burger's vector required to close the Burger's circuit shown in Fig. 1a (the red arrow), and is given by n0 =2[thin space (1/6-em)]sin(θ0′M/2)(cos(θ0L/2)e1 − sin(θ0L/2)e2), where, θ0′M = θ0M for 0 < θ0M ≤ 30°, while θ0′M = 60° − θ0M for 30° < θ0M ≤ 60°, and, e1,2 are the unit vector parallel and perpendicular to v1, respectively (ESI Section S1). Let the energy of this GB be γ0 = γ(θ0M,θ0L). Now consider a GB near this configuration, with the perturbed misorientation and line angle given by (θ0M + δθM,θ0L + δθL). The perturbation in the density of interfacial Burger's vector is given by δn = (∂n0/∂θ0M)δθM + (∂n0/∂θ0L)δθ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 R0°,60°,120°e1; thus δn = δn1e1 + δn2R60°e1 + δn3R120°e1, giving
δn1e1 + δn2R60°e1 + δn3R120°e1 = (∂n0/∂θ0M)δθM + (∂n0/∂θ0L)δθL, (1)
where δn2 is the perturbation in the density of Burger's vector due to dislocations in the R60°e1 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 δn
image file: c6ra07584c-t1.tif(2)
where G is the shear modulus, μ is the Poisson's ratio, b is the Burger's vector, and ci's are dimensionless constants for a given configuration (θ0M,θ0L) 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 ci's are known, then the minimization of the perturbed energy in eqn (2) with respect to the dislocation density δni can be performed analytically. However, the coefficients ci's are unknown functions of (θ0M,θ0L), 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).38 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 b1,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 v1 (or equivalently to the lattice Burger's vector b1, see Fig. 2c). Thus, the energy minimizing configuration results when all the (1,0) cores are aligned with v1 (equivalently, b1). Aligning the core with R60°e1 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 v1 (equivalently b1, Fig. 2d), it is energetically not beneficial to have a dislocation Burger's vector aligned with v1 (or b1). Hence, the high angle GBs have dislocations with Burger's vectors aligned with b2 and b3 only.
image file: c6ra07584c-f2.tif
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:39 first at the Σ7(θM = 21.78°, θL = 0°) GB, and the second at Σ13(θ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).23 This local straining leads to significant out of plane buckling near the dislocation, as seen in Fig. 3d.38 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 Σ7 = (θ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 Σ13 = (θ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.

image file: c6ra07584c-f3.tif
Fig. 3 (a) The energy of all simulated symmetric GBs in units of eV Å−1. 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 Å−1 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.

image file: c6ra07584c-f4.tif
Fig. 4 (a) Energy of all simulated symmetric and asymmetric GBs in units of eV Å−1. 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 Å−1. (c) and (d) The basis functions, |h1(θM,θL,θMc) + h2(θM,θL,θMc)|, |h1(θM,θL,θMc)| + |h2(θM,θL,θMc)| used to fit the cusp singularity at θMc = 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:

image file: c6ra07584c-t2.tif(3)
where pi, ai are dimensionless fitting parameters. The overall factor of |sin[thin space (1/6-em)]3θM| gives the correct asymptotic form at the cusps at θM = 0°, 60°. The first term (pi's) fits the smooth variation in the energy function. The second term (ai's) fits the cusps at angles image file: c6ra07584c-t3.tif. Although any desired number of cusps can be included, we include the two prominent cusps at the Σ7,13 CSL GBs, thus nc = 2, and image file: c6ra07584c-t4.tif. 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 p2 = −3.70 × 10−2, p3 = 6.18 × 10−3, p4 = 1.99 × 10−2, a1 = 8.91 × 10−2, a2 = 2.00 × 10−1, 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 Å−1. 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,58 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[thin space (1/6-em)]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 (θ0M,0) to its periodic counterpart at (60° − θ0M,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)
image file: c6ra07584c-t5.tif(4)
where pij, ati, asi are fitting parameters, Iij ≡ (1 + (−1)i+j)/2 is an indicator function that is 1 if both i, j are even or odd, and zero otherwise, image file: c6ra07584c-t6.tif, and image file: c6ra07584c-t7.tif. 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 p00 = p10 = p01 = p11 = 0 because the constant term is not needed, and the other harmonics are contained in his. The ati terms model ridges that turn, while asi 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 Å−1, 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 Å−1 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 Σ13 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 Σ7 GB have basically the same structure as Σ7, 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 Σ13 GB. It is remarkable that the GB generation algorithm is able to capture all the features expected from well annealed graphene GBs.

image file: c6ra07584c-f5.tif
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 Σ7 GB. (b)–(d) Perturbations about the Σ7 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 Σ13 GB. (f)–(h) Perturbations about the Σ7 GB with (δθM,δθL) = (−2.2°,0°), =(1.8°,0°), and =(0°,6.4°), respectively.

We have found that all of the approximately 79[thin space (1/6-em)]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 boundaries28,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.70 The ubiquity of pentagon–heptagon pairs in graphene grain boundaries is consistent with our HRTEM study of 176 boundaries,37 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.32 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.62 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 LAMMPS61 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 Å−1. 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 Å−1 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) = (EtotalnatomsEbulk)/lGB, where Etotal is the net energy of the configuration, natoms is the number of atoms in the configuration, Ebulk is the energy per-atom in the reference crystal (=−7.81 eV for the AIREBO potential), and lGB 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.


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.


  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 CrossRef CAS.
  2. A. K. Geim and K. S. Novoselov, Nat. Mater., 2007, 6, 183–191 CrossRef CAS PubMed.
  3. K. Novoselov, A. K. Geim, S. Morozov, D. Jiang, M. Katsnelson, I. Grigorieva, S. Dubonos and A. Firsov, Nature, 2005, 438, 197–200 CrossRef CAS PubMed.
  4. A. A. Balandin, S. Ghosh, W. Bao, I. Calizo, D. Teweldebrhan, F. Miao and C. N. Lau, Nano Lett., 2008, 8, 902–907 CrossRef CAS PubMed.
  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 CrossRef CAS PubMed.
  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 CrossRef CAS PubMed.
  7. Y. Liu, X. Dong and P. Chen, Chem. Soc. Rev., 2012, 41, 2283–2307 RSC.
  8. F. Bonaccorso, Z. Sun, T. Hasan and A. Ferrari, Nat. Photonics, 2010, 4, 611–622 CrossRef CAS.
  9. X. Sun, Z. Liu, K. Welsher, J. Robinson, A. Goodwin, S. Zaric and H. Dai, Nano Res., 2008, 1, 203–212 CrossRef CAS PubMed.
  10. W. Ren and H.-M. Cheng, Nat. Nanotechnol., 2014, 9, 726–730 CrossRef CAS PubMed.
  11. O. V. Yazyev and Y. P. Chen, Nat. Nanotechnol., 2014, 9, 755–767 CrossRef CAS PubMed.
  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 CrossRef CAS PubMed.
  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 CrossRef CAS PubMed.
  14. K. Kim, Z. Lee, W. Regan, C. Kisielowski, M. Crommie and A. Zettl, ACS Nano, 2011, 5, 2142–2146 CrossRef CAS PubMed.
  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 CrossRef CAS PubMed.
  16. H. I. Rasool, C. Ophus, Z. Zhang, M. F. Crommie, B. I. Yakobson and A. Zettl, Nano Lett., 2014, 14, 7057–7063 CrossRef CAS PubMed.
  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 CrossRef CAS PubMed.
  18. O. V. Yazyev and S. G. Louie, Nat. Mater., 2010, 9, 806–809 CrossRef CAS PubMed.
  19. A. Bagri, S.-P. Kim, R. S. Ruoff and V. B. Shenoy, Nano Lett., 2011, 11, 3917–3921 CrossRef CAS PubMed.
  20. A. Y. Serov, Z.-Y. Ong and E. Pop, Appl. Phys. Lett., 2013, 102, 033104 CrossRef.
  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 CrossRef.
  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 CrossRef PubMed.
  23. Y. Wei, J. Wu, H. Yin, X. Shi, R. Yang and M. Dresselhaus, Nat. Mater., 2012, 11, 759–763 CrossRef CAS PubMed.
  24. L. Yi, Z. Yin, Y. Zhang and T. Chang, Carbon, 2013, 51, 373–380 CrossRef CAS.
  25. J. Zhang, J. Zhao and J. Lu, ACS Nano, 2012, 6, 2704–2711 CrossRef CAS PubMed.
  26. R. Grantab, V. B. Shenoy and R. S. Ruoff, Science, 2010, 330, 946–948 CrossRef CAS PubMed.
  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 CrossRef CAS PubMed.
  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 CrossRef PubMed.
  29. C. Ma, H. Sun, H. Du, J. Wang, A. Zhao, Q. Li, B. Wang and J. Hou, Nanoscale, 2015, 7, 3055–3059 RSC.
  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 CrossRef CAS PubMed.
  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 CrossRef CAS PubMed.
  32. Y. Mishin, M. Asta and J. Li, Acta Mater., 2010, 58, 1117–1151 CrossRef CAS.
  33. A. Rollett, D. Srolovitz and M. Anderson, Acta Metall., 1989, 37, 1227–1240 CrossRef CAS.
  34. C. Herring, Phys. Rev., 1951, 82, 87–93 CrossRef CAS.
  35. G. Hasson, J.-Y. Boos, I. Herbeuval, M. Biscondi and C. Goux, Surf. Sci., 1972, 31, 115–137 CrossRef CAS.
  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 CrossRef CAS.
  37. C. Ophus, A. Shekhawat, H. Rasool and A. Zettl, Phys. Rev. B: Condens. Matter Mater. Phys., 2015, 92, 205402 CrossRef.
  38. O. V. Yazyev and S. G. Louie, Phys. Rev. B: Condens. Matter Mater. Phys., 2010, 81, 195420 CrossRef.
  39. T.-H. Liu, G. Gajewski, C.-W. Pao and C.-C. Chang, Carbon, 2011, 49, 2306–2317 CrossRef CAS.
  40. J. M. Carlsson, L. M. Ghiringhelli and A. Fasolino, Phys. Rev. B: Condens. Matter Mater. Phys., 2011, 84, 165423 CrossRef.
  41. Y. Liu and B. I. Yakobson, Nano Lett., 2010, 10, 2178–2183 CrossRef CAS PubMed.
  42. S. Malola, H. Häkkinen and P. Koskinen, Phys. Rev. B: Condens. Matter Mater. Phys., 2010, 81, 165447 CrossRef.
  43. J. Zhang and J. Zhao, Carbon, 2013, 55, 151–159 CrossRef CAS.
  44. A. Romanov, A. Kolesnikova, T. Orlova, I. Hussainova, V. Bougrov and R. Valiev, Carbon, 2015, 81, 223–231 CrossRef CAS.
  45. Z. Zhang, Y. Yang, F. Xu, L. Wang and B. I. Yakobson, Adv. Funct. Mater., 2015, 25, 367–373 CrossRef CAS.
  46. X. Zhang, Z. Xu, Q. Yuan, J. Xin and F. Ding, Nanoscale, 2015, 7, 20082–20088 RSC.
  47. J. Han, S. Ryu, D. Sohn and S. Im, Carbon, 2014, 68, 250–257 CrossRef CAS.
  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 RSC.
  49. W. Read and W. Shockley, Phys. Rev., 1950, 78, 275 CrossRef CAS.
  50. A. P. Sutton and R. W. Balluffi, Interfaces in crystalline materials, Clarendon Press, 1995 Search PubMed.
  51. W. Gui-Jin and V. Vitek, Acta Metall., 1986, 34, 951–960 CrossRef.
  52. H. S. Seung and D. R. Nelson, Phys. Rev. A: At., Mol., Opt. Phys., 1988, 38, 1005–1018 CrossRef.
  53. S. Chen and D. Chrzan, Phys. Rev. B: Condens. Matter Mater. Phys., 2011, 84, 214103 CrossRef.
  54. F. Frank, Acta Metall., 1953, 1, 15–21 CrossRef CAS.
  55. B. Bilby, R. Bullough and E. Smith, Proc. R. Soc. London, Ser. A, 1955, 263–273 CrossRef CAS.
  56. D. Hull and D. J. Bacon, Introduction to dislocations, Pergamon Press Oxford, 1984 Search PubMed.
  57. J. Hirth, R. Pond, R. Hoagland, X.-Y. Liu and J. Wang, Prog. Mater. Sci., 2013, 58, 749–823 CrossRef.
  58. V. R. Coffman and J. P. Sethna, Phys. Rev. B: Condens. Matter Mater. Phys., 2008, 77, 144111 CrossRef.
  59. S. J. Stuart, A. B. Tutein and J. A. Harrison, J. Chem. Phys., 2000, 112, 6472–6486 CrossRef CAS.
  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 CrossRef CAS.
  61. S. Plimpton, J. Comput. Phys., 1995, 117, 1–19 CrossRef CAS.
  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 CrossRef CAS.
  64. D. Wolf, Acta Metall., 1990, 38, 781–790 CrossRef CAS.
  65. D. M. Duffy and P. W. Tasker, Philos. Mag. A, 1983, 47, 817–825 CrossRef CAS.
  66. S. Tsurekawa, T. Tanaka and H. Yoshinaga, Mater. Sci. Eng., A, 1994, 176, 341–348 CrossRef CAS.
  67. V. V. Bulatov, B. W. Reed and M. Kumar, Acta Mater., 2014, 65, 161–175 CrossRef CAS.
  68. J. Lahiri, Y. Lin, P. Bozkurt, I. Oleynik and M. Batzill, Nat. Nanotechnol., 2010, 5, 326–329 CrossRef CAS PubMed.
  69. B. Yang, H. Xu, J. Lu and K. Loh, J. Am. Chem. Soc., 2014, 136, 12041–12046 CrossRef CAS PubMed.
  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 CrossRef.


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

This journal is © The Royal Society of Chemistry 2016