Energetics of atomic scale structure changes in graphene

The presence of defects in graphene has an essential influence on its physical and chemical properties. The formation, behaviour and healing of defects are determined by energetic characteristics of atomic scale structure changes. In this article, we review recent studies devoted to atomic scale reactions during thermally activated and irradiation-induced processes in graphene. The formation energies of vacancies, adatoms and topological defects are discussed. Defect formation, healing and migration are quantified in terms of activation energies (barriers) for thermally activated processes and by threshold energies for processes occurring under electron irradiation. The energetics of defects in the graphene interior and at the edge is analysed. The eﬀects of applied strain and a close proximity of the edge on the energetics of atomic scale reactions are overviewed. Particular attention is given to problems where further studies are required.


Introduction
Since the discovery of one-and two-dimensional carbon nanostructures, nanotubes 1 and graphene 2 respectively, a variety of interesting fundamental properties and applications of these nanostructures have been found. Contrary to bulk matter, oneand two-dimensional nanostructures do not correspond to the ground state of an infinite 3D system. Thus, studies of the kinetics of atomic scale processes are especially important for understanding atomistic formation mechanisms of these nanostructures. Such kinetics are determined by the energetics of bond realignment, emission and insertion of carbon atoms.
Transformation of carbon nanostructures under heating 3-6 or electron irradiation 3,4,[7][8][9][10][11][12][13][14] leads to a production of entirely new species (see also ref. 15  transformation processes under electron irradiation). Among the most striking examples of such processes are the formation of graphene nanoribbons from different organic molecules inside carbon nanotubes under heating 3,4 and electron irradiation, 3,4 the formation of flat nanotubes from bilayer graphene nanoribbons, 7 fullerenes from initially flat graphene flakes 8 and carbon chains from graphene 13,14 under electron irradiation. Others include the transformation of polyhedral graphitic nanoparticles into quasi-spherical onions under heating 5 and electron irradiation, 9,10 the formation of double-walled nanotubes from single-walled nanotubes filled with fullerenes under heating 6 and electron irradiation, 11 and the formation of a trilobate structure from three La@C 82 endofullerenes inside a carbon nanotube under electron irradiation. 12 Mesoscopic processes such as the shrinking of nanotube diameter under electron irradiation, 16 migration and merging of large holes in the outer wall of nanotubes at high temperature, 17 superelongation of carbon nanotubes at high temperature, 18,19 and the migration and annihilation of graphene grain boundaries under electron irradiation 20 have also been observed. These large transformations in nanostructures take place via bond realignments complemented by atom emission, and are determined by energetic characteristics of the corresponding atomic scale reactions.
There are currently several reviews devoted to the properties of carbon nanostructures with irradiation defects 97,98 and properties of graphene with defects. 99,100 A very recent review 101 covers the general study of graphene with TEM, while the various structures of interest in graphene under TEM are detailed in ref. 72. However, the energetics of rate determining bond reorganisations in graphene has not yet gained considerable attention within a single review. In the present review we provide an overview of the energetic characteristics of various thermally activated and irradiation-induced reactions in the graphene interior and at graphene edges. In addition to static characteristics, such as formation energies of defects relative to perfect structures, we consider activation energies and threshold energies under irradiation for various bond reorganisations and atom emission reactions. The relationship between atomic scale energetic characteristics and experimentally observable macroscopic processes is discussed.
As most energetic characteristics related to atomic scale structure changes in graphene cannot be determined from experiment, the majority of the data presented in the review is obtained by chemical calculations. However, the activation energies for thermally activated processes and the threshold energies for processes occurring under electron irradiation are

Andrey M. Popov
Andrey Popov is a senior scientific researcher at the Institute of Spectroscopy of Russian Academy of Science, Moscow. He received his MSci degree (1989) from the Moscow Institute of Physics and Technology, and his PhD in Physics (2007) from the Institute of Spectroscopy of Russian Academy of Science. Andrey is an author of 80 refereed publications including 4 reviews concerned with the energetics, structure, formation, transformation, thermodynamic and tribological properties of nanostructures and related nanotechnology problems, particularly applications of carbon nanotubes and graphene in nanoelectromechanical systems.

Elena Bichoutskaia
Elena Bichoutskaia is an Associate Professor in Theoretical and Computational Chemistry at the University of Nottingham, UK. Her research interests include the development of theoretical and computational approaches to prediction of materials properties; computational modelling of the behaviour, properties and manipulation of carbon nanomaterials; electrostatic interactions at the nanoscale; gas storage and interactions in porous solids. Elena is a recipient of EPSRC Career Acceleration Fellowship rather sensitive to the level of theory used in calculations. Section 2 details some common computational approaches to calculating the key parameters of thermally-and irradiationinduced structure transformations. Section 3 presents the structure of defects in graphene as well as their formation energies and barriers for atomic scale reactions in thermally activated processes. The threshold energies for reactions under irradiation are described in Section 4. Conclusions based on published results are discussed in Section 5, together with some problems which remain unsolved.

Methods for calculating interatomic interactions in carbon nanostructures
Depending on the complexity of the problem under consideration, different methods should be chosen for the description of atomic interactions in carbon nanostructures. The most accurate results are obtained by quantum chemical methods that do not require a priori knowledge of any specific characteristics of the considered system. The optimal balance between accuracy and efficiency is achieved by density functional theory (DFT), 102,103 in which the total energy of the system is considered as a functional of the electron density and the description of the system is reduced to a set of self-consistent one-particle Schrödinger equations. In the simplest local density approximation (LDA), 104 the exchange-correlation energy depends only on the local electron density, which is sufficient for a description of periodic solids but can lead to significant errors for molecular systems and surfaces. Introduction of the dependence on the gradient of electron density in the generalized gradient approximation (GGA) improves results considerably. Perdew-Wang 105,106 (PW91) and Perdew-Burke-Ernzerhof 107,108 (PBE) GGA functionals are successfully used for calculations of defects in carbon nanostructures. Admixing of a fraction of exact exchange in hybrid functionals (PBE0, 109 B3LYP 110,111 ) or inclusion of higher derivatives of the electron density in metaGGA functionals (such as M06-L 112,113 ) further improves DFT performance, while dispersion corrections (vdW-DF, 114 116 ) are important for accurate assessments of geometrical and energetic characteristics of defects in few-layer carbon nanostructures. The account of spin is critical for the analysis of defects with unpaired electrons, such as monovacancies, or graphene edges.
Though DFT offers an efficient framework for full-scale quantum modelling, the size of the considered systems is still limited to hundreds of atoms. Tight-binding models [117][118][119][120][121][122][123][124] represent a cheap alternative to first-principles methods, where quantum effects are captured by a direct, albeit simplified, description of electronic structure. Tight-binding parameters determining single-electron energies and a classical repulsive potential energy term are fitted to reference datasets or results of first-principles calculations. Tight-binding models can therefore be considered as a bridge between classical methods and full modelling of the quantum nature of chemical bonding. The incorporation of self-consistent calculations of the charge distribution 117,118 further improves the accuracy of tightbinding methods.
Classical force fields enable simulations with increased computational efficiency and system size, using analytical potentials to describe the interactions between atoms. These potentials are constructed and fitted to describe the particular behaviour of atoms in the system, observed in experiment or in accurate quantum mechanical calculations. Alongside simple analytical forms of pair potentials, such as the Morse potential, 125 more complex reactive potentials have been developed, capable of describing covalently bonded crystal structures such as graphene. In traditional forms of potentials all bonds are defined explicitly, which makes them unable to model chemical reactions due to the requirement of breaking and forming bonds. Reactive forms of potentials, however, deliberately avoid using explicit bonds in favour of bond orders, thus allowing for continuous bond formation and breaking. The development of reactive potentials is based on the decomposition of chemical bonding into individual contributions (angular, stretching, dihedral terms etc.) to the binding energy, while introducing cross-terms as a penalty for over-and/or under-coordination. Thirty years ago, Abell introduced a general expression for the binding energy as a sum of nearest neighbour pair interactions moderated by the local atomic environment. 126 In the 1990s, Tersoff 127 129 and non-bonded and dihedral-angle interactions were included (AIREBO). 130 Extensions of these potentials to metals 87,131-133 make simulations of complex catalytically activated processes possible. Van Duin, Goddard and co-workers developed a new formulation of the reactive potential, ReaxFF, based on accurate benchmarking density functional theory studies and extended its use to various materials, including hydrocarbon reactions, 134 transition metals, 135 silicon 136 and other materials such as polymers and ceramics.
In spite of this significant progress in the development of reactive empirical potentials for carbon systems, they are still not sufficiently accurate for studies of some problems related to graphene energetics. As mentioned below, the reactive bond-order potentials (initially developed for hydrocarbons) fail to reproduce the adequate ground state for monovacancies in graphene. Further improvement of existing potentials for the description of specific problems is possible, by fitting their parameters to certain properties that are important for the phenomena under consideration. For example, the parameters of the Brenner potential have recently been fitted to energies of zigzag and armchair graphene edges and elastic energies of the C 60 and C 70 fullerenes. 133 This modification of the potential improved the description of the balance between the fullerene elastic energy and graphene edge energies, which is important for the energetics of graphene formation or its transformation to other carbon nanostructures.

Analysis of thermally activated reactions of defects
A significant contribution to the understanding of the formation and transformation of defects in carbon nanostructures is provided by the analysis of their free energy surface. To reduce computational cost, these reactions are considered in simulation cells with periodic boundary conditions (PBC) or for finite-size flake models. The reliability of such calculations is determined not only by the accuracy of the description of the interatomic interactions, but also by the convergence of the results with respect to the size of the simulation cell (so that periodic images of the defects do not interact) or size of the finite model (to eliminate edge effects) and other simulation parameters. The saddle point in reactions of defects can be sought by restricting some degrees of freedom of the system or by modified geometry optimisation methods. 137,138 A more efficient nudged elastic band method, 139,140 which has been widely used in recent years, is based on the optimisation of the total energy of a string of system images connecting known initial and final states. Another useful method for searching for saddle points is action-derived molecular dynamics, which generates dynamical trajectories with fixed pre-assigned initial and final boundary conditions. 141 Such a method was used, for instance, to find a path with a set of saddle points for vacancies approaching and coalescing in graphene. 142 Using a simple Arrhenius formula, the reaction barriers obtained by these methods can be used to estimate the reaction rates at experimental conditions. The lifetime t s of a state can be calculated in this way as where n is the characteristic frequency, DE t is the barrier between states, k is the Boltzmann constant, and T is the temperature. Molecular dynamics (MD) 143,144 and Monte Carlo 143,145 methods allow computational modelling of the evolution of carbon nanostructures with time. In Monte Carlo methods, one of the possible elementary processes at each simulation step is chosen on the basis of their probabilities, obtained according to their energetics. A list of possible elementary processes and their parameters can be known in advance or determined ''on-the-fly'' by probing the free energy surface. Molecular dynamics, in which trajectories of atoms are obtained by the numerical integration of Newton's equations of motion for nuclei, with the forces evaluated according to the classical force fields or ab initio methods, does not require any information on possible elementary processes or analysis of the free energy surface. However, even for classical systems, this method is typically limited by the time scales of microseconds. Limitations of ab initio molecular dynamics methods are much more severe and the corresponding numerical experiments are three orders of magnitude shorter in time. Practically all processes of formation and transformation of defects are activated and would not fit into this timeframe under experimental conditions (e.g., at room temperature). A common approach to overcome this difficulty is to accelerate kinetics by increasing the simulation temperature. In this way, one can access valuable information on the process mechanism and gain insight into the kinetics of the process. However, extraction of kinetic parameters requires extensive statistics. Furthermore, care should be taken in the analysis of the results. Overheating the system in the numerical experiment distorts the free energy surface and affects the kinetics of various reactions in different ways, which may complicate the interpolation of the data to experimental conditions. Therefore, static analysis of the potential energy surface at zero temperature is still helpful for supplementing and verifying molecular dynamics results.

