Structural, chemical and dynamical trends in graphene grain boundaries

Grain boundaries are topological defects that often have a disordered character. Disorder implies that understanding general trends is more important than accurate investigations of individual grain boundaries. Here we present trends in the grain boundaries of graphene. We use density-functional tight-binding method to calculate trends in energy, atomic structure (polygon composition), chemical reactivity (dangling bond density), corrugation heights (inflection angles), and dynamical properties (vibrations), as a function of lattice orientation mismatch. The observed trends and their mutual interrelations are plausibly explained by structure, and supported by past experiments.


I. INTRODUCTION
Real graphene has always defects, edges, point defects, chemical impurities-and grain boundaries. Grain boundaries (GBs) are topological defects, trails of disorder that separate two pieces of pristine hexagonal graphene sheets. They reside practically in any graphitic material, graphite 1-3 , soot 4-6 , single-and multi-layer graphene 7,8 , fullerenes 9,10 , carbon nanotubes 11,12 , with varying abundance. In graphene fabrication, for instance, chemically synthesized samples tend to have more GBs than mechanically cleaved samples 13 . For electron mobility high purity is important and GBs are best avoided, while other properties may tolerate defects better. We refer here to strictly two-dimensional GBs (albeit not necessarily planar) which should not be confused with genuinely three-dimensional GBs in graphite.
But not are defects, impurities and GBs always something bad; they are interesting also in their own righteven useful 13,14 . After all, being extended, GBs modify graphene more than point defects. They affect graphene's magnetic, electronic, structural, and mechanical properties 13,15 .
Considering the relevance and prevalence of graphene GBs, they ought to deserve more attention. Most related studies are on GBs in graphite and moiré patterns therein 3,16-18 , or on intramolecular junctions in carbon nanotubes [19][20][21] ; GBs in single-or few-layer graphene, whether experimental 2,7,22 or theoretical 23 , are still scarce. Often graphitic grain boundaries are idealized by pentagon-heptagon pairs 18,23 , where the pair distances determine lattice mismatch. While this ideal model works in certain occasions 18,23 , it does not work for GBs that are rough, corrugated, have ridges 16 , andmost importantly-are decorated by dangling bonds or reactive sites alike 2 .
In this paper we investigate graphene GBs beyond simple models. We explore computationally the trends in * Corresponding author structural properties, chemical reactivity and vibrational properties for an ensemble of GBs that are not ideal, but have certain roughness. What this roughness really means will be clarified in Sec.III.

II. COMPUTATIONAL METHOD
The electronic structure method we use to simulate the GBs is density-functional tight-binding (DFTB) [24][25][26] , and the hotbit code 27 . DFTB models well the covalent bonding in carbon [28][29][30] , and suits fine for our simulations that concentrates on trends. We leave the details of the method to Ref. 26 and only mention-this will be used later-that the total energy in DFTB, including the band-structure energy, the Coulomb energy, and the short-range repulsion, can be expressed as a sum of atomic binding contributions, Here N is the atom count and i atom i's contribution to cohesion energy; in pristine graphene i ≡ gr = −9.6 eV for all atoms (DFTB has some overbinding). The quantity i allows us to calculate energies in a local fashion: for a group of carbon atoms in a given zone S (set of atoms), the sum measures how much the energy of zone S is larger than the energy of equal-sized piece of graphene. In addition, since gr = −9.6 eV comes from 3 equivalent bonds, the value i = 2/3 · gr ≈ −6.4 eV infers a dangling bond for atom i. Hence, albeit less accurate than densityfunctional theory, DFTB enables straightforward analysis of local quantities, in addition to fast calculations and extensive sample statistics.

2θ 2θ
FIG. 1: (color online) Grain boundary sample, one out of 48, made by merging two graphene edges with the same chirality, but different orientation and random translation along the boundary; our samples are hence characterized by chirality. In this figure chiral indices are (11,6), chiral angle θ = 20.4 • . Dashed vertical lines stand for periodicity across the horizontal direction.

