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

Abstract


Introduction
The topic of point defects formation and migration in graphite has been a long subject of debate in carbon science 1 as, for instance, it is of key importance in understanding radiation damages in nuclear graphites 2 or structure-property relationships in carbon/carbon composites, especially when exposed to high temperatures. 3Even 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. 4In 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. 1 Later on, and before its actual experimental realization, electronic structure calculations (ESC) were performed on single graphene sheets used as model graphite. 5,6Kaxiras and Pandey obtained a formation energy of 7.6 eV for a bare SV 5 and a migration energy E a of 1.6 eV, involving a saddle point with a fourfold atom in the center of the vacancy (55-66 configuration). 5These 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. 6The 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 7 and 1.4 8,9 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 2 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 6

making
2][13][14][15][16] Using ac-TEM Hashimoto et al. were the first to create and image in situ some atomic defects in graphene. 11More precisey 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 2 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, [12][13][14][15][16][17] to delocalized and pseudo-amorphous defect areas. 12,13These 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.[20] Using tight-binding molecular dynamics (TBMD) at elevated temperatures Lee et al. 18 have shown that the 5-8-5 reconstruction (noted here 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. 20 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., 21 then in graphene by kotakoski et al. 13 Finally, another stable form, the 5555-6-7777 (form C) was recently observed in a few TEM experiments. 13,14,19It 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. 19All 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.2][33][34][35][36][37][38][39][40] In most of these works, the Brenner reactive empirical bond order (REBO) potentials were employed.It includes the first 41 (REBOI) and second 42 (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. 43In 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.
[30] 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 44 and LAMMPS 45 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 charectirizing the equilibrium distribution of the different reconstructions of DVs and their dynamical properties.These simulations also show the occasionnal 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.

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, 46 for comparison.The adapative interatomic reactive empirical bond order (AIREBO) potential of Stuart et al. 43 is used to describe carbon-carbon interactions in all the calculations reported in this paper.
As in former work, 26 the Lennard-Jones component of the potential is discared 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 simmulated annealing scheme using a hybrid Monte-Carlo (HMC) simulation porocedure. 26In these simulations, pure canonical HMC 47 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. 26Maximum 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 intiated at 500 K.
Formation energies were computed according to 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 (≈ 10 × 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 (≈ 4 × 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 decreased by a factor 0.75 every 2000 steps and stopped at 10 −2 K.This procedure is similar to the one adopted by Kim et al. with DFT calculations. 20

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 48 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 49 with a timestep of 0.2 fs.The collision frequency of the thermostat is set to a low value, typically 1 ps −1 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  and morphology with time.For SV defects this can simply be achieved by counting the number of non sp 2 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 50 ) 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: 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 ... 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: In these calaculations 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 equilibruium and dynamic property calculations, the initial parts of the trajectories (up to at least five mobility events) are discared.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 7 or in metal supported graphene 51,52 can differ significantly from the latter.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.94][55] 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-ofplane deformation during the geometry optimization has also been shown to have an impact on the obtained energy. 55The 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. 5,6Also 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 6 ).However, as observed by Krasheninnikov et al. for the REBOII potential, 29,30 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). 6ooking 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  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).

Energetics of vacancy defects in graphene
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.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 2 sheets in the MD simulations.Very little is said in the litterature regarding buckling in the 5-8-5 defect.El-Barbary et al. 6 mentionend 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, [18][19][20]24 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 or 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 13 respect to the larger one.5555-6-7777 (C) 6.6 6.0 6.0 19 5.7 19 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.53,54,59 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 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-  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 reasonnable 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 contant and T the temperature.

Equilibrium distributions and dynamical properties
While the AIREBO potential poorly describes the geometry of the SV defect, it behaves much better for its dynamics.7][8][9] As can be seen in Figure 6, and as expected from both experiments ans 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. 7In 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 2 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 α and P β of being in state (configuration) α or β should follow ln where H i and S i are respectively the enthalpy and entropy of form i. Measured ln The evolution with temperature of the average residence times (or lifetimes) τ where τ is the time lap between the apprition 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: τ 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 encoutered but survives for longer times when it is formed.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 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.

Lifetimes values for form B have intermediate values between these two. Remembering that the
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 recostructions of the DV defect  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 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).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 posseses 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.

Transition mechanisms
In what follows the transformation mechanisms are further investigated by looking at some se-  form J into form I (also taken from an MD ran at 3250 K).Finally in Figure 12e is shown a concerted rotation mechanism of two bonds, observed at 4000 K, transforming form C into form L in less than 2ps.
Figure 13: Three reactive sequences of configurations taken from a 92 ps long part of a DV trajectory obtained at 3500 K (the color code is redefined at the beginning of each sequence).
In Figure 13 we show a sequence of snapshots charecterizing 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 snaphots of a particularly disordered and reactive sequence 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 form 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 occsionaly 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 2 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. 60In 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 relaxtion through delocalization.To conclude this discussion we would like to add a few comments on strutures counting fourmember rings.Four-member rings have been often observed in recent TEM experiments, essentialy under the form of reconstructed vacancies V n where the number of vacancies n ≥ 4. 13,16 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 structures (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 obtained 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 62 or graphitic carbons 63,64 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.
three commonly observed DV reconctructions (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 4 :
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 ).

Figure 5 :
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 2 graphene sheets.

5
reconstruction after 228 ps (Figure 5f) by the successive closure of the four dangling bonds (an intermediate 5-db 2 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 aquire 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 betwene 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.

Figure 6 :
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.

P
α data are shown in the inset of Figure7for 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 evolves much in the temperature range considered here.

Figure 8 :
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 unsufficient data.The escape rates (1/τ) are plotted with respect to 1/T in the inset together with their Arrhenius fit.

Figure 9 :
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).

Figure 10 shows
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 τ (including those shown in Figure 8 for forms A, B and C) are given in Table4.Tables for the number of

Figure 11 :
Figure 11: Formation energies (AIREBO) of full reonstructions 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).
quences 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 Figure12.In the sequence of Figure12a, 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, Figure12billustrates 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 (Figure12c) 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 Figure12dare shown a single atom (yellow) migration mechanism transforming form B into form J, immediately followed by a bond rotation (orange-yellow) mechanism transforming

Figure 12 :
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).

Figure 14 :
Figure 14: Snaphots 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).
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. 12 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 61 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 roation of a bond joining a pentagon of a SW defect to a hexagon.

Table 1 :
Computational details of MD simulations.

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.

Table 3 :
Migration energies E m of DV through bond rotations (the case of SW defect formation in graphene is given for comparison).

Table 4 :
Average lifetimes τ (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).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 τ = 304 ps at this temperature.Form H is also a rather long-lived defect with τ = 88 ps (to be compared with 108 ps for the 5-8-5 defect at this temeprature).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).