# A Large-Scale Molecular Dynamics Study of the Divacancy Defect in Graphene

Jean-Marc Leyssale, Gérard Vignoles

► To cite this version:

Jean-Marc Leyssale, Gérard Vignoles. A Large-Scale Molecular Dynamics Study of the Divacancy Defect in Graphene. Journal of Physical Chemistry C, 2014, 118, pp.8200-8216. 10.1021/jp501028n. hal-01761328

HAL Id: hal-01761328

https://hal.science/hal-01761328

Submitted on 8 Apr 2018

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
A large scale molecular dynamics study of the single
and di-vacancy defects in graphene

Jean-Marc Leyssale*,† and Gérard L. Vignoles‡

CNRS, Laboratoire des Composites Thermostructuraux, UMR 5801:

CNRS-Herakles-CEA-Université de Bordeaux, 3 allée de la Boetie, 33600 Pessac, France, and
Université de Bordeaux, Laboratoire des Composites Thermostructuraux, UMR 5801:

CNRS-Herakles-CEA-Université de Bordeaux, 3 allée de la Boetie, 33600 Pessac, France

E-mail: leyssale@lcts.u-bordeaux1.fr

<footnote>*To whom correspondence should be addressed</footnote>

<footnote>†CNRS, Laboratoire des Composites Thermostructuraux, UMR 5801: CNRS-Herakles-CEA-Université de Bordeaux, 3 allée de la Boetie, 33600 Pessac, France</footnote>

<footnote>‡Université de Bordeaux, Laboratoire des Composites Thermostructuraux, UMR 5801: CNRS-Herakles-CEA-Université de Bordeaux, 3 allée de la Boetie, 33600 Pessac, France</footnote>
**Introduction**

The topic of point defects formation and migration in graphite has been a long subject of debate
in carbon science¹ as, for instance, it is of key importance in understanding radiation damages
in nuclear graphites² or structure-property relationships in carbon/carbon composites, especially
when exposed to high temperatures.³ Even if their concentration may be very low in single crystals
under equilibrium conditions, a particular attention has been given to single vacancies (SV) due to
their possibly higher concentrations in more disordered or irradiated materials.⁴ In 1978, reviewing
earlier experimental works, Thrower and Mayer have proposed values of 7.0 ± 0.5 and 3.1 ± 0.2
eV for respectively the formation and the migration (in the basal plane) energies of SVs.¹ Later
on, and before its actual experimental realization, electronic structure calculations (ESC) were
performed on single graphene sheets used as model graphite.⁵,⁶ Kaxiras and Pandey obtained a
formation energy of 7.6 eV for a bare SV⁵ and a migration energy *E*<sub>a</sub> of 1.6 eV, involving a saddle
point with a fourfold atom in the center of the vacancy (55-66 configuration).⁵ These results were
revisited by El-Barbary *et al*. showing that SVs reconstruct by closure of two of the three dangling
bonds to form a pentagon (noted 5-db), lowering the formation energy by 0.2 eV.⁶ The migration
energy reported in this work was of 1.7 eV with a similar saddle point similar yet slightly less
symmetric. The calculated migration energies (1.6-1.7 eV), significantly lower than previously
hypothesized by Thrower and Mayer indicates that self-diffusion of defects in graphite cannot be
accounted for by SVs only (more recent calculations indicates values of around 1.1⁷ and 1.4⁸,⁹ eV
for the vacancy migration in respectively graphite and graphene). In addition, these authors also
studied the case of divacancies (DVs), finding a reconstruction in a fully sp² defect made of two
pentagons and a central octogon (585) by closure of the four dangling bonds and with a formation
energy of 8.7 eV. The migration energy of this defect was estimated to be of around 7 eV⁶ making
DVs almost immobile at low temperatures.

The facile synthesis of single graphene sheets¹⁰ together with the development of atomic resolution microscopy techniques like STM or abberation corrected TEM (ac-TEM) have opened a new era in the characterization of point defects in carbon materials.¹¹⁻¹⁶ Using ac-TEM Hashimoto *et al.* were the first to create and image *in situ* some atomic defects in graphene.¹¹ More precisely they were able to detect a bare SV, a bare DV, a reconstructed DV (585) and an adatom on a very tiny surface of a few tens of nm² only. In the following years many kinds of defects were imaged ranging from point defects like SVs, DVs and other topologies composed of 4 to 8-member rings,¹²⁻¹⁷ to delocalized and pseudo-amorphous defect areas.¹²,¹³ These experimental results have considerably renewed the interest in DVs defects by showing that, in addition to being abundant, DVs can show a significant mobility under an electron beam. DFT calculations have shown that mobility is achieved through series of Stone-Wales like bond rotations, involving several defect topologies.¹³,¹⁴,¹⁸⁻²⁰

Using tight-binding molecular dynamics (TBMD) at elevated temperatures Lee *et al.*¹⁸ have shown that the 5-8-5 reconstruction (noted here as form A) spontaneously transform into a more stable form composed of three pentagon/heptagon pairs, the 555-777 (noted form B) defect. LDA DFT calculations presented in the same paper have shown that the latter has a formation energy 0.9 eV lower than the former. Later on, Kim *et al.*²⁰ have reported that the energy barrier (GGA DFT), of the rotation transforming form A into form B is of around 5.3 eV. In the same paper, an another reconstruction made of two pentagon/heptagon pairs, the 55-66-77 defect (form D) was considered as a stable intermediate for the translation or rotation of the 5-8-5 defect. Its formation energy is however much higher than the two formers (4.37 eV above form B). From the experimental viewpoint, form D was first observed in reduced graphene oxide by Gómez-Navarro *et al.*,²¹ then in graphene by kotakoski *et al.*¹³ Finally, another stable form, the 5555-6-7777 (form C) was recently observed in a few TEM experiments.¹³,¹⁴,¹⁹ It consists of four pentagons, four heptagons and a central hexagon rotated of 30 degrees with respect to the graphene lattice and is formed by a bond rotation from the 555-777 defect. Calculations by Cretu *et al.* have shown that
its formation energy lies in between those of forms A and B and the rotation barrier is around 6 eV.¹⁹ All these migration steps are certainly of great interest in understanding the behaviour of defects in graphene but the picture is still probably quite incomplete.

Owing to the extremely high computational cost associated to ESC, these techniques have been limited so far to the study of small time and length scales (typically on the order of one or two nm and for a few ps). Molecular dynamics (MD) simulations with empirical interatomic potentials (EIP), although supposedly less accurate, allows bypassing these limitations. Such methods have been used in many studies regarding the structure and energetics of defect-containing graphene or nanotubes²²⁻³⁰ as well as their mechanical and thermal properties.³¹⁻⁴⁰ In most of these works, the Brenner reactive empirical bond order (REBO) potentials were employed. It includes the first⁴¹ (REBOI) and second⁴² (REBOII) generations as well as the adaptive intermolecular REBO potential (AIREBO) complementing REBOII with an adaptive treatment of van der Waals interactions and a torsional potential for single bonds.⁴³ In many of these studies, the results obtained with such force fields are found in excellent agreement with experiments or results from DFT or other ESC. Especially, in the case we are interested here of vacancy defects in graphene (or nanotubes), it was shown that these potentials give a rather fair account of the formation energies of the SW defect and of the vacancy defects, including SV and A, B and C reconstructions of the DV defect.²⁴,²⁸⁻³⁰ It was also shown that the second generation REBOs (including REBOII and AIREBO) predicts the 55-66 reconstruction to be the ground state of the SV defect instead of the 5-db defect.²⁸⁻³⁰ Although the REBOI potential is somehow more accurate for the SV defect, REBOII and AIREBO are considered as more advanced potentials and are generally used. They also benefit from being implemented in some widely-used simulation softwares like GULP⁴⁴ and LAMMPS⁴⁵ codes.

In this paper, we use MD simulations with the AIREBO potential to study the behaviour of single SV and DV defects. After a short methodological section, we first thoroughly evaluate the validity of the AIREBO potential to describe SV and DV defects in graphene by comparing their formation and migration energies to ESC values from the litterature and give an account of possible size effects. We then present results from MD simulations at elevated temperatures
allowing characterizing the equilibrium distribution of the different reconstructions of DVs and their dynamical properties. These simulations also show the occasional formation of many other reconstructions. In the last part of the paper we discuss the possibility of having some of these new structures observed in a real experiment.

# Systems and Methods

## General information

Calculations are performed on single graphene sheets containing a unique point defect, either SV or DV, as well as on pristine graphene and on a graphene sheet with a Stone-Wales defect,⁴⁶ for comparison. The adaptive interatomic reactive empirical bond order (AIREBO) potential of Stuart et al.⁴³ is used to describe carbon-carbon interactions in all the calculations reported in this paper. As in former work,²⁶ the Lennard-Jones component of the potential is discarded as we are only concerned with single isolated carbon sheets. For practical reasons we have been working with orthorhombic cells with periodic boundary conditions (PBCs). A large cell dimension (z axis) is taken perpendicular to the sheets (xy plane) to avoid artificial interactions due to PBCs. The number of unit cells in the x and y directions are chosen to obtain graphene sheets with geometries as close as possible to a square. The different kinds of calculations performed in this work are described in what follows.

## Static energy calculations