III. SAMPLE CONSTRUCTION
To investigate trends, we constructed an ensemble of 48 graphene GBs using a procedure we now describe. (i) Cut two ribbons out of perfect graphene in a given chirality direction, with random offset. Chirality index (n, m) means we cut graphene along the vector C = na 1 + ma 2 , where a i are the primitive lattice vectors. The cut direction can be expressed, equivalently, vy the chiral angle θ = tan −1 [ √ 3m/(m + 2n)]. (ii) Passivate ribbons' other edges with hydrogen, leaving carbon atoms on the other edges bare. (iii) Thermalize the ribbons, still separated, to 1500 K using Langevin molecular dynamics (MD). (iv) Perform reflection operation on the other ribbon with respect to a plane along ribbon axis and perpendicular to ribbon plane. (v) Perform random translations along the ribbon direction. (vi) Merge the two ribbons' bare carbon edges with MD at 1500 K. (vii) Cool the system to room temperature within ∼ 1 ps. (vii) Optimize the geometry using the FIRE method 31 , and finally obtain structures like the one in Fig.1. The total process takes about 3 ps.
Our ensemble of 48 GB samples are chosen to get representatives for all chiral angles from θ = 0 • to θ = 30 • . Our limit for the width in cutting the ribbons is 10Å. Ribbons have to be wide enough (> 9Å) for faithful description of a semi-infinite graphene, and, by limiting the total number of atoms below 170, also the periodicity of GB becomes limited below 40Å. (Edges with θ ∼ 0 • and θ ∼ 30 • , that would have short periodicity, are scarce.) Our procedure to make GBs is neither the only nor the best one, but is suitable for trend-hunting.
Sure enough, it's a downside that the procedure excludes GBs from ribbons with different edge chirality, such as merged zigzag and armchair edges. Different chirality, unfortunately, would mean different periodicity and practical problems. To have a single number, θ, to identify GBs is vital for finding the trends we concentrate on. For GBs identified by two numbers, θ 1 and θ 2 , trend-hunting would be harder; those GBs can be investigated afterwards, using the insights we learn here. Further, while the temperature 1500 K, motivated partly by experiments 16 , allows rearrangements in merging, the heat of fusion renders initial temperatures irrelevant. We use canonical MD in the merging to allow GB formation process to be steered by intrinsic driving forces and energetics.
The time scale of the process is limited, as usual, in part by computational constraints. However, prolonging the process would cause no fundamental structural changes as the merging process itself is nearly instantaneous; there is no diffusion after merging and saturation of the dangling bonds, there is only annealing of the worst local defects and stress-release that causes buckling (Sec.VI). The usual time scales in graphene growth are hence not relevant to our GB formation process. Only intrinsic roughness remains in the GB samples-roughness related to chiral angle and randomness in cutting offsets and translations.
Usually all atoms in ribbons were also bound to GBs, but in few structures (near armchair chirality) the merging process squeezed dimers out, and either inflected them out of plane or detached them completely. Since ribbons found these own ways to optimize geometriesways not anticipated beforehand-it suggests that the design of our construction process is not dominant.

IV. PRELUDE: ENERGIES FOR FREE GRAPHENE EDGES
Before going into GBs themselves, let us look at the free bare edges of graphene. Fig.2 shows the edge energies as a function of the chiral angle, θ = 0 • meaning zigzag and θ = 30 • armchair edges. To ignore the effect of hydrogen passivated edge, carbons bound to hydrogen are neglected in edge energy calculation, The edge energy between zigzag and armchair varies linearly; fluctuations in energy are due to random offset in the cut, occasionally producing pentagons. The edge energy in zigzag is high due to strong and unhappy dangling bonds; in armchair the dangling bonds are weakened by the formation of triple bonds in the armrest parts 30,32 .
It was recently predicted theoretically 32 and later confirmed experimentally 33 that the zigzag edge is actually metastable, and prefers reconstruction into pentagons and heptagons at the edge, forming a so-called reczag ac zz FIG. 2: (color online) Energies of the free edges of graphene as a function of chiral angle; zigzag edge has θ = 0 • and armchair θ = 30 • . Free zigzag edge segments are metastable and prefer local reconstruction (two adjacent hexagons reconstruct into a pentagon and a heptagon); these so-called reczag edges are shown with gray connected symbols. The reczag segments are, however, important only for free edges and irrelevant for grain boundaries studied here.
edge. As it turned out, reczag is energetically even better than armchair. Hence, if we reconstruct the zigzag segments in edges with small θ, we usually lower the edge energy, as seen in Fig.2 where reconstructions are added by hand. Reczag edges are, however, irrelevant for our GBs where edges are not free, and are ignored because we want the dangling bonds to spontaneously find contact from the other merging edge. The edge energies were investigated and presented here for comparison with DFT calculations 32 . The accuracy in edge energies is better than 10 %, and we expect same accuracy in GB energetics.

