Inclusion of radiation damage dynamics in high-resolution transmission electron microscopy image simulations : The example of graphene

Adriano Santana,1 Alberto Zobelli,2 Jani Kotakoski,3,4 Andrey Chuvilin,5,6,* and Elena Bichoutskaia1,† 1School of Chemistry, University of Nottingham, University Park, Nottingham, NG7 2RD, United Kingdom 2Labaratorie de Physique des Solides, Univ. Paris Sud, CNRS, UMR 8502, F-91405 Orsay Cedex, France 3Department of Physics, University of Vienna, Boltzmanngasse 5, A-1190 Vienna, Austria 4Department of Physics, University of Helsinki, P. O. Box 43, F-00014 Helsinki, Finland 5CIC nanoGUNE Consolider, Avenida de Tolosa 76, E-20018 Donostia-San Sebastian, Spain 6IKERBASQUE, Basque Foundation for Science, E-48011 Bilbao, Spain (Received 1 May 2012; revised manuscript received 11 September 2012; published 21 March 2013)


I. INTRODUCTION
Recent advances in electron microscopy, in particular implementation of aberration correction of electromagnetic lenses, have stimulated unprecedented growth of interest in low-voltage high-resolution transmission electron microscopy (HRTEM) capable of spatial resolution at the atomic level when using energy of the imaging electrons near or below the ejection threshold.][3][4][5][6][7][8][9][10][11][12][13][14] The ability of HRTEM to observe the dynamics of individual atoms under the controlled influence of the electron beam (e-beam) brings another dimension to the method potentially providing a tool for direct measurements of diffusion coefficients, cross sections of the damage events, chemical constants, and other characteristics of the dynamic processes that take place at the atomic scale.These advances in experimental HRTEM techniques require theoretical solutions capable of accurate treatment of the dynamic evolution of structures under the e-beam combined with subsequent image simulation at exactly the same e-beam conditions.
6][17][18][19][20] The first step in conventional image simulation includes modeling of the dynamic scattering of electrons by a sample, which can be described by a variety of well established theoretical methods such as the Bloch wave method (see, for example, Ref. 21), different variations of the multislice method, 22,23 the Born series as developed by Fujiwara, 24 the scattering-matrix method, 25,26 or the path integral method. 27The influence of electron optics of a microscope is next introduced through the instrumental optical transfer function.This approach accounts for perturbations of the e-beam caused by a sample and the imaging system, while the fact that the e-beam in turn affects the sample is usually omitted.However, it is well recognized, especially in the low-dose electron microscopy of beam sensitive biological materials, 28,29 that the e-beam may cause significant structural transformations of a sample.In general, the radiation effects of high-energy electrons on nanomaterials may include a large variety of processes (see, for example, a recent review in Ref 30) such as radiolysis triggered by ionization, local heating, direct atom ejection damage, atom rearrangements, and chemical reactions with residual gases, to mention a few.
In the majority of samples used in materials science, including small-band-gap semiconductors, metals, and graphitic materials, radiolysis and heating are suppressed due to the high density of delocalized electrons and high thermal conductivity.2][33][34][35] The atom ejection damage can be reduced or completely avoided by decreasing the energy of the e-beam below the threshold.On the other hand, the e-beam energy determines the resolution and contrast of the image obtained in HRTEM, and a careful selection of the e-beam energy is required to achieve a fine balance between reduction of the sample damage and quality of the acquired image.The resolution of HRTEM images is also limited by the aberrations of the electron optics and by the signal-to-noise ratio 32 defined by the accumulated electron dose.Thus the quality of HRTEM images can be related to sample damage by means of two experimental parameters: the e-beam energy and the electron dose.In addition to atom ejection damage some other radiation-induced processes such as back adsorption, surface diffusion of the atoms, and chemical etching may also contribute to the structure transformation.However, substantial experimental efforts are still required to determine the cross sections of these more subtle processes.
In this paper we introduce a computational framework for simulation of the dynamics of an atomic structure under the ebeam and subsequent HRTEM image simulation.The structure evolution is described by a sequence of externally initiated discrete damage events with a frequency determined by the event cross sections, which in turn depend on the energy of the e-beam.The image simulation part inherits the same energy of the e-beam in the calculation of the contrast and the same frequency of the damage events in the account for electron statistics.In order to demonstrate the simplicity and feasibility of the proposed computational approach, we study the effect of atom ejection damage on the structural transformation of a finite graphene flake into a fullerene cage under the 80 keV e-beam recently observed in experiment. 5We show that the proposed approach reproduces the experimentally observed transformation path, transformation rate, and the signal-tonoise ratio.Furthermore, we evaluate the consequences of the use of higher electron energies for imaging of the same process.We conclude that at the 200 keV e-beam the flake-to-fullerene transformation will also take place but it would be hardly observable due to the fast transformation time and low signalto-noise ratio.