Geometry optimizations including cell relaxation have been performed with a simulated annealing scheme using a hybrid Monte-Carlo (HMC) simulation procedure.²⁶ In these simulations, pure canonical HMC⁴⁷ moves (with a 70 % probability) were complemented by attempts to increase or decrease the cell parameters in x or y (30 % probability) with an acceptance probability derived to fix the diagonal elements of the stress tensor to zero.²⁶ Maximum volume changes at-
tempts as well as the lengths of the molecular dynamics runs in terms of step length and number of steps (in the range of 5 to 20 steps) are adjusted along the runs to get around 50 % acceptances. Starting from a high temperature the system is slowly quenched down to $10^{-3}$ K, a sufficient value for the convergence of energy, forces and geometries, by multiplying the temperature by 0.9 every $10^5$ steps. Most of the calculations are performed with an initial temperature of 2000 K apart for some unstable structures (bare SV od DV defects essentially) that would undergo a quick reconstruction at such a high temperature and for which the simulations were initiated at 500 K.

Formation energies were computed according to

$$E_f = E(N) - Ne_0 \quad (1)$$

where $E(N)$ is the energy of the optimized sheet of $N$ atoms and $e_0$ is the energy per atom in a perfect graphene sheet (optimized in the same way). Most of the static calculations reported in this paper have been performed on sheets based on a 3840-atom ($\approx 10 \times 10$ nm$^2$) graphene. In order to look for possible finite size effects when smaller sheets are considered, some calculations were also performed on 640-atom graphene ($\approx 4 \times 4$ nm$^2$). This size can be considered as an upper limit for ESC calculations are most of them, reported so far, deal with systems of around 100 atoms. For the same concern, calculations on the 5-8-5 reconstruction of the DV defect have been performed with the petagon-octagon-pentagon axis parallel to a simulation cell edge and rotated of 30° with respect to the latter as PBCs can induce interactions between the defect and its periodic images when small system sizes are used; these interactions, if any, should depend on the defect orientation.

Static energy calculations were also used to characterize the rotation barriers connecting different forms of reconstructed divacancies as well as the Stone-Wales defect and graphene. These barriers were computed by series of simulated annealing runs at increasing (or decreasing) values, by step of 2.5°, of the orientation in the $xy$ plane of a rotating dimer. In each of these simulations the dimer orientation was held fixed in plane and move attempts were performed to increase or
decrease its length (15 %), to move the other N-2 atoms using HMC (70%) or to change the cell
dimensions (15 %). In these runs the simulations were initiated at 10 K, the temperature was de-
creased by a factor 0.75 every 2000 steps and stopped at 10<sup>-2</sup> K. This procedure is similar to the
one adopted by Kim et al. with DFT calculations.²⁰

**Molecular dynamics simulations**

Series of molecular dynamics simulations were performed to study both the mobility (diffusion) and the relative importance of the different possible forms of SV and DV defects, at different elevated temperatures. The general setup of the simulations is shown in Figure 1. In each case the simulation is initiated with a unique bare defect (either SV or DV) located in the center of a large, almost square, ideal graphene sheet (the lattice parameters of the sheet are obtained from an HMC simulation on a pristine graphene sheet at the suited temperature). In a circled area surrounding the defect, of radius close to half the square edge, atoms are allowed to move using unperturbed newtonian equations. On the square corners, on the contrary, some atoms are fixed to their lattice positions avoiding global translational motions of the sheet, while in the intermediate area a simple Andersen thermostat⁴⁸ is applied to fix the working temperature.

Details of the molecular dynamics simulations are as follows: initial velocities on moving
atoms are drawn from a Maxwell-Boltzmann distribution. Equations of motion are integretad using
a Velocity-Verlet algorithm⁴⁹ with a timestep of 0.2 fs. The collision frequency of the thermostat is
set to a low value, typically 1 ps⁻¹ per thermostated atom, in order to perturb as less as possible the
newtonian dynamics. Due to their different mobilities, SV and DV defects are respectively studied
from 2000 to 3000 K and from 3000 K to 4000 K. Different number of runs, run length and
system size are used for different defect type/temperature couples. These details are summarized
in Table 1. During most of these runs, configurations were stored every 2 ps for analysis. In some
cases and for shorter times, configurations were stored every 0.02 ps for analysis of short-time
scale dynamics.

The principal element of the analysis consists in determining the evolution of the defect location
<figure><img src="image_2.png" /><figcaption>Figure 1: General setup for the MD simulations; large spheres: free atoms; medium-sized spheres: thermostated atoms; tiny spheres: fixed atoms.</figcaption></figure>

<table><caption>Table 1: Computational details of MD simulations.</caption><thead><tr><th>Defect</th><th>Temperature</th><th># of runs</th><th>length (ns)</th><th># of atoms</th></tr></thead><tbody><tr><td>SV</td><td>2000 K</td><td>2</td><td>48</td><td>2719</td></tr><tr><td>SV</td><td>3000 K</td><td>1</td><td>20</td><td>2719</td></tr><tr><td>DV</td><td>3000 K</td><td>4</td><td>150</td><td>958</td></tr><tr><td>DV</td><td>3000 K</td><td>4</td><td>60</td><td>2718</td></tr><tr><td>DV</td><td>3250 K</td><td>5</td><td>90</td><td>2718</td></tr><tr><td>DV</td><td>3500 K</td><td>4</td><td>80</td><td>2718</td></tr><tr><td>DV</td><td>3750 K</td><td>3</td><td>40</td><td>2718</td></tr><tr><td>DV</td><td>4000 K</td><td>3</td><td>24</td><td>2718</td></tr></tbody></table>
and morphology with time. For SV defects this can simply be achieved by counting the number of non sp² or non hexagonal atoms. The vacancy location is then defined as the center of mass of groups of such atoms. For DV defects, such a programmable analysis is more difficult to achieve, especially at the highest temperatures where there can be some overlap, due to thermal vibrations, between first and second neighbour distances. To avoid such problems, DV are identified by visual inspection (using the VMD software⁵⁰) of every configuration and their location defined by the center of mass of the four closest atoms to the vacancy center.

From these analysed trajectories, equilibrium distributions of the different forms of SV and DV defects are simply obtained by averaging over the different trajectories performed at a given temperature. In order to characterize their mobility, the average mean square "in-plane" displacements (MSD) of the defects are computed according to:

$$MSD(t) = \left\langle \left(x_I(t-t_0) - x_I(t_0)\right)^2 + \left(y_I(t-t_0) - y_I(t_0)\right)^2 \right\rangle_{t_0} \quad (2)$$

where $x_I(t)$ and $y_I(t)$ are the positions of the vacancy at time $t$ on the x and y axes respectively and $\langle \dots \rangle_{t_0}$ indicates averaging on different time origins. In our calculations each configuration is used as a time origin and the MSD(t) are computed up to different times ranging from 500 ps to 20 ns, depending on the defect mobility. The defects "in-plane" self-diffusivities ($D_s$) are then computed from MSD(t) according to Einstein's equation:

$$D_s = \frac{1}{4} \lim_{t \to \infty} \frac{d(MSD(t))}{dt} \quad (3)$$

In these calculations only the second half part of the MSD(t) curves (the high $t$ linear part) is considered. It is also worth noting that prior to these equilibriums and dynamic property calculations, the initial parts of the trajectories (up to at least five mobility events) are discarded.
# Results

## Energetics of vacancy defects in graphene

<figure><img src="image_3.png" /><figcaption>Figure 2: Optimized geometries of SW, SV (bare, 5-db and 55-66 configurations) and DV (bare, 5-db$_2$, 5-8-5, 555-777, 5555-6-7777 and 55-66-77 configurations) defects in a 3840-atom graphene sheet as computed with the AIREBO potential. (the color scale indicates the out-of-plane positions from -1.5Å(red) to 1.5Å(blue)).</figcaption></figure>

The geometries of several vacancy-defects, as optimized using the AIREBO potential in an ≈ 10× 10 nm² graphene sheet are shown in Figure 2. Are included the bare SV and DV defects and their known reconstructed forms: 5-db and 55-66 configurations for SV; 5-8-5 (A) 555-777 (B), 5555-6-7777 (C) and 55-66-77 (D) configurations for DV. The hypothetical partial reconstruction 5-db$_2$ of DV, the Stone-Wales (SW) defect and graphene are also given for comparison. Their formation energies with respect to graphene, as computed using Eq. (1) are compared to ESC
values from the literature, when available, in Table 2. In this table only values regarding pristine free standing graphene are reported as it was shown that the corresponding values in 3D graphite⁷ or in metal supported graphene⁵¹,⁵² can differ significantly from the latter.