Calculations of threshold energies
The threshold energy of an irradiation-induced process is the minimum energy that must be transferred from an electron of the electron beam (e-beam) to an atom in graphene in order for the process to take place without immediate reversal. Although this can in principle be determined by experiment, due to the development of atomic resolution microscopy, currently there are few examples of experimentally determined thresholds. 146,267 The threshold energy of a process of interest is therefore usually obtained computationally, via a series of MD simulations. At the beginning of each simulation, a kinetic energy corresponding to the energy transferred from the beam electron is assigned to the impacted atom. The structural evolution of the entire system over time is then computed using MD in order to determine whether or not the process occurs. This is repeated over a range of initial kinetic energies, and occasionally over a range of impact angles, in order to determine the threshold energy to the desired accuracy. The binary search algorithm has previously been used 225 to choose the value of kinetic energy to apply in order to determine the threshold energy. However, if the computational cost of an individual MD simulation is high while the number of threshold energies to be determined is low (for example when calculating a single threshold energy at a very high level of theory), manual inspection of the MD trajectories and estimation of a suitable kinetic energy based on the results can be more efficient.
In addition to computationally inexpensive classical MD, density functional tight-binding (DFTB) has also been widely used to calculate threshold energies in graphene. 147,151,225,265,268 Its low computational cost relative to DFT enables calculations on systems with large numbers of atoms, or allows a large series of MD simulations to be run in order to obtain precise threshold energies. Apart from calculating threshold energies, DFTB can be generally very useful when modelling the effect of electron impacts on graphene structures, for example in the case of Stone-Wales rearrangements at a variety of impact angles, 147 in which the 27 000 DFTB-MD simulations performed would have been unfeasible at a higher level of theory. Despite the low computational cost, these calculations typically give comparable results to DFT for carbon atoms in graphene, although caution must be used when employing DFTB; emission threshold energies of zigzag edges in graphene for example have been shown to be overestimated by approximately 2.5 eV compared to the DFT case, ascribed to stronger local bonding in the DFTB relaxed structures (shown in Fig. 13). 148 Disparities for other systems have also been shown, such as in hexagonal boron nitride (h-BN), for which DFTB simulations predicted similar threshold energies for the emission of boron and nitrogen 149 while DFT simulations 150 later agreed with the experimentally shown asymmetry. This discrepancy was assigned to an inadequate description of charge transfer in the DFTB model.
DFT-MD provides an increase in accuracy and versatility compared to DFTB, at the cost of increased calculation time. Calculations using LDA and GGA exchange-correlation functionals are commonly employed in calculations of threshold energies. 148,150,266,269,271 Higher levels of theory, such as hybrid exchange-correlation functionals that include a degree of exact exchange (most notably B3LYP), can be used to obtain more accurate values, but the computational cost of these methods tends to limit the possible system size.
Infinite graphene can be emulated computationally in one of two ways: by using PBC and by using large graphene flakes or nanoribbons. Supercells used in calculations using PBC must be large enough that the atom hit by the electron beam is not influenced by its periodic image, while graphene flakes must be large enough to negate edge effects, shown to extend at least 10 Å from the edge. 67 Simulations using PBC are far more commonly used for calculating threshold energies than large flakes, due to the requirement of a much larger number of atoms to neutralize the large effects of the extended edges compared to a point defect. The atomic processes of structure change studied in graphene in this way typically take of the order of 100 fs to occur when the energies close to the threshold energy are applied, and a 1 fs time step is usually used. Zero-point velocities or thermal vibrations derived from the Maxwell-Boltzmann distribution are used for the initial velocities of other atoms in the system.

The CompuTEM algorithm
The CompuTEM algorithm 15,151 is a method for linking molecular dynamics and high-resolution transmission electron microscopy (HRTEM) image simulations. The link is established by incorporating structural information from the MD simulations into multislice image simulations, as well as using realistic estimations of the signal-to-noise ratio from the electron dose rate and detector limitations. The effect of the e-beam on a sample is described by incorporating structural changes caused by the electron irradiation from the MD simulations, at realistic rates determined from the electron dose rate and threshold energies. Fig. 1 illustrates this algorithm, in which a random irradiation-induced event occurring in a sequence of such events is described as follows: (1) the nanostructure is equilibrated at a temperature corresponding to experimental conditions in HRTEM, (2) each atom in the nanostructure is classified with respect to the number and strength of its chemical bonds, (3) the probability of an irradiation-induced event (such as atom removal and/or changes to the local atomic structure) is assigned to each atom in accordance with the atom type determined at step 2, so that the sum of the probabilities over all atoms and all considered types of events for each atom is equal to unity, (4) a single random irradiation-induced event is introduced, (5) MD simulation at a temperature corresponding to experimental conditions for a duration sufficient for bond reorganisation or atom removal, (6) MD simulation at the elevated temperature taking into account the structure relaxation after the irradiationinduced event. The introduction of the irradiation-induced event can be described by one of two approaches, an in depth discussion of which can be found in ref. 15. 3. Structure and energetics of defects in graphene 3.1 Vacancies 3.1.1 Effects of vacancies on mechanical, electronic and magnetic properties. The influence of vacancies on the physical properties of graphene has been widely studied. The Young modulus, [23][24][25][26][27][28] Poisson ratio, 25 and tensile strength 22,27,31 of graphene have been shown to decrease with an increase in vacancy concentration. The dependence of fracture strain on the concentration is non-monotonic; there is a decrease in strain at low concentrations and an increase at higher concentrations. 22 A set of calculations demonstrated that the presence of vacancies induces magnetism in graphene. 21,24,35,[45][46][47][48][49][50][51][52][53][54] According to other calculations the presence of divacancies leads to a decrease of graphene conductivity. 33,34 A decrease in the conductivity of graphene nanoribbons is also observed upon the incorporation of vacancies. 41 Vacancies, as all local defects, do not affect the universal quantization of lowtemperature thermal conductance, while they lead to a significant decrease in the thermal conductivity of graphene 26,27,58 and graphene nanoribbons 56,59,60 at room temperature. The thermal conductance of narrow graphene nanoribbons is found to be very sensitive to the position of the vacancy relative to the edge. 56 3.1.2 Structure and energetics of monovacancies. Although pioneer works consider the energetics of the monovacancy with D 3h symmetry and dangling bonds, 119,152 subsequent calculations showed that the local three-fold symmetry breaks down to C 2v symmetry due to Jahn-Teller distortion. 21,25,47,52,[153][154][155][156][157][158][159][160][161][162] Most of these calculations show that the reconstructed vacancy has a 5/9 structure 21,47,52,67,119,[152][153][154][155][156][157][158][159][160][161][162] where a new bond with length b v arises between a pair of the three atoms with dangling bonds, forming a pentagon with displacement d v of the third atom out of graphene plane (Fig. 2a). The structure of a reconstructed monovacancy with two new equivalent bonds and without any out of plane displacement of atoms has also been reported. 48 The formation energy of a defect is defined as where E d and E p are total energy of defective and perfect systems, respectively, m is the chemical potential estimated by calculations as the total energy per atom in graphene, and n is number of removed atoms for vacancies (positive value) or added adatoms (negative value). The values of formation energies E v and structural parameters b v and d v of the monovacancy in graphene calculated using different methods 21,25,35,47,48,52,54,67,119,[152][153][154][155][156][157][158][159][160][161][162][163][164] are listed in Table 1. To summarise the results presented in Table 1, recent DFT calculations for infinite graphene give vacancy formation energies in the range E v = 7.6-7.9 eV and the length of the newly formed bond b v = 1.8-2.0 Å. However, the out of plane distortions caused by defects in graphene are not fully understood. A reliable value for the out of plane displacement d v has not yet been obtained, as these deformations are highly sensitive to the system size and edge structure (for flake geometries). Note that the calculated value for the length of the newly formed bond coincides with the experimental one b v = 1.9 Å, measured by an aberrationcorrected TEM study. 159 TEM studies cannot directly measure out of plane displacement, however the observation of the projections of two bonds of the third under-coordinated atom of the vacancy being compressed to 1.37 and 1.21 Å (ref. 159) can be considered as an indirect argument that out of plane displacement d v takes place.
In graphite, the ground state of a single vacancy occurs when it is positioned over the hexagon centre of a neighbouring layer. This configuration has a lower formation energy than a vacancy positioned over an atom of a neighbour layer, by 0.03-0.04 eV for bulk graphite 47,160,165 and 0.03 eV for a graphite surface. 165 The newly formed bonds of this vacancy in graphite are 0.01 Å (ref. 47) and 0.03-0.04 Å (ref. 160 and 165) longer than when positioned over an atom of a neighbour layer. The formation energy of the vacancy on a graphite surface is found to be 0.7 eV lower than for bulk graphite. 165 Analogously to graphite, in AB-stacked bilayer graphene, a monovacancy positioned over a hexagon centre is preferred by up to 0.04 eV compared to a vacancy over a carbon atom. 35,54 The formation energy of a vacancy in twisted bilayer graphene is 0.01-0.07 eV larger than in AB-stacked bilayer graphene, depending on the relative positions of the layers. 35 The energetic and structural parameters of the vacancy in bulk graphite, on a graphite surface and in bilayer graphene, calculated using different methods, 35,47,54,152,160,165,166 are also listed in Table 1. The map of bond lengths around the vacancy in graphite at its ground and transition states has also been presented. 165 According to DFT calculations, the barrier for transition between the three possible equivalent states of the 5/9 vacancy is only DE t = 0.13 eV. This leads to a small lifetime of t s = 15 ps and t s = 32 ms at room and liquid helium temperature, respectively, 153 estimated by the Arrhenius formula (eqn (1)) with a characteristic frequency of n B 10 À13 s À1 . In the case of such a small lifetime, the superposition of three equivalent states should be observed in TEM observations, which have imaging timescales of the order of one second. However, both the reconstructed 5/9 vacancy 159,167,168 (observed for at least 90 s, ref. 159) and the symmetrical vacancy 159 have been found by HRTEM studies. The authors of this study proposed that the under-coordinated carbon atom is functionalised by a low mass contaminant such as hydrogen, precluding the oscillations between three equivalent states. However, the 80 keV electron beam used in this study should be able to easily remove hydrogen atoms. Note also that these rather old calculations giving the small barrier for transition between three equivalent states of the 5/9 vacancy were performed on a small graphene flake C 120 H 17 where the considered vacancy is close to the edge even in the middle of the flake. 153 However as discussed below, the proximity of the edge has drastic effects on the energetic characteristics of the vacancy. 56,67,68 The barrier DE t for transitions between equivalent states of the vacancy should correlate with values of the structural parameters b v and d v of the vacancy, which characterise the distortion of the symmetrical vacancy. Table 1 shows the wide spread of these parameters calculated using different methods. Note also the barrier DE t cannot exceed the relaxation energy (the energy difference between vacancies with reconstructed 5/9 structure and the vacancy with all bond lengths as in pristine graphene and three-fold symmetry). The relaxation energy of the vacancy with D 3h symmetry into a 5/9 vacancy can be considered as upper estimate of the barrier DE t . DFT-calculated values of this energy are also considerably scattered: 0.9 eV (ref. 156) and 0.29 eV (ref. 53) for graphene, 0.11-0.15 eV for bulk graphite, 165 and 0.14-0.16 eV for graphite surface. 165 Recent studies using density-functional methods show that the vacancy in graphene forms a dynamical JT centre in graphene (where the nuclei tunnel between the minima) owing to the small quantum mechanical barrier for nuclear tunnelling. 169 Note also that calculations for infinite graphene based on the tight-binding (TB) model give a considerably greater barrier for the transition between equivalent states of vacancy, DE t = 0.78 eV, 142 than the DFT calculations for the small graphene flake. Thus further experimental and theoretical studies are necessary to clarify the ground state structure of the vacancy in graphene.
3.1.3 Structure and energetics of divacancies. Three different reconstructed states of the divacancy in graphene are commonly observed on graphene under electron irradiation in HRTEM: 65,147,167,170 the V 2 (5-8-5), V 2 (555-777) and V 2 (5555-6-7777) states, see Fig. 3. Divacancies with the 5-8-5 structure were also found by scanning tunnelling microscopy in Ar + irradiated graphene. 171 60 kV HRTEM studies of monolayer graphene following bombardment with energetic gold particles allowed observation of the unreconstructed divacancy. 172 The transitions between the three reconstructed divacancy states under electron irradiation were observed in HRTEM. 65,147,167 The formation energies of the V 2 (5-8-5) and V 2 (555-777) divacancies Table 1 Energetic and structural parameters of the reconstructed monovacancy in graphene (see Fig. 2), graphite surface and bulk graphite: formation energy E v (in eV), barrier D E m for vacancy migration (in eV), bond length b v of new formed bond and displacement d v of the atom colored in dark blue out of graphene plane (both in Å). The system considered (periodic boundary conditions (PBC) or flake), number of atoms in the computational cell per layer without the vacancy and calculation method (SP and NSP stand for ''spin-polarised'' and ''non-spin-polarised'' calculations, respectively) are indicated

