Vacancy diffusion and coalescence in graphene directed by defect strain fields

Please note that technical editing may introduce minor changes to the text and/or graphics, which may alter content. The journal’s standard Terms & Conditions and the Ethical guidelines still apply. In no event shall the Royal Society of Chemistry be held responsible for any errors or omissions in this Accepted Manuscript or any consequences arising from the use of any information it contains. Accepted Manuscript Nanoscale


Introduction
Atomic scale defects in graphene and other graphitic structures can have a profound inuence on the overall behaviour of these materials, either as the root cause of degradation that needs to be predicted and controlled, or as the source of new properties to be exploited.The 'defect engineering' of graphite, graphene and other carbon nanostructures both during fabrication and post-growth with energetic particle irradiation is now considered a key tool to modify and control important macroscopic material properties. 1,2The effect that irradiation induced point defects will have on the electrical, chemical and mechanical properties of a graphene layer will depend upon the structures and properties of aggregates and complexes that are formed if individual point defects are free to diffuse and react via thermal activation. 3,4rguably the most important point defect in graphene is the single vacancy.This defect is formed when a single carbon atom is removed from its lattice position, oen as a result of a collision with an energetic particle when under irradiation, leaving three dangling bonds on the adjacent atoms.This structure undergoes a spontaneous Jahn-Teller distortion where two of the neighbouring atoms bond weakly to form a 5 fold ring, leaving one dangling bond remaining (labelled as the 5-9 structure). 5,6This asymmetric structure can switch between the three equivalent orientations via as m a l l( ca.0.1 eV) barrier.Extensive quantum chemical calculations employing various methods have determined that the activation energy for vacancy diffusion between neighbouring lattice sites is in the range 1.1-1.3eV, [7][8][9][10][11][12][13] which means that the defect will begin to become mobile at slightly above room temperature -this has been observed and conrmed directly in real time on the graphite (HOPG) (0001) surface with atomic resolution scanning probe microscope imaging. 14When vacancies in graphene become mobile, they will react with each other and coalesce to form vacancy pairs and higher order complexes, whichcantaketheformofavarietyofdifferent morphologies.T h en a t u r eo ft h e s ed e f e c ts t r u ctures will substantially affect the properties of the sheet, 3,4,15 and predicting these morphologies is crucial for optimizing the performance of graphene based materials.Different overall defect morpholo g i e sw i l ll e a dt ov e r yd i fferent material property changes.Aggregation into a line of missing atoms creates a slot which can 'heal' leaving dislocations or stacking faults and causing a contraction in the plane.Aggregates with aspect ratios closer to one may lead either to 'holes' which will reduce elastic moduli and allow molecules to pass through the plane, or other more complex topological defects involving rotated bonds and grain boundaries.
When two individual diffusing mono-vacancies meet they will coalesce into a nearest-neighbour di-vacancy defect (equivalent to removing a carbon dimer from the lattice).This process releases approximately 8 eV (ref.7][18][19] The process of di-vacancy coalescence has been studied by several authors, 7,10,11 and found to occur readily at moderate temperatures with the highest activation barrier to be crossed in the range 1.2-2.1 eV.This di-vacancy structure is immobile, with a high barrier to thermal diffusion, but it can transform to the even more energetically stable (by approximately 1 eV) 555-777 haeckelite structure via a bond rotation, a morphology which is also observed in electron beam damaged graphene. 12,16,19The energy barrier to this transformation is greater than 5 eV, 7,10,11 and therefore this state is considered to be inaccessible from thermal activation alone, but can evidently be accessed by sub-threshold excitations from an electron beam. 20An additional bond rotation can transform this structure into the 5555-6-7777 structure 16 which is a grain boundary loop ('ower') defect. 21nce formed, the immobile 5-8-5 di-vacancy will act as a nucleation site for further coalescence of additional mobile mono-vacancies to form higher order multi-vacancy structures.The structures arising from this coalescence will be more complex, and it is not immediately clear which morphologies will result.HRTEM images of graphene reveal a number of different structures formed from electron beam damage: these oen consist of lines of vacancies running in the two low index directions, as well as extended haeckelite type structures forming grain boundaries. 16,17Vacancy lines may heal (i.e.3][24] Furthermore, a vacancy line in the armchair direction has also been observed in Scanning Tunnelling Microscopy experiments on HOPG irradiated with carbon ions. 25n general even numbered vacancy complexes will be more stable than odd numbered ones, 26 since in the former there will be an even number of dangling bonds that can mutually bind.But with just four vacancies coalescing there are several different possible stable morphologies that do not involve any bond rotations, that would require the crossing of large (>5 eV) barriers to form.7][28][29] They have all found that for a four vacancy aggregate (V 4 ) the lowest energy structure is the 'loop' or 'hole' type morphology -corresponding to a single C atom and its three nearest neighbours being removed from the lattice.This trend continues as the size of the multi-vacancy increases, with the hole-type structure lowest in energy.The two types of line isomers, running in the two main crystallographic directions (a straight line of atoms removed in the [10  10] direction and a zig-zag line of atoms in the [11  20] direction), are found to be higher in energy -even though these are the morphologies typically observed in atomic resolution images of electron beam damaged graphene. 12lthough the 'hole' type morphology, or structures incorporating bond rotations, may be lower in energy, the structures that actually result from the coalescence of vacancies at nite temperatures will be determined by the mechanism of their formation, and the kinetic pathways and barriers to the formation of each morphology or isomer.Although higher energy structures are technically metastable, the barriers that separate them from the global minimum will be substantially larger than the barrier for monovacancy diffusion.
][32][33] This is caused by the reconstruction induced by the pair-wise removal of dangling bonds.Another recent study has also determined that the activation energy for the diffusion of a single vacancy in graphene is extremely sensitive to in-plane uni-axial strain, both parallel and perpendicular to the direction of the hopping motion of a single vacancy, with the barrier doubled or halved with respectively tensile or compressive strain of just 2-3%. 34,35In this report we show that these two phenomena cooperate to pull and channel mobile vacancies into particular directions surrounding a multi-vacancy complex, through the operation of the inhomogeneous strain eld induced by the complex strongly biasing the motion of diffusing mono-vacancies.This in turn leads to the growing multivacancy adopting particular morphologies that do not necessarily correspond to the global or thermodynamic minimum energy structure.
To achieve this, we have mapped out the potential energy surface for the reaction of a mono-vacancy with both the V 2 and V 3 multi-vacancy defects, employing chemically accurate density functional theory (DFT) calculations.This, along with the potential energy surface for vacancy dimerisation, provides the complete accessible energy surface for the diffusion limited coalescence of mono-vacancies up to V 4 .We use this energy surface to simulate the coalescence process with the kinetic Monte Carlo algorithm which gives both quantitative rates for the coalescence and also denitive probabilities of different morphologies forming at nite temperatures.