<table><caption>**Table 2: Formation energies of SW, SV and DV defects in graphene. Results from simulated annealing calculations using the AIREBO potential on a 3840-atom graphene sheet are compared to ESC from the literature.**</caption><thead><tr><th>Vacancies</th><th>Name</th><th>Configuration</th><th>E<sup>AIREBO</sup><sub>f</sub> (eV)</th><th>E<sup>ESC</sup><sub>f</sub> (eV)</th></tr></thead><tbody><tr><td>-</td><td>SW</td><td>55-77</td><td>5.38</td><td>4.8-5.9<sup>53-56</sup></td></tr><tr><td>SV</td><td>bare</td><td>db<sub>3</sub></td><td>7.55</td><td>7.6<sup>5,6</sup></td></tr><tr><td>SV</td><td>reconstructed</td><td>5-db</td><td>7.02</td><td>7.4<sup>6</sup></td></tr><tr><td>SV</td><td>reconstructed</td><td>55-66</td><td>5.34</td><td>7.6<sup>6</sup></td></tr><tr><td>DV</td><td>bare</td><td>db<sub>4</sub></td><td>10.08</td><td></td></tr><tr><td>DV</td><td>1/2 reconstr.</td><td>5-db<sub>2</sub></td><td>9.47</td><td></td></tr><tr><td>DV</td><td>A</td><td>5-8-5 (C<sub>2v</sub>)</td><td>7.35</td><td>7.5,<sup>24</sup> 8.25,<sup>56</sup> 8.7<sup>6</sup></td></tr><tr><td>DV</td><td>A</td><td>5-8-5 (C<sub>2h</sub>)</td><td>7.70</td><td></td></tr><tr><td>DV</td><td>B</td><td>555-777</td><td>7.29</td><td>6.8,<sup>24</sup> 7.25-7.8<sup>57</sup></td></tr><tr><td>DV</td><td>C</td><td>5555-6-7777</td><td>7.87</td><td>7.5,<sup>24</sup> 7.65-8.2<sup>58</sup></td></tr><tr><td>DV</td><td>D</td><td>55-66-77</td><td>9.78</td><td>10.8-11.9<sup>20</sup></td></tr></tbody></table>

With the AIREBO potential, the energy per atom in a relaxed graphene sheet, $e_0$ in Eq. (1), is -7.408 eV, corresponding to a CC bond length of 1.398 Å. The formation energy of the Stone-Wales defect, 5.38 eV, is in excellent agreement with the values of 4.6-5.9 eV obtained with DFT calculations.⁵³⁻⁵⁵ The large range of values reported in the literature for the SW defect stems from the use of different different DFT formalisms and system sizes. Allowing (or not) out-of-plane deformation during the geometry optimization has also been shown to have an impact on the obtained energy.⁵⁵ The formation energy of a bare SV was used in the fitting process of REBO potentials so that the obtained value of 7.55 eV is almost equal to the one (7.6 eV) found with ESC.⁵,⁶ Also in agreement with DFT, the formation energy of the 5-db SV defect is slightly lower than the one of the bare SV (0.5 eV lower with AIREBO against 0.2 with DFT⁶). However, as observed by Krasheninnikov *et al.* for the REBOII potential,²⁹,³⁰ the AIREBO potential fails in reproducing the energetics of the 55-66 SV. Indeed, with $E_f=5.34$ eV, it is found as the most stable form of SV with AIREBO while it should be a transition state in between equivalent 5-db configurations according to DFT ($E_f=7.6$ eV).⁶
Looking now at DV defects, we can see that reconstructions A, B and C show $E_f$ values around 2-3 eV below the bare ($db_4$) defect while form D is almost iso-enegertic to the latter (we could not find any ESC value for $db_4$ and 5-$db_2$ defects). Formation energies obtained with AIREBO for the three commonly observed DV reconstructions (A, B and C) are found in the range of DFT values published in the literature (they are also very close to those obtained with the REBOI potential$^{24}$); $E_f$ of form D is slightly underestimated with AIREBO. In agreement with DFT, form B (555-777) is the most stable.

<figure><img src="image_4.png" /><figcaption>Figure 3: Optimized geometry of a 5-8-5 defect in a 3838-atom sheet; (a) stable ($C_{2v}$) form; (b) metastable ($C_{2h}$) form (same color scale as is Figure 2).</figcaption></figure>

As can be seen in Table 2, two $E_f$ values are reported for the 5-8-5 defect. Indeed, doing our simulated annealing calculation starting from a perfectly flat structure, the $C_{2h}$ geometry displayed in Figure 3b was obtained with $E_f$ = 7.70 eV. Repeating the calculation on a much smaller system (638 atom), a more buckled $C_{2v}$ form came out with $E_f$ = 7.26 eV; the $C_{2h}$ geometry being only recovered in this small system by initiating the calculation at a much lower temperature (10 K). The formation energy of the latter (7.68 eV) was close, though slightly lower than the one of the corresponding form for the large system. By initiating the calculation on the large system with a curvature compatible with the $C_{2v}$ geometry we were able to obtain the structure shown in
Figure 3a, which with $E_f = 7.36$ eV is the AIREBO groundstate for the 5-8-5 defect.

<figure><img src="image_5.png" /><figcaption>Figure 4: Root mean square deviation (RMSD) of the CC bond length with respect to pristine graphene as a function of the radial distance to the 5-8-5 defect center (circles: $C_{2v}$; squares: $C_{2h}$).</figcaption></figure>

The evolutions of the deviation of the CC bondlength with respect to its value in pristine graphene as a function of the distance from the defect center are shown in Figure 4 for the two geometries of Figure 3. It shows that buckling in the $C_{2v}$ structure allows reducing the strain on the CC bond length, thus reducing elastic energy. Considering elastic energy also explains the slight increase of $E_f$ with the system size for a given geometry: the smaller the size, the smaller the number of strained bonds and so the lower the elastic energy (more data regarding size effects are given in the supporting information (SI)). Finally, in Figure 4 we can see that the RMSD of the CC bond length reaches a plateau at distances from the defect core above 3 nm, justifying the use of 8 × 8 nm² sheets in the MD simulations. Very little is said in the literature regarding buckling in the 5-8-5 defect. El-Barbary *et al.*⁶ mention a flat defect yet the formation energy reported by these authors (8.7 eV) is the highest among reported values for this defect. In agreement with former studies,¹⁸⁻²⁰,²⁴ defects B and C were found to have rather flat geometries for both system sizes (638 and 3838 atoms) with little size effects on the energies as well ($E_f$ on 638-atom systems were found 0.05 and 0.14 eV lower than in 3838-atom systems for respectively B and C reconstructions). As for form A, form D is a highly buckled defect. Our calculations have converged to a $S_i$ symmetry structure for both system sizes with the smaller one stabilized by 0.3 eV with
respect to the larger one.

<table><caption>**Table 3:** Migration energies $E_m$ of DV through bond rotations (the case of SW defect formation in graphene is given for comparison).</caption><thead><tr><th>α</th><th>β</th><th>E<sup>AIREBO</sup><sub>m,αβ</sub> (eV)</th><th>E<sup>AIREBO</sup><sub>m,βα</sub> (eV)</th><th>E<sup>DFT</sup><sub>m,αβ</sub> (eV)</th><th>E<sup>DFT</sup><sub>m,βα</sub> (eV)</th></tr></thead><tbody><tr><td>graphene</td><td>Stone-Wales</td><td>9.4</td><td>4.0</td><td>9.2<sup>53</sup> 10<sup>59</sup> 10.3-11<sup>54</sup></td><td>4.4<sup>53</sup></td></tr><tr><td>5-8-5 (A:C<sub>2h</sub>)</td><td>555-777 (B)</td><td>5.8</td><td>6.2</td><td>5.3<sup>20</sup></td><td>6.2<sup>20</sup></td></tr><tr><td>5-8-5 (A:C<sub>2h</sub>)</td><td>55-66-77 (D)</td><td>5.0</td><td>2.6</td><td>5.3<sup>20</sup></td><td>1.8<sup>20</sup></td></tr><tr><td>555-777 (B)</td><td>5555-6-7777 (C)</td><td>6.6</td><td>6.0</td><td>6.0<sup>19</sup></td><td>5.7<sup>19</sup></td></tr></tbody></table>

In Table 3 we compare the energy barriers of the SW transformations connecting different reconstructions to DFT data from the literature. The $C_{2h}$ form of the 5-8-5 defect was chosen as it is the one with closest $E_f$ to DFT calculations. Again, AIREBO performs almost equally with DFT, especially when looking at the large range of values obtained by different authors for the SW defect.<sup>53,54,59</sup> Note that because we fix the rotating dimer in the xy plane in our barrier calculations the difference in the energy barriers from reactive to product and from product to reactive is not strictly equal to the difference in the formation energies reported in Table 2. This is especially the case for curved structures like forms A and D for which $E_f$ can significantly be affected (a few tenths of eV) by such a constraint.

## Molecular Dynamics

Having shown in the previous section that the AIREBO potential can fairly well reproduce the energetics of vacancy defects in graphene, we present in what follows the dynamical behaviour of isolated SV and DV defects in graphene simulated at high temperatures.

## Description of trajectories

Figure 5 shows some snapshots characterizing the early steps of MD simulations initiated with respectively a bare SV and a bare DV defects. Because of their different migration energies, SV and DV were simulated at respectively 2000 and 3000 K. At 2000 K, the SV defect (Figure 5a) reconstructs in the 5-db defect (Figure 5b) after 56 ps, then quickly evolves at t = 64 ps towards
<figure><img src="image_6.png" /><figcaption>Figure 5: Initial steps of MD simulations of a SV defect at 2000 K (a-d) and a DV defect at 3000 K (e-h) in 8 × 8 nm² graphene sheets.</figcaption></figure>