II. METHODOLOGY
A computational methodology has been implemented that takes into account the effect of an external damage event on the atomic structure deformation during observation and imaging in HRTEM.It is assumed that evolution of the structure is driven by a sequence of externally stimulated damage events with different probabilities (or cross sections).After each damage event, the structure undergoes relaxation to the equilibrium within a few nanoseconds, and this process is too fast to be observed in a microscope.Between the damage events, the structure thermally oscillates around the equilibrium and does not change its chemical configuration; it is assumed that no changes occur in the structure without an external stimulus.
At the initial stage of the simulation cycle a sample is prepared and optimized at a temperature which corresponds to experimental conditions.A subsequent evolution of the sample structure under a continuous flow of high-energy electrons is described using five basic iteration steps: (1) topological analysis of a sample structure, classification of every atom with respect to the probability of an external damage event, and an assignment of the corresponding cross sections; (2) simulation of a damage event based on the probability defined by the event cross section; (3) molecular dynamics (MD) simulation of the structure evolution at the elevated temperature for the annealing to the equilibrium configuration; (4) molecular dynamics simulations at 300 K to provide "frozen phonons" configurations at the equilibrium for subsequent image simulations; (5) image simulation step which utilizes the obtained atomic configurations and applies electron statistics followed by the time upscaling.
An intrinsic assumption of the adopted approach is that evolution of the structure is driven by a sequence of discrete externally stimulated damage events each having different probabilities and, hence, cross sections.The frequency of these events is determined not only by the event cross section but also by the electron current density (flux).The total time of structure evolution can be therefore expressed as a sum of the time periods between the subsequent damage events.The time period between the events, t ev , is defined as the inverse of the product of the overall cross section corresponding to all considered damage events, σ = i N i σ i , and the electron current density, j , where σ i is the cross section of a particular type of external damage event, and N i is the number of atoms in the structure that can experience this type of damage depending on the local atomic configuration and chemical state.Equation (1)  gives the rate of the structure evolution under the influence of the e-beam, and it allows a direct comparison of the simulated process with the experimentally observed dynamics.It should be noted that Eq. ( 1) does not assume a priori any particular type of damage event and thus it has general applicability.In this study, however, we make an emphasis on the ejection damage events caused by a direct collision of electrons of the imaging beam with sample nuclei.
In HRTEM, a typical electron flux on a sample is of the order of 10 7 electrons/(s nm 2 ).In graphene, for example, there are approximately 40 carbon atoms per 1 nm 2 of the network.If ejection of carbon atoms due to direct collision with the imaging electrons is considered to be the only damage event causing structural transformation then for a typical value of the cross section of 10 barn (10 −9 nm 2 ), an interval between subsequent damage events (ejection in the considered case) taking place in the area of 1 nm 2 is 2.5 s.This implies that after each ejection, once a graphene-based system has undergone structural relaxation into a stable state (typically within a few nanoseconds) and the kinetic energy gained in the collision is fully dissipated, it will reside in the proximity of the minimum (within kT energy interval) until the next ejection occurs.The adopted event driven molecular dynamics approach is related to exact stochastic simulation algorithms, 36,37 which also use a single random number per simulation event.The assumption that radiation damage occurs merely by way of removal of atoms by a direct knock-on collision between electrons of the e-beam and nuclei of a sample is generally applicable to materials with good electrical and heat conducting properties.In electrical conductors, electron vacancies created in collisions between electrons of the imaging beam and electrons of the sample (process known as ionization) are readily filled with electrons of the conduction band.It happens on a time scale, which is much faster than that of any structural transformations, and thus in these materials ionization does not contribute directly into the structure evolution.Good heat conductivity allows a "hot spot," created in the structure by a colliding electron, to be fully dissipated before the next collision takes place, and thus accumulation of thermal energy, which could in principle initiate structural transformations, does not occur.On the contrary, in soft (bio-) materials due to the low thermal and electrical conductivity, ionization and heating are the dominant radiation damage mechanisms.The cross sections for these processes can be very high, approaching 10 4 barn, and thus the time period [Eq.(1)] between the subsequent damage events is significantly shorter.As long as the main assumption of the proposed simulation approach holds, i.e., the damage events are considered to be independent, it can be applied to describe structure evolution of soft (bio-) materials; however, simulation over an extended time period will require significant computational effort.
Through implementation of the e-beam driven morphological changes into image simulations we study the recently observed process of structural transformation of a small graphene flake into a fullerene cage. 5At the initial stage of each computation cycle, topological analysis of the structure of a small graphene flake has been performed, and every carbon atom has been classified with respect to the probability of its removal from the flake.We select three local atomic configurations within the flake with different values of the threshold displacement energy: a two-coordinated atom at the edge or an atom with a dangling bond, an atom in a pentagon site, and an atom in a bulk environment.The cross sections corresponding to ejection of atoms from the flake due to head-on collisions with the imaging electrons can be either precomputed or estimated directly from experiment. 38he displacement threshold energies have been calculated by employing the density functional theory-based tight binding model 39 in molecular dynamics simulations similar to those presented in Refs.40-42.An initial kinetic energy has been assigned to a selected target nucleus and the time evolution has been followed in order to deduce whether the given energy was sufficient to eject the target atom from the structure.This procedure has been repeated until the limiting kinetic energy, which corresponds to the displacement threshold, has been found.The threshold energies for a two-coordinated atom at the edge, an atom in a pentagon site, and an atom in a bulk environment were found to be 13.4,16.9, and 23.0 eV, respectively.The corresponding ejection cross sections were calculated from the obtained threshold energies using the McKinley-Feshbach approximation 43 taking into account the role of lattice vibrations. 38,44This approach yields for the 80 keV e-beam the cross section of 10 −6 barn for a three-coordinated bulk atom residing in a hexagon site, 0.7 barn for a three-coordinated atom residing in a pentagon site, and 14.6 barn for a two-coordinated atom at the edge; for the 200 keV e-beam the cross sections are 11.8 barn for bulk atoms in a hexagon site, 22.9 barn for bulk atoms in a pentagon site, and 34.4 barn for the edge and dangling atoms.
The initial structure of the flake has been optimized using the B3LYP/6-31G * level of theory as implemented in the Q-CHEM quantum chemistry package. 45The GULP program 46 has been employed for the constant temperature molecular dy-namics (MD) simulations using empirical Brenner potential. 47he simulation temperature was taken to be 2500 K.Each simulation step was produced with the equilibration time of 3.5 ps and the production time of 100 ps.Once the local thermodynamical equilibrium was reached the structure was quenched to 300 K, and further MD simulations produced to provide "frozen phonons" configurations at the equilibrium.A single molecular dynamics simulation run of a "flake-tofullerene" structural transformation containing 120 simulation steps takes approximately 12 h to complete in a serial run.Calculations of the frozen phonon configurations using molecular dynamics take a similar amount of time.In order to estimate fullerene yield under the e-beam radiation and acquire good statistical distribution of possible pathways simulation runs have been repeated 300 times.
The atomic configurations obtained at every step of the system evolution are subsequently used for the HRTEM image simulations performed using a standard multislice code. 20tomic thermal vibrations are considered using a frozen phonons approach 48 rather than the standard Debye-Waller factor, which allows taking into account anisotropic vibrations and fast conformational changes of the system.It is especially important for weakly connected (almost free) fragments such as carbon chains and loops.The images of 20 different "phonon" configurations have been averaged at every step of the structure evolution.The obtained images have been next superimposed by the electron statistics; namely, for a given electron flux, j , and exposition time, t exp , the intensity in every pixel has been calculated as where I sim (x,y) is the image intensity resulting from the multislice calculation and normalized to unity, and x y is a pixel size.The experimentally measured 49 modulation transfer function of a CCD camera (CCD MTF) was applied to the image.This approach gives a realistic estimation of the signal-to-noise ratio of the image, which can be obtained at given experimental conditions.As a final result, an image sequence has been produced showing the dynamic process of the atomic structure evolution under the e-beam with a realistic time scale, and a correct account for the atomic motion and the signal-to-noise ratio.The structural calculation part is linked to the image simulation part by means of two experimental parameters: (1) the energy of the e-beam, which defines the cross sections of damage events and thus the path of the structure evolution, as well as the scattering and imaging parameters in the image simulation part; and (2) the electron flux, which defines the time scaling factor for the MD calculation [Eq.(1)] and the signal-to-noise ratio in the image simulations.The described simulation approach has been implemented in the in-house software package COMPUTEM.