Results and discussion
To simulate the graphene system with both chemical accuracy and computational efficiency, we employ the density functional theory (DFT) method with both the Local Density Approximation (LDA) 36 and Generalized Gradient Approximation (GGA). 37hese calculations were performed using the AIMPRO simulation code, 38 employing atomic basis functions 39 and HGH psudopotentials. 40The graphene defects were represented in orthorhombic supercells, containing 8 Â 4 orthorhombic unit cells (128 atoms), 12 Â 6 unit cells (288 atoms) and 14 Â 7 unit cells (392 atoms) to ensure the convergence of the results with system size.A Monkhorst 5 Pack k-point mesh of 4 Â 4 Â 1 points was determined to give total energy convergence for each system size.A vacuum gap of 1.3 nm between periodic images in the direction perpendicular to the plane is used throughout, ensuring there is no interaction.Potential energy minima are determined by optimizing the position of all atoms with respect to the total energy until the maximum force on any atom falls below 0.1 eV nm À1 .The minimum energy path between adjacent minima is calculated using the climbing image nudged elastic band method (NEB) 41,42 and chains consisting of 9 images are employed.Activation energies for diffusion are determined by the height of the saddle-point relative to the energy minimum.
Fig. 1a shows the optimized structure of the 5-8-5 di-vacancy structure in the 12 Â 6 supercell, where the total energy of this defect relative to two separated mono-vacancies (the coalescence energy) is 8.4 eV with the LDA and 8.2 eV with the GGA.The sp 2 bonds in this gure are coloured according to the bond length change in order to visualize the strain imposed upon the surrounding lattice by the reconstruction.The surrounding lattice is under tension across the defect (perpendicular to the line formed by the two missing atoms) and under compression above and below it (parallel to the line): this is a direct result of the formation of the two bonds creating the two 5-fold rings.
The V 3 , or tri-vacancy, complex will be formed as another mono-vacancy diffuses to and reacts with the di-vacancy.Far from the di-vacancy, in the undisturbed graphene lattice, this diffusion will occur via a jump between equivalent neighbouring lattice positions over the saddle point depicted in Fig. 1c.][9][10][11][12][13] In the vicinity of the di-vacancy, the distortion caused by the strain may change the height of this barrier between states and/or change the absolute energies of the minima dened by the relative position of the mono-vacancy.
To illustrate the effect that the strain created by the divacancy has on the potential energy surface for the monovacancy, we calculated the minimum energy paths using the NEB method for a mono-vacancy to move between lattice positions 1 and 2, and between positions 3 and 4 in Fig. 1a.It should be noted that the asymmetric strain from the interaction with the di-vacancy forces the mono-vacancy to adopt a single reconstruction (as opposed to the three different orientations it can assume in perfect graphene 29 ).The barrier for moving the vacancy from position 1 to 2 is signicantly higher than for the perfect lattice, at 1.8 eV (with both LDA and GGA), even though the energy difference between the two cong u r a t i o n si sl e s st h a n0 .1e V .F o rt h et r a n s i t i o no ft h e vacancy from lattice position 3 to 4 the barrier is reduced by more than half to 0.4 eV (with both LDA and GGA), even though again, the two neighbouring sites are very close in e n e r g y ,w i t ho n l ya0 .1e Vd i fference between them.These dramatic changes are a result of the sensitivity of the saddle point energy to the bond strains experienced at these different locations: the bond linking positions 1 and 2 is stretched by 1.6%andlinkingpositions3and4compressedby1.2%bythedi-vacancy reconstruction.The barrier is increased in the direction of tension and decreased in the direction of compression.
Based on the strain depicted in Fig. 1 and these two selected transitions, intuition may suggest that vacancies would diffuse towards the areas above and below the di-vacancy that are under overall compressive strain parallel to y-axis.However to denitively analyse the coalescence process it is necessary to systematically map out the entire potential energy surface for the mobile mono-vacancy by determining each possible different non-equivalent potential energy minimum for the vacancy in relation to the di-vacancy and all of the transition paths that separate them, including the nal transitions to the ground state (using the NEB method).The results of this endeavour are shown in the schematic in Fig. 2a, where the potential energy surface for the mono-vacancy is represented on the graphene lattice.The circles on lattice sites represent potential energy minima for the vacancy located at that site (i.e. the atom at that position missing) and the rectangles between them (the 'bonds') represent the saddle points on the minimum energy paths connecting the minima.Minima and saddle points are coloured according to the total energy of the system relative to the nal V 3 ground state.The dashed area represents both the nearest neighbour positions and un-stable (inaccessible) 2nd nearest neighbour positions.The potential energy surface for the mono-vacancy is determined up until the 9 th nearest neighbour positions (from the di-vacancy), where the barriers to move away from the di-vacancy differ by less than 0.05 eV from that of the unstrained graphene.It can be seen from the schematic that barriers are increased in the areas either side of the di-vacancy and decreased above and below it, corresponding to the areas of tension and compression in Fig. 1a.However, there is no direct relationship between individual bond strain and change in barrier due to the inhomogeneous nature of the strain and the three different hopping directions.The barrier for a particular mono-vacancy transition is not solely dependent on the corresponding bond length in that direction, but also on the surrounding strain: i.e. the corresponding bond lengths to the other nearest-neighbour positions.
All of the transitions that constitute the movement of the mono-vacancy consist of single atom hops of one of the three neighbouring atoms into the vacancy position (as in Fig. 1b-d), up until the mono-vacancy gets to a third nearest neighbour position from the di-vacancy.At these points, the nal transitions into the ground state occur via a more complex movement involving a cooperative motion of the 1 st and 2 nd nearest neighbour atoms, as depicted by the three non-equivalent transitions (from positions 1, 2 and 3) depicted in Fig. 2b-d.
As can be seen from Fig. 2, there is a single stable morphology resulting from the coalescence of a mono-vacancy with the di-vacancy, even though there are different pathways to its formation.This consists of two reconstructed 5-membered rings and a 10-membered ring with a single dangling bond, and is shown in Fig. 3, coloured with the corresponding bond strains.Here it can be seen that, as in the case of the di-vacancy, the reconstruction causes a tensile strain across the defect as it is pulled in from the bond reconstructions forming the two pentagons, and compressive strain emanating from the vertices of the two pentagons.This structure has been observed directly in HRTEM images of irradiated graphene. 43he formation of a V 4 defect will result from the reaction of another mobile vacancy with the V 3 structure.As in the case for the di-vacancy, we systematically mapped out the potential energy surface for a mono-vacancy reacting with the V 3 defect, including exploring all possibilities for structures formed as a result of the coalescence.The results of this are summarized in Fig. 4a, where the schematic has the same meaning as before, and the colouring scheme for the potential energy minima and barrier heights is relative to the global minimum of the V 4 (the 'hole' structure).This system has only one mirror plane, reducing the symmetry, however the effect of the strain eld converges more quickly in this case, only requiring he determination of the potential energy surface up to the 8 th nearest neighbour positions.As in the case of the formation of V 3 , all transitions are single atom hops except for the nal transitions to coalesced stable structures which consist of more complex cooperative motions.Three of the nal transitions result in two of the stable structures and are depicted in Fig. 4b-d along with their activation energies.The nal transitions to the other two structures -the global minimum 'hole' morphology and a higher energy 'U' shaped morphology -are not shown as barriers of more than 3 eV must be crossed to access them and are therefore considered thermally inaccessible at low temperatures (below about 1000 K).Unlike in the case of the formation of the V 3 , there are several (four) distinct stable structures that are the result of the coalescence occurring via different pathways on this potential energy surface.Each of these structures are stable in that all barriers to further structural change are high, and that for the technically meta-stable states, the barriers for transitions to move into the global minimum are thermally inaccessible.
The 'hole' conguration that would result from removing the under-coordinated atom from the V 3 (Fig. 3) is the global minimum for the V 4 .This coalescence process releases 8.4 eV in total (8.3 eV with GGA).The next lowest energy conguration is the arm-chair direction line (shown in Fig. 4b), which appears as two adjacent 5-8-5 di-vacancy structures (with the two coincident pentagons forming a rectangle -i.e.5-8-4-8-5), and is 1.7 eV higher in energy than the hole (1.5 eV with GGA).The formation of this conguration occurs via a more complex transition that involves the concerted motion of three atoms, nevertheless it has a relatively low barrier.
The zig-zag line morphology can result from two nonequivalent transitions, depicted in Fig. 4c and d.The transition in Fig. 4c consists of the same type of movement as in Fig. 2c (and the trans-3rd nearest neighbour collapse of the divacancy 7 ), whereas the transition in Fig. 4d consists of the same type of movement as in Fig. 2d (and the cis-3rd nearest neighbour collapse of the di-vacancy 7 ).