the 55-66 form (Figure 5c), the groundstate for the SV defect with the AIREBO potential (see Table 2). This configuration survives for another 436 ps before the defect migrates to an equivalent configuration at t = 500 ps (Figure 5d) with its center of mass translated to the red atom instead of the blue one (see also the destruction of the bond between blue and green atoms and the creation of a bond between red and green atoms). At 3000 K, a bare DV (Figure 5e) evolves to a 5-8-5 reconstruction after 228 ps (Figure 5f) by the successive closure of the four dangling bonds (an intermediate 5-db₂ structure, not shown, is observed between t = 224 and t = 228 ps). It then remains unchanged for more than 3.5 ns before the 555-777 (Figure 5g) reconstruction gets formed at t=3746 ps, survive for ≈ 2.7 ns and transforms into the 5555-6-7777 (Figure 5h) defect. This succession of 5-8-5, 555-777 and 5555-6-7777 configurations was observed in a recent TEM experiment by Kotakoski et al. (see Fig. 2 in Ref. 13). The two sequences of configurations displayed in Figure 5 are typical of what we have observed in our simulations: SVs acquire mobility through translation of the 55-66 defect, DVs through series of bond rotations involving different fully reconstructed forms. In all these simulations the transient times (time required to actually perform the transition between two configurations) can be estimated to be of the order of 2 ps for DV migration and even lower for SV migration. No dependence on the temperature could
be estimated for these transient times. Another worth mentioning point is that at 3000 K, lattice
vibrations are such that no distinction can be made anymore between the $C_{2v}$ and $C_{2h}$ forms of the
5-8-5 defect.

Equilibrium distributions and dynamical properties