V. TRENDS IN ENERGY AND STRUCTURE
The GB energy per unit length is measuring how much GB costs energy relative to the same number of carbon atoms in graphene. The energies for the whole ensemble of GBs are shown in Fig.3, as a function of the chiral angle. The value ε GB = 0.0 eV/Å appears only with θ = 0 • and θ = 30 • -meaning pristine graphene.
Most GB energies are less than or equal to the edge energies of the free graphene edges, albeit with variation. This means that GBs regain, on average, other ribbon's edge energy during the merging, and the energy of fusion ranges from 1 . . . 2 eV/Å; on average half of the free edges' dangling bonds get passivated. The picture is not this simple, however, as MD simulation creates different polygons that cause strain. The randomness of the polygon formation is manifested by the energy variations in Fig.3. Compared to ideal GBs with pentagons and heptagons only, our ensemble reveals the full complexity that rough GBs have.
GB energies are the highest around θ ∼ 15 • , which can be understood from structural analysis, shown in Fig.4a as polygon distribution within the GB zone. What we count into GB zone are all the polygons that were not part of the ribbons prior to merging. Polygons larger than nonagons, which appear more like vacancies instead of polygons, are omitted here. It is around θ ∼ 15 • where the abundance of hexagons is at minimum and GBs are invaded by other polygons. The abundance of pentagons and heptagons is as high as the abundance of hexagons, but also squares, octagons and nonagons are found. Close to θ ∼ 0 • and θ ∼ 30 • hexagons prevail. From this we may conclude that around θ ∼ 15 • the edge geometries have the largest mismatch, resulting in various polygons, consequent strains, and high energy.
To understand how much different polygons cost energy, we used the quantity measuring how much atoms, on average, cost more in ngon relative to atoms in pristine graphene; this quantity is for illustration only and considers polygons as separate items-the sum of ε n 's for given GB is not the total energy. Fig.4b shows ε n 's for the polygons within GB zones as a function of the area of the polygon (determined by triangulation). For simple geometrical reasons for small polygons the area distribution is narrow; large polygons have more freedom to change their shape. For small polygons the distributions in ε n , on the contrary, are wider; this is partly due to smaller n in Eq.(5). Clearly, the cost ε n of any polygon depends on its environment, just as it also depends for hexagons, for which ε 6 = 0 . . . 1 eV. But note that Fig.4b already contains the effect of the polygon environment and all potential cross-correlations (such as pentagons often neighboring heptagons). It would be interesting to investigate polygon statistics also from transmission electron microscopy, now that aberration-corrected measurements can achieve atom accuracy 34 .