Ref.
System  [119][120][121] and DFT methods, are listed in Table 2. To summarise the results presented in Table 2, recent DFT calculations for infinite graphene give formation energies of the V 2 (5-8-5) divacancy in the range E 2v = 7.5-8.0 eV and energy gains for the transition from the V 2 (5-8-5) to the V 2 (555-777) state in the range DE 2 = 0.8-1.0 eV. The energy gain due to the transition from the V 2 (5-8-5) to V 2 (5555-6-7777) state was calculated to be 0.33 eV. 162 Maps of bond length changes around the divacancy have been presented. 25,161 The calculations of divacancy energetics in graphene and molecular dynamics simulations of its dynamical behaviour at high temperature (3000-4000 K) have also been performed 174 using the AIREBO potential. These calculations show that this potential gives the correct ground state of the divacancy but the wrong ground state of the monovacancy. Although 18 other possible states of the divacancy have been found during these simulations, these can only be considered as a qualitative description of possible divacancy metastable states. The barriers and transition states for transformation between the V 2 (5-8-5) and V 2 (555-777) states and between the V 2 (555-777) and V 2 (5555-6-7777) states have been found by DFT calculations 162 and are shown in Fig. 3. Close values of the barrier for the V 2 (5-8-5) to V 2 (555-777) transformation were found by other DFT studies, 5.27 eV, 170 5.17 eV (ref. 173) and 5.1 eV. 160 All three states of divacancy are therefore stable at room temperature.
3.1.4 Structure and energetics of multi-vacancies. A set of studies is devoted to the structure and energetics of trivacancies and multi-vacancies in graphene, obtained by structure reconstruction after the removal of neighbour atoms without bond rearrangements. 49,50,156 The trivacancy with this structure was obtained by bombardment with energetic gold particles. 172 For the vacancies with structures in which an odd number of carbon atoms are missing, at least one dangling bond persists, Fig. 3 Schematic representation of the structure and energetics for three states of a divacancy in graphene: V 2 (5-8-5) state, V 2 (555-777) state and V 2 (5555-6-7777) state, and transition states for transition between V 2 (5-8-5) and V 2 (555-777) states and between V 2 (555-777) and V 2 (5555-6-7777) states. Reprinted (adapted) with permission from ref. 162. Copyright (2013) American Chemical Society. Table 2 Energetic and structural parameters of reconstructed divacancies in graphene (see Fig. 3): formation energy of the V 2 (5-8-5) divacancy E 2v (in eV), energy difference between the V 2 (5-8-5) and V 2 (555-777) states of divacancy DE 2 (in eV) and bond length d v of newly formed bond of pentagon of V 2 (5-8-5) divacancy (in Å). The system considered (periodic boundary conditions (PBC) or flake), number of atoms in the computational cell per layer without the vacancy and calculation method (SP and NSP stand for ''spin-polarised'' and ''non-spin-polarised'' calculations, respectively) are indicated whereas for structures missing an even number of atoms (2, 4 and 6 missing atoms) more stable structures lacking dangling bonds are possible. 49,156 Other interesting types of multi-vacancies form in graphene as a result of vacancy coalescence and then subsequent bond rearrangement under electron irradiation in HRTEM. 65,66,167 In the case of multi-vacancies with an even number of missing atoms, the formation of lines of two or three V 2 (5-8-5) vacancies is observed in the graphene membrane under electron irradiation in HRTEM. 65,66,167 In such pairs or lines, adjacent divacancies have a common pentagon when aligned in the zigzag direction, while a tetragon is formed when aligned in the armchair direction, due to the overlap of two pentagons. The DFT-calculated energy gain is 1.32 and 2.01 eV per divacancy pair, as compared to isolated divacancies, for alignment in the zigzag and armchair directions, respectively. 167 The multi-vacancy created by combining two and three V 2 (5555-6-7777) divacancies has been considered. 167 The DFT-calculation-based comparison of formation energies for several types of multi-vacancies in graphene with even missing atoms has also been performed. 50,175,176 For the case of multi-vacancies with odd numbers of missing atoms, the formation of a bridging atom stabilizing the structure is observed in HRTEM. 66 The DFT-calculated energy gain is 1.55 eV for a trivacancy with a bridging atom, as compared to a trivacancy with a dangling bond. 66 3.1.5 Vacancy migration. Recent DFT calculations of formation energies of mono-and divacancies in infinite graphene 50,156,157,[160][161][162][163] show that the V 2 (5-8-5) divacancy has a 7.2-8.7 eV lower energy than a pair of separated vacancies. This means that the coalescence of two monovacancies into one divacancy is an exothermic process, which can take place via the thermally activated diffusion of vacancies or under electron and ion irradiation. In a similar manner, trivacancies have a lower energy than a separated monovacancy and divacancy, and so on. 49,156 The coalescence of large numbers of vacancies can lead to the formation of amorphous graphene structure. 146 The directed motion of a vacancy to the edge, due to a decrease in total energy at each step towards the edge, has been proposed based on DFT calculations 67,68 and molecular dynamics simulations. 67 The rates of the above processes are determined by the barriers for vacancy motion. The structures of the initial, transition and final states corresponding to one step of monovacancy motion is shown in Fig. 2. The transition state is a so-called spiro structure, where the moving atom is of equal distance to four nearest neighbours. 161,165 Note that due to the high symmetry of the spiro structure it should correspond to the minimum, maximum or saddle point of the potential energy of the system. The described step of vacancy migration was directly observed by a TEM study. 159 It is interesting to note that during the time between two subsequent image exposures, this TEM study observed the migration of a vacancy by two steps significantly more frequently than the migration by one step. Calculated values of the barrier DE m of migration of a vacancy in graphene, 23,24,152,153,155,157,160,164,173 bulk graphite 152,160,165 and on the graphite surface 165 are listed in Table 1. Some DFT calculations give the value of the barrier DE m of migration of a vacancy in graphite as being 0.15-0.26 eV lower than in graphene, 160 whereas DFT calculations with an empirical van der Waals correction give the opposite result; the barrier DE m is greater for graphite by 0.4 eV compared with the graphite surface. 165 Only two experimental estimates of the barrier DE m are currently available for bulk graphite and its surface. The first, DE m = 1.8 AE 0.3 eV, was obtained by Raman measurements of disorder relaxation in He + irradiated graphite. 177 Recent STM studies of vacancy aggregation on vacancydecorated graphite surfaces at different temperatures give the value DE m = 0.9-1 eV. 178 The calculated values of DE m listed in Table 1 range between these two experimental estimates and are only in tentative agreement with few experiments. Further studies are necessary to obtain accurate values of the barriers for migration of vacancy in graphene, bulk graphite and graphite surface.
Contrary to monovacancy migration, the barrier for divacancy migration is considerably large. The estimation of this barrier for a small graphene flake (C 120 H 17 ) gives the value of about 7 eV. 153 An estimation based on migration in two steps, via the dissociation and merging of two vacancies, gives the value 7.49 eV. 157 HRTEM studies reveal the migration of the V 2 (5-8-5) divacancy under the action of electron irradiation through the formation of an intermediate metastable dislocation dipole (the V 2 (55-77) state). 167,170 Divacancy migration under electron irradiation in HRTEM has also been observed as a result of transitions between the V 2 (5-8-5), V 2 (555-777) and V 2 (5555-6-7777) divacancy states. 147,170 It has been found using DFT calculations in combination with careful analysis of HRTEM images that one step of the migration of the V 2 (5-8-5) divacancy occurs through two intermediate V 2  states which have a higher total energy by 3.44 eV (other DFT calculations 3.72 eV (ref. 167)). The barrier for transformation from the V 2 (5-8-5) to V 2 (55-77) state is 5.27 eV, whereas the migration of the V 2 (555-777) divacancy occurs through a transformation to the V 2 (5-8-5) state with a barrier of 6.20 eV. 170 Thus, thermally activated migration of divacancies at room temperature can be excluded.
3.1.6 Effects of edges and deformation on the energetics of vacancies. The structure and energetics of the monovacancy near the edge of the graphene layer have been studied on examples of narrow graphene nanoribbons with the vacancy in the middle, 52,56 graphene flakes with dangling bonds at the edge 67 and graphene flakes with the edge terminated by hydrogen atoms 68,69 with different positions of the vacancy. Note that the case of an edge with dangling bonds is realised under the action of electron irradiation in TEM and is important when considering irradiation-induced processes such as the graphene flake to fullerene transformation. 8 It is evident that the structural parameters and formation energy of vacancies should change near the edge. For armchair graphene nanoribbons with widths of 7 and 5 hexagons, a stronger new bond is formed for the reconstructed 5/9 vacancy with the length b v of 1.9 and 1.8 Å, respectively, slightly less than the value b v = 1.96 Å obtained for bulk graphene using the same parameters of spin-polarised DFT calculations. 52 DFT-calculated formation energies of reconstructed 5/9 vacancies in the middle of graphene flakes (ranging in size from C 52 to C 116 ) with dangling bonds at the edge 67 are listed in Table 3. Table 3 also presents the difference DE in the total energy of the flake with the vacancy in the middle and at the edge. A considerable difference of DE = 6.6-8.9 eV in the total energies has been found for all flakes. The total energies corresponding to the subsequent steps of vacancy motion from the middle to the edge have been calculated for the C 116 67 and C 41 H 16 68 graphene flakes.
Moreover, upon decreasing the distance between the vacancy and the graphene layer edge, a change of energetic and structural parameters of the reconstructed 5/9 vacancy is predicted, as well as a drastic change of the whole vacancy structure. 56,67,68 DFT calculations show that the vacancy in the middle of an armchair graphene nanoribbon with the width of 5 hexagons has a spiro ground state with an sp 3 atom with four equivalent bonds 1.64 Å in length 56 in the centre of the vacancy (the same bond lengths as in the spiro state monovacancy in the middle of the C 116 graphene flake 67 ). The coexistence of the 5/9 and spiro states of the monovacancy near the graphene edge has been obtained by DFT calculations for graphene flakes with both dangling bonds 67 and hydrogen atoms 68 at the edge. It was also found that the spiro state of the vacancy is more stable than the 5/9 state in the middle of a graphene flake with hydrogen atoms at the edge when the distance between the vacancy and the edge is less than 7 Å. 68 Furthermore, these calculations give the stable spiro state of the vacancy even at a distance of 15.6 Å from the edge, 68 contrary to in infinite graphene where it corresponds to a saddle point. 161 Note that calculations using the old version of the Brenner potential gives the 5/9 structure for the ground state of the vacancy whereas the new version of this potential gives the spiro structure. 154 The influence of the vacancy on the surrounding structure is more significant near the edge than in infinite graphene. For flakes consisting of approximately one hundred carbon atoms, the presence of a vacancy in the centre resulted in considerable bending of the entire structure, for both flakes with dangling bonds at the edge 67 and flakes terminated at the edge by hydrogen atoms. 69,153 Vacancies in the spiro state in the middle of narrow graphene nanoribbons causes them to coil into spiral helices. 68 The vicinity of the edge has an even more drastic effect on vacancy migration than on its structure and energetics. 67,68 As both 5/9 and spiro states of the vacancy near the graphene edge correspond to local potential energy minima, the only vacancy migration step near the edge is the transition between the 5/9 and spiro states. 67 The scheme of possible vacancy migration steps near the edge is presented in Fig. 4. A transition from the initial 5/9A state to two spiro states SA and SB is possible. Four transitions are possible from each spiro state SA and SB, which can be achieved by breaking one of the four bonds: a return to the initial 5/9A state, transition to another 5/9 state of Table 3 The DFT-calculated energy characteristics of flakes with a monovacancy located at distance d (in Å) from the edge: E v (in eV) is the energy of vacancy formation, and DE (in eV) is the difference in the total energy between a flake with the vacancy in the middle and at the edge.  View Article Online the same vacancy (5/9B and 5/9F) and migration of the vacancy (5/9C, 5/9D, 5/9E, and 5/9G). 67 The energies of local minima and transition states for vacancy migration in the middle of the C 116 67 and from the middle to the edge of the C 41 H 16 68 graphene flakes have been calculated. These calculations reveal that the barrier for vacancy motion towards the edge is considerably lower than the barrier for motion away from the edge. Thus, both the difference in the total energy of graphene with the vacancy at the edge and at distance d from the edge (see Table 3) and the barrier to vacancy migration decrease with the decrease of distance d. This leads to the conclusion that thermally activated directional motion of the vacancy towards the edge should take place, 67,68 as confirmed by molecular dynamics simulations. 67 As these barriers for the vacancy motion towards the edge are smaller than barriers for vacancy motion in infinite graphene, the motion to the edge will occur too fast to be observed by STM or TEM studies at room temperature. According to DFT calculations, a 3% compression of graphene leads to a change in the ground state of the vacancy to the spiro state. 51 The deformation of graphene also exerts considerable influence on the barriers DE m for vacancy migration. DFT calculations show that the value for pristine graphene DE m = 1.17 eV changes to 3.15 and 0.75 eV for elongation by 5% in the armchair and zigzag directions, respectively, and to 0.12 and 1.79 eV for compression by 5% in armchair and zigzag directions, respectively. 23 Thus the possibility of directional motion of vacancies caused by applied graphene layer deformation has been proposed. 23 3.1.7 Interaction and coalescence of vacancies. A set of studies is devoted to the interaction of vacancies and their coalescence into divacancies. 50,142,157,160,162,173 This set includes a TB molecular dynamics study of the convergence and coalescence of vacancies, 173 action-derived molecular dynamics for energetics of certain paths of vacancies migrating towards each other, 142 DFT calculations of energies of the subsequent states at the convergence of vacancies and barriers between these states, 160,162 and barriers of the final step of coalescence of pairs of vacancies. 157,160,162,173 The formation energy of vacancies in graphene as a function of their concentration has also been studied. 50 Thorough DFT calculations give the energies of pairs of vacancies for different neighbour locations, as well as barriers to transition between these locations. 160 The barriers for the final step of coalescence of a pair of vacancies into a V 2 (5-8-5) divacancy are found in this work to be 2.10 and 0.82 eV, for different final neighbour locations 160 162 Two examples of energy changes for the migration towards and coalescence of vacancies, based on a TB model, are shown in Fig. 5. 142 These calculations give the values of the barrier for the final step of coalescence as 1.3 and 1.9 eV for different approaching pathways with V 2 (555-777) divacancy formation. Detailed DFT-calculated maps of energies of vacancies located near the V 2 (5-8-5) divacancy and a trivacancy, as well as transition barriers between neighbouring vacancy positions, have been used to study the dynamical evolution of the coalescence of monovacancies with divacancies and trivacancies by kinetic Monte Carlo simulations. 161 Analogously to the coalescence of a pair of monovacancies, the barriers for the final stage of coalescence depend on the pathway of the vacancy approaching the di-and trivacancies. These barriers lie in the range 0.9 to 4.7 eV and 0.9 to 1.5 eV for different neighbour positions for the V 2 (5-8-5) divacancy and trivacancy, respectively. 161 Examples of the final steps of coalescence of the monovacancy and V 2 (5-8-5) divacancy are shown in Fig. 6.