<figure><img src="image_7.png" /><figcaption>Figure 6: Arrhenius plot of the in-plane self-diffusivity coefficients with respect to temperature for SV (small filled squares) and DV (large empty squares; plain line: linear fit for DV; dashed line: linear fit for DV excluding the 3250 K data; dotted line: linear fit for SV.</figcaption></figure>

We comment now on average properties computed from the different trajectories performed
in this work (see Table 1). The self-diffusivity coefficients $D_s$, computed from the simulations,
are plotted with respect to the temperature in Figure 6. Data collected at 3000 K for DV defects
were not sufficient to obtain a reasonable estimate of $D_s$ at this temperature. Linear fit to the data
were computed according to the Arrhenius equation $\ln(D_s) = \ln(k_0) + E_a/k_B T$ where $k_0$ is the
pre-exponential factor, $E_a$ the activation energy, $k_B$ the Boltzmann constant and $T$ the temperature.
While the AIREBO potential poorly describes the geometry of the SV defect, it behaves much
better for its dynamics. Indeed, the activation energy $E_a$ of 1.7 eV is very close to the migration
barriers of 1.2-1.7 eV obtained by DFT calculations.⁶⁻⁹ As can be seen in Figure 6, and as expected
from both experiments and calculations, the mobility of the DV defect is much lower - the $D_s$ value for the DV at 4000 K being even slightly lower than the value for SV at 2000 K - and the error bars larger. Fitting the complete set of data we found an $E_a$ value of 4.55 eV which is significantly lower than what is expected from the calculated migration barrier (around 6 eV). This difference is possibly due to large error bars on the $D_s$ values, especially at the lower temperatures. For instance repeating the fit while excluding the 3000 K (the less accurate) data leads to a more satisfactory value of 5.17 eV. In all cases the kinetic prefactors $k_0$ are in satisfactory agreement (i. e. within two orders of magnitude) with the commonly assumed values derived from the Debye frequency of Graphite.⁷

<figure><img src="image_8.png" /><figcaption>Figure 7: Distribution of the relative probabilities of the three main DV reconstructions: A (circles), B (squares), C (lozenges). Data computed from MD simulations performed on a 8 × 8 nm² sheet for T ≥ 3250 K and a 5 × 5 nm² sheet for T = 3000 K. The inset shows $\ln(P_β/P_α)$ vs. $1/T$ (symbols) and $(H_α - H_β)/(k_B T) + \text{cst}$ plots for different (α,β) couples; circles & full lines (A,B); squares & dashed lines (C,B); triangles & dotted lines: (C,A).</figcaption></figure>

In Figure 7 we plot the relative fractions of time spent by a DV defect in respectively the 5-8-5 (A), the 555-777 (B) and the 5555-6-7777 (C) reconstructions. For simplicity the bare DV and 5-db₂ defects as well as a few other structures with dangling bonds obtained by the breaking of
a bond in the 5-8-5 form were counted as part of the latter (they amount for around 2 % of the time spent in 5-8-5 form at 3000 K up to 15 % at 4000 K, where they are certainly stabilized by entropic effects). In agreement with the formation energies reported in Table 2, the 555-777 defect is the most encountered configuration with around 60 % of the simulated time spent in this form. The 5-8-5 and 5555-6-7777 forms both account for around 20 % of the time. This, however, is in apparent contradiction with the formation energies: form A should have a significantly larger weight than form C. To further investigate this point we have computed the relative enthalpies of the three forms (assuming the pressure is zero) by averaging the potential energy obtained for each forms over the 60 ns of simulations performed at 3000 K on the large system. We found out that at this temperature forms A and C have very close enthalpies of respectively 0.55 ± 0.23 eV and 0.57 ± 0.45 eV above the one of form B. According to basic statistical mechanics, the equilibrium respective probabilities $P_\alpha$ and $P_\beta$ of being in state (configuration) $\alpha$ or $\beta$ should follow $\ln(P_\beta/P_\alpha) = \frac{H_\alpha-H_\beta}{k_B T} - \frac{S_\alpha-S_\beta}{k_B}$ where $H_i$ and $S_i$ are respectively the enthalpy and entropy of form $i$. Measured $\ln(P_\beta/P_\alpha)$ data are shown in the inset of Figure 7 for AB, CB and CA couples together with the corresponding predictions by considering the above mentioned enthalpies (the entropic terms being adjusted). The agreement is very good, keeping in mind the width of the error bars. This Figure also shows that these probabilities do not evolve much in the temperature range considered here.

The evolution with temperature of the average residence times (or lifetimes) $\langle \tau \rangle$ where $\tau$ is the time lap between the apprision of the defect and its annihilation are given in Figure 8 for reconstructions A, B and C. These data quantifies two facts that were clearly observed in our trajectories. First, the lifetimes of the different forms diminish of around two order of magnitudes when the temperature is increased from 3000 to 4000 K. Second, although having similar equilibrium distributions, forms A and C show completely different lifetimes: $\langle \tau \rangle$ for the latter being one order of magnitude larger than for the former. In other words, form A is often encountered yet for short times while form C is rarely encountered but survives for longer times when it is formed. Lifetimes values for form B have intermediate values between these two. Remembering that the
<figure><img src="image_9.png" /><figcaption>Figure 8: Average residence time with respect to temperature for the three main reconstructions of the DV defect. Error bars for form C at T = 3000 K and 3250 K were not computed due to insufficient data. The escape rates (1/τ) are plotted with respect to 1/T in the inset together with their Arrhenius fit.</figcaption></figure>
lifetime $\tau$ of a given form is actually the inverse of the escape rate from this form, we show in the
inset of Figure 8 an Arrhenius plot of the escape rates from these three defect forms.

<figure><img src="image_10.png" /><figcaption>Figure 9: Schematic of the potential energy surface (left) compared to the apparent enthalpy surface (right) deduced from MD simulations (side numbers correspond to the barrier heights in eV).</figcaption></figure>

From the data collected in the MD simulations we can thus propose an apparent enthalpy
surface (AES) for the DV defect that we compare in Figure 9 to the calculated potential energy
surface (PES). To build this diagram we have considered the enthalpies computed at 3000 K for
the three forms. The barriers from one form to another were based on the activation energies E_a
obtained in fitting the escape rates (see the inset in Figure 8). Indeed if the previously described
bond rotation schemes are the dominant migration mechanism, the activation energy obtained for
the disparition of forms A and C can be identified in both cases as the enthalpy barrier to form B.
Similarly the E_a value for the disparation of form B can be attributed to the formation of form
A, as it should have a higher rate than the formation of form C. As can be seen in Figure 9 the AES
is in qualitative agreement with the PES, however, the barriers are underestimated of around 1 eV
with respect to the PES, which is certainly too large a difference to be attributed to the thermal
energy only. We will come back to that point later on, in the discussion section.

In addition to forms A, B, and C the MD trajectories gathered in this work have also revealed
the occasional apparition of the already described form D and of not less than 17 additional fully
reconstructed forms of DV. In what follows we briefly enumerates these structures in a first place and highlight some observed reactive mechanisms in a second place.

## 17 more full reconstructions of the DV defect

Figure 10 shows snapshots of DV reconstructions E-U after geometry optimization. Their formation energies with respect to their number of atoms in non-hexagonal rings, as a measure of the defect size, are plotted in Figure 11 (numerical values are given in the SI). As can be seen in Figure 10, DV reconstructions can contain up to 42 atoms (form U) in non-hexagonal rings, against 14 for form A. Apart from form Q ($E_f$ = 15.06 eV), all their formation energies are in the range [7:13] eV. If globally the formation energy tends to increase with the defect size, the data appear very scattered. Based on Figure 11 the reconstructions can be classified into two groups: a low energy group consisting of forms A, B, C and E, and a high energy group with all the other forms. All these structures but two (L and Q) can be transformed into at least another one by a single bond rotation. Numerical values of the rotation barriers (migration energies) for all possible SW transformations between these structures are given in the SI. Their average lifetimes $\langle\tau\rangle$ (including those shown in Figure 8 for forms A, B and C) are given in Table 4. Tables for the number of occurrences and the longest-lived event of each of these forms are also given in the SI. We now briefly enumerate these different defects.

Form D was obtained in simulations ran at 3250 K and higher. It is a particularly unstable defect with lifetimes of only 1-4 ps due to a very low barrier (2.6 eV) for the migration back to the 5-8-5 defect (form A). At the same temperature were also obtained forms I and J. Among these 17 structures, form I is the one with the most numerous occurrences and with non-negligible lifetimes (up to 86 ps at 3250 K) while form J was extremely rarely encountered and never survived more than a few ps. Form I was most of the time obtained from form B, form J from form I. Increasing the temperature to 3500 K we were able to observe forms F and G with lifetimes up to 180 ps for both structures. These two defects are particularly interesting as they are both obtained by a true SW rotation, transforming four hexagons into two pentagons and two adjacent heptagons,
<figure><img src="image_11.png" /><figcaption>Figure 10: Seventeen other full reconstructions of the DV defect obtained in MD simulations. The geometries were optimized using simulated annealing with the AIREBO potential in a 3838-atom sheet (same color code as in Figure 2).</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Figure 11: Formation energies (AIREBO) of full reconstructions of the DV defect in a 3838-atom sheet as a function of the number of atoms in non-hexagonal rings (Form Q, having a much higher energy, is not displayed).</figcaption></figure>

in the immediate neighborhood of a 5-8-5 (A) defect. The rotation barriers of 7.4 and 7.6 eV to respectively form defects F and G from the 5-8-5 defect, 2 eV lower than the barrier to SW defect formation in graphene, clearly show that the presence of strain close to the 5-8-5 defect strongly facilitates the bond rotation mechanism. F and G defects are only two of the many possibilities to extend the size of the 5-8-5 defect by a true SW transformation. In the SI we consider all the other possibilities involving rotations around the first and second neighbour bonds. Among 11 possible such structures, forms N and O were also obtained in MD simulations performed at 4000 K. Interestingly, forms F, G, N and O are those combining rather low formation energies (below 12 eV) and low rotation barriers (from 7.4 to 7.7 eV). In addition to be connected to form A, form F is also connected to form C, form G to form I. Three structures containing a four-member ring, P, R and S, were also observed at 3500 K although with very short lifetimes of 2-3 ps. According to our calculations, both static and dynamic, form P is a metastable intermediate for the rotation of form D, obtained from the latter by the 30° rotation and destruction of the central bond. This structure
possesses two fourfold atoms in an almost planar geometry. As for the 55-66 reconstruction of the SV defect, its appearance is certainly an artefact arising from the AIREBO potential. Forms R and S are obtained by a SW rotation of a bond shared by a pentagon and a hexagon in respectively the 555-777 and the 5-8-5 defects. Form M is another short-lived structure obtained at 3750 K from form A.

<table><caption>**Table 4:** Average lifetimes $\langle \tau \rangle$ (ns) of the different reconstructions of the DV defect as obtained from MD simulations on 3838-atom system (error bars are given in parentheses when sufficient data were available for their evaluation).</caption><thead><tr><th>Temperature (K)</th><th>3000</th><th>3250</th><th>3500</th><th>3750</th><th>4000</th></tr></thead><tbody><tr><td>A</td><td>4.94 (1.45)</td><td>1.82 (0.40)</td><td>0.694 (0.026)</td><td>0.255 (0.045)</td><td>0.108 (0.021)</td></tr><tr><td>B</td><td>21.75 (6.64)</td><td>4.40 (1.15)</td><td>1.81 (0.20)</td><td>0.568 (0.038)</td><td>0.222 (0.030)</td></tr><tr><td>C</td><td>52.66</td><td>13.68</td><td>3.12 (1.64)</td><td>1.24 (0.23)</td><td>0.364 (0.075)</td></tr><tr><td>D</td><td></td><td>0.005</td><td>0.004</td><td>0.004</td><td>0.003</td></tr><tr><td>E</td><td></td><td></td><td></td><td></td><td>0.304</td></tr><tr><td>F</td><td></td><td></td><td>0.186</td><td>0.010</td><td>0.024</td></tr><tr><td>G</td><td></td><td></td><td>0.089</td><td>0.026</td><td>0.110</td></tr><tr><td>H</td><td></td><td></td><td></td><td></td><td>0.088</td></tr><tr><td>I</td><td></td><td>0.047</td><td>0.043</td><td>0.011</td><td>0.014</td></tr><tr><td>J</td><td></td><td>0.008</td><td></td><td></td><td>0.002</td></tr><tr><td>K</td><td></td><td></td><td></td><td></td><td>0.030</td></tr><tr><td>L</td><td></td><td></td><td></td><td></td><td>0.040</td></tr><tr><td>M</td><td></td><td></td><td></td><td>0.002</td><td></td></tr><tr><td>N</td><td></td><td></td><td></td><td></td><td>0.002</td></tr><tr><td>O</td><td></td><td></td><td></td><td></td><td>0.009</td></tr><tr><td>P</td><td></td><td></td><td>0.002</td><td></td><td></td></tr><tr><td>Q</td><td></td><td></td><td></td><td></td><td>0.002</td></tr><tr><td>R</td><td></td><td></td><td>0.002</td><td></td><td></td></tr><tr><td>S</td><td></td><td></td><td>0.003</td><td></td><td>0.002</td></tr><tr><td>T</td><td></td><td></td><td></td><td></td><td>0.002</td></tr><tr><td>U</td><td></td><td></td><td></td><td></td><td>0.010</td></tr></tbody></table>

Rising the temperature to 4000 K, seven more structures were identified in addition to the already presented forms N and O. Form E, belongs to the low energy structures in Figure 11, despite its 38 atoms in non-hexagonal rings. It is formed from form C by exactly the same process by which form C is obtained from form B. As for form C, it is a very long-lived structure with $\langle \tau \rangle = 304$ ps at this temperature. Form H is also a rather long-lived defect with $\langle \tau \rangle = 88$ ps (to be compared with 108 ps for the 5-8-5 defect at this temperature). It can either be formed by a true
SW rotation from form B, or by a bond rotation from form E. As for forms E and H, forms K, L and U are three of the largest full reconstructions observed in this work (U even being the largest of all with 42 atoms in non-hexagonal rings). Their lifetimes range from 10 (U) to 40 (L) ps. Forms K and U were both obtained from form E: by a bond rotation for form K and by a true SW rotation for form U. The formation of form L involves two simultaneous bond rotations from form C (see below). Finally, forms Q and T are two short-lived four-membering containing structures. With a formation energy of 15.06 eV, form Q is by far the reconstruction of highest energy observed in this work. Form T, in addition to a four member-ring, also contain a nine-member ring. This defect is obtained by a bond rotation from A. An unsaturated variant of this form, with the bond connecting the 4- and 9-member rings broken, was often encountered in transient states between fully reconstructed structures (see below).

## Transition mechanisms

In what follows the transformation mechanisms are further investigated by looking at some sequences of configurations extracted from the MD trajectories. It will actually show that the Stone-Wales like bond rotation mechanism considered so far is far from being the only one involved in the dynamics of the DV defect. Five typical kinds of mechanisms are shown in Figure 12. In the sequence of Figure 12a, taken from a MD simulation ran at 3000 K, a 5-8-5 defect is transformed into a 555-777 defect in around 2 ps by the usual bond rotation mechanism. Unlike, Figure 12b illustrates the rotation of a 555-777 defect at 3000 K. As can be seen, this transition is obtained through a completely different mechanism involving the successive rotation of three atoms around the central atom. The third mechanism (Figure 12c) shows the rotation and translation of a 5-8-5 defect at 3250 K. It first involves the formation of the 55-66-77 defect (D) through the rotation of the red bond; the isomerization of the 55-66-77 defect through the rotation of the bond between a red and a green atoms; and eventually the rotation of the bond between the upper red and green atoms. In Figure 12d are shown a single atom (yellow) migration mechanism transforming form B into form J, immediately followed by a bond rotation (orange-yellow) mechanism transforming
<figure><img src="image_13.png" /><figcaption>Figure 12: Typical migration mechanisms of the DV defect observed in MD simulations. a: A → B transition with a bond rotation mechanism (3000 K); b: rotation of the 555-777 defect (B) with a concerted three-atom rotation mechanism (3000 K); c: migration-rotation of the 5-8-5 (A) defect (3250 K); d: B → J transition with an atom translation mechanism immediately followed by a J → I transition using a bond rotation (3250 K); e: C → L transition using two concerted bond rotations (4000 K). (colored atoms and bonds are guides to the eyes).</figcaption></figure>
form J into form I (also taken from an MD run at 3250 K). Finally in Figure 12e is shown a con-
certed rotation mechanism of two bonds, observed at 4000 K, transforming form C into form L in
less than 2ps.

<figure><img src="image_14.png" /><figcaption>Figure 13: Three reactive sequences of configurations taken from a 92 ps long part of a DV trajec-
tory obtained at 3500 K (the color code is redefined at the beginning of each sequence).</figcaption></figure>

In Figure 13 we show a sequence of snapshots characterizing 92 ps of MD performed at 3500 K. Starting from a 5555-6-7777 configuration (C), a 4-555-77777 (R) configuration is obtained after 4 ps by a three-atom rotation mechanism. As for the rotation of the 555-777 defect (see Figure 12b) it proceeds by the successive counter-clockwise rotation of the red, green and orange atoms around the central blue atom. This structure quickly evolves towards form G (t = 8 ps) by the rotation of the bond between the blue and orange atoms. Form G survives for around 30 ps until a long sequence of configurations with unsaturated bonds is observed. As can be seen, this sequence corresponds to the 180° rotation of the bond between a red atom and an orange atom; defect G is recovered at t = 52 ps with only these two atoms inverted. At t = 90 ps this configuration then evolves to form I (t = 92 ps) by the rotation of the red bond.

Finally, in Figure 14 we show some snapshots of a particularly disordered and reactive sequence
<figure><img src="image_15.png" /><figcaption>Figure 14: Snapshots of a MD simulation of a DV defect at 4000 K (the simulation time is indicated in the bottom right corner of the images; the color code is a guide to the eyes).</figcaption></figure>

taken from a MD simulation ran at 4000 K. Starting from a 5555-6-7777 (C) defect a structure containing four dangling bonds is obtained at t = 6 ps by i) the insertion of a blue atom in between the two orange atoms (t = 4 ps); and ii) the rotation of the bond between a green and an orange atoms (t = 6 ps). This structure is actually an open form of the 4-55-77-8 (S) defect obtained by creating two bonds between the four unsaturated atoms (compare to the structural relationship between the 5-8-5 and the bare DV defects). It survives for more than 10 ps until the bond shared by the remaining five- and seven-member rings breaks up (t = 18 ps). At this point starts a complex series of transformaitons involving either a four-member ring bonded to a twofold atom (see the snapshots at t = 18 or t = 26 ps) or a chain of two twofold atoms (t = 20 or t = 28 ps) which eventually evolves toward a T defect at t = 34 ps. Form T does not survive for more than 2 ps and is transformed to form S (t = 46 ps) through another complex series of disordered unsaturated configurations. Among them, one contains a planar fourfold atom (at t = 36 ps), certainly overstabilised by the AIREBO potential. More interestingly, the structure observed on the snapshots at t = 42 and 48 ps differs only from form S (t = 46 ps) by the breaking of the bond shared by the four- and eight-member rings in the latter, this time in close analogy with the 5-db₂ partial reconstruction of the DV defect. This structure eventually gets transformed into the 5-8-5 (A) defect at t = 54 ps.
Discussion