III. RESULTS AND DISCUSSION
In large graphene sheets, a loss of carbon atoms at the edge and its reconstruction do not cause any significant structural changes within the sheet. 1,50However, in small fragments of graphene a loss of even one carbon atom at the edge can enable a sequence of thermodynamically driven transformations which trigger the curving of the fragment and its subsequent closure into a fullerene molecule. 5Static ab initio calculations 5 of the initial steps of the flake-tofullerene evolution have demonstrated a viable route for a "top-down" mechanism of fullerene formation in HRTEM.However, in order to provide a direct link between the experimentally observed process and simulations a detailed knowledge of the dynamics and kinetics of this intriguing process is required.In this paper, we simulate the entire process of the structural transformation of a flat finite graphene flake into a fullerene cage under the e-beam radiation.We study this process at two values of the energy of the imaging e-beam, 80 V and 200 keV, and we make a comparison of these results with a focus on both the structure evolution and imaging in HRTEM.
Under the 80 keV e-beam we observe an expected etching of the two-coordinated edge atoms by the e-beam and subsequent immediate relaxation of the flake into a curved structure containing pentagons near the edge.When an atom is sputtered from the edge, the flake bends so that the edges follow the direction of displacement away from an underlying graphene substrate; thus the curving is additionally enhanced by the e-beam.As the ejection cross section of three-coordinated atoms at 80 keV is extremely low (10 −6 barn), no defects are initially generated in the internal middle part of the flake, and pentagons, which are formed at the edge, do not drift readily inside the flake.However, we do observe a rapid diffusion of these pentagons along the edge with a very small energy barrier.The edge pentagons distort the flake structure but these distortions are not sufficient to promote the closure into a fullerene cage so that the flake continues being etched away until it vanishes completely.Our MD simulations show that the probability of the flake-to-fullerene transformation is negligible in this case.In order for the curving process to begin it should be seeded, and we find that a vacancy created inside the flake near the edge where all bonds are slightly stretched can serve as such a seed.Multiple MD simulations show that once seeded, the transformation process results in an approximately 50% success rate; i.e., half of the flakes evolve into the closed cages.
Unlike the case of the 80 keV e-beam, the flake-to-fullerene transformation process takes place readily if driven by the 200 keV electrons, as in this case the ejection cross section for the three-coordinated atoms residing in hexagon sites is high (11.8barn), and defects inside the flake are easily generated.As shown in Fig. 1 the ratio of the number of two-coordinated (edge) to three-coordinated (bulk) atoms (the 2/3 ratio), which reflects the degree of closure during the transformation, remains almost the same at 80 keV (dashed line 1) and 200 keV (solid line 2).The 2/3 ratio for a flake, which remains flat, has been plotted for comparison (dotted line 3).In the size range of 80 to 50 atoms the equilibrated structures fall into the "island of stability," where the closed structure is thermodynamically favorable (the 2/3 ratio is below that of a flat flake).In the flakes containing less than 50 atoms the curvature cannot be maintained any longer, and the cage begins to open again exposing new edge atoms.Remarkably, the clusters containing fewer than 40 atoms in both cases evolve into structures less ordered than that of a flat flake, and at the latest transformation stage the weakly connected fragments of monocarbon linear chains and loops can be clearly seen.FIG. 1.The the number of (edge) atoms three-coordinated (bulk) atoms as a of the total number of carbon atoms in the model flake structure obtained for the transformation under the 80 keV e-beam (dashed line 1) and under the 200 keV e-beam (solid line 2).The dotted line 3 corresponds to a flake, which remains flat during the transformation.In the range of 100 to 45 atoms, the number of edge atoms in the relaxed structures, which lead to the closure into a fullerene cage, is lower than that of a flat flake reflecting the fact that the structure gets stabilized during the process of closure.For flakes containing fewer than 50 atoms, a curved/closed configuration is no longer energetically favorable, and cages begin to open again while the number of edge atoms grows.The structure with less than 40 atoms is no longer a cage but remains a disordered cluster with a number of edge atoms exceeding that of an ideal flake.Decrease of the number of atoms from left to right corresponds to the direction of time evolution of the structure.
We next consider the kinetics of the flake-to-fullerene transformation process determined by the cross sections of ejection damage events.Figure 2 shows the time evolution of the structures calculated using Eq. ( 1) and assuming an electron flux rate to be j = 4.1 × 10 6 electron/(s nm 2 ) as estimated directly from experimental images. 5At the beginning of the transformation process, the time evolution curve corresponding to the closure at 80 keV (dashed line 1) follows the time evolution curve of a flat flake (dotted line 3).Once the structure begins to close into a fullerene cage the value of the overall ejection cross section decreases as the contribution from the three-coordinated carbon atoms (having a negligibly small value of the ejection cross section of 10 −6 barn) increases, the structure becomes robust and stable under the influence of the e-beam and the time evolution curve 1 acquires a long plateau.This plateau extends from 80 to 40 atoms; this is the region in which the ratio of the number of two-coordinated to three-coordinated atoms is low as shown in Fig. 1.The stabilization effect is even more pronounced and the number of created perfect closed cages is higher for flakes containing about 80 atoms.It should be emphasized that the observed stability does not correspond directly to the thermodynamic stability, but it rather reflects a robustness of the structure to the influence of external stimuli such as the electron beam.When the fullerene is formed the etching process slows down significantly but still continues due to the appreciable value of the ejection cross sections for atoms residing in pentagon rings.The structure remains until it the size of about 40 atoms, at which point the cage opens again, the evolution speeds up, and the cluster vanishes very fast.It can be clearly seen from Fig. 2 that the flake-to-fullerene transformation driven by the 200 keV electrons (solid line 2) does not show any signs of structure stabilization and leads to a complete etching of the flake in a few minutes time.The proposed method for estimation of the process rate predicts the total flake-to-fullerene transformation time, which is in excellent agreement with the transformation time observed experimentally.In experiment, the closure of the cage took approximately 6 min, while the calculated transformation time is 8-10 min.The remaining discrepancy can be attributed to the fact that currently only direct atom ejection damage of a sample has been considered.
From the multiple simulations we have found out that the initial size of graphene flake dictates the size of the fullerene cage formed from it.The flakes with larger initial size close up into larger cages.For example, flakes containing about 300 atoms typically get closed into cages ranging in size between C 140 and C 160 ; these larger cages have also been found in HRTEM. 5The simulated yield of larger fullerene cages is higher than that of C 60 , and once formed large cages remain closed but continue decreasing in size through the loss of carbon atoms in pentagon rings and subsequent bond rearrangements.A similar mechanism has been proposed by Irle and Morokuma in the "shrinking hot giant" scheme for fullerene production, [51][52][53] and by Yakobson for the observed process of shrinkage of giant fullerenes, 54 in which smaller fullerenes have been created in TEM by removal of carbon atoms from a giant fullerene.If the flake contains several thousands of carbon atoms the curving step carries a large energetic penalty due to a significant van der Waals interaction with the underlying graphene sheet.Therefore, the edges of larger flakes continue being etched until they reach a size suitable for the fullerene formation, typically between 110 and 350 carbon atoms.On the other hand, the transformation of very small flakes, less than 60 carbon atoms, into fullerenes is suppressed by the excessive strain on the carbon-carbon bond imposed by the high curvature of small fullerene cages.Indeed, our simulation experiments confirm that fullerene cages formed from graphene flakes in TEM fall into a relatively narrow range of sizes averaged at approximately 1 nm in diameter.Typically, the structure of a newly formed fullerene does not correspond to the most stable isomer and contains seven-member rings as well as abutting pentagons.Under the continued high-temperature exposure, such imperfect cages could anneal via a series of Stone-Wales rearrangements 41,55 and other annealing schemes 56 giving rise to perfect fullerenes, which obey the isopentagon rule.
The number of electrons arriving per unit area of a sample before a damage event takes place can be expressed as n el = 1/σ , where σ (having the dimension of area, barn) is the overall cross section of damage events used in Eq. ( 1).The value of n el therefore defines the electron dose (number of electrons per unit area) available in the electron microscope to create an image of a sample at a given step of the sample evolution.Note that we assume that structural relaxation after each event happens within a few nanoseconds, which makes this process too fast to be observed in the microscope.An equilibrium structure adopted between the damage events stays unchanged for a few seconds and can be therefore imaged.Between the damage events the system undergoes only thermal oscillations around the equilibrium and does not change its chemical configuration.The knowledge of the electron dose available for imaging allows estimating a signal-to-noise ratio (SNR) achievable while imaging this evolution step.In such estimation the dose should be related to an area corresponding to the expected resolution, i.e., n el should be multiplied by (r/2) 2 , where r is the resolution.The signal, S, is defined as the image intensity (number of electrons) multiplied by the contrast of the atom, C, Similarly, the statistical noise, K, attributed to the resolution, r, is defined as Therefore for a given damage process with the overall cross section σ , initiated by the imaging e-beam, the maximum SNR achievable in the individual imaging step can be defined by the ratio of the signal, S, to the statistical noise, K, given by Eqs. ( 3) and ( 4): This relation, in different forms, is well recognized in low-dose electron microscopy, 31,57,58 and it has been recently derived for the atomic imaging under low-dose conditions. 59n addition to the concept of "critical dose" previously considered in the literature, our approach attributes configurational The values of the image contrast and the instrumental resolution are taken to be C 200 keV = 0.035, r 200 keV = 0.1 nm, C 80 keV = 0.025, r 80 keV = 0.14 nm, as described in the text.Under the 200 keV e-beam, the signal-to-noise ratio is close to the critical value of 5 during the whole transformation process due to the large values of cross sections for a damage event.At the 80 keV e-beam, the signal-to-noise ratio remains above 10 at all transformation steps, and it reaches the value of 25 for the stable closed fullerene cage.For a small system the SNR value increases in both cases due to decrease in the overall damage cross section.changes of the object to a number electrons available for imaging per object configuration, and thus Eq. ( 5) allows estimation of the SNR for every step over the entire duration of structure evolution, rather than an integral value for the whole process.Figure 3 shows the changes in the SNR value during the studied process of flake-to-fullerene transformation at 80 and 200 keV.The following values of the image contrast, C, and the instrumental resolution, r, have been used: C 200 keV = 0.035, r 200 keV = 0.1 nm, C 80 keV = 0.025, r 80 keV = 0.14 nm.These values have been obtained from the simulated images before application of the electron dose and CCD MTF, so that the curves shown in Fig. 3 represent the reference values assuming an ideal image detector.The SNR for the flake-to-fullerene transformation at 200 keV is very low for the entire duration of the transformation, and it is close to a critical value of 5 known as "Rose detection criterion." 60As Eq. ( 5) shows this is due to the large values of the total cross sections predicted for this case.We can therefore conclude that at 200 keV the imaging of the studied process is limited by the image noise rather than by the instrument optics.For the same process driven by the 80 keV imaging electrons, the value of the SNR is higher than 10 at all steps of the transformation, and it reaches 25 for the "island of stability" corresponding to the closed fullerene cage.Equation ( 5) gives merely an upper estimation of the achievable SNR.There exists a number of further details related to image acquisition such as movements of an object, detector modulation transfer function (MTF), detector dark noise, fixed exposition time, dead time for data transfer, limited detector frame rate, etc. which could also be accounted for in order to achieve a true likeness with experimental images.
The proposed link between molecular dynamics and image simulation allows reproducing very closely the existing experimental data as well as providing realistic snapshots of the flake-to-fullerene transformation at different imaging conditions.To illustrate this modeling capability we simulate the same transformation pathway predicted by molecular dynamics simulations in three different HRTEM imaging conditions: (1) 80-kV accelerating voltage and 4.1 × 10 6 electron/(s nm 2 ) electron flux; these imaging conditions have been used during the observation of the top-down fullerene formation in HRTEM; 5 (2) 200-kV accelerating voltage and 4.1 × 10 6 electron/(s nm 2 ) electron flux; (3) 200-kV accelerating voltage and 0.25 × 10 6 electron/(s nm 2 ) electron flux.Imaging condition (3) leads to the formation of fullerene in approximately the same time (1500 s) as case (1) since these conditions provide the same rate of transformation.For all three cases the quality of the imaging system was set to be the same; namely, correction of coherent aberrations provided a flat phase plate up to 20 mrad and incoherent aberrations resulted in a 4-nm focus spread.Among a variety of possible corrector settings [61][62][63] we have selected those used in experimental study: 5 positive value of spherical aberration (Cs = 10 μm) and defocus close to Scherzer 64 value (df = −8 nm for 80 kV and df = −6 nm for 200 kV).At these imaging conditions, resolution is determined by incoherent envelope, and the contrast of atoms does not change substantially at slight variation of Cs and df .Following classical image simulations, the electron dose has been emulated on the images as described by Eq. (2), 1 s exposition time has been assumed, and the MTF of the CCD camera for corresponding beam energy has been applied.
Figure 4 represents the snapshots of the simulated video obtained for the three sets of imaging conditions.The full video is available in Supplemental Material 65 (CompuTEM.avi).Due to a correct account for MTF of CCD 66 images at 80 keV reproduce well the contrast of carbon atoms and the SNR obtained experimentally.As mentioned above, the time scale of the transformation in the simulated video at 80 keV is in close agreement with the experimentally observed process.Image simulations for the presented video took about 7 h on a Core i7 1.6GHz 4Gb RAM notebook running Windows 7 64-bit (10 s per one image of process step times 20 phonon configurations per image times 120 steps).