Kinetic Monte Carlo simulations
The transitions to the nal states for the formation of both V 3 and V 4 depicted in Fig. 2 and 4 respectively can give an indication as to which pathways and therefore structures are inaccessible due to prohibitively high barriers.However the full sequence of transitions for the diffusion of a mono-vacancy from a large separation needs to be considered to properly understand the complete coalescence process.Due to the countless different possible trajectories for the mobile vacancy to follow on these potential energy landscapes, the real-time dynamical evolution of the system must be simulated and sampled with the Kinetic Monte Carlo (KMC) method.In this approach, since the activation barriers separating each accessible state are known, we can approximate the transition rate between neighbouring states using the Arrhenius expression, choosing a pre-exponential attempt frequency equal to the Debye frequency of graphene ($10 14 Hz).The KMC method [44][45][46] employs a sequence of random numbers to choose a transition away from a particular state based on a list of the rates of all the possible transitions, and also to determine the residence time in that state.By repeating this as the system moves from state to state, the method can simulate the dynamical evolution of the system at a given temperature -provided that all possible transition mechanisms are known in advance.Assuming the rates are accurate and the system follows the Markov property (un-correlated jumps) the method is exact.The method is computationally inexpensive and therefore enables us to evaluate statistically converged behaviours over many individual trajectories.This is in contrast to a direct dynamical simulation method such as molecular dynamics, which although does not require an exhaustive search for accessible states in advance, is very limited in its sampling of the system.
The initial conguration is with the mobile mono-vacancy far from the multi-vacancy complex so that they do not interact.
The KMC simulation is then run and the system follows a particular trajectory for the diffusing mono-vacancy until it falls into a kinetically stable state (i.e. one of the nal congurations in Fig. 2 and 4).In the following simulations, the initial congurations are created by placing the mono-vacancy at a sequence of positions surrounding the complex, at the 20 th nearest neighbour away from either the V 2 or V 3 .The system is periodic in two dimensions (for the mobile vacancies) with the complex at the centre of a 7.3 nm Â 4.2 nm cell.
The animation included as ESI S1 † shows the mono-vacancy coalescing with the 5-8-5 di-vacancy at a temperature of 600 K.Each frame of this animation represents one KMC step, the duration of which in real time varies according to the calculated residence time in that state.Where the mono-vacancy is initially far enough away from the di-vacancy to be not interacting with it, it undergoes a standard random walk.When the vacancy diffuses to within 9 nearest neighbour positions of the divacancy, the diffusion ceases to be random and is biased (i.e.dri).In the case of the trajectory depicted in the animation (and in Fig. 5 a) the vacancy is repelled from the side of the divacancy until it diffuses to the region above or below where it is rapidly channelled (or 'pulled') into the position 3 (Fig. 2) -and then the structure collapses to the V 3 ground state.This animation, and the trajectory depicted in Fig. 6a, are from particular starting positions: to analyse the average behaviour we repeat this simulation 800 times with different initial congurations each followed by a different sequence of random numbers to build up a statistical ensemble.
The schematic in Fig. 5a depicts how this ensemble of trajectories accesses different states on the potential energy Here, the intensity of the shade of red at a particular lattice position is proportional to the number of trajectories from the ensemble where the vacancy has visited that position at least once.The white lattice positions indicate that these states are never occupied and are therefore completely inaccessible.It is clear from this schematic that 100% of the trajectories access the ground state V 3 via position 3 (Fig. 2) -which is a direct result of the biasing of the energy surface from the di-vacancy induced strain.Here it can be seen that the areas under tension across the di-vacancy in Fig. 1a, leading to increased barriers in Fig. 2a, are inaccessible to the diffusing mono-vacancy.At this temperature (600 K), the average time taken for the monovacancy to coalesce with the di-vacancy from the 20 th nearest neighbour separation is 0.5 s.
In the case of the formation of the V 4 from the 3 and monovacancy, the process is more intricate, and as implied by the potential energy surface (Fig. 4), there are several different possible outcomes.The animation included in ESI S2, † and the sequence depicted with the black arrows in Fig. 5b, show examples of KMC trajectories at 600 K where the mobile vacancy (positioned initially 20 nearest neighbour positions from the V 3 as before) diffuses to position 3 (Fig. 4), and then collapses into the zig-zag line V 4 morphology.The schematic in Fig. 5b depicts the kinetic accessibility of the potential energy surface in the same way as before, based on 800 separate KMC trajectories.As before, large regions of the lattice, which are under tensile strain, are completely inaccessible to the mobile vacancy.It is also clear that a substantial majority of the trajectories end up coalescing via the transition in Fig. 4d.A small proportion, however, end up diffusing to position 1 (Fig. 4), and then collapsing into the arm-chair line structure: an example of one of these trajectories is overlaid on Fig. 5b with blue arrows, and also shown in the animation included as ESI S3. † At 600 K, approximately 10% of trajectories in the ensemble result in the arm-chair morphology and the remaining 90% result in the zig-zag morphology, all via the transition in Fig. 4d.As the temperature is increased, the proportion of trajectories resulting in the armchair morphology increases, reaching approximately 2 out of 10 at 900 K.
The structures of the zig-zag and arm-chair V 4 morphologies, and the corresponding bond strain resulting from their reconstruction, are shown in Fig. 6.There are very strong similarities in the form of the zig-zag V 4 strain eld with the V 3 strain eld due to the similar reconstruction.
Again for the arm-chair V 4 , the form of the strain eld is the same as for the V 2 , but in this case with a longer range tensile strain across the defect due to the formation of the four-fold ring in the centre of the defect 24 (this will only enhance the inaccessibility of the regions either side of the defect).We should therefore expect that additional mobile vacancies diffusing to these defects will be channelled by the strain-distorted potential energy surfaces in the same way.Following the same accessibility of the potential energy surface for the V 2 +V and V 3 + V, additional mobile vacancies will coalesce with the ends of the lines, extending them in the seeded direction.This analysis strictly applies to diffusion limited coalescence, however in the case of a higher concentration of vacancies, two or more mono-vacancies may simultaneously come into the eld of inuence of a multi-vacancy defect.We expect that mono-vacancies will be channelled to the compressive regions in the same way.Due to the fact that the residence times of individual vacancies in these regions will be lowered, this in turn will reduce the chance of them coalescing with each other as opposed to coalescing with the multi-vacancy defect creating the strain eld.This will have the effect of promoting the growth of the initial multi-vacancy defect at the expense of the nucleation of additional di-vacancies.
Even though our calculations suggest that the zig-zag line will be the dominant morphology, it is higher in energy than the arm-chair line due to the increased separation of the undercoordinated atoms in its structure.However it may be energetically favourable for the two unsaturated atoms to bond, creating a 5-7-7-5 defect.This would substantially increase the strain in the surrounding lattice and require very large simulation cells to accommodate it, however it has been observed in tight binding calculations. 47As the length of the line increases, it will certainly become energetically favourable to heal, creating a pair of 5-7 dislocations 23 (see for example Fig. 7 which shows the structure of a 'healed' V 10 zig-zag line).It is clear that the increased strain due to the zig-zag line healing will dramatically enhance the channelling of mobile mono-vacancies to the ends of the line, effectively moving the 5-7 dislocation pair apart by climb as vacancies are added to the ends of the line.In principle the lowest energy state for a dipole of 60 edge dislocations is where they are separated with the line connecting them making an angle of 45 with their glide plane.It has already been shown that the 60 dislocation in graphene is immobile when it is in the 5-7 structure arising from an even number of vacancies, but much more mobile (barrier of $2 eV) in the 6-8 structure for an odd number of vacancies. 48Since the coalescence process considered here are single vacancy additions, this makes it likely that above 300 C, the constituent dipoles can move apart and migrate to the minimum energy state via thermal activation.