VI. TRENDS IN BUCKLING (INFLECTION ANGLES)
In the generation process GBs are free to deform. The optimized GBs will typically end up having inflection angles, as the inset in Fig.5 illustrates. The data points in Fig.5 show the inflection angles for the GB samples as a function of chirality. Flat GBs occur with θ = 0 • and θ = 30 • , that is, with pristine graphene alone.
The notable trend in this scattered plot is the systematically small inflection angles, meaning flat GBs, around θ ∼ 30 • ; near θ ∼ 0 • inflection angles are more scattered meaning geometries that vary from flat to sharply kinked GBs. This trend can be understood by edge profiles: the stronger dangling bonds at zigzag edge, when brought into contact with another edge, can induce larger distortions than inert armchair edge (dangling bonds at armchair are partly quenched by triple-bond formation). The steps in edge profiles near θ ∼ 0 • , moreover, are 2.1Å, whereas the steps in edge profiles near θ ∼ 30 • are only 1.2Å. This means that, to make bonds, the edge atoms near θ ∼ 0 • need more pulling than edge atoms near θ ∼ 30 • , causing the buckling.
It is clear that the numbers in Fig.5 have no direct relevance as such, because our geometrically optimized GB samples are in vacuum, and measure only ∼ 20Å across the GB. We argue, however, that the inflection angle measures how much given GB would buckle in an experiment-large inflection angle meaning tendency to stick out (sticking out would be bound by geometric constraints and by surface adhesion). For example, a GB with θ ∼ 30 • on a support will always remain smooth, whereas a GB with θ ∼ 0 • will either remain smooth or buckle, to be seen as a ridge or a mountain range in scanning probe experiments. 16 Indeed, this buckling tendency was seen in an experiment byČervenka at al., measuring corrugation heights up to 3Å within GB regions 2 . Small corrugation heights and inflection angles are natural in samples on flat substrates, but also large inflection angles are realistic, for example in soot particles. 6

VII. TRENDS IN CHEMICAL PROPERTIES
Now we pose the general question: how reactive are GBs? We approach this question by examining hydrogen adsorption energies and dangling bonds (DB). To this end, we first develop a connection between hydrogen adsorption energy and the electronic structure given by DFTB. Fig.6b shows the hydrogen adsorption energies for carbon atoms in 6 representative GB samples, from total 350 hydrogen adsorption calculations. It appears that adsorption to given carbon atom i is strong if atom's cohesion − i decreases. Carbon atoms part of regular hexagons have adsorption around −2 eV ( in Fig.6a) but adsorption increases if atom is surrounded by other polygons and bonding angles deviate from 120 • (, , The error bars remind of technical ambiguities in angle determination. The main trend here is that zigzag edges cause both large and small inflection, whereas armchair edges cause only small inflections. Inset: inflection angle illustrated. and in Fig.6a). The common denominator in these examples is that the change in hydrogen adsorption energy is caused by the strain in bond angles and in bond lengths.
The strongest adsorptions around −5 . . . − 6.5 eV, in turn, are caused by dangling bonds that have i > −7 eV ( in Fig.6); such strong adsorption never occurs with three-coordinated atoms. (We gave the argument about DB energetics already in Sec.II.) We note that DFTB hydrogen adsorption energy to zigzag (5.8 eV), for example, agrees reasonably with DFT energy (5.4 eV) 32 . Therefore, we can characterize the reactivity directly by the DFTB electronic structure using the quantities i , without any adsorption calculations. Surely, dangling bonds can be identified from geometry by defining criteria for coordination numbers, but this approach is prone to errors, especially in disordered regions that GBs are.
Using i > −7 eV as a criterion for a dangling bond, we then analyzed the reactivity for the whole GB ensemble. Fig.6c shows the average number of dangling bonds per unit length in a GB with a given chiral angle. The highest density, one DB per ∼ 14Å, occurs with θ ∼ 20 • . The lowest density, one DB per ∼ 30Å, occurs with θ ∼ 15 • ; pristine graphene with no DBs becomes probable with θ ∼ 0 • and θ ∼ 30 • .
By analyzing the histogram in Fig.6c another way, 30 % of the samples have no DBs, 54 % of the samples have one DB every 7.5 . . . 20Å, and 16 % of the samples have one DB every 20 . . . 50Å. These numbers agree with a recent scanning tunneling microscopy (STM) experiment. Namely, dangling bonds are highly localized states just below the Fermi-level, and hence seen as a bump in constant-current STM with low bias. The STM a b c FIG. 6: (color online) Characterizing reactivity in grain boundaries. a) Grain boundary sample with (9, 5) chirality, including various polygons. b) Characterizing chemical reactivity from the DFTB electronic structure directly: we calculated hydrogen adsorption energy for selected 350 atoms as adsorption sites and plot for each atom adsorption energy and i. Given the correlation between the adsorption energy and i, we access reactivity directly from the electronic structure, without additional calculations. The adsorption site numbers refer to panel a); atom , having three ∼ 120 • angles has a small adsorption energy, atom has a dangling bond, and other atoms have strained environment. c) Using the correlation established in panels a) and b), we calculate the averaged density of dangling bonds per unit length along grain boundary for all samples, as a function of θ. 20Å periodicity for 47 % of the samples. Even if these periodicities should be caused by adsorbed impurities and not from bonds that dangle, it is still likely to be result from reactive sites within GBs-and answers the original question we posed in this section. The agreement with experiment gives confidence in the objectivity of the GB construction process. Understanding the trends in defects and reactivity will hopefully help in the design of functionalized graphene compounds 13 .