Adatoms
The healing of various defects of graphene by the migration of adatoms has been observed by HRTEM. 64,65 Possible magnetism of adatoms on graphene has been discussed, 179,180 and the structure and energetics of the adatom on graphene have been considered in a set of papers. 157,163,[179][180][181][182][183] A DFT study gives three stable positions of the adatom on graphene: a bridge position over the centre of a bond, a dumbbell arrangement where the adatom has bonds with three atoms of graphene, and an off-top structure where the adatom and one atom of graphene each have four bonds and the same position relative the graphene plane. 181 The bridge position is found to be a ground state, 166,179,181 with the dumbbell arrangement and offtop structure having 0.37 and 0.22 eV higher energies, respectively. 181 The bridge position of an adatom on graphene was observed by HRTEM 184 and scanning transmission electron microscopy (STEM). 185 182,183 The migration of adatoms can therefore easily take place at room temperature. DFT calculations predict an energy gain of up to 0.25 eV for close positions of neighbour adatoms, with a decrease of the barrier to migration to only 0.225 eV, 180 meaning that the agglomeration of adatoms at low temperatures is possible. The structure and energetics of adatom pairs on graphene (i.e. the case where chemical bonds exist between adatoms) have been also considered by DFT. 181 The most stable pair is found to be the so-called 7-5-5-7 defect that is comprised of two pentagons and two heptagons and introduces a local elevation of 2 Å out of the graphene plane.
As monovacancies can easily migrate and merge into divacancies, whereas thermally activated migration of divacancies is not possible, the interaction of adatoms with divacancies is an important process observed by HRTEM. 184 DFT energy calculations of states of an adatom and divacancy approaching one another, as well as the barriers between these states, have been performed. 162,163 For the V 2 (5-8-5) divacancy, the coalescence of the divacancy and adatom takes place with three barriers during the approach of the adatom (0.91, 0.73 and 1.36 eV) before the coalescence and a total energy release of 5.79 eV, 162 whereas a different study gives even higher barriers of 1.55 and 2.49 eV between neighbouring bridge positions near the V 2 (5-8-5) divacancy. 163 Thus a state in which an adatom near a V 2 (5-8-5) divacancy is stable at room temperature is possible, in agreement with HRTEM observations. 184 However, the adatom does not coalesce with the V 2 (555-777) divacancy -a dumbbelllike configuration is formed where the adatom is adsorbed over the central atom of the V 2 (555-777) divacancy. 162 The motion of the adatom to this final location occurs with the subsequent barriers 0.43, 1.00, 0.78, 0.39 and 0.56 eV and the total energy release is 1.8 eV. The creation of an inverse Stone-Wales defect on the periphery of a divacancy, by one adatom hopping within the vicinity of another, was proposed. 163