IV. CONCLUSIONS
We have developed a general computational approach for the inclusion of the dynamics of atomic rearrangements under the imaging electron beam into transmission electron microscopy image simulations.The effect of the electron beam on atomic structure is described by the event driven molecular dynamics, in which events of the electron collision with a sample are represented by the corresponding displacement cross sections specific to the energy of the electron beam and the type of atom interacting with the imaging electron.The influence of the sample on the imaging electron wave is also included and described using a standard multislice approach utilizing the snapshots of structure transformation and thermal vibrations obtained from molecular dynamics simulations.A software implementation of this approach, COMPUTEM, provides a sequence of images of sample evolution with time during its observation in HRTEM for given experimental conditions.Two main computational parts, molecular dynamics and image simulations, are linked by the experimental value of the electron flux, which allows the upscaling of MD simulation time to the experimental time and determines the signal-to-noise ratio of the simulated images.
The recently observed process of the structural transformation of a flat finite graphene flake into a fullerene cage in HRTEM 5 has been selected for a case study.The simulated image series showing fullerene formation under the 80 keV e-beam closely reproduces experimental HRTEM image series in regard to the structure evolution route, evolution rate, and signal-to-noise ratio.It has been shown that under the increased e-beam energy of 200 keV similar observations will be obscured by high damage rate and/or low signal-to-noise ratio.
The proposed approach adds another dimension to the existing image simulation concept as it includes robustness of a sample as one of the input parameters.Thus, the HRTEM imaging conditions can be evaluated and optimized not only with respect to instrumental resolution and contrast, but also with respect to the electron dose.This approach will be especially beneficial for the-beam sensitive samples, where resolution is defined by the signal-to-noise ratio rather than by the quality of the instrument.

