## Abstract

Atomic engineering is envisioned to involve selectively inducing the desired dynamics of single atoms and combining these steps for larger-scale assemblies. Here, we focus on the first part by surveying the single-step dynamics of graphene dopants, primarily phosphorus, caused by electron irradiation both in experiment and simulation, and develop a theory for describing the probabilities of competing configurational outcomes depending on the postcollision momentum vector of the primary knock-on atom. The predicted branching ratio of configurational transformations agrees well with our atomically resolved experiments. This suggests a way for biasing the dynamics toward desired outcomes, paving the road for designing and further upscaling atomic engineering using electron irradiation.

## INTRODUCTION

Controlling the exact atomic structure of materials is an ultimate form of engineering (*1*, *2*). Atomic manipulation and atom-by-atom assembly can create functional structures that are hard to synthesize chemically (*3*–*6*), e.g., exactly positioning atomic dopants to modify the properties of carbon nanotubes and graphene (*7*). Nitrogen (N) or phosphorus (P) dopants might be useful in quantum informatics due to nonzero nuclear spin (*8*).

Successful atomic engineering requires understanding of two parts: (i) how the desirable local configurational changes can be induced to increase the speed and success rate of control and (ii) how to scale up basic unit processes into feasible structural assemblies of 1 to 1000 atoms to produce the desired functionality. Historically, scanning tunneling microscopies (*9*) have demonstrated good stepwise control of single atoms, leading to physicochemical insights and technological advances (*10*). However, their scalability and throughput are severely limited by the mechanical probe movements. Recently, aberration-corrected scanning transmission electron microscopy (STEM) has emerged as a versatile tool for characterizing the precise atomic structure of materials (*11*–*18*). Despite its very early stage of development, STEM also shows great promise as a tool for atomic manipulation: In two-dimensional (2D) graphene, Si dopants are found to be stepwise controllable (*19*–*21*), and iterating these basic steps enables long-range movement with a high throughput (*22*), whereas in a 3D silicon crystal, the projected location of Bi dopants was also manipulated (*23*). However, imprecise understanding of the dynamics of the basic steps, which involves relativistic electron-nucleus collisions, electronic excitation and relaxation, dynamic ion trajectories, momenta dephasing, and heat conduction, add uncertainties to this technique. While the traditional theory of radiation damage provides a basis for understanding, instead of trying to minimize beam effects, atomic engineering seeks to utilize them to achieve desired configurational changes. Concepts like the displacement threshold energy *E*_{d}, which is, in most cases, approximated as scalar, turn out to be too crude to guide the design of the precise cross-sections of different configurational outcomes (*24*, *25*).

Here, we use STEM to both drive and identify the atomic motion of individual phosphorus (P) dopants in graphene and construct a theoretical scheme for evaluating their relative probabilities with respect to electron energy and momentum direction. We have categorized the dynamics into four types: (A) direct exchange and (B) Stone-Wales (SW) transition, which conserve the atoms, and (C) knockout of a carbon neighbor and (D) replacement of the dopant atom by C, which do not conserve the local composition. We choose to use *c* = 1.3377 × 10^{8} m/s) to minimize (C) and (D) while maximizing the rates of direct exchange and SW transition. Hereafter, the variables with tilde (~) on top represent quantities before atom-electron collision, and variables without tilde represent those after collision. Furthermore, instead of aiming the electron beam directly at the dopant itself, it has been established that dynamics can be more effectively induced when the electron beam is aimed at a carbon neighbor of the dopant (*3*, *21*).