Topological defects and bond-realignment reactions
An important class of reactions in carbon nanostructures is related to bond rotations that lead to the formation or transformation of topological defects, while not being accompanied by a loss or addition of atoms. The simplest example of such a reaction is the 90 degree rotation of a bond in a perfect hexagonal network, leading to the formation of a Stone-Wales (SW) defect 72,168 comprised of two pentagons and two heptagons (Fig. 7a).
3.3.1 Effects of SW defects on mechanical, electronic and magnetic properties. SW defects only very weakly affect the elastic properties of graphene such as the Young modulus [26][27][28][29][30] and Poisson ratio, 30 although this effect becomes more pronounced with increased defect concentrations. [26][27][28] Numerous studies show that pre-existing SW defects serve as nucleation centres for fracture 22,27 and substantially decrease the failure strain and tensile strength of graphene. 22,27,29,31,32 The accumulation of SW defects in graphene degrades the ultimate tensile strength and failure strain to saturated levels that are 30-50% lower than in pristine graphene. 22,27 The SW transformations are also direct participants of plastic deformation in graphene, which is discussed below.
The analysis of the density of states in graphene showed that p z orbitals of the carbon atoms of the rotated bond in SW defects give rise to a defect band B0.5 eV above the Fermi level, with the width and height dependent on the defect concentration. [33][34][35][36][37][38][39][40] The conductivity of defected graphene exhibits systematic degradations around defect resonance energies. 33,34 Electron backscattering regions are also introduced in transmission spectra of graphene nanoribbons with SW defects. 41-44 SW defects generally do not carry a magnetic moment [33][34][35]40 and are neutral, 40 although some charge transfer was detected between layers of defected bilayer graphene 35 and asymmetric arrangements of SW defects in graphene nanoribbons can affect their magnetic moment. 55 The effect of SW defects on the thermal conductivity of graphene 26,27,57,58 and graphene nanoribbons [59][60][61] is similar to that of other local defects, such as vacancies.
3.3.2 Structure and energetics of SW defects. The combined DFT 35,36,54,70,164,166,[186][187][188] and TB data 164,189 show that the formation energy of a SW defect in a graphene layer is 4.5-5.0 eV (Table 4). A slightly larger formation energy (by 0.3 eV) is predicted on the basis of quantum Monte Carlo (QMC) simulations. 186 Significant variations 35,36,53,54,70,91,119,152,164,166,173,[186][187][188][189][190][191][192][193][194] are observed depending on the supercell used in the calculations (Table 4), due to interactions of periodic images of SW defects through generated long-range stress and strain fields. 190,191 In the calculations for finite graphene flakes, the formation energy [195][196][197][198] of SW defects might be affected by the vicinity of the edges. While the results of TB methods 119,121-123 are very close to those of DFT, 188 the semiempirical potentials [127][128][129] and continuum models 199 are not very accurate in exact numerical values for energies. 70,71,138,[200][201][202] However, they give the correct trends and qualitative behaviour for bond realignment reactions. 70,91,188,[202][203][204] It is seen that although all atoms in a SW defect have the same number of bonds as in the perfect hexagonal network, it has a significant formation energy. The reason is that upon the SW rotation, many bonds get compressed or stretched 36,152,164,166,186,198,205,206 (Table 5), and the angles between the bonds deviate from 1201. The angle of the pentagons at the atoms forming the rotated bond was calculated to be 115.51, 36 115.01 206 and 118.21. 205 The most deformed bond is the central rotated bond (separating two heptagons in the SW defect), which experiences compression by almost 10% (Table 5 and Fig. 8a). 152,164,166,186,198,205,206 Some compression is also observed for the bonds between pentagons and hexagons, while the bonds between heptagons and hexagons and the bonds between heptagons and pentagons are stretched 152,164,166,205,206 (Fig. 8a). This behaviour is in agreement with general observations that pentagons and heptagons are centres of compression and tensile stress, respectively. [207][208][209][210]   Geometrical optimisation of a SW defect starting with all atoms in the same plane yields a flat SW defect structure. 119,152,166,190,191,205 However, this structure has two imaginary frequencies, and the true minimum is characterised by the carbon atoms forming the rotated bond at the defect core buckling out-of-plane above and below the graphene plane 35,36,186,194 (Fig. 7b). The energy gain from out-of-plane buckling of the core atoms is 60 -300 meV depending on the method used and the supercell considered, 35,36,186,194 and comes from slight elongation of the compressed rotated bond at the defect core. 186 The maximum displacement between atoms perpendicular to the graphene plane in the buckled energy minimum configuration is 0.8-1.6 Å (ref. 35, 36 and 186) (Fig. 7b). Table 4 Formation energy E SW , activation energy E a,SW , energy gain from out-of-plane buckling dE b (in eV) and maximum difference in out-of-plane coordinate of atoms h (in Å) for SW defects in single-layer and bilayer graphene and graphite. The system considered (periodic boundary conditions (PBC) or flake), number of atoms in the computational cell per layer and calculation method are indicated The analysis of correlation effects in the formation energy of two SW defects has revealed that both repulsion and attraction between two SW defects is possible, depending on their relative position and orientation. 138,201 Attraction valleys correspond to a 'diagonal' location of defects with respect to each other for both aligned 138 and misaligned defects. 201 The most energetically favourable configuration is immediately adjacent but not overlapping.
The formation of SW defects in bilayer graphene with the AB stacking of graphene layers (AB-BLG) 35,54 and with twisted graphene layers (TBLG) 35 was investigated. It was shown that in bilayer graphene, SW defects are also stabilised by out-ofplane buckling, but the magnitude of this distortion is reduced to 0.7 Å compared to 1.2 Å in single-layer graphene (SLG), 35 i.e. the presence of the second layer inhibits the distortion. The stabilization energy from buckling is also decreased from 220 meV in SLG to 40 meV in AB-BLG and 60 meV in TBLG. 35 The defect formation energy in TBLG was found to be almost the same as in AB-BLG, while there was only slight variation (by tens of meV) for different positions of the SW defect due to the differences in the interlayer coupling. 35  Typical DFT values for the activation energy of the SW transformation lie in the range 9-11 eV 70,91,119,164,166,188,[193][194][195][196] (Table 4 and Fig. 8b). The barriers calculated with semiempirical potentials can be substantially underestimated, 70,71,138,197,200,202 although they are still large enough to reflect that two carbon bonds are broken at the same time in this transformation. The transition state for the SW transformation corresponds to a rotation of the carbon bond by 45-551 (Fig. 8b). 119,138,152,166,194,195,200 As in the final SW defect, compressive stress in the transition state caused by the shortened rotating bond (1.29 Å) 194 and other bonds of the central carbon atoms (1.36 Å) 194 is relieved by outof-plane buckling 138,194,195,200 with an energy gain of about 0.9 eV, 194  According to the Tersoff-Brenner potential, 128 the reaction path S À+ is the most probable, followed by S + , S ++ and finally S 0 . 138,200 However, DFT studies 194,195 predict that the true transition state with only one imaginary frequency corresponds to the S + structure with the central bond tilted by 151 195 or one of the central atoms buckling out by 0.95 Å. 194 The maximum displacement between atoms perpendicular to the graphene Table 5 Geometry of SW defect in single-layer graphene: lengths (in Å) of bonds between two heptagons (r 77 ), between a heptagon and a pentagon (r 57 ), between a pentagon and a hexagon (r 56 ) and between a heptagon and a hexagon (r 76 ). The system considered (periodic boundary conditions (PBC) or flake), cell size, number of atoms in the computational cell and calculation method are indicated  plane in the buckled transition state configuration is of the order of that in the final SW defect, at 1.42 Å. 194 It has been shown that the barrier to the SW transformation can be considerably reduced in the presence of carbon adatoms, 193,196,198 hydrogen and hydroxyl groups 194 or transition metals. 91 Participation of a nickel atom reduces the activation energy of the SW transformation in graphene by 1.8 eV. 91 One hydrogen atom, two hydrogen atoms and two hydroxyl groups reduce the activation energy by 2.5 eV, 5.6 eV and 3.4 eV respectively. 194 The effect of a carbon adatom is even more pronounced as the corresponding mechanism involves only one bond breaking at a time and the activation energy is decreased by 5-7 eV (Table 4). 193,196,198 3.3.3 Modification of SW energetics by strain. It is seen that the activation energy for a SW transformation is very high, both for the formation of non-hexagonal rings and their annihilation. Therefore, SW defects are unlikely to be formed in pristine graphene even at high temperatures. However, additional manipulations such as electron irradiation or mechanical deformation can induce the formation of SW defects. As the rotated bond of the SW defect is highly compressed, stretching a graphene layer along the direction of the rotated bond leads to a decrease in the formation energy of SW defects 138,188,193,195,[199][200][201][202][203][204]206,211 and activation energy of SW transformations. 31,138,188,193,195,200,202,211 The formation of SW defects can therefore be facilitated by strain, and it has been considered as the first step in the plastic deformation of carbon nanotubes and graphene. The energy gain associated with the rotation of a bond depends on its orientation, and increases with a reduction in the angle of the rotated bond with respect to the loading direction. 31,138,188,195,[199][200][201]206,211 This implies that relieving lateral strain in graphene by SW rotations is most efficient when the layer is stretched along the zigzag direction. The critical strain at which formation of SW defects in graphene becomes thermodynamically favourable was calculated to be 6-9% 188,193,199,202 and 12-17% 188,199 when the strain is applied in the zigzag and armchair directions, respectively. The annihilation of unfavourably oriented SW defects under strain has been demonstrated with MD simulations. 212 Although both the formation energy of SW defects and activation energy of SW transformations are reduced by several eV at the critical strain, such transformations are still kinetically limited 138,188,193,195,200,202,211 and greater strains should be applied to observe them experimentally. The Arrhenius equation was used to estimate that the yield strain at which the formation of SW defects should take place during experiments at room temperature is on the order of 17-22% 195 for different loading directions in graphene or nanotubes of different chirality, in agreement with the experimentally measured yield strains for carbon nanotubes up to 10% 213 and 17%. 214 The minimum activation energy was reached for an intermediate angle of the rotated bond between 01 and 301 with respect to the loading direction. Calculations 200 on the basis of the Tersoff-Brenner potential 128 predicted that this angle is 111, while DFT calculations 195 for nanotubes gave an estimate of 221. The approximate formation energy of SW defects 138,201 and activation energy for SW transformations 195,200 have been proposed as functions of strain magnitude and angle of the rotated bond with respect to the loading direction.
3.3.4 Dislocations. Nucleation of a SW defect ''unlocks'' carbon nanostructures for further plastic relaxation under strain: 215 (1) plastic flow 138,203,204,211,216,217 and (2) the generation of multiple SW defects 138,203,204 and large rings, such as octagons. 138,216 The first scenario that explains experimentally observed superelongation of carbon nanotubes under strain 18,19 is realised by a separation of pentagon-heptagon (57) pairs 184,218,219,220 (Fig. 10) that represent crystal dislocations 138,190,191,203,204,211,216 (a SW defect is a dislocation dipole in which two edge dislocations with opposite Burgers vectors are displaced by one lattice unit). The glide of 57 pairs, i.e. displacement by one Burgers vector at a time, is possible by a SW rotation of the bond between the heptagon and the hexagon adjacent to the pentagon (Fig. 10A, B, D, E). In unstrained graphene, 138,190,191,211,218 separating 57 pairs with one intervening hexagon costs 1-6 eV, although the energy required decreases with further separation of the 57 pairs 211,216,218 (Fig. 11). As with the initial formation of a SW defect, separation of 57 pairs is stabilised by the application of an appropriate load. 138,203,204,211,216 Therefore, the long range strain fields favour gliding of 57 dislocation cores in opposite directions, while the local curvature energy of two oppositely directed 57 pairs is minimised when they are adjacent and form a SW defect. 138,190,191,203,204,211,216,218 The glide of 57 pairs can also be facilitated by the formation of chains of dislocation dipoles (so-called ''dislocation worms''), which effectively screen stress fields of the separated 57 pairs. 203,204 In addition to the glide of dislocation cores by SW rotations, another type of motion, a climb step, is possible. Emission of a pair of carbon atoms from the edge between the pentagon and the hexagon adjacent to the heptagon 217 (Fig. 10B, C, E, F) results in motion along the defect axis in the direction from the heptagon to the pentagon. This was shown to be responsible for nanotube sublimation that preserves perfect hexagonal structure even upon significant loss of mass. 224 A combination of the two mechanisms, glide and climb, explains the motion of 57 pairs along any arbitrary trajectory. 217 A combined motion of 57 pairs in graphene by both glide and climb 218,225 (Fig. 10) as well as by glide, assisted by the formation of ''dislocation worms'', 225 was demonstrated for graphene experimentally and is in agreement with HRTEM observations 219,226 for carbon nanotubes.
Unlike in nanotubes, in graphene the formation of a second SW defect is preferable to the dissociation of the first SW defect into 57 pairs, both under strain and in its absence. 138,203,204 Simple dissociation of SW defects into 57 pairs in graphene is therefore thermodynamically unfavourable, and in most cases the plastic deformation of a graphene layer should start through the formation of additional defects, multiple SW defects or larger rings such as octagons. Once the defects are present in graphene, subsequent SW rotations become much easier. For example, the activation energy of the transformation of the V 2 (5-8-5) structure to V 2 (555-777) was found to be 5.74 eV, according to the TB calculation 173 (5.17 eV according to LDA calculation), which is almost half the 10.2 eV activation energy for the formation of a SW defect in perfect graphene, obtained in the same paper. It is important to note, however, that plastic deformation through the creation of multiple defects can also finally lead to the formation of isolated 57 pairs. For instance, chains of SW defects are known to rearrange giving 57 pairs screened by ''dislocation worms''. 203,204 Deformation of graphene under strain, nevertheless, is hardly the main source of lone separated 57 pairs so frequently observed in HRTEM. 218,227,228 Creation of dislocations under HRTEM conditions can be explained by the reconstruction of vacancies and adatoms introduced by electron irradiation.  Tight-binding molecular dynamics simulations for carbon nanotubes with multi-vacancies demonstrate that such nanotubes are sewed up almost perfectly with only two separated 57 pairs left and a reduced diameter. 92,93 The transformation of graphene with a chain of missing atoms to a nearly ideal layer with only two separated 57 pairs was shown to be preferred over the formation of a local haeckelite structure composed of collective 555-777 defects when the multivacancy size exceeds ten, 176 in agreement with the experimental data. 225 The formation of separated 57 pairs was found to be favourable over other possible reconstructions for multivacancy sizes above 4 in narrow graphene nanoribbons of width within 1 nm, and above 26 in wide graphene nanoribbons of width greater than 10 nm. 175 According to the experimental data, 225,228 however, even scattered vacancies tend to reconstruct and form a nearly ideal hexagonal network with only two separated 57 pairs. Divacancies can dissociate into 57 pairs by glide and pseudoclimb similar to SW defects. 217,227 However, 57 pairs formed in this way turn out to be unstable, as seen in tight-binding molecular dynamics simulations. 228 Stable dislocations are rather produced by agglomeration and collective reconstruction of several vacancies assisted by carbon adatoms helping to convert non-hexagonal rings into hexagons, as follows both from computer simulations 228,229 and experimental observations. 228 The geometry of stable dislocations depends critically on whether they are formed from adatoms or vacancies. 225,228 The dislocation cores formed from adatoms (i.e. pointing towards each other) have only shallow minima at nanometre distances associated with the antisymmetric hillock-basin configuration, while at short distances they prefer the symmetric hillockhillock configuration and tend to migrate towards one another by SW rotations. 225 The dislocation cores formed by multi-vacancy reconstruction, on the other hand, repel at short distances and have only one minimum at the distance determined by the number of vacancy units, in which the antisymmetric hillockbasin configuration is slightly preferred over the symmetric hillock-hillock one. 228 In addition to 57 pairs, double pentagon (octagon) rings have been recently identified as dislocation cores in graphene. 230 As with separated 57 pairs, these defects are formed by the reconstruction of vacancy constellations and move by glides involving one bond rotation at each step. 230 3.3.5 Graphene grain boundaries. The non-hexagonal rings in polycrystalline graphene play an important structural role as grain boundaries 20,231-236 between arbitrarily oriented graphene flakes. Dislocations characterised by any Burgers vector as well as grain boundaries, covering the whole range of possible misorientation angles, can be constructed on the basis of coupled ((1,0) dislocation) and dissociated ((1,1) dislocation) 57 pairs separated by hexagons [221][222][223]237 (Fig. 12). The stress fields of isolated dislocation cores in the limit of misorientation angles y -01 and y -601, which correspond to small-angle armchair and zigzag regimes, respectively, induce strong buckling of graphene with protrusion heights on the order of 2-3 Å. 221,222 The possibility of such buckling in free-standing graphene results in a considerable reduction of the grain boundary energy. 203,204,221,238,239 These theoretical results are supported by experimental STM images Buckling is also predicted to take place in disordered grain boundaries. 240 The energies of the most favourable straight grain boundaries covering all angle ranges lie within 0.5 eV Å À1 (Fig. 12e), [221][222][223]238,239,241 i.e. they are smaller than typical values for bare graphene edges (1.0-1.2 eV Å À1 ). 96,164,[242][243][244][245][246][247][248] Especially stable large-angle grain boundaries are formed by close packing of (1,0) (LAGBI) and (1,1) (LAGBII) dislocation cores (Fig. 12c-e). 221 These grain-boundaries with misorientation angles 21.81 and 32.21 have formation energies of only 0.338 eV Å À1 and 0.284 eV Å À1 , respectively. 221 Though they contain the maximum number of dislocation cores, their stress fields mutually cancel each other providing the lowest formation energy and maximum mechanical strength. 208,221,249,250 Unlike in small-angle grain boundaries, these structures are flat, 220,221 as confirmed by the experimental data. 233 The formation energy of LAGBII per 57 pair is almost half that of SW defects (1.3 eV vs. 2.5 eV according to ref. 251), leading to its abundance in experimental observations [233][234][235][236] and realisation of grain boundaries at other large angles as pieces of LAGB II separated by kinked sites. 235 Similar results were obtained for rotational grain boundaries, where the structure formed by a loop of six alternating pentagons and heptagons (the flower defect) was identified as the topological defect with the lowest formation energy per dislocation core (1.2 eV), 251 explaining why its formation is frequently observed via coalescence of mobile dislocations or SW defects. 20,225 The glide and climb of 57 pairs constituting grain boundaries 222,252 make their migration possible. 20 The driving force for this migration depends only on the local in-plane boundary curvature and does not depend on the atomistic structure, providing that curved grain boundaries tend to become straight, while grain boundary loops annihilate upon annealing. 20 Tightbinding molecular dynamics simulations 252 of grain boundary motion showed that C 2 dimer emission, corresponding to a climb of the 57 pairs, is preceded by the formation of an adatom.
3.3.6 Effects of edges on the formation of SW defects. The close vicinity of graphene edges allows for a more efficient relaxation of SW defects, providing a decrease in both the formation energy of SW defects and the activation energy for this process. 70,71 The formation energy of SW defects was shown to reduce from 4.4 eV inside a large graphene flake to 2.3 eV at the zigzag flake edge 70,71 according to the Tersoff-Brenner potential 128 and from 5.2 eV to 3.1 eV according to DFT calculations. 70 The corresponding changes in the activation energy were from 7.9 eV to 6.7 eV according to the semiempirical potential 70,71 and from 9.9 eV to 6.8 eV according to DFT calculations. 70 The formation energy and barrier for SW defects at the armchair edge of the flake were calculated to be 2.8 eV and 7.4 eV, respectively. 71 3.3.7 Edge reconstructions. In addition to the effect on the energetics of reactions inside graphene layers, edges enable new reactions with the participation of under-coordinated atoms at the very edge. As these atoms are already destabilised compared to ones in the graphene interior, the corresponding reactions have considerably reduced barriers and reaction energies 70,71 This means that such processes can occur at experimentally accessible times at high temperatures even without the assistance of catalysts, strain or irradiation, as opposed to SW transformations inside the layers.
In particular, reconstruction of the zigzag (ZZ) edge via the transformation of pairs of adjacent hexagons to pentagon-heptagon pairs (Fig. 13a) is energetically favourable 77,96,148,164,[242][243][244][245][246][247][248]253 and is often observed experimentally. 13,75,254,255 This reconstruction (ZZ(57)) decreases the zigzag edge energy by 0.1-0.4 eV Å À1 , making it even slightly more stable than the armchair (AC) edge (Table 6). These results can be understood by considering edge geometries. The pristine zigzag edge is known to be much more expensive energetically than the armchair one, as twocoordinated atoms at the zigzag edge have dangling bonds, while two-coordinated atoms at the armchair edge form triple bonds, as seen by short bond lengths between such atoms (1.22-1.25 Å, 164,242,243,245 close to the 1.20 Å triple bond in acetylene (Table 7)). The same triple bonds of 1.23-1.26 Å length 78,164,242,243,245,247,253 are observed between two-coordinated atoms in the reconstructed ZZ(57) zigzag edge, guaranteeing its stability ( Table 7). The initial compressive stress of the pristine zigzag edge of À0.5 eV Å À1 , 247 À2.05 eV Å À1 , 76 and À2.09 eV Å À1 248 is completely reversed to a tensile stress of 1.5 eV Å À1 , 247 0.24 eV Å À1 , 76 and 0.14 eV Å À1 248 in the reconstructed ZZ (57) edge that provides its planarity. The 100% reconstructed edge also becomes non-spin-polarised. 73,247,248 Reconstruction of the armchair edge through the transformation of three adjacent hexagons to two heptagons and a pentagon (an incomplete SW defect, Fig. 13a, AC(677)) was shown to be energetically unfavourable although still preferred over the unreconstructed zigzag edge (Table 6). 243,247 Nevertheless, it was also predicted that partial reconstruction of the armchair edge can further stabilise it. 247 The minimum edge energy of 0.98 eV Å À1 , very close to the edge energy of the reconstructed ZZ(57) zigzag edge, was calculated for about 25% of hexagons transformed into 6757 fragments. 247 Such a structure was also shown to have a small edge stress, indicating that it is stable against deformations, while the fully reconstructed armchair edge has a considerable tensile stress of 4 eV Å À1 . 247 It should be emphasised, however, that these considerations are only valid at very low hydrogen pressures, while at ambient conditions the edges are passivated and no dangling bonds are left. 246 In this case, reconstruction of the edges has a considerable energy cost, while the edges consisting only of hexagons and with two-coordinated atoms bonded to hydrogen are the most preferred. 192,246,253,256 The energetics and warping of various hydrogen-passivated reconstructed edges was analysed. 192,246,256 The reconstruction of the zigzag edge into a sequence of alternating pentagons and heptagons is also unfavourable on transition metal surfaces, while other new types of reconstructions with an increased number of low-coordinated carbon atoms become possible. 257 Decoration of graphene edges with silicon 258 and other adatoms 259 is also known to affect their structure.
The activation energy for simultaneous and complete reconstruction of the zigzag edge is proportional to the edge length, and was calculated to be 0.6 eV per hexagon pair. 243 The simulations using a semiempirical potential 245 gave a rather close value of 0.7 eV at the low temperature limit. 245 Though it was shown to have non-monotonic behaviour at high temperatures (increasing up until 700 K and slowly decreasing at higher temperatures 245 ), the cost for this simultaneous transformation remains too high and the transformation actually occurs stepby-step via the formation of pentagon-heptagon pairs. Therefore, sequential steps of the process have to be analysed.
3.3.8 Bond rotation reactions at graphene edges. The formation energy of a 57 pair at the zigzag edge of a 12-atom square hole in graphene was found to be 1.35 eV by firstprinciples calculations, much smaller than 3.83 eV for the SW transformation inside the layer. 187 The formation energy of a 57 pair at the pristine zigzag edge of a graphene flake was found to be about 1 eV according to the semiempirical potential, 128 3.4 eV smaller than for a SW defect inside the graphene flake. 70,71 The activation energy of this reaction was found to be around 3 eV, which is less than half that of the SW transformation inside the layer. 70,71 Such a decrease was explained by the fact that bond switching at the edge requires only one bond to be broken at a time, compared to two in the perfect graphene layer. A consecutive energy decrease by 1.51 eV and 1.78 eV was observed in successive steps of the generation of 57 pairs 96 according to TB and unpolarised LDA calculations, which are known to underestimate the stability of the pristine Table 6 Energies (in eV Å À1 ) of pristine zigzag (ZZ), pristine armchair (AC), reconstructed zigzag (ZZ(57)) and reconstructed armchair (AC(677)) edges of graphene. The calculation method (SP and NSP stand for ''spinpolarised'' and ''non-spin-polarised'' calculations, respectively) is indicated     77,245,247,248 Activation energies of just 0.69 eV and 0.44 eV for each respective step were obtained, indicating that such a transformation should take place spontaneously at room temperature. 96 A similar free energy barrier of 0.83 eV was obtained at room temperature using a semiempirical potential. 245 The barrier was shown to increase up to 1 eV with increasing temperature. 245 Recent spin-polarised DFT calculations have, however, refined these values. The energy increase by about 0.3 eV was found for the formation of the first 57 pair, followed by a total energy gain of 0.6 eV from the generation of the second one 77 (Fig. 14). An activation energy of 1.80 eV (ref. 77) or 1.12 eV (ref. 73) was obtained for the first step, and 1.11 eV for the second one 77 (Fig. 14), demonstrating that the pristine zigzag edge is stable at room temperature, as confirmed by experiments. 72,74,260 An energy decrease of 0.21 eV was predicted for the formation of the first 57 pair in ref. 79 and the activation energy for this reaction was found to be 1.61 eV. It was also suggested that the formation of these pairs can occur through the recombination of an adatom and a vacancy, with an activation energy of just 0.89 eV. 79 The barriers for the migration of an adatom and vacancy on the zigzag edge were calculated in the same paper to be 0.76 eV and 2.12 eV, respectively. Tensile stress was shown to increase energies of the pristine and defected zigzag edges, decreasing the activation energy for the firststep transformation with 0.21 eV per % rate. 73 Therefore, a uniaxial strain of 5% is sufficient to provoke reconstruction of the zigzag edge at room temperature. 73,77 Fracture of a monolayer graphene was shown to be governed by the competition between bond breaking and bond rotation at a crack tip. 261 Using atomistic reaction pathway calculations, a kinetically favourable fracture path was identified as an alternating sequence of bond rotation at the graphene edge and bond breaking. 261 The generation of 57 defects was also shown to play a crucial role in graphene evaporation, as the emission of carbon atoms was found to proceed preferentially from heptagons. 96 In particular, the energy gain from the formation of several 57 pairs at the zigzag edge was revealed to be sufficient to compensate for atom evaporation. 96 Similar energetics are observed for other defects and graphene edges. The formation of a 757 fragment on the armchair edge of a hole in graphene was calculated to be 2.01 eV, which is significantly reduced as compared to the formation energy of a SW defect. 187 The energetics of various defects at zigzag and armchair edges of graphene flakes were investigated. 70,71 Similar barriers on the order of 3 eV were obtained for all reactions proceeding through the breaking of only one bond, including diffusion of pentagons and heptagons along the edges. Reactions forming carbon chains by breaking a bond between adjacent rings were recognised as important for all types of edges at high temperatures. Though the potential energy change in such reactions is relatively high, on the order of 2 eV, the formation of carbon chains becomes favourable at high temperatures due to an increase in entropy, related to their low-frequency vibrations. 70,71 The formation of carbon chains and their reconnection followed by edge reconstruction was shown to be the principal mechanism of transformation of a graphene flake to a fullerene 70,71,87 and formation of nanotube caps from open nanotube ends. 88 The transformation of narrow graphene nanoribbons into carbon chains was observed experimentally under electron irradiation. 13,14 The formation of carbon chains also took place in simulations of strained graphene nanoribbons at high temperature. 262,263