3 FIG. 2 .
FIG.2.The time evolution of the flake-to-fullerene transformation under the 80 keV e-beam (dashed line 1) and under the 200 keV e-beam (solid line 2).The dotted line 3 corresponds to a flake, which remains flat during the transformation; in this case the ejection cross sections corresponding to the 80 keV e-beam have been used, and the edge atoms have been removed one by one.A pronounced stabilization of the structure was observed at 80 keV due to the formation of a closed fullerene cage.The stability plateau ranges from 80 to 40 atoms.All curves are calculated for the electron flux of 4.1 × 10 6 electron/(s nm 2 ).

2 FIG. 3 .
FIG.3.The signal-to-noise ratio as a function of the total number of carbon atoms in the model flake structure obtained for the transformation under the 80 keV e-beam (dashed line 1) and under the 200 keV e-beam (solid line 2).The curves represent data averaged over ten simulation runs.Decrease of the number of atoms from left to right corresponds to the direction of time evolution of the structure.The values of the image contrast and the instrumental resolution are taken to be C 200 keV = 0.035, r 200 keV = 0.1 nm, C 80 keV = 0.025, r 80 keV = 0.14 nm, as described in the text.Under the 200 keV e-beam, the signal-to-noise ratio is close to the critical value of 5 during the whole transformation process due to the large values of cross sections for a damage event.At the 80 keV e-beam, the signal-to-noise ratio remains above 10 at all transformation steps, and it reaches the value of 25 for the stable closed fullerene cage.For a small system the SNR value increases in both cases due to decrease in the overall damage cross section.

FIG. 4 .
FIG. 4. Simulations of the flake-to-fullerene transformation for three sets of experimental conditions in HRTEM.Top panel: 80 keV accelerating voltage and 4.1 × 10 6 electron/(s nm 2 ) electron flux, these values correspond to the experimental conditions used during the observation of fullerene formation; 5 middle panel: 200 keV accelerating voltage and 4.1 × 10 6 electron/(s nm 2 ) electron flux; bottom panel: 200 keV accelerating voltage a nd 0.25 × 10 6 electron/(s nm 2) electron flux, these experimental setting determines the same rate of the flake-to-fullerene transformation as in the case shown in the top panel.For all three cases the quality of imaging system was considered to be the same and described in the text.