Atomic structure and energetics of large vacancies in graphene

We present a computational study on the topology, energetics and structural deformations for a large number of experimentally observed defect configurations in graphene. We find that both the number of lost hexagonal carbon rings and introduced non-hexagonal rings increase linearly as a function of the vacancy order (number of missing atoms). The formation energies of the defects increase by about 2.2 eV per missing atom after an initial offset, establishing these defects as the lowest energy vacancy configurations studied in graphene to date. In addition, we find that even small point defects, which have been until now assumed flat, cause graphene to bend out of plane when not restricted into prohibitively confined geometries. This effect reaches to relative long distances even for some of the smallest defects, significantly reducing the stress otherwise imposed on the surrounding lattice.


I. INTRODUCTION
Two-dimensional materials, such as graphene and single layers of hexagonal boron nitride and transition metal dichalcogenides, have provided an unprecedent opportunity to directly image the atomic structure of various lattice imperfections.For example, high resolution transmission electron microscopy (HR-TEM) studies have revealed the atomic structure of various graphene defects [1][2][3] and grain boundary structures [4,5], while scanning tunneling microscopy (STM) has provided a complementary view [6][7][8].The energetic electrons, which are used for the imaging within TEM, have also been utilized to intentionally create defects [9][10][11][12][13] and to drive their dynamics [9,[14][15][16].However, with both TEM and STM, one remains typically agnostic about any small height variations that may occur around the defects.This is because changes in the sub-nm range are difficult to discern based on the focus in a TEM experiment, whereas STM measurements for supported graphene are more sensitive to height variations of the substrate than those of graphene, and interactions between the STM tip and a freestanding graphene can significantly alter its shape during the measurement.[17] From the theoretical point-of-view, defects in graphene have offered an intriguing playing field allowing for a multitude of different atomic configurations via incorporation of various non-hexagonal carbon rings into the lattice.An overview on research of structural defects in graphene until 2011 is provided in Ref. [18].Since then, further studies have been conducted on the atomic structure of point defects [11][12][13], dislocations [3,16] and grain boundaries [19].However, few attempts [20,21] have * jani.kotakoski@iki.fibeen made to understand the properties of realistic larger vacancy complexes in graphene.Further, only rarely in studies of defects in graphene (with a few notable exceptions [16,22,23]), have out-of-plane variations been adequately addressed.Instead, the structures have been typically assumed flat, and the atomic projections in the TEM images have been interpreted as corresponding to a flat structure under significant negative strain [13,19].However, by bending out from the flat configuration, a thin membrane can locally reduce the stress around defects.
In this work, we study a set of defect configurations in graphene, produced and observed during a TEM experiment under a 100 kV electron beam.The formation of these defects has been previously described in Refs.[9,10].TEM images of the defects can be found in the supplementary video S4 of Ref. [10].We find that, on average, 1.6 hexagons will be transformed into 1.0 nonhexagonal carbon rings per one missing atom.All of the studied structures are observed to deform graphene in the out-of-plane direction, which significantly reduces the local stresses around the defects.The formation energies increase linearly as a function of the vacancy-order (by ∼ 2.19 eV per missing atom after an initial offset), establishing this set of defects as having the lowest formation energies for any vacancy complexes hitherto studied in graphene.

II. METHODS AND RESULTS
We start the analysis by looking at the topology of the defects (see Fig. 1a for four example structures).In Fig. 1b, we plot the number of introduced non-hexagonal polygons and the number of removed hexagons (with respect to the pristine lattice) as a function of the va-arXiv:1404.5396v2[cond-mat.mtrl-sci]22 May 2014 cancy order (number of missing atoms).Out of all nonhexagonal carbon rings, 0.8% were tetragons, 55.3% pentagons, 38.2% heptagons and 5.8% octagons.The number of removed hexagons increases (on average) linearly with the number of removed atoms with a rate of ca.1.6 per atom.This is comparable to value of 2 for the V 2 (5-8-5), but somewhat lower than those for the V 1 (5)(6)(7)(8)(9) or the energetically more stable divacancies [18], namely, the V 2 (555-777) and the V 2 (5555-6-7777), which have values of 3, 3.5 and 4.5 per missing atom, respectively.The number of non-hexagonal rings, on the other hand, increases at a rate of about 1.0 polygons/missing atom.Comparing to the average value of ∼ 0.6, for all of the above-mentioned divacancies, the difference between removed hexagons and introduced other polygons is exactly 0.5 per missing atom (1.0 for the single vacancy).We stress that all of the above analysis is based on the average properties of all of the defect structures, which means that for any specific defect, the values can deviate from the ones listed.All of the analyzed atomic structures are available through Ref. [24].
The largest defect structure considered in this study has 44 missing atoms and involves nearly 80 lost hexagons with up to 40 other polygons introduced.Thus, in order to contain each of the structures within a pristine lattice of the same size, we had to create a relatively large simulation system into which the defects were then (separately) introduced.As a consequence, also the small wellstudied defects were in this study relaxed in an atypically large system.We carried out structural optimization for all of the defect structures using the conjugate gradient energy minimization scheme employing two different analytical interaction models to describe the energetics and inter-atomic forces in the system, namely the PII parametrization by Brenner from Ref. [25] and the improved version of the potential (AIREBO) by Stuart et al. [26].Both potentials reproduced all of the trends reported within this work.However, since the bond lengths and formation energies predicted by AIREBO are closer to density functional theory (DFT) values, as will be to some extent described below, all of the results shown were obtained with AIREBO.These simulations were carried out with the LAMMPS code [27,28].The DFT results were calculated with VASP [29,30] using projector augmented wave potentials [31].Plane wave cut-off was set to 300 eV, and the exchange and correlation were described with the generalized gradient approximation parametrized by Perdew, Burke and Ernzerhof [32].Only one k-point was used (Γ).These parameters were dictated by the relative large system size (up to 880 atoms for the DFT simulations).
During initial relaxation of the structures with both analytical potentials, we noticed that the defects caused the lattice to bend out from the flat configuration.In order to check that this observation was not caused by a simulation artefact, we selected the V 2 (5-8-5) divacancy to study this effect.This defect has been extensively studied in the past (see, f.ex., Ref. [33]).Curiously, it is known to lead to a local curvature change in carbon nanotubes [33,34], but for graphene, as far as we know, only flat configurations have hitherto been reported.As can be seen in Fig. 2a, the formation energy of this defect decreases with increasing system size up to about 1,000 atoms due to repulsion between the defect and its mirror images over the periodic boundaries at shorter interdefect distances.At system sizes exceeding 100 atoms, the flat configuration becomes less favorable than the buckled one, and remains so up to the largest studied systems.From this data, we can conlude that defectto-defect interactions reach up to at least 15 nm (corresponding to the system size for 20,000 atoms in this study) even in the case of small point defects in graphene.Interestingly, we find two possible buckled structures for this defect with an anti-parallel and a parallel symmetry with respect to the graphene plane [see Fig. 2b].Both of these configurations are favored over the flat structure.
We point out that no restriction was required to keep the structure flat during the relaxation at any system size.According to the analytical potentials, this configuration presents thus a metastable defect state even when it is less favored than the nonflat structure.However, one can reasonably assume that in any practically relevant situation, the graphene membrane will tend towards the energetically favored configuration due to external influences (f.ex., thermal vibrations).We assume that these reasons have until now hidden the corrugated nature of graphene with point defects in theoretical studies.In contrast to rather large system sizes required to reveal the out-of-plane bending of the divacancy, earlier quantum monte carlo and DFT results have shown that the Stone-Wales defect bends the graphene sheet already at much smaller system sizes (down to just a few tens of atoms) [22].
Also in Fig. 2a (second y-axis), we plot the maximum distance for any two atoms in the out-of-plane direction for the optimized locally buckled configurations (∆z).The local height of the membrane around the defect increases continuously indicating that the further decrease in strain, mediated by the out-of-plane deformations, does not significantly alter the energetics of the defect and that the defects become sufficiently spatially isolated at the largest system sizes.
We further looked into how the out-of-plane relaxation correlates with the release of the local negative strain.The length of the shortest bond in each relaxed structure was found to increase monotonously with the defect height for the buckled structures (from 2.4% shrinkage with respect to ideal graphene up to 1.7% shrinkage over the studied range), whereas for the flat structure the shortest bond remained at 2.4% shrinkage regardless of the system size.The number of bonds with more than 0.5% shrinkage saturated at 26 for the non-flat structures and at 36 for the flat ones at the largest system sizes.These numbers clearly show that the out-of-plane relaxation is an important way for the graphene lattice to relieve negative strain caused by the formation of de- fects.We also checked how the deformation reacts to strain in the lattice.It turns out that 0.5% is enough to make the divacancy structure flat for a 20,000-atom system.However, strains up to several % are required to flatten any of the larger defects.A more in-depth study of these deformations is out of the scope of the present work.
Based on the above analysis, we established 20,000 atoms as a reasonable system size for the AIREBO simulations (we tested system sizes up to 180,000 to check that there are no significant changes at formation energies even for the largest defects).However, structures this large are clearly beyond what can be typically modeled with any first principles method, ∼ 1, 000 atoms being more typical.A comparison of our AIREBO results for systems of 20,000 atoms and 880 atoms showed that the formation energies of small defects (∼ 6 missing atoms) were already accurate within 0.5 eV for the small system with a few exceptions being accurate within 0.1 eV, while larger defects exhibited deviations in the range of 2 − 3 eV.Additionally, as was shown above, the buckled structures appear energetically favored over the flat ones already for the 880-atom system.(Although, the preferred mode of bending for an isolated defect may remain hidden.) In order to qualify our AIREBO simulations, we repeated some of the V 2 (5-8-5) simulations with DFT for system sizes of 72, 200 and 880 atoms.For the 880-atomic system, we obtained formation energy of 7.13 eV for the buckled structure (anti-parallel) and 7.21 eV for the flat one, confirming our AIREBO results with a similar energy difference (0.04 eV for AIREBO for the same system size).The value for the flat structure is in line with the results reported in the literature [18].(The AIREBO overestimates E f for this defect by about 1 eV.)The maximum atom-to-atom distance in the out-of-plane direction (∆z), as obtained from DFT simulations is ca.0.76 Å, indicating that AIREBO slightly overestimates the magnitude of the deformations.For both of the smaller systems (72 and 200 atoms), DFT simulations converged into the flat structure, showing that they remain too small to accompany the out-of-plane relaxation.
After establishing the reliability of the AIREBO results regarding energetics and structural deformations, we expanded the formation energy study for all of the defect structures obtained from TEM images.The results are shown in Fig. 3 (notice that there are several different defects with same number of missing atoms, N ).The formation energies scale linearly with the increasing vacancy order (E f /N ∝ 1/N ), increasing by ca.2.19 eV per missing atom, with an initial offset of about 5.35 eV.This energy penalty is associated with the local deformation of bonds around the locations of removed atoms.The linear dependency is similar to that found by Jeong et al. [20] for dislocation lines and local haeckelite structures [composed of merged V 2 (555-777) defects] up to 12 missing atoms, although the formation energies of the vacancy structures considered here are substantially lower (by about 0.5 eV per missing atom).Introduction of the first defect into the pristine lattice is associated with the highest energy penalty (the initial offset in E f ), because the growing defect can easier accumulate the non-ideal bonds due to effects like overlapping strain fields with opposing signs [34], which leads to lower energy penalty for the removal of the additional atoms.This behavior is the physical reason why graphene can be turned into an amorphous 2D carbon glass via introduction of a growing number of defects [9,35].Finally, we turn to look at the out-of-plane deformation caused by the larger vacancies.As an example, four height maps for defects with different N are presented in Fig. 4a.Also the ∆z and an estimation of the in-plane size of the out-of-plane deformation for all defects are plotted in Fig. 4b,c.This estimation was made by plotting height maps, similar to those in Fig. 4a but for a larger area, and by drawing a circle around each defect so that all areas with a height deviation more than an arbitrarily selected ∆z 0 were contained inside the circles.The diameter of these circles was then used as the measure for the size.This approach resulted in estimation of deformation length in the in-plane direction, which is linearly dependent on ∆z, showing that each of the buckled structures has approximately the same local curvature.∆z itself appears to increase for defects with N < 15 proportional to √ N , which agrees with the assumption that the area of the defect grows linearly with N .At N ≈ 15, however, ∆z saturates to approximately a constant value of ca.4.3 Å.This is most likely due to the fact that at this point the growth of the defect leads to buckling at the defect itself (increasing the frequency of the out-of-plane deformations around the defect) without further contribution to the buckling amplitude of the membrane.This can be seen as the increased number of ripples for the largest defect in Fig. 4a as compared to the smaller defects.These trends, and even absolute values, remained almost completely unchanged also for the largest studied system sizes (up to 180,000 atoms).The actual interaction distances reach up to tens of nm, as was already discussed above in the case of V 2 (5-8-5).

III. CONCLUSIONS
As a conclusion, we have presented a computational study on the topology, energetics and atomic structure of experimentally observed defect structures in graphene.Unlike what has been previously assumed, all of the defects, including the smallest ones, cause local buckling of graphene to lower the negative strain otherwise imposed on the lattice.The number of removed hexagonal carbon rings and introduced other polygons are shown to systematically increase as a function of the number of missing atoms.The ratios of the other polygons remain approximately constant.Formation energies of the defects increase linearly with the increasing defect size by ca.2.2 eV per removed atom after an initial onset of about 5.4 eV.Finally, the height of the buckling is shown to increase as a square root of the number of removed atoms until 15 missing atoms before saturating to up to 4 Å for the largest defects.These deformations result in an interaction length between defects in graphene in the regime of tens of nanometers.

Figure 1 .ideal 6 −
Figure 1.(Color online) (a) Four examples of the 44 vacancy structures used in this study with varying number of missing atoms (N ).x and y are the cartesian coordinates in the in-plane direction.(b) Number of missing hexagonal carbon rings (N ideal 6 − N6) and number of all other polygons ( i =6 Ni) as a function of N .Note that for many N , there are several different configurations.Lines are fits to the data.

Figure 2 .
Figure 2. (Color online) Effect of system size and out-ofplane corrugations on the formation energy of a V2(5-8-5) divacancy structure in graphene.(a) Formation energy E f as a function of the system size (Natoms, number of atoms) for the flat and the two non-flat configurations (with parallel and anti-parallel symmetries) and the maximum distance of two atoms within the relaxed non-flat structures in the outof-plane direction (∆z).(b) Out-of-plane corrugations for the two non-flat configurations.The lines are guides to the eye.x and y are the cartesian coordinates in the in-plane direction.

Figure 3 .
Figure 3. Formation energies per number of missing atoms E f /N for vacancy structures, as calculated with the AIREBO potential for a system of 20,000 atoms.The dashed curve is a fit to the data.

Figure 4 .
Figure 4. (Color online) Out-of-plane buckling.(a) Height maps for four example defects with different number of missing atoms (N = 4, 8, 16 and 28, from left to right).x and y are the cartesian coordinates in the in-plane direction.(b) Maximum difference in the out-of-plane direction (∆z) for each of the defects as a function of the number of missing atoms, N .The dashed and dotted lines are fits to the data for N < 15 and N > 15, respectively.(c) Relative increase of the in-plane size of the local out-of-plane deformation as compared to the largest one as a function of ∆z.