Reactions under irradiation
While the important values when discussing thermally activated reactions are the activation and formation energies, threshold energies are key for irradiation-induced reactions. The threshold energy of a process is the minimum energy required for that process to occur without immediate reversal, considering a specific direction of atom motion induced by the electron beam. The difference between the activation energy and threshold energy for the same reaction can be explained as follows. Firstly, as the reaction pathway for an irradiation-induced reaction is determined by the direction of momentum transferred to a given atom from the incident electron, it does not necessarily coincide with the optimal reaction way for the thermally activated reaction. Secondly, even for optimal directions of transferred momentum (at the minimum threshold energy), the collective motion of multiple atoms will usually not coincide with the optimal reaction way for thermally activated reactions. The threshold energy of an irradiation-induced process therefore tends to be higher than the activation energy of the equivalent thermally induced process.
Under the electron beam, graphene is constantly exposed to an external stimulus. In these non-equilibrium conditions, low formation energies do not therefore imply the structures most likely to form and survive. In fact, the presence of possible routes induced by the beam towards a structure and the subsequent stability of the structure against the electron irradiation are key. The preferred configuration of graphene edges in TEM is a good example of the importance of these considerations. Two main edge configurations are observed: zigzag (ZZ) and armchair (AC). Although early observations and calculations 260 concluded that the ZZ configuration would be the most stable, the consideration of dynamic effects following electron impacts demonstrated the very high stability of AC edges with respect to irradiation. 148 For another demonstration of this effect consider the structure of the most commonly observed tetravacancy, the extended linear structure. 264 It was proposed that the high occurrence of the most common V 4 structures in TEM images is due to the specifics of the irradiation-induced pathways towards their creation, but an explanation for the long lifetime of the linear structure over the more frequently modelled structures is perhaps due to the shared structural features between this vacancy and the V 2 (5555-6-7777) divacancy (or 'butterfly' vacancy). It has been shown that the central atoms of this divacancy have an even larger emission threshold than pristine graphene, 147 implying a very strong stability against electron irradiation. It would therefore be expected that this is also true of the very similar tetravacancy structure.