VIII. TRENDS IN VIBRATIONAL PROPERTIES
The structure of a GB, given its constituent polygons, affects directly on its vibrational spectrum, and gives experimentally complementary information. Fig.7 shows the projected vibrational density of states (PVDOS) for the GBs as a function of chiral angle and wave number. PVDOS was calculated by solving vibrations for the whole GB, by projecting the eigenmodes to atoms within GB area, and by renormalizing the modesfor all the 48 GB samples. Hence Fig.7 shows a lot of data in a compact form.
Spectra show three main bands, two at low energies ∼ 750 cm −1 and ∼ 500 cm −1 , and one-the socalled G-band in Raman spectroscopy-at high energy ∼ 1600 cm −1 . The band at ∼ 500 cm −1 is steady across all θ's and can not be used for structural identification. The G-band, in turn, loses intensity and comes down some 100 cm −1 in energy around θ ∼ 15 • . This can be explained by structure. As discussed in Sec.V, around θ ∼ 15 • hexagons are at relative minimum, implying a non-uniform and less rigid structure, and causing floppiness in high-energy modes. The high-energy modes have bond stretching between nearest-neighbors and are hence sensitive to local structural changes. Low-energy modes, again, are more collective and hence insensitive to local changes in structure, as long as they remain somehow "graphitic". If Fig.7 and Fig.4 are compared carefully, one finds that the local intensity maxima of the G-band occur precisely for θ with local maxima in the abundance of hexagons. The intensity increase around ∼ 750 cm −1 and θ ∼ 15 • is mainly due to renormalization of PVDOS.
The observations above are qualitatively supported by earlier experiments. Using Raman spectroscopy, a Gband shift from 1600 cm −1 to 1510 cm −1 was reported by Ferrari and Robertson in Ref. 35; the shift was identified to be a consequence of structural change from nanocrystalline graphite to amorphous phase. While still being mainly planar (apart from inflection), GB zone around θ ∼ 15 • , such as Fig.6a with θ = 20.6 • , indeed appears amorphous.

IX. CONCLUDING REMARKS
Some defects, like singular Stone-Wales defect or reczag edge of graphene, have their own name because they are well defined. Grain boundaries, on the contrary, are like snowflakes-there is no flake like another, but it's enough to know that they are roughly hexagonal, flat, small and cold. For the same reason it's enough to know how grain boundaries usually look and feel. Knowing trends is valuable.
We investigated GBs from different viewpoints, discovering trends with complementary information, measurable also experimentally. We investigated (i) geometry (polygons and inflection angles) measurable with transmission electron microscope or atomic probe microscope; (ii) energy, that is manifested in geometry and thereby measurable; (iii) reactivity, measurable with adsorption experiments; (iv) vibrational properties, measurable with Raman spectroscopy. We are confident that the trends are genuine, part because of agreements with earlier experiments, part because the trends all make intuitive sense: energy, inflection angles, reactivity and vibration trends make sense given the structure, and the structure trends make sense given the graphene edge mismatch. We hope the trends help to explain experiments-and also help simulations and experiments to design graphene structures for given functions.