Silicon-carbon bond inversions driven by 60 keV electrons in graphene

We demonstrate that 60 keV electron irradiation drives the diffusion of threefold coordinated Si dopants in graphene by one lattice site at a time. First principles simulations reveal that each step is caused by an electron impact on a C atom next to the dopant. Although the atomic motion happens below our experimental time resolution, stochastic analysis of 38 such lattice jumps reveals a probability for their occurrence in a good agreement with the simulations. Conversions from three- to fourfold coordinated dopant structures and the subsequent reverse process are significantly less likely than the direct bond inversion. Our results thus provide a model of non-destructive and atomically precise structural modification and detection for two-dimensional materials.


Recent breakthrough developments in imaging and spectroscopy in (scanning) transmission electron microscopy [(S)TEM] have enabled the study of structural modifications that
occur very literally at the atomic scale. Due to their low dimensionality, materials such as carbon nanotubes, and especially graphene, have proven ideal for these investigations [1][2][3][4][5][6]. At the same time, (S)TEM instruments can also be turned into nano-sculpting tools: for example, graphene ribbons with specific geometries [7], or perforations of controlled sizes [8,9], can be fabricated via adjustments of the local chemistry, electron beam energy and density. Heteroatom doping is another way to tailor the properties of graphene [10], which is otherwise ill suited for many applications due to its lack of an electronic band gap [11]. Exchanging some of the carbon atoms by boron or nitrogen can result in an opening of the gap [12,13], while localised enhancements of plasmon resonances can be created around single silicon substitutions, which then act as atomic antennae [14]. The ability to directly observe the effect of single dopants is thus of the utmost importance in the further development of nano-engineering.
However, doping changes the effect an electron beam has on the atomic structure of graphene, as we have recently shown for nitrogen substitutions [15]. In that case, the slightly higher mass of nitrogen as compared to carbon leads to an increased likelihood to knock out the carbon atoms next to the dopant than the dopant itself or C in pristine areas. We also predicted that damage would be negligible at a primary beam energy of 60 keV, which was recently confirmed by atomic resolution imaging and electron energy loss spectroscopy (EELS) [5,6]. Recent studies established that silicon dopants-significantly heavier and larger in covalent radius than either carbon or nitrogen-can bond in two distinct ways within the lattice: a non-planar, threefold coordinated configuration (denoted Si-C 3 ) where the Si atom replaces a single carbon atom, binding to three neighbouring C and buckles out of the plane; and a planar, fourfold coordinated configuration (Si-C 4 ) where the Si atom is bonded to four C atoms and occupies a divacancy in the lattice [3,4]. Although beam damage was occasionally observed, apart from a study on the dynamics of Si 6 clusters in a graphene pore [16], the effects of electron irradiation have not been reported in detail.
In this Letter, we show that structural changes in silicon-doped graphene (Si-graphene) drastically differ from those in nitrogen-doped graphene. For Si-graphene, they predominantly take the form of a random walk by the Si atoms through the lattice, with no other changes in the structure. Through first principles molecular dynamics simulations, we show that each step is a result of an electron impact on one of the C atoms neighbouring the Si, and how the non-planarity of the structure plays a crucial role. The probability calculated for this process agrees well with an estimate obtained through a stochastic analysis of the experimental data. We further discuss the few observed events that lead to the conversion of Si-C 3 sites into Si-C 4 ones, and show that they are well accounted for as knock-on damage.
The significantly greater probability of the non-destructive reorganisation coupled with its directionality should allow the motion of the Si atoms to be controlled with atomic precision.
Our graphene samples were synthesised by chemical vapour deposition [17], and typically contain a low concentration of silicon incorporated as individual dopants in the lattice [3,4,18]. We observed the samples using a Nion UltraSTEM TM 100 electron microscope equipped with a cold field emission gun operated at a 60 keV primary beam energy in near-ultrahigh vacuum (2×10 −7 Pa) [4,22]. The Si atoms can be directly identified within STEM images by their brighter contrast with respect to the graphene lattice C atoms due to their higher atomic number [19]. The identification has also been verified by atomically resolved EELS [3,4].
In Fig. 1, we illustrate the predominant beam-induced process, whereby the Si atom is seen to move from one lattice site to the next during continuous imaging with an estimated dose rate of 2.2×10 7 e − /Å 2 s. Since atomic motion happens at sub-picosecond timescales, electron microscopy can only capture static snapshots of structures that are in effect fully relaxed. Although processes happening below our experimental time resolution of 88 ms cannot be ruled out, our observations give us a high degree of confidence that we capture the relevant dynamics. The observed events are non-destructive reorganisations of the structure, similar to the process of Stone-Wales transformations that were earlier shown to be due to sub-threshold electron impacts [20], rather than thermally driven bond rotations [21].
The process is also not limited to perfect hexagonal arrangements (see the Supplemental Material [22]).
To gain insight into the dynamics of such processes at the atomic scale, we used density functional theory molecular dynamics (DFT/MD) calculations as described in more detail in Refs. [15,20,23]. Although DFT describes the electronic ground state, the atomic dynamics take place over tens to hundreds of femtoseconds, whereas the relevant electron dynamics occur on sub-fs timescales in a metallic system [24]. Ionisation effects were also explicitly ruled out by experiments with 12 C and 13 C graphene [25]. Thus the ground state approximation is valid to a good degree of accuracy. In high-energy irradiation, the displacement threshold T d is defined as the minimum kinetic energy required by an atom to be removed from its position in a material. We estimated it by increasing the starting kinetic energy of a target atom until it escaped the structure during the course of an MD simulation. For atomic rearrangements, the same procedure was used to establish threshold limits for a particular reorganisation.
The calculations were performed using the grid-based projector-augmented wave code (GPAW [26,27]; for details on the computational parameters see Ref. [22]). To speed up the calculations, we used a double-zeta linear combination of atomic orbitals (LCAO) basis. However, we directly compared the calculated knock-on thresholds for C in pristine graphene and the Si atom in Si-graphene with our earlier methodology [28,29], and established agreement within our computational accuracy. Furthermore, we double-checked the Si jump threshold by full accuracy finite-differences calculations.
The non-planar Si-C 3 configuration can have two distinct positions with respect to the beam direction: the Si atom either protrudes 'below' (in the instrument geometry) the plane towards the incoming beam, or 'above' it along the direction of the beam. These present different T d and cannot be experimentally distinguished from a 2-dimensional projection of the lattice recorded with a normally incident beam. However, we suspected that the threshold for flipping a Si atom from 'below' to 'above' should not be very high even for 60 keV electrons, as no bonds need to be broken for the transformation. Indeed, a nudged elastic band [30] calculation yields a barrier of only ca. 1.1 eV for this process (see Ref. [22]).
Thus all Si configurations are effectively 'above' the graphene plane under observation, and the beam-induced alignment is expected to be thermally stable.
Our simulations yield a T d of approximately 13.25 eV for the Si atom bonded in the Si-C 3 configuration. However, due to the large mass of Si (28 amu), the probability for 60 keV electrons to transfer this much energy to the dopant is low, resulting in a cross section of 6.6×10 −7 barn when out-of-plane lattice vibrations with a Debye temperature of 1287 K are taken into account [25]. This agrees well with the observation that Si atoms are rarely lost [3,4,16]. The knock-on threshold for one of the three C neighbours to Si is higher, about 16.875 eV, but leads to a much larger cross section (∼ 0.022 barn) due to the lower mass of C (unless otherwise stated, we refer to 12 C). This process leads to a conversion Si-C 3 → Si-C 4 , which can occur at high irradiation doses.
More interesting are impacts below T d , which nevertheless lead to local changes in the structure. For example, energies between 15.0 and 16.25 eV result in the C atom being ejected from the lattice, but with a trajectory curving first slightly away and then towards the Si due to their mutual interaction. The Si simultaneously relaxes towards the vacated 5 lattice site, as illustrated for the 15 eV-case in Fig. 2 (a movie is available through Ref. [22]). to estimate a cross section for the process (see Ref. [22] for details).
We obtained an expectation value of 2.57×10 8 e − /Å 2 for the event dose, with a 95%  Figure 3c shows a histogram of the event doses, which-as expected for a Poisson process-are found to be well described by an exponential with the fitted Poisson mean (apart from an excess of one event in the third bin, and a deficit of one event in the fifth).
Our data also contains eight events where Si-C 3 is transformed into Si-C 4 (see Fig. 4).
Taking into account the total experimental dose on trivalent sites (9.71×10 9 e − /Å 2 ), we get an estimated cross section of 0.08 barn for this process, in good agreement with the calculated value of 0.066 barn [34]. Interestingly, the knock-on threshold for the four C neighbours of the Si-C 4 was calculated to be 17.125 eV, i.e., slightly higher than that for the threefold site, suggesting that Si-C 4 is more stable towards knock-on damage than Si-C 3 .
Indeed, a fourfold site was only once observed to damage further by the loss of atoms.
Each of the eight Si-C 4 sites converted back into Si-C 3 with the addition of a carbon atom (see Fig. 4 sites with a single C adatom initially bonded to C-C bridge sites 1-4 bonds away from the Si. We found the total energy of the system with the C adatom at the closest site was ∼2.3 eV lower than when the C was three or four bonds away (see Ref. [22] for details). This suggests that there is an attractive force drawing in mobile adatoms into the Si-C 4 , possibly contributing to driving the observed recombinations.
Although our analysis relied on the stochastic nature of events when all the area around the Si dopants was irradiated, it should be stressed that due to the ca. 1Å beam diameter in modern 60 keV aberration-corrected instruments, it is possible to restrict the irradiation to a chosen carbon atom. For example, we estimate that a 1.1Å beam with a Gaussian profile, centred on a selected C neighbour of the dopant, would deposit < 0.3% of the irradiation dose on the other two C neighbours. Thus the motion of silicon atoms in the lattice can in principle be controlled with atomic precision. To explore this idea, we calculated the total energies of systems with two Si atoms separated by 1-4 lattice sites, and found that the energies were lower for closer separations. Remarkably, we once also experimentally observed two Si dopants moving under electron irradiation from two sites' separation to become nearest neighbours (see Ref. [22]).

COMPUTATIONAL UNIT CELL
The unit cell for the molecular dynamics calculations needs to be large enough so that atoms at the edge of the cell do not appreciably move before the ejected atom has either escaped or reached a turning point in its trajectory, as this would cause an overestimation of the restoring forces. In our case a 8×6 graphene supercell of 96 atoms was sufficient.

TOTAL ENERGY CALCULATIONS
All calculations described below were performed in a larger 10×8 unit cell of 160 atoms (+ possible adatoms) to avoid finite size effects, using a 3×3×1 k-point grid. The structures were allowed to fully relax without constraints so that the maximum forces were lower than 0.01 eV/Å.

Relative stabilities
The relative stabilities of the Si-C 3 and Si-C 4 configurations were calculated as where E f orm denote the formation energies of the two configurations, E tot are their fully relaxed total energies, E gra tot is the total energy of pristine graphene system, E C is the chemical potential for carbon in graphene (calculated at -9.2232 eV), and E Si is the energy of an isolated Si atom (note that this cancels out in the calculation). This results in a relative energy of -1.023 eV, establishing Si-C 3 as the more stable configuration.

Si-C 3 flipping barrier
The assess the energy barrier for flipping the Si dopant bonded in the Si-C 3 configuration from 'below' to 'above' the graphene plane, we performed a simple nudged elastic band Since the relaxed starting and ending configurations are symmetric with respect to the graphene plane, the trajectory is a simple line with the maximum energy at the midpoint where the Si atom is in the plane of the lattice. The energy barrier given by the calculation is 1.08 eV ( Figure S1), and we also confirmed by a DFT/MD calculation that already a 2 eV starting kinetic energy was sufficient to dynamically induce the flip.

Si-C 4 → Si-C 3 reconversion
Since each of the Si-C 4 we observed spontaneously reconverted back into the more energetically stable Si-C 3 configuration-presumably by adatom diffusion and recombination-we studied the energetics of the Si-C 4 defect and a single C adatom (Fig. S2). In the initial configurations, the C adatom was placed on a C-C bridge site 1. . . 4 bonds away from the defect site, and then the structures were fully relaxed. During geometry optimisation, we observed that adatoms placed 1 or 2 bonds away spontaneously moved closer to the defect site, bonding directly with the Si atom. Thus although there seems to be an energy barrier for the recombination that our geometry optimisation could not overcome, there is a clear attractive interaction between the Si-C 4 defect and any diffusing C adatoms in the system.

Two nearby Si atoms
To explore the possibility of moving several Si dopants close to each other, we calculated the energetics of systems with two Si dopants incorporated into the graphene lattice separated by 1. . . 4 lattice sites ( Figure S3). We found that the energies of the configurations with both Si atoms on either one (unidirectional) or on different (corrugated) sides of the graphene lattice were almost equal, apart from the single lattice site separation, where the corrugated configuration was more stable. Furthermore, the energies were lowered when the two Si atoms were brought closer together, with a shallow energy minimum for the unidirectional configuration found at a two lattice sites' separation.
In one particular instance, we experimentally observed two Si dopants near each other in the lattice (see Figure S4). In the beginning of observation, the Si atoms occupied next nearest neighbour lattice positions. After a certain time of continuous imaging, one of the Si atoms moved one lattice site closer, making the two Si atoms nearest neighbours in the lattice (with one carbon atom possibly removed) This offers direct experimental support for the calculations shown in Figure S3, demonstrating that at least two Si atoms can be brought arbitrarily close.

STATISTICAL ANALYSES
The doses for the 38 jump events observed in our experimental data can be found in Table   I. An illustration of the Monte Carlo integration used to estimate the precise irradiation doses on each of the carbon atoms neighbouring the Si dopant the scan frame was centred on can be seen in Figure S5. If the integration resulted in a proportion of dose larger than 1 on the C atom that jumped (that is, if that C received a higher dose than the Si due to the position of the frame), the precise event dose was simply equal to the experimental dose.
Assuming the jump data are stochastic, the waiting times (or, equivalently, the doses) should arise from a Poisson process with mean λ. Thus the probability to find k events in a given time interval follows the Poisson distribution To estimate the Poisson expectation value, the cumulative doses (calculated by summing the doses incurred by those specific carbon atoms that were observed to jump) given in Table I were combined into a single dataset, which was divided into bins of width w, and the number of bins with 0, 1, 2. . . occurrences were counted. The binwidth was optimised by whereλ is the estimated mean and Z 2.5 is the normal distribution single tail cumulative probability corresponding to a confidence level of (100 − α) = 95%, equal to 1.96.
The statistical analyses were conducted using the Wolfram Mathematica software, and the Computable Document Format Mathematica script used is included as a Supplemental File. For the entire dataset, a binwidth of 1.06×10 8 e − /Å 2 was found to be optimal, yielding an estimated mean dose of 2.57 × 10 8 e − /Å 2 . However, note that when the experimental cumulative dose data is plotted along with 20 simulated random Poisson processes with the fitted mean, we see that the agreement between the data and the fit is poor due to the influence of high dose events (ie., the outliers) showing as high steps in the cumulative data (see Figure S6).
We then calculated the probabilities that each event was a random event corresponding to a Poisson process with the fitted mean, also shown Table I. This reveals that the three events with the highest doses have probabilities lower than 10 −3 to be a result of the same process as the rest of the cases. Reanalysis of the data without these outliers results in a revised expectation value of 1.63×10 8 e − /Å 2 (optimal binwidth 1.215×10 8 e − /Å 2 ). Although the Poisson fit is seemingly slightly less good than with the full dataset, a much better fit between the experimental cumulative doses and the simulated random Poisson processes with the fitted mean can be seen ( Figure S7). The expectation value was found to not be very sensitive to the value of the binwidth, provided the chosen width resulted in a good fit. Figure S8 shows a Si atom moving by one lattice site along a graphene grain boundary under electron irradiation.