Atom emission
4.1.1 From pristine graphene. The emission of carbon atoms from pristine graphene in TEM has been extensively studied, both theoretically and experimentally. It has been observed that in order to eliminate atom emission and the creation of holes and defects, an accelerating voltage of 80 kV or lower must be used. This corresponds to a maximum transmittable energy to a stationary carbon atom of 15.8 eV, implying that the emission threshold energy of carbon must be slightly above this. However, Meyer et al. 146,267 demonstrated the importance of including the effects of the thermal motions of the graphene lattice when considering the transfer of momentum from beam electron to sample atom, explaining the apparent overestimation of theoretical emission thresholds compared to experiment. This inconsistency had previously been attributed to the absence of any beam-induced electronic excitations in calculations, 150,265 but was revealed to be an effect of the increased transmittable energy from the beam due to thermal motion of the atoms; the threshold energies are generally accurate when taking this into consideration (Table 8). An experimentally derived rate of atom emission from graphene was found at various accelerating voltages and dose rates, and used to calculate the threshold energy of 23.6 eV including these lattice vibration effects. 146,267 It is worth noting that this is higher than any published theoretical value, by at least 0.6 eV. A comparison 266 of theoretical threshold energies calculated by DFT to excited state time-dependent DFT (TDDFT) calculations that explicitly included the electron dynamics concluded that the ground state Born-Oppenheimer dynamics universally used when calculating threshold energies is indeed justified. Qualitative trends were very similar between DFT and TDDFT, despite small quantitative differences; DFT resulted in the transfer of slightly more energy to the surrounding graphene lattice, while TDDFT showed similar kinetic behaviour of the emitted atom and the graphene lattice at a delay of 10 fs.
Although areas of pristine graphene are generally considered to be protected from irradiation damage at or below 80 kV, the use of very high electron dose rates (410 8 e nm À2 s À1 ) at 80 kV has been shown to create defects. 65 Similar defect creation is not witnessed under lower dose rates (o10 6 e nm À2 s À1 ) and long irradiation times, where large total numbers of electrons are passed through the sample (B10 10 e nm À2 ), 146,267 meaning that a dose-rate dependent effect is either lowering the threshold energy or increasing the energy that can be transferred by increasing the out-of-plane atomic vibrations. This lowered threshold energy is calculated to be 19.7 eV, and weakening of the bonds due to ionization or plasmon excitations was suggested as a possible cause. Alternatively, while beam-induced heating effects were shown to be not responsible, a flexural phonon mode induced by inelastic collisions could result in the required increase in out-of-plane vibrations. This electron Table 8 Threshold energies E d (in eV) calculated for various processes in graphene. The system considered (periodic boundary conditions (PBC), flake or nanoribbon), number of atoms in the computational cell and calculation method are indicated  65,159,264 The threshold energies given for atom emission only correspond to the case when the transmission of energy from beam electron to graphene atom occurs perfectly perpendicular to the graphene plane. This is the angle at which a minimum energy will be required to emit the atom, while more energy will be needed if transferred at a shallower angle. When energy is transferred at non-perpendicular angles, there will also be a variation in threshold energy with the azimuthal angle j; the in-plane direction of energy transfer. In the typical setup of the electron beam being perpendicular to the graphene sheet, this minimum threshold energy coincides with the angle at which most energy can be transferred (y = p, where y is the electron scattering angle). This means that for emission processes, the single value of the threshold energy is highly informative, despite the anisotropic dependence of threshold energy on angles of momentum transfer. For other processes, such as the SW rearrangement discussed later, this overlap of the minimum threshold energy E d with the angle at which the maximum energy can be transferred from the beam does not hold, and so mapping of the anisotropy is vital. Mapping of the anisotropy of the threshold energy of emission from pristine graphene, shown in Fig. 15, 265 confirmed the very sharp increase in E d at y o p, rising from 23 eV to 43-780 eV at y ¼ p 2 , depending on j.

From vacancy structures.
Following the emission of a single atom from pristine graphene, a symmetrical vacancy or reconstructed 5-9 structure is formed. In any structure with an odd number of missing atoms, there will always be at least one under-coordinated atom bonded to fewer than three atoms. These atoms will be more easily lost via collisions with electrons from the e-beam due to a decreased emission threshold, calculated as 12 eV for the unreconstructed vacancy 266 and 14.7 eV in the case of the reconstructed monovacancy. 148 As this energy is easily accessible to even an 80 kV electron beam, correlated sputtering of two atoms from areas of pristine graphene is normally observed. Even-numbered vacancies are fully coordinated with three bonds to each atom and subsequently have larger threshold energies, although the value varies depending on the ring size and to a smaller extent on the specific local structure. For example, atoms in pentagons typically have threshold energies around 16.9 eV, 151 depending on the local atomic structure. This generally leads to the formation of separate divacancies on the graphene lattice, as once a divacancy is formed there is a much smaller preference as to which atom will be emitted next, and the typically low concentration of defect sites means that additional emissions are likely to occur on pristine graphene.
This effect is increased with the existence of the two commonly observed reconstructed divacancy structures shown in Fig. 3. As the central atoms of these structures have been shown to have a larger emission threshold than pristine graphene, 147 the probability of direct trivacancy formation is further reduced.
In dislocation dipoles, the emission of two atoms leads to a climb step, increasing the distance between the two dislocations (shown in Fig. 10). 218 Emission of atoms from the area surrounding the dislocation has been observed to occur more frequently than in pristine graphene, confirmed by calculations of the threshold energies for all surrounding atoms. 225 The threshold energies ranged from 18.64 to 24.14 eV, compared to 21.34 eV calculated for emission from pristine graphene using the same method. This study also demonstrated the importance of the beam direction in defects which introduce out-of-plane topology. As the buckled structure of the dislocation is not flat, there are two distinct directions in which it can be oriented with respect to the beam, up or down. Depending on the direction of the e-beam in the z-axis, the emission threshold energy of an individual atom varied by up to 2.9 eV (Fig. 16).
4.1.3 From graphene edges. A study of the emission threshold energies of various simple and reconstructed edge configurations used graphene nanoribbons of various widths. 148 DFTB calculations were undertaken for all edge configurations shown in Fig. 13. Termination of various edge configurations by hydrogen atoms was shown to have no effect on the emission threshold energies of edge atoms. The threshold energies calculated by this method are shown, ranging from B11 eV for the AC(56) edge to B22 eV for the AC(677) edge. Also shown are the results of more computationally demanding DFT simulations, carried out for the unreconstructed AC (19.0 eV), ZZ (12.0 eV) and KL (12.9 eV) edges. This set of calculations showed that the specific local structure of edges has a very large impact on the irradiation stability, and demonstrated the lack of correlation between irradiation stability and low formation energy; the AC and KL edges exhibit very flexible behaviour following impacts from the electron beam, which increases their resistance to atom emission. The discrepancy between these predictions of the stability of AC configurations over ZZ, compared to experimental observations of preferential ZZ formation in holes in graphene, 260 was attributed to the differences between an approximately circular hole and extended edges. The predicted stability of the KL edge under the e-beam was very recently confirmed by experimental HRTEM images showing their existence in graphene nanoribbons and at the edge of bulk graphene. 270 4.2 Bond rotations 4.2.1 In pristine graphene. Other than emission, the only major process of irradiation-induced atomic structure change in graphene is the Stone-Wales (SW) rearrangement, a 901 rotation of a single carbon-carbon bond. In pristine graphene a SW rearrangement leads to a SW defect, composed of two 5and two 7-membered rings. This has been observed in TEM and is seen to heal quickly, on the order of the electron dose per image used (B10 7 e nm À2 ). 147 Theoretical calculations confirmed the creation of the SW defect via impacts with electrons from the e-beam, revealing two mechanisms of formation. The energy required for SW defect formation was calculated as being typically 1 eV lower than the emission threshold (calculated as 22.5 eV), with 2 eV below this resulting in SW formation at certain space angles of transferred momentum.
The threshold energies for the formation and healing of a SW defect were calculated as being 22.5 eV and 14 eV respectively, confirming the reversible nature of the defect. 271 The same study revealed the ability of substrates to alter the energy required for each process; placing the graphene on hexagonal boron-nitride increased them to 31 eV and 17 eV respectively. An additional study 266 (using a much larger, converged supercell) gives comparable numbers for the pristine graphene case, although giving slightly lower values of 19 eV and 13 eV for formation and healing of the SW defect. The anisotropy of the energy required to heal the SW defect at space angles (y,j) was mapped, 266 revealing the strong dependence on the direction of momentum transfer, as previously discussed for atom emission. In addition, this paper revealed the ability of a SW defect to migrate due to a single electron impact, provided one of two combinations of space angles and energy (18 or 20 eV) are achieved.
4.2.2 Vacancy structure rearrangements. Divacancy structures have been seen to convert between different stable states under electron irradiation. 147 These transformations can be understood in terms of individual or multiple SW rearrangements on bonds in the vacancy structures. A very in-depth study 266 mapped the threshold energies required for singleimpact structure changes at different space angles for each of the V 2 (5-8-5), V 2 (555-777), V 2 (555-6-777) and V 2 (55-77) divacancy structures. It showed that the amount of energy required for a SW rearrangement is highly dependent on the angles of transferred momentum, ranging from 14 to 22 eV, and that different structures can result from impacts on the same atom, depending on the direction of momentum transfer. Fig. 17 shows a summary of these findings, giving the threshold energy for each SW rearrangement indicated by an arrow. It is worth noting the inclusion of the reconstructed V 2 (55-77) divacancy, as typically only the previous three structures are discussed. 99 This structure can be achieved via the rotation of a different bond in the V 2 (5-8-5) structure, however with a threshold energy of 21 eV compared to the 16 eV required for the V 2 (5-8-5) to V 2 (555-777) transformation. In addition, the reverse process of V 2 (55-77) to V 2 (5-8-5) has the lowest threshold energy of any divacancy rearrangement at 14 eV, and so while this structure should be able to form under irradiation, it would be expected to have a shorter lifetime compared to the other reconstructed divacancies. Indeed, while it is very rarely seen in TEM, an example of this structure forming from a V 2 (5-8-5) vacancy, before quickly undergoing the reverse transformation, has been observed. 167 This structure also provides a low energy route to the formation of trivacancies, with only 14 eV required to remove an atom and form the V 3 (5-10-5) structure.
Although the SW rearrangement has been shown to be responsible for a wide range of dynamical behaviour in larger vacancy structures, the need to map the anisotropy of the threshold energy, combined with an increasing number of inequivalent atoms and potential structures in larger vacancies, makes the calculation of threshold energies prohibitive. While several key structures of particular stability were found for the tetravacancy, 264 a large number of other structures were observed experimentally, each living for seconds at a time. Fig. 18 shows the 17 most frequently observed tetravacancy structures together with the SW rearrangements required to convert between them. Tight-binding MD simulations were carried out in order to demonstrate the evolution of the tetravacancy under electron irradiation in this manner, 264 however a full search of the impact angle space in order to determine threshold energies for defect structures larger than divacancies has not been performed.
4.2.3 Vacancy migration. In addition to conversion between different stable structures, vacancy structures have been observed to migrate across the graphene lattice. 167 As the barrier to the thermal migration of monovacancies is in the region of 1-1.5 eV (Table 1), it is assumed that the monovacancy migration is primarily thermally driven. Impacts from the electron beam will typically transfer much more energy than this (at 80 kV the maximum transferrable energy to a carbon atom is 15.8 eV), and so electron impacts on monovacancies are generally assumed to result in the emission of the under-coordinated carbon atom, as discussed above. However, the barriers for thermal migration of larger vacancy structures are unreasonably large, and so the effects of the e-beam are assumed to be responsible for the experimentally witnessed migration. Two consecutive SW rearrangements are proposed to be responsible for a migration of the V 2 (5-8-5) vacancy by one lattice unit via the V 2 (55-77) structure. 167 The theoretical study described above 266 provided confirmation of the possibility of this proposed two-step mechanism; it confirms the ability of SW rearrangements to convert between the two divacancy structures, and the low irradiation stability of the V 2 (55-77) structure (with the lowest threshold energy calculated, 14 eV) provides an explanation for both rearrangements happening in quick succession. An additional two step mechanism was proposed 167 for the V 2 (5-8-5) vacancy rotating 601 around a pentagon, via an intermediate V 2 (555-777) structure. The calculated threshold energies of these two steps, 16 eV and 19 eV, suggest that while this mechanism will be less likely to occur than the migration via the V 2 (55-77) structure, it will be possible given an electron beam of sufficient energy. During the mapping of the threshold energies of various divacancy structures, 266 a previously unpredicted single-impact mechanism of V 2 (5-8-5) migration was discovered, via bond breaking and reforming rather than through a direct bond rotation. The threshold energy of this process is 19 eV, and due to the requirement of only one electron impact it would be expected that this mechanism of migration would be especially prevalent at sufficiently high e-beam energies.
4.2.4 At graphene edges. The threshold energies for the bond rotation between the zigzag (ZZ) and reconstructed zigzag pentagonheptagon ZZ(57) edges in graphene were calculated from experimental TEM images of torn graphene. 255 These calculations assumed that every scattering event that transferred energy Fig. 17 A summary of the threshold energies calculated for irradiation-induced processes in ref. 266. Emission processes are labelled in red and bond rotations (SW rearrangements) in black. Rings are coloured according to their size: 5, yellow; 6, white; 7, blue; 8, purple; 9, green; 10, red. Threshold energies are also given for the direct migration of structures across the graphene lattice due to a single impact from the e-beam, only witnessed for a small number of combinations of angle space and transferred energy.
above the threshold energy would result in a transformation; they ignored the anisotropic dependence of the threshold energy. By observing the rates of conversion between the two structures, it was shown that the threshold for the ZZ to ZZ(57) rotation is 1.3 eV, while the ZZ(57) to ZZ rotation is 2.5 eV. These values are far lower than theoretical calculations of thresholds for SW rearrangements in the graphene interior, consistent with the smaller number of bonds that must be broken at the edge in order for a SW rearrangement to take place. It would be expected that in general, bond rotations will be far more prevalent at the edge than in the graphene interior.
In addition to conversions between the ZZ and ZZ(57) structures via a bond rotation, other SW rearrangements have been observed occurring during simulations of electron impacts below the emission threshold at various edge configurations. 148 Conversions between the AC and AC(57) edges and the AC and AC(667) edges were witnessed, while the AC to AC(56) transformation was always seen to occur following the emission of an atom. Threshold energies of these and other similar processes have not yet been calculated, but the knowledge of rates of conversion between structures could provide valuable insights into the dynamical behaviour of graphene edges under electron irradiation.