To achieve atomic configurational change, the postelectron collisional energy of the primary knock-on atom (PKA; here, it is carbon), *E*, needs to be on the order of 10 eV. This requires the penetrating electron to pass very close to the PKA nucleus (impact parameter *b < b*_{c} ~ 10^{−14} m), with corresponding collisional cross section on the order of barns (σ ~ 10^{−28} m^{2}). Such elastic collision and large energy transfer occur mainly within 0.1 zeptosecond time scale (τ_{c} ~ 10^{−22} s), inducing a postcollisional PKA momentum labeled by a vector (θ, φ, E). With a total beam current of *I* ~ 50 pA, this amounts to about one relativistic electron penetrating the graphene every 3 ns (τ_{p} ≡ *e*/*I* ~ 3 ns), and one can focus the e-beam to a spot with a full width at half maximum (FWHM) of ~1 Å (which provides a sufficient description of the scanning beam). The collisional probability (defined as imparting the PKA with *E* ~ 10 eV energy, which may cause “immediate” configurational change within picoseconds) is thus only ~σ/FWHM^{2} ~ 10^{−8} per penetration event or σ/FWHM^{2}/τ_{p} ~ 10 per s (0.1 per s for events like direct exchange with cross section of 0.01 barn); the rest of the penetration events cause electronic excitation and small ionic rattling but not immediate local configurational change.

Regardless of whether a penetrating electron gets within *b*_{c} or not, a penetration event will cause electronic excitation, occurring with an attosecond time scale τ_{e} ~ 3.4 Å/^{−18} s (3.4 Å being the graphene thickness), which, however, in the case of graphene, will relax collectively on the femtosecond time scale (τ_{E} ~ 10^{−15} s) to the electronic ground state (*26*). Thus, after τ_{e} + τ_{E}, the electronic subsystem falls back to electronic equilibrium, one may use the Born-Oppenheimer (BO) approximation to describe the ion dynamics, which can achieve either one of the (A) to (D) configurational changes (labeled by *i* = 1 to 4) or remain unchanged (*i* = 0) on the BO surface, within a few picoseconds (τ_{I} ~ 10^{−12} s). Since τ_{I} > > τ_{e} + τ_{E}, it is justified here to apply ground-state density functional theory (DFT) to track the main portion of the ion dynamics and to obtain the probability of success, *P _{i}*, of atomic dynamics that lead to configurational outcome

*i*.

Throughout *t* = τ_{c}, τ_{E}, τ_{I}, the PKA momentum history needs to be tracked; thus, we build a theoretical scheme called primary knock-on space (PKS) for estimating the relative scattering cross sections of different electron-induced dynamics due to either sample or electron beam tilt and for selectively activating the desired outcomes. We further provide experimental verification of our calculations, thus opening new avenues for atomic engineering using focused electron irradiation.

## RESULTS

We find that the P dopant in graphene can serve as a good example for covering many categories of electron-induced dynamics. With highly collimated and focused (e-beam FWHM ~1 Å) electron irradiation on a carbon atom neighboring the phosphorus dopant, we occasionally create a single energetic PKA, with rate ~*d*σ/FWHM^{2}/τ_{p}, where *d*σ is the differential cross section corresponding to a particular postcollisional PKA differential momentum volume. To clarify, the term PKA here exclusively refers to the energetic carbon neighbor of the phosphorus dopant, so “PKA” and “C neighbor hit by an electron of the beam” are equivalent throughout this paper. This energetic carbon atom then drives a short burst of atomic motions nearby within τ_{I} ~ picoseconds.

In Fig. 1, four types of dynamics are shown and categorized into two groups: atom-conserving dynamics (which is desirable) and atom-nonconserving dynamics (which is often not desirable). Atom-conserving dynamics include (A) direct exchange between P and C [(Fig. 1A), earlier dubbed “bond inversion” in the context of Si (*21*)] and (B) SW transition (Fig. 1B), i.e., 90° rotation of a P─C bond (*25*). Atom-nonconserving dynamics include (C) knockout (Fig. 1C), where the PKA is knocked out by the electron beam (P turns from three fold-coordinated to four fold-coordinated, after which we found that it is no longer possible to further manipulate the configuration with 60-keV e-beam) and (D) replacement (Fig. 1D), where a diffusing carbon adatom that happens to be nearby receives energy from a penetrating electron and replaces the dopant atom, which then diffuses away quickly. These wandering C adatoms are always present on graphene surfaces (*21*, *27*), but they diffuse too quickly to be imaged. In the above experiments, we scanned the beam over a square area covering the dopant atom so that the configurational changes could also be captured in frames (often as a broken “transit” frame, where part of the scanned image is discontinuous with the rest of the image that is scanned later).

In Fig. 1A, three consecutive frames of direct exchange including a transition frame are presented. As a result, the P dopant atom exchanges site with the PKA, while the e-beam is scanning from left to right across the PKA (white dashed line; note that, at each pixel, most of the electron dose is distributed within an ~Å-sized area surrounding it according to the beam intensity profile). In Fig. 1B, a SW transition is preceded by a direct exchange. After the direct exchange (frames 1 to 2), the P─C bond is rotated by 90° (frames 2 to 3), and the hexagonal lattice is locally transformed into two pairs of five- and seven-membered rings (55-77 structure hereafter). The 55-77 structure is only stable for less than 0.2 s before reverting back to hexagons (frames 3 to 4) due to the subsequent electron irradiation. In Fig. 1C, a threefold-coordinated P (frame 1) turns into fourfold-coordinated P (frame 2) when the PKA is knocked out by the electron beam. Once this happens, we find that the P can no longer be manipulated. In Fig. 1D, P is replaced by C, and then diffuses away, which is the most commonly observed outcome of P impurities—in stark contrast to Si, which are almost never removed or replaced. It should be noted that we never observed a phosphorus being simply knocked out leaving a vacancy behind, consistent with the prediction that its knockout cross section is several orders of magnitude smaller than that of the lighter C atoms.

As a basic test of controlling the P dopant for atomic engineering, a direct exchange is intentionally initialized by targeting the highly focused e-beam at a neighboring C atom. Since the out-of-plane dynamics of the energetic C neighbor are responsible for the change in the structure (*21*), the outcome of the exchange can be controlled by selecting the PKA among the three possible carbon neighbors. The initial position of the P dopant is shown in Fig. 1E. The yellow crosses indicate where the electron beam is parked for 10 s, and afterward, a second frame is immediately captured, as shown in Fig. 1F. As a result, the P atom hops to the site as directed, but this occurred only after 68 ineffective 10-s spot irradiations (another P jumped after 12 10-s spot irradiations). Compared to Si impurities, P is much harder to induce direct exchange with (Fig. 1A): Irradiating the neighbor C site typically triggers the replacement process (Fig. 1D) instead. We tried to manipulate 10 P impurities, 2 of which exchanged (outcome A or B), 1 lost a C neighbor (outcome C), and 7 were replaced by freely diffusing carbon adatom (outcome D) after, on average, 22 ± 5 (mean ± SE) 10-s spot irradiations. Note that because the electron beam is not scanning continuously, we cannot resolve SW transition, which was actually found to have the largest cross-section (Table 1) when we later scanned the beam, since the 55-77 metastable configuration is quite beam-sensitive.

To reduce the replacement of the dopant by freely diffusing C adatoms, we also used double-layer graphene (fig. S1), where atom diffusion on one side is suppressed. It is interesting to observe that the phosphorus dopant in a double-layer graphene is much less likely to be replaced than in single-layer graphene. With a similar dose rate, the P atom was not replaced during our observation (~12 min), which is more than four times longer than in single-layer graphene (~3 min). It should be noted that the difficulty of manipulating P atoms represents a generic challenge in atomic engineering, where a desired configurational outcome is disrupted by other unwanted ones. Our work is specifically focusing on dealing with this issue.

To explain these processes, we have performed extensive ab initio molecular dynamics (abMD) and climbing-image nudged elastic band (cNEB) calculations. With a clear separation of time scales, in particular, τ_{E} < <τ_{I}, it is a reasonable approximation to simulate the configurational change processes on the BO surface, assuming each dynamic step evolves according to the Hellman-Feynman forces calculated on the basis of the electronic ground state.

The distribution of various types of dynamics is shown in Fig. 2 (A to C), which corresponds to initial postcollision kinetic energies of the PKA at *E* = 15, 16, and 17 eV, with the angular space sampled with an interval of 15° for the azimuthal angle φ and 5° for the polar angle θ (up to 25°). Figure 2 (D to G) shows four examples representing different dynamical processes, shown in the order of SW transition, knockout, direct exchange, and unchanged structure. All of these beam-induced dynamics of P dopants are initiated by an out-of-plane momentum imparted on PKA by the backscattering of a single electron, which occurs stochastically with a small probability. The definitions of spherical coordinates θ and φ (momentum direction of the PKA whose energy is *E*) are plotted in the first frame of Fig. 2G, along with an example of an unchanged structure (θ = 25°, φ = 285°, with the kinetic energy *E* = 15.0 eV). If the initial momentum is not strictly upward, but tilted at an angle (θ = 20°, φ = 75°, *E* = 15 eV in this case), then a SW transition occurs (Fig. 2D) (*25*). As an example of knockout in Fig. 2E, the initial momentum of PKA is tilted toward θ = 20°, φ = 180°, with *E* increased to 17.0 eV. Last, in Fig. 2F, an initial PKA momentum perpendicular to the plane (θ = 0°) yields a direct exchange when *E* = 17 eV.

From these plots, several conclusions can be drawn for the phosphorus dopant: (i) A SW transition can be initiated with a lower PKA energy *E* (starting from 15 eV) than direct exchange. (ii) Increasing from 15 to 17 eV, direct exchange gradually becomes the dominant dynamical process. (iii) When *E* reaches around 17 eV, knockouts begin to occur. (iv) Somewhat counterintuitively, direct exchange is easier when the PKA momentum is pointing away from the target P atom (φ = 180°), instead of pointing toward it (φ = 0°). As we shall see, these polar plot features are derivable from the PKS theory, from which the relative scattering cross sections of each configurational outcome can be estimated.

The replacement dynamics (Fig. 1D) are due to the freely diffusing C adatoms on graphene surface. In Fig. 2H, our calculation shows that C adatoms can bond stably on a C─C bridge close to the underside of a P site (shown as the initial state). By performing a cNEB calculation, we see that, to transit from this initial state to a final state where the P has been replaced by C, the system only needs to cross a 0.4-eV barrier, available thermally or from the 60-keV electron beam (*28*), and subsequently reducing the potential energy of the system by 4.5 eV.

Comparing different graphene dopants, we found that P hops less actively in experiment than what has been reported for Si (*22*). To explain this, we compare the PKS-predicted energy range of direct exchange for Si, P, and Al cases when assuming a head-on collision (θ = 0°; Fig. 3A). We find that the Si case covers the greatest energy range, resulting in larger probability of direct exchange than for P. The displacement threshold of the C neighbor of an Al dopant is much lower than those for Si and P, so knockout of the PKA is a more likely outcome for Al dopants. We have observed an Al dopant and its surrounding atoms being knocked out together by a 60-keV electron beam (Fig. 3B), while we never observe this process for Si or P at the same electron energy. This implies that a lower acceleration voltage (

On the contrary, a SW transition is more likely to be observed for a P dopant compared with Si. Related cNEB calculations are shown in Fig. 3C. As a broader comparison, we compute six elements, including C, N, B, P, Si, and Al, any of which theoretically could experience a SW transition. To be able to observe the SW transition in STEM, the 55-77 structure must be sufficiently stable to capture an image frame. Its stability is proportional to the depth of the potential energy well of the 55-77 structure (energy barrier between the highest energy transition state and the 55-77 structure), which is given as the activation energy *E*_{a} (the original cNEB curves can be found in fig. S4). We note that the cNEB calculation can only provide qualitative ranking, not quantitative characterization of the beam-induced dynamical process, since the electron-imparted momentum is localized on the PKA and does not necessarily exactly follow the collective reaction pathway of the minimum energy path. The stability of 55-77 structures follows the order C > N > B > P > Si > Al. Among all the dopants we observed, we indeed find that N has the most stable 55-77 structure [Fig. 3D; the single-atom electron energy-loss spectroscopy (EELS) characterization of this particular N dopant can be found in (*29*)]. Purely thermally, for a preexponential factor of 2 × 10^{12}/s estimated from harmonic transition-state theory analysis in (*30*), the Si 55-77 structure back-transformation rate at 300 K is 0.073/s, making such defects (and all the dopants with a higher energy barrier), in principle, STEM observable if they are created.

### PKS formalism

Predicting and comparing the scattering cross sections of different dynamic processes within a unified framework is essential for atomic engineering, so we have developed the PKS formalism. Illustrated on the polar plots in Fig. 2 (A to C), the azimuthal angle φ and polar angle θ correspond to the direction of the momentum of postcollisional PKA (Fig. 4A), and the radius of the polar plot represents its kinetic energy *E* (Fig. 4B). Every point in PKS describes the momentum status of the PKA in terms of its momentum direction and kinetic energy right after collision (*t* = τ_{c}), all of which lead to a dynamic outcome that correspond to the points in Fig. 2 (A to C). In Fig. 4C, these outcomes are grouped to differently colored blocks represented in three dimensions in PKS. The momentum distribution of PKA after an electron collision has an ovoid profile, whose shape changes with respect to the energy and direction of an incoming electron and the precollisional momentum of the atom. We conceptualize this momentum distribution as a “Doppler amplification effect” because small changes in precollisional momentum can lead to a much greater change of the postcollision momentum, as illustrated in Fig. 4B. While the theory presented in Materials and Methods is more general, for the illustration here, only atoms vibrating perpendicular to the graphene plane are considered. The Doppler amplification effect is essential here because our calculation shows that, if there were no precollisional kinetic energy of the atom, then there would not be a chance of direct exchange, SW transition, or knockout of a carbon neighbor in the experiments (Fig. 1). In Fig. 4C, the intersection of the colored regions and the ovoid of a vibrating carbon atom (we use *24*, *31*). A detailed formulation of the PKS theory can be found in Materials and Methods.

For atomic engineering, we use a decision tree to show the possible paths of evolution (Fig. 4E), with its root node indicating the initial structure, and the child nodes indicating the next possible structures with different branching probabilities. In this example, the root trifurcates into three different paths, corresponding to the three different dynamic processes. By using the PKS formalism and taking the mean squared velocity of the carbon atom to be 3.17 × 10^{5} m^{2} s^{−2} (*32*) or *5*, *33*), partly due to possible inaccuracies in the description of bond breaking in standard (semi-) local DFT (*34*), and the scheme we have constructed here does not consider electronic excitations. For these reasons, we considered only relative probabilities in the PKS above. Further improvements in theoretical modeling are needed before quantitative predictions are possible for impurity sites in graphene (*5*) and to account for the in-plane vibration components of the C atoms. However, we want to stress here that PKS correctly predicts that the SW transition has a much higher probability than direct exchange, whereas static in-plane transition paths calculated with cNEB cannot rationalize the branching ratio of the two processes.

To further experimentally test our theory, we tilted a Si-doped sample so that the electron beam was incident at a specific angle (

## DISCUSSION

The long-term vision of atomic engineering is to precisely position individual atoms with desired internal states, including the nuclear spin, image and control atomic assemblies from 1 to 1000 atoms, and use precisely controlled atoms and their electronic/nuclear states for devices such as atomic clock and memory. Successful atomic engineering comes from understanding two parts: (i) how the desirable local configurational changes can be induced to increase the speed and success rate of control and (ii) how to scale up these basic unit processes into feasible structural assemblies of 1 to 1000 atoms to produce the desired functionality. Here, we focused on the first part by surveying the single-step dynamics of graphene dopants, primarily phosphorus, caused by electron irradiation both in experiment and simulation, and developed a theoretical scheme for describing the probabilities of competing configurational outcomes through the postcollisional momentum vector of the PKA. However, a brief description of the second part is also warranted.

What one would want is to arrive at a predesigned configurational state *i*_{final} ≡ {*r** _{n}*} of the atoms as quickly as possible, through a series of collisions with focused electrons, which are known to have enough energy to displace atoms in the radiation damage context (

*35*) but exploited here in a controlled fashion to bias the configurational evolution, some of which may conserve mass locally and some of which may not. We start with an initial configurational state

*i*

_{initial}that is precisely imaged and wish to travel across intermediate configurations … →

*i*→

*i*′ →

*i″*→ v… and finally arrive at

*i*

_{final}, similar to playing Rubik’s cube but with probabilities. One obviously must balance “risk” against “speed” in playing this game, since there could exist trap states {

*i*

_{trap}} that severely delay the arrival to

*i*

_{final}or even make achieving

*i*

_{final}impossible (for example, fourfold-coordinated P is a trap state with a 60-keV e-beam, since we found that it is no longer possible to further alter the configuration once it is reached). Through the PKS formalism, we see that ideally we can affect

*i*through the following control variables: (a) choosing the PKA atom and the e-beam spot center, (b) choosing the FWHM of the e-beam and the beam dose, (c) choosing

*36*) and artificial intelligence to understanding the unit processes, as well as the assembly process, in the future. However, first-principles theory has, at this stage, been demonstrated to be tremendously helpful.

Specifically, here, we have categorized four types of electron-induced dynamics of atomic dopants on graphene and constructed a scheme for controlling them. By explaining the mechanisms for each process by first-principles calculations, we provide a convenient categorization for generic dopant dynamics. We have demonstrated the possibility of electron-beam manipulation of P and selectively induced directional SW transitions of Si. A vector-space theory (PKS) is proposed for calculating the relative ratio of scattering cross sections between different configurational outcomes, corresponding to branching probabilities in a decision tree. The two main ingredients of this theory, the outcome functions and the momentum-resolved differential cross sections (see Materials and Methods), are assessed and numerically computed by abMD and analytical relativistic electron collision kinematics, respectively. A “Doppler amplification effect” is discussed whereby a small change in precollisional PKA momentum results in a larger change in PKS momentum due to momentum conservation.

The PKS theory is developed on the basis of the fortunate separations of time scales of relativistic electron collision (τ_{c} ~ 10^{−22} s), electronic excitation (a penetration event, τ_{e} ~ 10^{−18} s), collective postpenetration electronic relaxation (τ_{E} ~ 10^{−15} s), ionic trajectory (τ_{I} ~ 10^{−12} s) leading to configurational change, and the frequency of these penetrations (τ_{p} ~ 10^{−9} s). While momenta will eventually be dephased after τ_{I}, what the PKA momentum vector was and how this vector may evolve up to τ_{I} are essential, so this is a truly dynamical problem. Up to τ_{c}, we have a relativistic collision problem that is PKA nucleus dependent but crystal independent. The physics revealed and the computational/analytical framework developed in this paper are general and can further help develop techniques for controlling single-atom dynamics in 3D materials (*23*), and ultimately, upscaling manipulations of multiple atoms to assemble 1 to 1000 atoms with high speed and efficacy.

## MATERIALS AND METHODS

### Chemical vapor deposition

Phosphorus-doped graphene was synthesized using chemical vapor deposition. First, a 25-μm-thick Cu foil (Alfa Aesar, no. 13382) was washed in 5% HCl solution for 3 min and then rinsed in deionized water several times. After that, the Cu foil was dried by nitrogen and quickly loaded into a tube reactor (diameter, 1 inch; length, 1.5 m). A quartz boat container with about 100 mg of triphenylphosphine (C_{18}H_{15}P; Sigma-Aldrich) used as a sole precursor source was placed upstream from the sample. The system was evacuated to a vacuum lower than 1 × 10^{−3} Pa. The zone-1 of the furnace was first heated to 1050°C at a rate of 20°C/min in 25 and 100 sccm of H_{2} and Ar, respectively. After annealing for 20 min, the temperature was decreased to 1000°C. Then, zone-2 of the furnace was heated to 80°C at a rate of 5°C/min. The triphenylphosphine vapor was carried into zone-1 by the flowing H_{2} and Ar, initiating graphene growth on the Cu foil. After 20 min, the system was cooled to room temperature with a cooling rate of about 50°C/min by opening the furnace. During growth and cooling, the flux of H_{2} and Ar remained unchanged. The resulting P-doped graphene was then transferred onto Quantifoil Au TEM grids for electron microscope imaging.

### Electron microscopy

The atomic structure of the sample was observed using the aberration-corrected Nion UltraSTEM 100 at the Oak Ridge National Laboratory’s Center for Nanophase Materials Sciences User Facility operated at 60 kV, and the same instrument at the University of Vienna. In a standard treatment, the sample was baked in vacuum at 160°C for 8 hours before insertion into the microscope column. The vacuum level at the sample volume during the experiments was <3 × 10^{−9} mbar. The convergence semiangle of electron probe is 30 mrad, and the collection semiangle of the electron energy loss spectrum is 48 mrad. The medium angle annular dark field (MAADF) collection semiangle range is 54 to 200 mrad. Electron current was kept in between 50 and 60 pA during all of the imaging. To estimate the doses of the spot irradiation manipulation experiments, we used a model of the expected probe intensity profile as described in (*2*, *22*).

### Postprocessing of images

The postprocessing of images was done to enhance their contrast and clarity. Figures 1 (A and B) and 3 (B and D) were processed in ImageJ using a median filter with a kernel radius of 2 pixels. Figure 1 (C and D) was binned by 0.5 (ImageJ-Scale) and a slight Gaussian blur (radius = 0.5), whereas Fig. 1 (E and F) are raw figures averaged by aligning and stacking several frames on top of each other. The insets of Figs. 1 (E and F) and 4H were processed using a double Gaussian filter plugin in ImageJ with parameters sigma1 = 0.21, sigma2 = 0.15, and weight = 0.3 (*12*) and then colored using lookup table “cyan.” Apart from Fig. 1 (E and F), all other dynamics were recorded during scanning.

### First-principles calculations

abMD was performed using DFT within the generalized gradient approximation in the form of Perdew-Burke-Ernzerhof’s exchange-correlation functional (*37*). The time step was chosen to be 1 fs, after testing that this is sufficient for predicting the dynamics (we simulated the direct exchange of P using a time step of 0.1 fs, but found no difference within the precision of our calculation ±0.1 eV). All calculations were performed using the Vienna Ab Initio Simulation Package (*38*), and the plane-wave cutoff energy was chosen to be 300 eV with PREC set to Normal?.

### PKS formalism

Every point in PKS describes the momentum status of the PKA in terms of its momentum direction and kinetic energy immediately after collision, which can be identified by a triplet Γ ≡ (θ, φ, *E*). Similarly, the energy-momentum triplet of a precollision electron (*t* = 0^{−}) will be denoted by *t* = 0^{−}) will be denoted by *d*Γ = *E*^{2}*d*Ω*dE*, where *d*Ω is the solid angle of the postcollisional PKA momentum direction and has a unit of eV^{3} despite conveying momentum vector-space information (one can think of Γ-space as a transformed momentum space with easy-to-read labels in electron volts).

The PKS framework involves a two-step process: (i) electron scattering from the nuclear potential of the PKA, denoted by “electron → PKA” (a 0.1 zeptosecond–time scale interaction, τ_{c} ~10^{−22} s) and described by function *Q*, the PKA momentum resolved electron differential scattering cross section; (ii) the ensuing dynamics of the PKA, denoted by “PKA → configurational change” (a picosecond–time scale interaction, τ_{I} ~10^{−12} s) described by *P _{i}*, the probability that outcome

*i*will take place. For every energy-momentum triplet Γ in PKS, the outcome functions

*P*(Γ) describe the probability that such a scattering event leads to an outcome configuration of unchanged (

_{i}*i*= 0), direct exchange (

*i*= 1), SW transition (

*i*= 2), knockout (

*i*= 3), etc., which is crystal structure dependent, and with 0 ≤

*P*(Γ) ≤ 1, ∑

_{i}*(Γ) = 1. Thermal and quantum perturbations of the surrounding crystal structure can smear the branching rates and make*

_{i}P_{i}*P*(Γ) neither 1 nor 0, but because

_{i}*E*has a much larger magnitude than such surrounding fluctuations, there tends to be a dominant outcome

*c*(Γ) ≡ arg max

*(Γ) for every Γ (*

_{i}P_{i}*c*stands for “configurational outcome” denoted by different colors). For example, if direct exchange is the most probable outcome at Γ, then

*c*(Γ) = 1; if SW transition dominates at Γ, then

*c*(Γ) = 2, etc. We used

*c*(Γ) to partition the PKS into different color blocks in the 3D visualization scheme shown in Fig. 4C (we use the color blue for

*i*= 1, magenta for

*i*= 2, etc.). In addition,

*c*(Γ) = 0 for regions where recovering to the same configuration is the dominant outcome. Different total cross sections of dynamic processes can be calculated considering the following two consecutive processes.

### “Electron → PKA” process

We introduce an intermediate function ^{3}, to describe the probability that a single penetrating electron can eject the PKA into a particular differential PKS volume dΓ (units of eV^{3}) by impinging on the corresponding impact-parameter differential area *d*σ = *Q* is essentially a probability density distribution, partly due to the impact-parameter dependence of the electron-PKA collision and the probabilistic nature of *32*, *39*). *Q* can be computed as*t* = 0^{−}) (*32*, *40*) and *33*, *34*), which describes the scattering probabilities of PKA with respect to its outgoing angles and energy. The δ function in energy is the result of energy-momentum conservation and is independent of the details of the nuclear potential. The function *f* defines the energy contour of PKA with respect to the outgoing angles θ, φ given the status of incident electron *f* as shown in eq. S10 in the Supplementary Materials. The

### “PKA → configurational change” process

The total cross section of a dynamic process *i* can then be computed by integrating *Q* in Eq. 1 weighted by the outcome function *P _{i}* over the whole PKS

*d*

^{3}Γ ≡

*E*

^{2}sin θ

*dEd*θ

*d*φ is the PKS differential volume element for postcollisional PKA. The cross sections of different dynamic processes are functions of

In computer-controlled atomic engineering, in evaluating Eq. 3, although *P _{i}*(Γ), however, is crystal and material dependent, and needs to be precomputed with expensive ab initio calculations, and tabulated or machine learned (

*36*) for efficient evaluation of Eq. 3.

For simplicity, in the graphical illustrations in the main text, the “PKA → configurational change” dynamics are assumed to be deterministic, making *P _{i}*(Γ) either 0 or 1, without any smearing at the boundaries. This is reflected in Fig. 4C as the sharp boundaries of the PKS regions, where the probability of configurational outcome

*i*is 1 within the boundary and is 0 everywhere else. On the other hand, the contour of

*, can thus be visualized easily in this fictitious limit of no thermal or quantum uncertainties: The intersection areas between the ovoid and the*

_{i}*c*(Γ) =

*i*regions represent the part of PKS space that can induce certain configurational change

*i*, which is then integrated with

*d*σ/

*d*Ω to get the total cross section for each of them.

To complicate the picture slightly, however, for a quantitative description of the outcomes, it has been shown that the precollisional momentum *32*, *39*), due to what we may conceptualize as a “Doppler amplification effect” on Γ. To illustrate this with an approximate example (see section S6 for details), the outgoing velocity, v, of a PKA with precollisional vibrational velocity, *E* due to the second term _{E}, the PKS momentum distribution *P _{i}*(Γ), a crystal-dependent quantity that one can precompute with abMD that integrates the evolution of atom trajectories on the ground-state BO surface (since we are beyond τ

_{E}). The overlap of

*P*(Γ) (crystal structure–dependent transition probability) in PKS space gives the net rate of configuational change (→

_{i}*i*), after which the correlated atomic momenta dephase and the momenta correlation information is lost, leaving only heat (which subsequently conducts away or radiates out with yet another, longer timescale). All these happen (or not) long before the next electron penetrates the system, and the system configuration evolves (

*i*→

*i*′ →

*i*″ → …) without carrying the detailed phase information about atomic momenta, so an uncorrelated probability distribution function of PKA momentum

## SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/5/5/eaav2252/DC1

Section S1. Overview micrographs

Section S2. Energy transfer from a 60-keV electron to a moving carbon atom

Section S3. EELS characterization of P and Al dopants

Section S4. Comparison of cNEB curves of various elements

Section S5. Primary knock-on space

Section S6. Ovoid modification by atom vibration (Doppler amplification effect)

Section S7. Atomic engineering: Manipulation decision tree

Section S8. Method of calculating experimental cross section

Fig. S1. STEM image of P-doped graphene.

Fig. S2. Energy transfer to vibrating carbon atom.

Fig. S3. EELS of P and Al dopant.

Fig. S4. cNEB curves.

Fig. S5. Construction of PKS.

Fig. S6. Selective dynamics by tilting beam.

Fig. S7. Ovoid modification by vibration.

Fig. S8. Decision tree for atomic engineering.

Fig. S9. Selective dynamics from 55-77 structure back to honeycomb.

This is an open-access article distributed under the terms of the Creative Commons Attribution license, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

## REFERENCES AND NOTES

**Acknowledgments:**

**Funding:**J.L. and C.S. acknowledge support by NSF ECCS-1610806 and U.S. DOE Office of Nuclear Energy’s NEUP Program under grant no. DE-NE0008827. J.Kon. and C.S. acknowledge the support from U.S. Army Research Office (ARO) under grant no. W911NF-18-1-0431, and through the MIT Institute for Soldier Nanotechnologies (grant no. 023674). J.Kon. and H.W. acknowledge support by NSF DMR/ECCS-1509197. T.S. and M.T. acknowledge funding by the Austrian Science Fund (FWF) via project P28322-N36, and T.S. further acknowledges funding by the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement no. 756277-ATMEN). J.Kot. was supported by the FWF via project I3181-N36 and by the Wiener Wissenschafts-, Forschungs- und Technologiefonds (WWTF) project MA14-009. The electron microscopy experiments were conducted at the Center for Nanophase Materials Sciences, which is a DOE Office of Science User Facility (to J.-C.I.). M.D. acknowledges grants from the Danish Council for Independent Research, AUFF NOVA project from Aarhus Universitets Forskningsfond, and EU H2020 RISE 2016-MNR4SCell project. Q.-B.Y. and G.S. are supported, in part, by the National Key R&D Program of China (grant no. 2018YFA0305800), the NSFC (grant no.14474279), and the Strategic Priority Research Program of the Chinese Academy of Sciences (grant nos. XDB07010100 and XDB28000000). L.B. acknowledges support by Escuela Politécnica Nacional (EPN) via project PIJ-15-09.

**Author contributions:**C.S., T.S., and J.L. conceived the idea. J.-C.I., C.S., M.T., C.H., and T.S. participated in the electron microscopy experiment. Z.W. synthesized the doped graphene using the chemical vapor deposition method. C.S., H.W., C.H., and M.T. prepared the TEM samples. C.S. performed the ab initio MD calculations. Q.-B.Y., Z.Z., and T.S. performed cNEB calculations. C.S., M.T., and C.H. performed the beam dose analysis for each event. C.S. wrote the MATLAB code for visualizing the PKS and calculating the cross sections. L.B. performed the STEM image simulation. C.S. performed the EELS simulation. J.L., T.S., J.-C.I., J.K., J.K., J.C.M., M.D., and G.S. supervised the whole project. C.S., T.S., and J.L. drafted the manuscript, and all authors participated in its revision.

**Competing interests:**The authors declare that they have no competing interests.

**Data and materials availability:**All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

- Copyright © 2019 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution License 4.0 (CC BY).