In this work we have used the AIREBO potential to study the energetics and dynamics of isolated
single and di-vavancies in graphene. In the first part of the paper we have shown that, in agreement
with former studies,²⁹,³⁰ potentials based on REBOII considerably overestimate the stability of
the 55-66 reconstruction of the SV defect. This is also certainly the case for other structures
possessing a fourfold carbon atom in a close to planar geometry. Such structures have been indeed
very occasionnaly encountered in this work (form P of the DV defect for instance). However our
dynamical results have shown that the migration kinetics of the SV defect is pretty well reproduced
by the AIREBO potential with an activation energy of 1.7 eV, very close to estimations based
on DFT. Furthermore, concerning the DV defect and its different reconstructed forms, on which
we have paid much of the attention in this paper, AIREBO performs very well. For the four
reconstructions identified experimentally (A, B, C and D), the formation energies computed in this
work with AIREBO are close to the values reported earlier from DFT calculations and so are the
migration barriers through the well-known bond rotation mechanism.

In the second part of the paper we have carefully analyzed the results of long MD trajectories
of a graphene sheet with a DV defect simulated at high temperatures to enhance the mobility of
the defect. Doing so we have been able to reproduce the behaviour observed in recent TEM exper-
iments showing a series of transformations (bond rotations) between forms A (5-8-5), B (555-777)
and C (5555-6-7777) leading to the diffusion of the defect’s location. Quantitative equilibrium and
dynamical descriptors were computed from the MD data including the self diffusion coefficients
of the defect as well as the statistical importance of the different reconstructions, their average
lifetimes and escape rates. The results have shown that a DV defect spend around 60 % of the time
in form B against around 20 % for both forms A and C, whatever the temperature (in the range
[3000K:4000K]). Also, it was found that form C is the one with the lowest occurrences but the
longest lifetimes while form A, the shortest-lived of these three structures is formed much more
often than the former. The activation energy for the diffusion of the defect (around 5 eV), calcu-
lated from the Arrhenius plot of self-diffusion coefficients has appeared slightly lower that what
would have been expected from the rotation barriers (≈ 6 eV according to Table 3). Similarly, the
apparent enthalpy surface (Figure 9) derived from the computed enthalpies (at 3000 K) and escape
rates seems to present slightly shorter barriers than those computed in the static calculations. This
discrepancy between static calculations and MD calculations could stem from several sources of
errors. First, the rotation barriers have been computed with a "circle" method, similar to the one
adopted by Kim *et al*,²⁰ which does not ensure that the minimum energy path is followed between
the reactive and the product. Furthermore, in our calculations the rotating dimer was constrained
in the main plane of the sheet: this can have a significant effect on the energy. Note that it is
possible that the energy barriers published by Kim *et al*. are slightly overestimated for the same
reason. The DFT reference value for the B → C transition, by Cretu *et al*.¹⁹ should be more ac-
curate however as these authors used a nudge elastic band approach. Another possible explanation
for the discrepancy reported here could simply be the large error bars on the data inferred from
the MD simulations, due to poor statistical convergence (see for instance the change of 0.6 eV
obtained on the activation energy derived from the *D*<sub>s</sub> values (Figure 6) by excluding the 3250 K
data from the fitting procedure). Nevertheless, at least two other possible explanations can also
be considered: the first one is related to the thermal contribution to the enthalpy and the appear-
ance of dangling-bond containing structures like the bare DV and the 5-db<sub>2</sub> structures from the
5-8-5 defect; the second one is the possibility of other mechanisms with similar or slightly lower
energy barriers. Looking closely at the simulations these two scenarios cannot be neglected. The
formation of the 5-db<sub>2</sub> defect, whose contribution to the time spent in form A (where the 5-8-5,
5-db<sub>2</sub>, and 5-db<sub>4</sub> defects are considered as form A) has been found to be up to 15 % at 4000 K,
could diminish the rotation barrier of the A → B transition. Also, the rotation-translation of form
A through the mechanism described in Figure 12c, involding form D as an intermediate, has a
barrier almost 1 eV lower than the direct rotation mechanism. Another important mechanism is
the three-atom rotation mechanism of form B diplayed in Figure 12b. In their interesting paper
Kim *et al*.²⁰ have considered such a three-atom rotation mechanism for the rotation of the 555-777
defect. They have obtained a rather large barrier - 1.7 eV higher than for the indirect way involving
the successive formation of form A, D, A (rotated) and finally B - for the simultaneous rotation of
the three atoms, without giving further details on the calculation methodology. It is possible that
considering the successive rotation of the three atoms, as is observed in the MD simulations would
reduce the overall energy barrier. Unfortunately we do not dispose of the adequate numerical tools
to investigate the energy barrier of such a complex mechanism, neither with AIREBO, nor with
DFT.

In addition to the four experimentally observed DV reconstructions A, B, C and D, our MD
simulations have revealed the occasional apparition of 17 new fully reconstructed forms. Apart
from form P, certainly artificialy stabilized by the AIREBO potential, and form Q, having a much
higher formation energy than the others, all these reconstructions have formation energies in the
range [7:13 eV], i. e. significantly lower than two isolated SVs. Also, this work has shown that,
apart from L, they can all be obtained from another one by a bond rotation mechanism for which
the activation energy is significantly lower than the one required to form a SW defect in graphene.
Indeed, SW defects have been occisionaly observed away from the DVs in our simulations, yet
with a rather low frequency compared to the apparition of the DV reconstructions, espacially if we
take into account the many possible sites for the formation of a SW defect in a 8 × 8 nm² sheet.
Therefore, most of the new DV reconstructions reported in this paper could certainly be formed
under the electron beam of a TEM. Among these structures forms E, F, G and I, are often formed
in the MD simulations (see the occurence table in the SI) and have relatively high lifetimes. Form
E is certainly the best candidate to bet on for an experimental observation. As already said, defects
B, C and E have closely related structures and seem to represent a lower energy limit for the data
cloud of Figure 11. Similarly to form C, form E is rarely encountered with only eight occurrences
at 4000 K over approximately 400 events observed during the 72 ns simulated. However, after
forms A, B and C, form E is the fourth most encountered configuration at 4000 K (with 3.5 % of
the time spent in this configuration) due to a long average lifetime of 316 ps, larger than those of
forms A and B (86 and 256 ps respectively) and only slightly lower than the one of form C (450
ps). Looking at the MD trajectories performed at 4000 K, the three longest-lived events were a
555-777 form (1200 ps), a 5555-6-7777 form (1182 ps) and a 55555-66-77777 - or form E - (1144 ps).