Conclusions
In the present paper we review the recent studies devoted to the energetics of formation, transformation, migration, coalescence and healing of defects in graphene. The cases of various vacancy-type and topological defects created due to bond realignment reactions in the graphene interior and at the graphene edge are considered, with both thermally activated and irradiationinduced processes taken into account.
The structure and formation energies for the monovacancy and many types of multi-vacancies in the graphene interior have been studied in detail. It has also been found that adatoms and single vacancies have low barriers for migration, about 0.5 eV and 1.0-1.4 eV, respectively, whereas barriers for divacancy migration exceed 5 eV. Thermally activated processes related to multi-vacancy formation and healing at room temperature will therefore be determined by the migration of adatoms and single vacancies. However, contradictory results have been obtained for the values of out-of-plane displacement, the barrier for transition between equivalent states of the 5/9 vacancy, and influence of interlayer interaction on the barrier for 5/9 vacancy motion in graphite. DFT calculations demonstrate a drastic change in the structure and energetics of a single vacancy near the graphene edge, and predict the fast directional motion of the vacancy to the edge at room temperature. 67,68 However these results have only been obtained on a few examples of rather small graphene flakes and narrow nanoribbons. Further detailed studies of the energetics of the vacancy near armchair and zigzag edges of semi-infinite graphene, both with dangling bonds and terminated by hydrogen and other atoms are necessary. A set of multi-vacancies with various structures formed due to vacancy coalescence have been found for free standing graphene. As strain has a drastic effect on the barrier to vacancy migration, 23,161 local strain could be used to control vacancy coalescence. Studies investigating the possibility of such control remain an open problem.
First-principles calculations clearly demonstrate that bond realignment reactions in the graphene interior have large barriers and should not take place even at high temperatures. The effect of strain on these reactions has been extensively studied in recent years, shedding light on the mechanisms of plastic deformations of graphene and carbon nanotubes. The catalytic effects of metal, carbon and hydrogen adatoms have been considered. Similar to vacancies, topological defects have been shown to result in out of plane buckling of graphene layers thus providing stress relaxation. However, the magnitude and energetics of the out of plane distortions are very sensitive to the calculation method used and the size of the periodic cell (or a finite system) considered. Therefore, more accurate studies are still required to get reliable data for the out of plane deformations caused by defects in graphene.
Graphene edges are found to be much more reactive than the graphene interior, and in the absence of strain or catalysing agents high-temperature transformations in graphene are concentrated at the edges. Significant efforts have been made to describe the reconstruction of zigzag edges. However, other types of bond realignment reactions at the edge have mostly been studied using semiempirical potentials, which are only qualitatively correct. Moreover, the investigation of sequential steps of defect generation has been focused on free (unfunctionalised) edges. More detailed first-principles analysis of various bond realignment reactions and the inclusion of functionalised edges into consideration are therefore of interest. In particular, the formation of carbon chains at edges has been shown to play a leading role in high-temperature transformations of graphene flakes 70,71,87 and open-ended nanotubes, 88 and has been observed experimentally for graphene nanoribbons under electron irradiation. 13 More detailed investigation into chain formation reactions is required to understand these processes.
Quantitative studies of irradiation-induced processes of structure change in graphene overwhelmingly deal with the case of atom emission, predominantly from pristine graphene. The threshold energy for this process is very well known, with good agreement between experiment 146,267 and theoretical studies ( Table 8), but similar studies for interesting defective types of graphene structure are lacking. While simulations considering edges 148 and small vacancy structures (up to divacancies) 266 have produced useful insights into the dynamics of these areas under electron irradiation, these remain the minority of cases. Even among the various divacancy structures, the V 2 (55-77) dislocation dipole typically receives little attention, presumably due to it being a much rarer sight experimentally. However, the presence of low energy routes from this structure mean that it could be playing an important part in irradiationinduced vacancy dynamics. A threshold energy of 14 eV for the rearrangement back to the V 2 (5-8-5) structure provides an accessible mechanism of divacancy migration, and a threshold energy of only 14 eV for the emission of an atom provides a low energy pathway to trivacancies and therefore to the creation of larger vacancy structures. Indeed, the low number of experimental observations of this structure may merely indicate that it typically has a short lifetime under the e-beam due to the presence of these low threshold processes.
Larger vacancy structures and commonly observed structures remain to be examined at all in terms of threshold energies of irradiation-induced processes. An especially important example of this is the stability and migration of grain boundaries in graphene. Emission threshold energies of different grain boundary configurations could help to explain experimentally seen structures, while the HRTEM observation of the atom-by-atom migration of a grain boundary has been reported, 20 mediated by individual bond rotations.
Bond rotations have been shown to be responsible for a wide range of atomic structure transformations in graphene, providing accessible routes between different stable structures. However, the theoretical estimation of threshold energies for this kind of process remains difficult due to the requirement of mapping the energy over the range of impact angles. This is reflected in the severe lack of data available for SW rearrangements in non-pristine graphene, with only several key studies providing quantitative understanding of these processes. For example, despite knowledge of the irradiation stabilities of a variety of edge configurations with respect to emission, 148 one of the very few examples of a quantitative study of bond rotations at the edge only deals with the conversion between a single pair of edge configurations, 255 and the threshold energies are calculated from the experimentally observed rates of conversion.
As this makes the large assumption that any impact above the threshold energy would result in a bond rotation, theoretical calculations are necessary for comparison; mapping the anisotropy of this case as has been performed for divacancies 266 could provide accurate values of the threshold energy. Simulations have shown that similar conversions can occur between a variety of other edge reconstructions, 148 however even the relative threshold energies of these processes, and therefore the relative likelihood of their occurrence under the e-beam, are unknown.
Knowledge of threshold energies of irradiation-induced processes are key to understanding the dynamical behaviour of graphene under the e-beam, with multiple examples showing that this behaviour cannot be explained by simply considering the equilibrium energetics. Despite this, and with HRTEM increasingly becoming the de facto tool for the experimental study of the atomic structure of graphene, theoretical calculations of threshold energies and resulting cross-sections remain scarce compared to the wealth of data on thermal activation and formation energies.
In summary, there are several areas in which further calculations of the energetic characteristics of atomic scale structure changes in graphene are necessary: (1) the presence and magnitude of out of plane displacements in 5/9 monovacancies, SW defects, dislocation cores and other defects, (2) the barrier for the transition between equivalent states of the 5/9 vacancy, (3) detailed energetics for the vacancy migration to graphene edges, (4) the control of vacancy migration and coalescence by local strain, (5) detailed energetics for bond realignment reactions at graphene edges, beyond the formation of pentagon-heptagon pairs at the zigzag edge, (6) detailed energetics of bond realignment reactions at hydrogenated and other types of functionalised edges, (7) detailed energetics of the formation of carbon chains at different types of edges, (8) threshold energies of structure changes (emission and migration by bond rotation) at grain boundaries, (9) both qualitative (relative rates) and quantitative (threshold energies) details of conversions between different edge configurations, (10) threshold energies of processes involving vacancy structures larger than divacancies.