Conclusions
By combining a systematic exploration of potential energy surfaces using accurate DFT calculations with explicit and realtime KMC simulations, we have demonstrated how vacancies coalesce in graphene layers from thermal activation.That multivacancy complexes can induce substantial bond strain that extends into the surrounding lattice is well known.We have shown how this strain dramatically alters the potential energy surface for the diffusion of mono-vacancies, and that this results in growing vacancy complexes adopting particular morphologies -line defects.That this effect is so dramatic in this particular system is due to two factors: rstly, the extent of the strain elds, which may be explained by a combination of the force of the reconstructions (due to the energy of C-C bonds) and the 2-dimensional nature of the material.The second is the sensitivity of the activation energy for vacancy hopping to the strain -the reason for this is not entirely clear and warrants further investigation, but may be related to the remarkable in-plane stiffness of graphene, again combined with the lack of moderating interactions above and below the plane (which would be present in a 3-dimensional crystal).
The results of this investigation suggest that the thermally activated coalescence of mono-vacancies will not result in 'holes' in the graphene structure, only lines which will create 5-7 dislocation pairs and/or stacking faults and will cause inplane contraction of the material.Of course 'hole' type defects, and other types of topological defects (such as grain boundary loops) are regularly observed in atomic resolution images of electron beam damaged graphene, but we suggest here that these cannot arise from thermally activated mono-vacancy coalescence alone.In a high energy electron beam, more complex electronic effects, which we will not describe in detail here, can excite different types of transitions (e.g. bond rotations), and can enhance the ejection of low coordinated carbon atoms (which is a potential mechanism for growing 'hole' type defects). 32In spite of this, our results do shed light on what may be an important mechanism of dislocation formation and climb in graphene.The effect that strain has on the migration of defects may also play a critical role in the evolution and growth of grain boundaries and other types of extended defect.It may be the case that similar effects to those demonstrated here can have a substantial effect on the migration of other defects such as carbon adatoms or other adsorbates and also direct the morphology of growing aggregates in a similar way.This will however depend on the sensitivity of migration barriers to the bond strain and requires further investigation.
In the case of the irradiation of graphite with neutrons, the mechanism described here may explain the formation of nonbasal dislocations, which will result in a-axis contraction, from the coalescence of mono-vacancies produced in collision cascades.This is of critical importance when it comes to understanding the property changes to graphite that occur as a result of neutron irradiation when the material forms the moderator and other structural components in a nuclear reactor core. 49