It is also interesting to notice that some of the reconstructions discovered in this work are extremely extended defects compared to the bare DV or to its 5-8-5 reconstruction. For instance the five largest of them (H, K, L, E and U) count between 34 and 42 atoms in non-hexagonal rings against 14 in the 5-8-5 defect and from two to four bond rotations are necessary to obtain them from the latter. As we have seen however, their formation energies are rather moderate. These extended point defects deriving from localized vacancy defects have been observed long ago in metals and termed "relaxions" by Nachtrieb and Handler.⁶⁰ In this paper we have shown that graphene is a particularly favorable material for the formation of such relaxions, despite its 2D nature which should in principle limit the possibility of stress relaxation through delocalization.

Another important result from this paper concerns the lowering of the barrier to a true Stone-Wales rotation, transforming four hexagons into two pentagon/heptagon pairs, in the immediate neighborhood of a reconstructed DV. In this way, forms F, G, N and O were obtained from form A, form H from form B and from U from form C. Around six years ago by Meyer *et al*.¹² published one of the first papers on the direct imaging of topological defects in graphene. The two defect structures made of three pentagons and three heptagons shown in Figures 3i and 3k in Ref. 12 are particularly interesting. The first one is made by a (true) SW defect formed in the immediate neighborhood of an inverse Stone-Wales⁶¹ defect. The second one, closely resembling defect I in the present work (they differ from each other by the translation of a pentagon/heptagon pair by one lattice parameter) is obtained by the rotation of a bond joining a pentagon of a SW defect to a hexagon.

To conclude this discussion we would like to add a few comments on structures counting four-member rings. Four-member rings have been often observed in recent TEM experiments, essentially under the form of reconstructed vacancies $V_n$ where the number of vacancies $n \ge 4$.¹³,¹⁶ A typical structure of this type is the 5-8-4-8-5 reconstruction of the quadruple vacancy ($V_4$) defect formed by the coalescence of two aligned neighbouring DVs. In this work we have shown that
structures containig 4-member rings can possibly be formed from DVs. However, all these struc-
tures (Q, R, S and T) are found in the high energy side of the data cloud in Figure 11. Also,
the MD simulations have shown that some of these structures (S and T particularly) are easily
prone to ring opening (see Figure 14) during the dynamics. In the SI section we show that with
the AIREBO potential unreconstructed forms of the 5-8-4-8-5, S and T defects are very close in
energy to the fully reconstructed variants. In the case of the 5-8-4-8-5 defect, the structure ob-
tained by opening a bond joining a 4-member ring and an octagonal ring has a lower formation
energy than the reconstructed defect. Similarly, the partial reconstruction of form S observed in
Figure 14 is stabilized of around 0.2 eV with respect to the full reconstruction. These results agree
with the usual underestimation of the fraction of 4-member rings obtained in atomistic models
of amorphous carbons⁶² or graphitic carbons⁶³,⁶⁴ prepared using the liquid quench MD method
and EIPs based on the REBOII parametrization with respect to quenchs performed with ab-initio
(DFT) MD. Therefore one can expect that the statistical importance of 4-member ring containing
structures is slightly underestimated in this paper.

Acknowledgement

Access to the computer resources of the Mesocentre de Calcul Intensif en Aquitaine (MCIA) is
gratefully acknowledged. This work was supported by the Agence Nationale de la Recherche
through the PyroMaN (Pyrocabon Matrices at the Nanoscale) project (Grant ANR 2010 BLAN
0929 01) and by the Aerospace Valley Foundation.

**Supporting Information Available**

Available as supporting information are: 1) tables of the formation energies (AIREBO) of all the structures studied in this work and of all the rotation barriers; 2) discussion on the size effects on the formation energy of the 5-8-5 reconstruction; 3) snapshots of the optimized structures not shown in the main text; 4) tables for the number of occurences and longest-lived events of each DV reconstruction; and 5) discussion on the stability of 4-member ring containing structures with the
AIREBO potential. This material is available free of charge via the Internet at http://pubs.acs.org/.
Notes and References

(1) Thrower, P. A.; Mayer, R. M. Point Defects and Self-Diffusion in Graphite. *Phys. Stat. Sol.*
*(a) **1978**, *47*, 11–37.

(2) Telling, R. H.; Heggie, M. I. Radiation defects in graphite. *Philos. Mag.* **2007**, *87*, 4797–4846.

(3) Savage, G. *Carbon-carbon composites*; Chapman and Hall: London, 1993.

(4) Hahn, J. R.; Kang, H. Vacancy and interstitial defects at graphite surfaces: Scanning tunneling microscopic study of the structure, electronic property, and yield for ion-induced defect creation. *Phys. Rev. B* **1999**, *60*, 6007–6017.

(5) Kaxiras, E.; Pandey, K. C. Energetics of Defects and Diffusion Mechanisms in Graphite.
*Phys. Rev. Lett.* **1988**, *61*, 2693.

(6) El-Barbary, A. A.; Telling, R. H.; Ewels, C. P.; Heggie, M. I.; Briddon, P. R. Structure and energetics of the vacancy in graphite. *Phys. Rev. B* **2003**, *68*, 144107.

(7) Latham, C. D.; Heggie, M. I.; Alatalo, M.; Öberg, S.; Briddon, P. R. The contribution made by lattice vacancies to the Wigner effect in radiation-damaged graphite. *J. Phys.: Condens. Matt.* **2013**, *25*, 135403.

(8) Krasheninnikov, A. V.; Lehtinen, P. O.; Foster, A. S.; Nieminen, R. M. Bending the rules: Contrasting vacancy energetics and migration in graphite and carbon nanotubes. *Chem. Phys. Lett.* **2006**, *418*, 132–136.

(9) Zhang, H.; Zhao, M.; Yang, X.; Xia, H.; Liu, X.; Xia, Y. Diffusion and coalescence of vacancies and interstitials in graphite: A first-principles study. *Diamond Relat. Mater.* **2010**, *19*, 1240–1244.

(10) Novoselov, K. S.; Geim, A. K.; Morozov, S. V.; Jiang, D.; Zhang, Y.; Dubonos, S. V.; Grigorieva, I. V.; Firsov, A. A. Electric Field Effect in Atomically Thin Carbon Films. *Science* **2004**, *36*, 666–669.
(11) Hashimoto, A.; Suenaga, K.; Gloter, A.; Urita, K.; Iijima, S. Direct evidence for atomic defects in graphene layers. *Nature* **2004**, *430*, 870–873.

(12) Meyer, J. C.; Kisielowski, C.; Erni, R.; Rossell, M. D.; Crommie, M. F.; Zettl, A. Direct Imaging of Lattice Atoms and Topological Defects in Graphene Membranes. *Nano Lett.* **2008**, *8*, 3582–3586.

(13) Kotakoski, J.; Krasheninnikov, A. V.; Kaiser, U.; Meyer, J. C. From Point Defects in Graphene to Two-Dimensional Amorphous Carbon. *Phys. Rev. Lett.* **2011**, *106*, 105505.

(14) Kotakoski, J.; Meyer, J. C.; Kurasch, S.; Santos-Cottin, D.; Kaiser, U.; Krasheninnikov, A. V. Stone-Wales-type transformations in carbon nanostructures driven by electron irradiation. *Phys. Rev. B* **2011**, *83*, 245420.

(15) Banhart, F.; Kotakoski, J.; Krasheninnikov, A. V. Structural Defects in Graphene. *ACS Nano* **2011**, *5*, 26–41.

(16) Robertson, A. W.; Allen, C. S.; Wu, Y. A.; He, K.; Olivier, J.; Neethling, J.; Kirkland, A. I.; Warner, J. H. Spatial control of defect creation in graphene at the nanoscale. *Nat. Comm.* **2012**, *3*, 1144.

(17) Girit, C. O.; Meyer, J. C.; Erni, R.; Rossell, M. D.; Kisielowski, C.; Yang, L.; Park, C.-H.; Crommie, M. F.; Cohen, M. L.; Louie, S. G.; Zettl, A. Graphene at the Edge: Stability and Dynamics. *Science* **2009**, *323*, 1705–1708.

(18) Lee, G.-D.; Wang, C. Z.; Yoon, E.; Hwang, N.-M.; Kim, D.-Y.; Ho, K. M. Diffusion, Coalescence, and Reconstruction of Vacancy Defects in Graphene Layers. *Phys. Rev. Lett.* **2005**, *95*, 205501.

(19) Cretu, O.; Krasheninnikov, A. V.; Rodríguez-Manzo, J. A.; Sun, L.; Nieminen, R. M.; Banhart, F. Migration and Localization of Metal Atoms on Strained Graphene. *Phys. Rev. Lett.* **2010**, *105*, 196102.
(20) Kim, Y.; Ihm, J.; Yoon, E.; Lee, G.-D. Dynamics and stability of divacancy defects in graphene. *Phys. Rev. B* **2011**, *84*, 075445.

(21) Gómez-Navarro, C.; Meyer, J. C.; Sundaram, R. S.; Chuvilin, A.; Kurasch, S.; Burghard, M.; Kern, K.; Kaiser, U. Atomic Structure of Reduced Graphene Oxide. *Nano Lett.* **2010**, *10*, 1144–1148.

(22) Liu, Y.; Yakobson, B. I. Cones, Pringles, and Grain Boundary Landscapes in Graphene Topology. *Nano Lett.* **2010**, *10*, 2178–2183.