Fig. 1
Fig. 1 (a) Structure of the nearest-neighbour 5-8-5 di-vacancy structure, with interatomic bonds coloured according to strain from the equilibrium graphene bond length of 0.142 nm.(b-d) The transition pathway for the diffusion of a single vacancy in graphene, where (b) is the initial structure, (c) is the saddle point and (d) is the final structure.The energy plots depict the energy barrier in perfect un-strained graphene (1.2 eV), the barrier is increased when the direction of the jump is under tension (position 1 to 2 where the corresponding bond length is stretched by 1.6%) to 1.8 eV, and decreased when under compression (position 3 to 4 where the corresponding bond length is shortened by 1.2%) to 0.4 eV.

Fig. 2
Fig.2(a) Schematic of the potential energy surface for the interaction of a single vacancy with a di-vacancy structure, relative to the ground state energy of the V 3 aggregate mapped onto the undistorted graphene lattice.The circles represent the potential energy minima for the vacancy located at particular lattice positions relative to the di-vacancy, and the rectangles the saddle-points of the minimum energy paths between them.(b-d) The final transition mechanisms and barriers to the V 3 ground state from the three non-equivalent 3rd nearest neighbour positions, where the two atoms involved in each transition have been coloured to track their movement.

Fig. 3
Fig.3Structure of the tri-vacancy ground state, with interatomic bonds coloured according to strain from the equilibrium graphene bond length of 0.142 nm.

Fig. 5
Fig. 5 (a) Accessibility of the potential energy surface for mono-vacancy diffusion in the vicinity of the 5-8-5 di-vacancy calculated from an ensemble of 800 separate KMC trajectories at 600 K.The lattice positions (states) are coloured according to the proportion of trajectories that visit the given state at least once during the trajectory.White indicates the state is not visited at all, and the darker the shade of red, the more trajectories in the ensemble have visited the state.An individual trajectory is overlaid, with arrows depicting the direction of individual jumps.(b) the same diagram for mono-vacancy diffusion in the vicinity of the V 3 .Overlaid is an individual trajectory resulting in a zig-zag line via transition 3 (black) and a trajectory resulting in an armchair line via transition 1 (blue).

Fig. 6
Fig.6(a) Structure of the zig-zag line V 4 complex, with the bonds coloured according to strain.(b) Structure of the arm-chair line V 4 complex, with the bonds coloured according to strain.