(23) Lehtinen, O.; Kotakoski, J.; Krasheninnikov, A. V.; Tolvanen, A.; Nordlund, K.; Keinonen, J. Effects of ion bombardment on a two-dimensional target: Atomistic simulations of graphene irradiation. *Phys. Rev. B* **2010**, *81*, 153401.

(24) Lee, S.; Hwang, G. S. Valence force field-based Monte Carlo bond-rotation method for the determination of sp 2-bonded carbon structures. *J. Appl. Phys.* **2011**, *110*, 093524.

(25) Liu, T.-H.; Gajewski, G.; Pao, C. W.; Chang, C. C. Structure, energy and structural transformations of graphene grain boundaries from atomistic simulations. *Carbon* **2011**, *49*, 2306–2317.

(26) Leyssale, J.-M.; Vignoles, G. L.; Villesuzanne, A. Rippled nanocarbons from periodic arrangements of reordered bivacancies in graphene or nanotubes. *J. Chem. Phys.* **2012**, *136*, 124705.

(27) Lehtinen, O.; Kurasch, S.; Krasheninnikov, A. V.; Kaiser, U. Atomic scale study of the life cycle of a dislocation in graphene from birth to annihilation. *Nat. Comm.* **2013**, *4*, 2098.

(28) Krasheninnikov, A. V.; Nordlund, K.; Sirviö, M.; Salonen, E.; Keinonen, J. Formation of ion-irradiation-induced atomic-scale defects on walls of carbon nanotubes. *Phys. Rev. B* **2001**, *63*, 245405.
(29) Krasheninnikov, A. V.; Nordlund, K.; Keinonen, J. Production of defects in supported carbon nanotubes under ion irradiation. *Phys. Rev. B* **2002**, *65*, 165423.

(30) Kotakoski, J.; Krasheninnikov, A. V.; Nordlund, K. Energetics, structure, and long-range interaction of vacancy-type defects in carbon nanotubes: Atomistic simulations. *Phys. Rev. B* **2006**, *74*, 245420.

(31) Grantab, R.; Shenoy, V. B.; Ruoff, R. S. Anomalous Strength Characteristics of Tilt Grain Boundaries in Graphene. *Science* **2010**, *330*, 946–948.

(32) Chen, S.; Ertekin, E.; Chrzan, D. C. Plasticity in carbon nanotubes: Cooperative conservative dislocation motions. *Phys. Rev. B* **2010**, *81*, 155417.

(33) Wei, Y.; Wu, J.; Yin, H.; Shi, X.; Yang, R.; Dresselhaus, M. The nature of strength enhancement and weakening by pentagon-heptagon defects in graphene. *Nat. Mater.* **2012**, *11*, 759–763.

(34) Pozrikidis, C. Effect of the Stone-Wales defect on the structure and mechanical properties of single-wall carbon nanotubes in axial stretch and twist. *Arch. Appl. Mech.* **2009**, *79*, 113–123.

(35) Bagri, A.; Kim, S.-P.; Ruoff, R. S.; Shenoy, V. B. Thermal Transport across Twin Grain Boundaries in Polycrystalline Graphene from Nonequilibrium Molecular Dynamics Simulations. *Nano Lett.* **2011**, *11*, 3917–3921.

(36) Haskins, J.; Kinaci, A.; Sevik, C.; Sevincli, H.; Cuniberti, G.; Cagin, T. Control of Thermal and Electronic Transport in Defect-Engineered Graphene Nanoribbons. *ACS Nano* **2011**, *5*, 3779–3787.

(37) Kotakoski, J.; Meyer, J. C. Mechanical properties of polycrystalline graphene based on a realistic atomistic model. *Phys. Rev. B* **2012**, *85*, 195447.

(38) Lherbier, A.; Roche, S.; Restrepo, O. A.; Niquet, Y.-M.; Delcorte, A.; Charlier, J.-C. Highly
defective graphene: A key prototype of two-dimensional Anderson insulators. *Nano Res.* **2013**, *6*, 326–334.

(39) Xu, L.; Wei, N.; Zheng, Y. Mechanical properties of highly defective graphene: from brittle rupture to ductile fracture. *Nanotechnology* **2013**, *24*, 505703.

(40) Han, J.; Ryu, S.; Sohn, D.; Im, S. Mechanical strength characteristics of asymmetric tilt grain boundaries in graphene. *Carbon* **2014**, *68*, 250–257.

(41) Brenner, D. W. Empirical potential for hydrocarbons for use in simulating the chemical vapor deposition of diamond films. *Phys. Rev. B* **1990**, *42*, 9458–9471.

(42) Brenner, D. W.; Shenderova, O. A.; Harrison, J. A.; Stuart, S. J.; Ni, B.; Sinnott, S. B. A second-generation reactive empirical bond order (REBO) potential energy expression for hydrocarbons. *J. Phys.: Condens. Matt.* **2002**, *14*, 783–802.

(43) Stuart, S. J.; Tutein, A. B.; Harrison, J. A. A reactive potential for hydrocarbons with intermolecular interactions. *J. Chem. Phys.* **2000**, *112*, 6472–6486.

(44) Gale, J. D.; Rohl, A. L. The General Utility Lattice Program. *Mol. Simul.* **2003**, *29*, 291–341.

(45) Plimpton, S. Fast parallel algorithms for short-range molecular dynamics. *J. Comput. Phys.* **1995**, *117*, 1–19.

(46) Stone, A. J.; Wales, D. J. Theoretical studies of icosahedral C60 and some related species. *Chem. Phys. Lett.* **1986**, *128*, 501–503.

(47) Mehlig, B.; Heermann, D. W.; Forrest, B. M. Hybrid Monte Carlo method for condensed-matter systems. *Phys. Rev. B* **1992**, *45*, 679–685.

(48) Andersen, H. C. Molecular dynamics simulations at constant pressure and/or temperature. *J. Chem. Phys.* **1980**, *72*, 2384.

(49) Allen, M. P.; Tildesley, D. J. *Computer Simulation of Liquids*; Oxford University Press, 1987.
(50) Nice, M. *Visual Molecular Dynamics*; Urbana Champain, 1901.

(51) Jacobson, P.; Stöger, B.; Garhofer, A.; Parkinson, G. S.; Schmid, M.; Caudillo, R.; Mitterdorfer, F.; Redinger, J.; Diebold, U. Disorder and Defect Healing in Graphene on Ni(111). *J. Phys. Chem. Lett.* **2012**, *3*, 136–139.

(52) Wang, L.; Zhang, X.; Chan, H. L. W.; Yan, F.; Ding, F. Formation and Healing of Vacancies in Graphene Chemical Vapor Deposition (CVD) Growth. *J. Am. Chem. Soc.* **2013**, *135*, 4476–4482.

(53) Li, L.; Reich, S.; Robertson, S. Defect energies of graphite: Density-functional calculations. *Phys. Rev. B* **2005**, *72*, 184109.

(54) Jensen, P.; Gale, J.; Blase, X. Catalysis of nanotube plasticity under tensile strain. *Phys. Rev. B* **2002**, *66*, 193403.

(55) Ma, J.; Alfè, D.; Michaelides, A.; Wang, E. Stone-Wales defects in graphene and other planar sp2-bonded materials. *Phys. Rev. B* **2009**, *80*, 033407.

(56) Zobelli, A.; Ivanovskaya, V.; Wagner, P.; Suarez-Martinez, I.; Yaya, A.; Ewels, C. P. A comparative study of density functional and density functional tight binding calculations of defects in graphene. *Phys. Stat. Sol. b* **2012**, *249*, 276–282.

(57) This range of values is obtained by considering the relative formation energy with respect to the 5-8-5 defect (-0.9 (Ref. 18) and -1 (Ref. 19) eV), and the different absolute values published for the 5-8-5 defect.⁶,⁵⁶

(58) Obtained by adding 0.4 eV to the formation energy of the 555-777 defect.¹⁹

(59) Dumitrica, T.; Yakobson, B. I. Strain-rate and temperature dependent plastic yield in carbon nanotubes from ab initio calculations. *App. Phys. Lett.* **2004**, *84*, 2775–2777.

(60) Nachtrieb, N. H.; Handler, G. S. A relaxed vacancy model for diffusion in crystalline metals. *Acta Metal.* **1954**, *2*, 797–802.
(61) Lusk, M. T.; Wu, D. T.; Carr, L. D. Graphene nanoengineering and the inverse Stone-Thrower-Wales defect. *Phys. Rev. B* **2010**, *81*, 155444.

(62) Marks, N. A.; Cooper, N. C.; McKenzie, D. R.; McCulloch, D. G.; Bath, P.; Russo, S. P. A comparison of density functional, tight-binding and empirical methods for the simulation of amorphous carbon. *Phys. Rev. B* **2002**, *65*, 075411.

(63) Leyssale, J.-M.; Da Costa, J.-P.; Germain, C.; Weisbecker, P.; Vignoles, G. L. Structural features of pyrocarbon atomistic models constructed from transmission electron microscopy images. *Carbon* **2012**, *50*, 4388-4400.

(64) Martin, N.; Pailler, R.; Leyssale, J.-M. Carbon - ceramic (AlN) interfaces from liquid quench ab initio molecular dynamics simulations. *Mol. Simul.* **2014**, *40*, 160-168.
Graphical TOC Entry

<figure><img src="image_16.png" /></figure>