Chemical sputtering from amorphous carbon under bombardment by deuterium atoms and molecules

We perform classical molecular dynamics simulations of the chemical sputtering of deuterated amorphous carbon surfaces by D and D2, at energies of 7.5–30 eV D−1. Particular attention is paid to the preparation of the target surfaces for varying impact projectile fluence, energy and species, to the vibrational state of D2 projectiles, as well as to the variation in sputtering yields with target surface and impact projectile. The methane and acetylene sputtering yields per deuteron, obtained with atomic and molecular projectiles, agree quantitatively with recent experimental values. We study the distribution of sputtered species, as well as their kinetic energy and angular spectra.


Introduction
Because plasma boundary physics encompasses some of the most important unresolved issues in future energy production, for both the International Thermonuclear Experimental Reactor (ITER) project and other fusion reactors, there is a strong interest in the fusion community for better understanding and characterization of plasma-wall interactions. Chemical and physical sputtering cause the erosion of the limiter/divertor plates and vacuum vessel walls, whether made of C, Be, W, or mixed materials, and degrades fusion performance by diluting the fusion fuel and excessively cooling the core. Hydrocarbon redeposition onto plasma-facing components can lead to long-term accumulation of large in-vessel tritium inventories. Although these processes are among the key issues for ITER development, they are currently poorly understood. Moreover, obtaining results from tokamak experiments in this area is stymied by the impossibility of direct measurement of the underlying processes. Fortunately, increased computational power, as well as improved understanding of atomic and plasma physics near the surfaces, provide conditions for a more robust computer simulation response to this challenge. In addition, progress in ion beam-surface experiments has provided much needed data under controlled conditions [1]- [7]. Our current work is devoted to the least studied, low-energy regime of chemical sputtering of amorphous carbon, 7.5-30 eV D −1 , by both deuterium atoms and molecules, focusing on the issues of target surface preparation and impact particle state, generation of a steady-state surface, analysis of the distributions of species, energies and angles of the sputtering ejecta, and understanding of the underlying physics. Use of deuterium projectiles, compatible with the d-t plasma environment, allows direct comparison of methane and acetylene sputtering yields with recent experimental results [5]- [7]. Apart from isotopic effects [3] (which are not a subject of the present study), the chemical sputtering process is expected to be equivalent to that with 1 H projectiles.
Chemical sputtering is a process where bombardment by ions, atoms or molecules induces a chemical reaction, producing a particle that is weakly bound to the surface and hence easily ejected or desorbed into the gas phase, as schematically shown in figure 1. Though many open questions on the nature of the chemical sputtering remain, some basic mechanisms are well understood for both high and low temperatures of the target. At higher temperatures, impacting ions weaken the bonds at the surface leading to thermal desorption [8]. However, at low temperatures of the carbon surface, it is hypothesized that incident ions break C-C bonds within the collision cascade, creating dangling bonds which are rapidly passivated by the hydrogen (deuterium, tritium) that is abundant in the fusion plasma environment [2]. This leads to the formation of hydrocarbon molecules underneath the surface, which desorb from the surface and are transported within the plasma, possibly dissociated and ionized, and re-deposited on plasma-facing materials. This mechanism of breaking the C-C bond is known as 'swift chemical sputtering' [9].
Molecular dynamics (MD) simulations are based on classical many-body dynamics using empirical potentials. When formation of a hydrocarbon is realized many atomic layers beneath the solid-plasma interface [8], the sputtering process is delayed by a period of time (milliseconds), defined by the time needed for the particle to diffuse to and desorb from the surface. Such a long time cannot be described within fully atomistic MD simulations (which typically can be performed up to ns). The low incident energies (and thus shallow penetration depths) considered here exclude the need for such extended ms simulations. Although the available time is not sufficient to model slow diffusive processes, the sputtered particles find their way to the surface 3 DEUTSCHE PHYSIKALISCHE GESELLSCHAFT Figure 1. Schematics of the chemical sputtering process with D and D 2 projectiles at a deuterated amorphous carbon surface; periodic boundary conditions were applied in the x-and y-directions. through the voids created by bombardment where they are desorbed kinetically, i.e. carrying a signature of a local temperature, which might be elevated by the impact cascade from the average surface temperature. We find that 5 ps after an impact allows for the ejection of the majority of hydrocarbons, enabling a complete study of chemical sputtering within the MD ps time frame. Moreover, we provide a scheme for preparing deuterated surfaces by cumulative bombardment which can be carried out within several nanoseconds.
Although carbon polytypes of interest for fusion range from crystalline to polycrystalline to amorphous, the target surface in our MD simulations is deuterated amorphous carbon (a-C : D). This is representative of the eroded and redeposited material on components facing a dominantly deuterium fusion plasma. The high flux (projected to reach up to 10 25 m −2 s −1 in ITER) of low energy plasma particles (primarily eV to tens of eV, extending to keV in edgelocalized mode bursts (ELM)) transforms any initial carbon surface into an amorphous one, after sufficient fluence. Simultaneously, carbon films redeposited from the plasma at the fusion reactor walls, originating from the erosion of carbon surfaces, are also amorphous. Finally, the deuterium-rich environment of the fusion plasma ensures a thorough deuterization of the amorphous carbon surface. The D/C ratio can reach 40% on average [9], although locally (in the surface layers) this value might be significantly elevated.
We compare our MD simulation results using D and D 2 impacts on a-C : D with the recently published experimental results [5]- [7] for sputtering of ATJ graphite surface by lowenergy beams of D + , D + 2 and D + 3 ions [7], which allows a more detailed benchmark of MD simulations. Comparison of simulated bombardment of a-C : D by neutrals with experimental data for bombardment of ATJ graphite by ions is justified because the sample is expected to become amorphous after significant fluence, and because the projectiles are predominantly neutralized 4 DEUTSCHE PHYSIKALISCHE GESELLSCHAFT before hitting the surface. These issues are addressed below. In order to mimic the experimental conditions as closely as possible, we pay particular attention to the preparation of the target surfaces, as discussed in section 2. This turns out to be a crucial parameter in determining the distributions of sputtered species and sputter yields. The detailed sputtering results are presented in section 3. Our conclusions are given in section 4.

Preparation of the target surface
The basic elements of an MD simulation of sputtering are the target surface microstructure, the interatomic potentials and the state of the impact species, all of which crucially influence the sputtering results. All simulations in this work were performed with the reactive empirical bond order (REBO) potential [12]. As a member of the classical bond-order family of potentials [13] (Tersoff-Brenner potentials) it provides a good, empirical description of the covalent bond for nonpolar systems, without incurring the enormous expense of a quantum mechanical calculation. We note that the most recent version of this potential [12] contains terms that specifically treat contributions to the bond order in radical species, and was fit directly to vibrational frequencies of molecular hydrocarbons and the elastic constants of bulk carbon-based materials, providing significant improvements in both cases over the widely used 1990 parameterization. Thus this model is quite well suited to describe both the bulk a-C : D substrate and the small molecular species that are ejected from the sample. More advanced bond-order potentials are available that include better descriptions of the torsional and van der Waals interactions (AIREBO [14]). These enhancements are not expected to contribute significantly for covalent network solids such as a-C : D, but would be important for crystalline graphite, describing the interaction between graphite layers [15]. One significant shortcoming in these potentials that is relevant for sputtering simulations is the repulsive core potential. While the newest version of the REBO potential contains a screened Coulomb repulsive potential [12], this repulsive potential was fit to reproduce the properties of hydrocarbon systems only within a few eV of covalent bonding minima, and little attention was paid to the high-energy repulsive wall at close encounters. Consequently, the repulsive part of these potentials for C-D and D-D interactions is correct up to approximately 100 eV, while for C-C it becomes inaccurate (in comparison to the repulsive, screened-Coulomb ZBL potential [16]) already above 30 eV. These observations limit the range of MD studies to less than 60 eV D −1 and applications at higher energies should be taken with caution. Of course, this is a conservative estimate, since deuterium impacts probe the highly repulsive regions of only the C-D and D-D interactions. Some authors, however, have used the unmodified REBO repulsive potential for impacts at 100 eV [17], conditions at which the scattering cross-sections may be poorly represented in the most energetic collisions. Work is currently underway to improve the short-range repulsive REBO potential, thus making it applicable to higher impact energies. Such improvements will be the subject of forthcoming articles; for the present study, all impacts were conducted at 30 eV D −1 or less, therefore in the 'safe' energy region. Impacting particles (D or D 2 ) were introduced at a distance of more than 20 Å above the initial surface position with a translational kinetic energy (KE) between 7.5 and 30 eV D −1 directed normal to the surface (the z-axis in figure 1). The position of the centre-of-mass of the incoming particle was chosen randomly within a x-y-plane parallel to the surface, with a uniform distribution. For molecular impacts, the initial orientation was also chosen randomly with a uniform distribution. The molecular projectile was given a rotational KE corresponding to 5 DEUTSCHE PHYSIKALISCHE GESELLSCHAFT a temperature of 1000 K, with the direction of the angular momentum vector chosen randomly from the unit sphere. This rotational temperature was chosen to qualitatively reflect the elevated, though not precisely known, rotational temperatures of the D 2 projectiles in the experimental work by Vergara et al [5]- [7]. In order to prevent net motion of the centre-of-mass of the system upon repeated impacts, the atoms in the 2 Å of the substrate furthest from the impact surface (bottom of the cell in figure 1) were held rigid in the direction of impact (orthogonal to the surface). These atoms were approximately 25 Å below the impact surface (figure 1).
The structure of amorphous carbon is strongly dependent on the method of preparation. In the hydrogen-rich environment of a fusion reactor the amorphous carbon is highly hydrogenated, leading to hydrocarbon sputtering from the carbon surfaces. For bulk material, the average H/C ratio can be as large as 0.4 in amorphous carbon surfaces at 300 K, decreasing with increasing temperature, and increasing under surface bombardment by energetic hydrogen (supersaturation) [9]. The depth of the hydrogenated layer can extend to hundreds of nm, depending on the incident energy of the hydrogen and the time available for the hydrogen diffusion away from the surface [18]. In MD simulations, this depth is limited to the size of the simulation cell and substantial diffusion cannot occur on the ns time scales of the simulations. Thus, our studies begin from an a-C : D substrate with a bulk-like D/C ratio of 0.4, while we rely on surface preparation by energetic impacts at high fluence to generate the correct supersaturation of D.
The surface preparation was performed in several phases. The initial sample consisted of a hydrogenated amorphous carbon sample with a density of 2.0 g cm −3 and hydrogen/carbon ratio of 0.4 (700 deuterium atoms and 1750 carbon atoms). This sample was prepared by heating a pure carbon (diamond) bulk sample to 10 4 K, allowing it to remain molten for several tens of ps, and quenching it back down to 300 K at a rate of 500 K ps −1 . Deuterium was then introduced into the simulation cell by substitution, at which point the system was again heated to a temperature of 10 4 K for 10 ps and requenched to 300 K over 62 ps. After this the periodic boundary conditions were removed from the simulation cell in the z-direction (see figure 1), generating a two-dimensional periodic slab. This surface was then relaxed at 300 K for 100 ps. In this way we generated an initial, 'virgin' a-C : D surface that was the starting point for preparation of individual surfaces by D and D 2 bombardment at different impact energies, for use in the sputtering studies. This resulting sample was a cube of approximately 26.5 Å in each direction.
Surface preparation was then performed by bombarding the 'virgin' surface with sequential impacts by projectiles of a particular energy, particle type (D or D 2 ), direction, and rovibrational state, with up to 2000 cumulative impacts on the surface. Thus, our calculation modeled an experiment with a well defined particle beam at a surface, at a fluence of particles large enough to result in substantial surface modification. The goal was to reach the same steady-state surface modification that is achieved at high fluences of > 10 22 D m −2 in the experiments. While we were only able to reach fluences of 2.8 × 10 20 D m −2 , we show below that this was sufficient to reach a steady state. It is not possible to reach a permanent steady state in simulations, because of the limited number of particles in the target; the D content will rise monotonically and the C content will decrease monotonically as the surface is progressively etched away by D impacts. But a steady out-flux of sputtered particles and steady structure of the subsurface region were reached at least temporarily after an initial transient period. The sputtering simulations were short enough to prevent erosion of more than 10% of the initial sample at highest energies of the considered range.

DEUTSCHE PHYSIKALISCHE GESELLSCHAFT
The fluence at which this steady state surface is reached and how long it is preserved depends on the impact conditions and the type of particle being sputtered. We find that for the sputtering of methane in the considered range of energies, surfaces prepared with a fluence of 200-700 impacting D 2 particles per cell, or approximately twice as many impacts for D (0.6-2 × 10 20 D m −2 ) reflect a temporary steady state. From the simulations of sequential impacts, we selected six surfaces in this range of fluences for each set of the impact particle initial conditions. All of the sputtering simulations described in the following section were done by averaging the results over these six surfaces, using the same energy, impacting particle and its vibrational state, if applicable, as were used to generate the surfaces. This is a key factor in the simulations presented here. All sputter yields were calculated using surfaces prepared at the same impact conditions as used to calculate the sputter yield. This is fully consistent with the experimental procedures, in both beam and reactor studies, but has been ignored in most previous simulation approaches [7,9,17]. Alman and Ruzic [19] have recognized the importance of the surface creation by cumulative bombardment, but in this procedure they used only one impact species (H) of one energy (20 eV) to prepare the surface for sputtering by C and various hydrocarbons in a range of energies. Träskelin et al [20] have significantly improved the surface preparation, using bombardment of the hydrogenated carbon target by H of the same impact energy as used in the sputtering simulation by mixed impact beam (H + an inert gas atom).
One consequence of the accumulation of deuterium at the surface is that it provides a possible pathway for neutralization of incident D + 2 ions, which are the typical projectiles in the beam-surface experiments and also relevant for plasma-facing surfaces. We have performed sputtering studies with various vibrational excitations of the D 2 projectiles, reported elsewhere in more detail [18]. Comparison of these results with experiment [5]- [7] serves as additional evidence for the dominant neutralization process of the incoming D + 2 above the surface. These studies supported the view that the leading mechanism of neutralization is resonant charge exchange between the incident D + 2 and a D atom in the hydrogenated surface, resulting in a vibrationally excited state of D 2 . Figure 2 illustrates this argument for the example of charge transfer in the collision of H + 2 and H. Figure 2(a) shows the diagram of the vibrational levels of initial (ground H + 2 (v = 0) + H(ls)) and final (H 2 (v f ) + H + ) states.
Here v f = 4 of H 2 is quasi-resonant with the (v = 0) ground state H + 2 (+H). As a consequence, the cross section for charge transfer in binary collisions to v f = 4, with an excitation energy of about 2 eV, is about two orders of magnitude larger than the cross-section for transitions to other neighbouring states, as can be seen in figure 2(b) [18]. For D + 2 projectiles, the interatomic potential is unchanged, but the higher nuclear mass compresses the vibrational spectrum, so that near-resonant charge transfer generates D 2 molecules primarily in the 5th and 6th excited states. These neutral molecules then impact the surface in that vibrationally excited state. Because their dissociation energy is thereby lowered by 2 eV, less of the impact energy is tied up in D 2 dissociation, leading to higher sputtering yields than with ground-state projectiles [22]. The situation is similar if the incident D + 2 projectiles are vibrationally excited, as might be the case with the projectiles originating from an ion beam [5]- [7]. As can be anticipated from figure 2(a), these would result in quasi-resonant charge transfer to D 2 molecules in states higher than the 5th and 6th vibrational state, increasing low energy sputtering yields, which would actually improve the agreement between simulation [22] and experiments [7] in the low-energy limit. An important outcome of these studies is that the rovibrational state has to be taken into account whenever it is likely that slowly impinging molecules are excited due to resonant charge transfer. Such conditions might be also met in a fusion divertor plasma at energies of several eV. Consequently we perform our surface preparation and sputtering simulations with vibrationally excited D 2 , at an energy of about 2 eV (i.e. with an initial classical bond length of λ = 1.3 Å) or higher. In order to illustrate the effect of the vibrational excitation of impact molecules we will also present results for ground-state D 2 (i.e. an initial bond length of λ = 0.764 Å). To address the possibility that the impacting D + 2 originates in an excited vibrational state in the molecular beam experiments of interest [7], we also use D 2 projectiles with a larger classical bond length. It is likely that D + 2 have a distribution of vibrationally excited states v, that peaks at v = 5 or 6 [23]. The molecular ion will again resonantly neutralize to a D 2 (v) where v now reaches 10, going close to the dissociative continuum edge. Since H-H covalent bonding in the REBO potential is short-ranged (not supporting bond lengths larger than approximately 2 Å), this high vibrational state of D 2 effectively means a dissociated D-D molecule. Thus we also perform calculations (and surface preparations) for just-dissociated D 2 molecules, having initial bond lengths of 2 Å.
During the surface preparation, incoming particles were generated at a flux of 1.4 × 10 29 D m −2 s −1 , or one impact every 2 ps on the 26.5 × 26.5 Å 2 surface. In the first ps following the implantation of the incoming particle, the dynamics were integrated freely, with no external forces or thermostats other than the constraint on the bottom layer. At times between 1 and 2 ps, a Langevin thermostat was applied to the entire simulation cell with a target temperature of 300 K and a time constant of 200 fs in order to remove the excess thermal energy that resulted from the impact. The simulated flux is larger, by many orders of magnitude, than that of the typical experiments (10 18 -10 20 D m −2 s −2 ) or even that of fusion applications (10 25 D m −2 s −1 for ITER). However, most of the conversion of kinetic to chemical energy occurs during the 1 ps following impact, and the thermostatting is designed to model the slow transport of the residual thermal kinetic energy from the system, as occurs in the beam experiment at much longer time scales of ms between local impacts. After the 1 ps impact and 1 ps equilibration, another impact is initiated. The surface is thus bombarded repeatedly, progressively accumulating more damage. A total of up to 2000 ps of simulation time (or twice as much for D impacts) was accumulated for the a-C : D surface for each set of initial conditions for the impact projectile (i.e. energy, type and internal state of particle). Thus, the simulations model fluences of up to 2.8 × 10 20 D m −2 (1000 D 2 /cell), which is sufficient to reach a steady state surface comparable to that achieved after tens of minutes in the experiment, when starting from a clean graphite interface [15].
Although considerable care was taken in constructing a 'virgin' carbon surface that is hydrogenated and amorphous, the surface composition evolves under bombardment, and the 'steady state'amorphous surface that is produced under sputtering conditions differs substantially from the initial surface. This makes our simulation approach substantially different from those used by Nordlund and Salonen [7,9] and by Marian et al [17], which use only the 'virgin' a-C : D surfaces for their sputtering simulations. These 'virgin' surfaces produce predominantly carbon atoms and dimers, as well as simple hydrocarbon radicals, as described in detail by Marian et al [17]. However, allowing creation of the surfaces by a fluence of the relevant impact particles, changes the structure of the sputtered species significantly, leading to heavier hydrocarbons, including methyl and methane, as will be discussed below.
As one method of characterizing the interface that is generated by cumulative bombardment, the density profile for deuterium and carbon atoms is plotted as a function of depth in figure 3. The initial density of deuterium is 40% that of carbon, by construction. Initially, after thermal equilibration but before any sputtering, the interface between a-C : D and vacuum was located at approx. z = 13 Å. Figure 3 shows the evolution of the C and D density under impacts by 40 eV D 2 Several features characterize the evolution of this surface under increasing fluence. First, the effect of the D 2 impacts is to gradually increase the D content at the interface, while progressively etching away the carbon. This is not surprising, and indeed is inevitable, given that the incoming particles consist solely of deuterium, while any carbon that escapes can never return. For fluences greater than 200 D 2 impacts per cell the surface expands by 5 Å, with a depletion of carbon atoms and an enrichment in deuterium down to 10 Å beneath the original interface. The carbon depletion and deuterium enrichment proceeds until fluences of ∼1 × 10 20 D m −2 (400 D 2 /cell), at which point the D/C ratio in the outermost 10 Å of the surface reaches a value of approx. 1. At that point, the composition of the surface changes little with additional damage, although the surface gradually migrates into the bulk as it is etched away. The kinetically formed steady-state surface is thus quite different from the equilibrium bulk material, with a D/C ratio substantially enhanced from the bulk value of 0.4 to a value of 1 at the interface. Similarly, Salonen et al [24] has reached at the interface a supersaturation ratio of 1, by bombarding the a-C : H surface by 1 eV H projectiles, explaining decrease of the sputtering with increase of the impact flux. Because of the significantly higher impact energy of D (and D 2 ) in the current study, the supersaturated penetration depth profile is noticeably different in our study and extends down to 10 Å beneath the interface. As we will see below, this modified surface generates very different sputtering products than does the initial surface, in addition to affecting the sputter yield. The higher D/C ratio in the modified surface produces more molecular and saturated products, while the equilibrium bulk surface generates ejecta that are unsaturated and carbon-rich. The D/C ratio, averaged over the simulation cell, steadily increases with fluence, as shown in figure 4. The rate of increase is larger at higher energies, due to the deeper penetration depth. The D/C ratio never reaches true saturation, in part because of the monotonic loss of carbons from the simulation cell and finite penetration depth of the projectiles.
The spectrum of sputtered hydrocarbon molecules is a quantity of substantial interest, both in beam-surface experiments and in fusion applications. This product distribution is in turn closely related to the hybridization and coordination state of the carbon atoms at the surface. Methyl radicals and methane molecules will result primarily from surface R-CD 3 moieties; more generally, saturated molecular products will result from cleavage at sp 3 -hybridized carbon atoms.
The depth-resolved variation in hybridization state together with the penetration depth of the incoming particles can be used to understand the detailed chemical changes induced by bombardment. Figure 5  of the penetration depth profile with increasing energy. The changes in the hybridization states, averaged over the six surfaces, closely mirror the structure of the penetration profile. Thus, the increase in sp 3 density occurs at generally the same depth as the final implantation position, or slightly deeper, indicating that the new sp 3 atoms are created primarily at the end of the impact cascade, when the nearly thermalized D atoms passivate unsaturated sp 2 atoms. This is also visible in the decrease of sp 2 density. The structure of this distribution is somewhat more complicated, with an increase at very shallow depths that likely results from passivation of sp atoms. This decrease in sp density is nearly independent of impact energy, indicating that the sp-hybridized atoms, located primarily at the interface, are highly reactive in the initial stages of the collision cascade.
The modifications in C atom hybridization at the surface allow some speculations regarding the sputtering yields from these steady-state surfaces. Cleavage of the C-C bond tethering an R-CD 3 moiety to the surface leads to the formation of a CD 3 radical, which can be deuterated, leading to a sputtered methane. It is expected that CD 3 will dominate, both from direct ejections, and because the short timescales accessible in the MD simulation do not allow for relatively slow D abstraction to take place before a loosely bound CD 3 species will be jostled from the surface by a subsequent impact. Likewise, R-CD 2 surface groups can easily produce CD 2 and CD 3 radicals, among which CD 2 will likely dominate.
Thus the increase in sp 3 carbon for the impact-modified surface would predict a steady increase in saturated sputtering products with increasing fluence as the surface approaches a steady state. However, not all of the sp 3 atoms will represent tethered methyl groups, of course. Thus, further insight regarding the methyl sputtering precursors can be obtained from figure 6, which displays the actual number of R-CD 3 , R-CD 2 , and R-CD terminal groups, averaged over the steady-state surfaces, and their variation with impact energy. These are evaluated separately over half of the system, capturing nearly all modifications due to bombardment (figure 6(a)) and for only the very outermost portion of the system ( figure 6(b)). In both cases, the R-CD 3 methyl precursors are most prevalent under bombardment at about 15 eV D −1 . Consequently, if this is the main source of creation of the CD 3 radicals (and CD 4 by passivation with an additional D atom), the CD 3 + CD 4 yield should be expected to show a maximum at about 15 eV D −1 impact energy, decreasing slowly at higher impact energies. Similarly, the sputtering yield of CD radicals should follow the behaviour of the R-CD terminal groups, which become monotonically more prevalent at higher impact energies and reach saturation at energies above 20 eV D −1 . In turn, the CD 2 yield may be expected to saturate at energies above about 15 eV D −1 . Likewise, it can 12 DEUTSCHE PHYSIKALISCHE GESELLSCHAFT be predicted that the CD 3 yield and CD 2 yield should be comparable, and both considerably larger than the CD yield. Although these predictions are simplistic, they are largely confirmed by MD simulations of sputtering, as described in detail in the following section. This indicates that, although sputtering is a dynamic process that depends on following the details of complicated collision cascades, a great deal of information can be obtained by careful study of the structural properties of the surface being bombarded.

Sputtering simulations
The sputtering yield from a particular surface is evaluated using a second type of simulation. In order to obtain a statistically meaningful estimate of the sputtering yields (which are typically less than 1%), more than 2000 statistically independent impacts of D 2 (4000 of D) are performed on the same surface with the same initial conditions and the same projectile as were used to prepare the surface (which was repeated for six surfaces at various fluences, resulting in 24 000 D-impacts for each case studied). For each of these statistically independent impacts, the position, orientation and direction of angular momentum of the incident particle are chosen randomly. These simulations evolve for 5 ps in the absence of any global thermostat. Only atoms in the bottom 2 Å of the simulation cell were again constrained in the z-direction while a Langevin thermostat at a temperature of 300 K is applied to those atoms in the x-y-plane (see figure 1). Ejected molecules are removed once they are more than 10 Å from the surface. Simulations were performed for the sputtering of a-C : D by atomic D and molecular D 2 at impact energies between 7.5 and 30 eV for D projectiles, and between 15 and 60 eV for D 2 projectiles. All molecular dynamics simulations were performed with the velocity Verlet integrator and a time-step of 0.2 fs. Figure 7 shows the species-specific sputtering yields for a number of different carboncontaining compounds, on a series of progressively damaged surfaces generated with impacts by 40 eV D 2 molecules (i.e. 20 eV D −1 ). The total yield of all carbon-containing ejecta is also presented. There are several noteworthy features in these species-specific sputtering yields, and their variation with increasing fluence. The total sputtering yield shows some variation with accumulation of damage to the surface, but this variation is relatively weak, and consistent with the statistical fluctuations. The species-specific sputtering yields, on the other hand, vary quite dramatically from the initial surface to the more damaged surfaces. These variations appear to be somewhat systematic, with atomic (C) and smaller unsaturated radical species ejected much more frequently from the undamaged amorphous surface. For example, the most likely species ejected from the undamaged surface is atomic carbon, with a sputtering yield of nearly 1%. Aside from acetylene, the only other species that are sputtered from the undamaged surface are the radical species CD and C 2 D. These highly unsaturated, radical species are typical of compounds that have been observed to be sputtered from carefully prepared, but not-bombarded amorphous carbon surfaces in previous studies [7,9,17,25].
After the surface has accumulated some damage, however, the distribution of ejected species takes on a very different character. Larger molecular species begin to dominate, with the largest contributions coming from acetylene, and from methane and methyl radicals. Complex polyatomic hydrocarbons begin to appear, with C 2 D x species appearing in significant amounts at fluences above 0.2 × 10 20 D m −2 and a substantial signal from C 3  sputtering yield from the 'virgin' surface. Likewise, the small atomic and radical species that were significantly sputtered from the undamaged surface gradually decrease in importance as the original surface becomes modified with increasing fluence. The most reactive of these species, such as C atoms and C 2 D radicals, contribute negligibly to the total carbon sputtering yield for large fluences, although yields from some unsaturated radical species, such as CD and CD 2 do persist. There is also an interesting transient period, at fluences of 0.2-1.0 × 10 20 D m −2 , during which the sputtering yields of some unsaturated species, such as acetylene and CD 2 , reach a peak before decaying to lower levels. Thus it appears that the surface has been 'worn in' by the sputtering process by the time a fluence close to 1 × 10 20 D m −2 has been achieved (roughly one impact per square angstrom), and appears to have reached a steady state. This does not guarantee that the sputtering yields (either total or species-specific) will not exhibit statistical variations with increasing fluence, however. The ejecta for independent impacts on a single surface at a single value of the fluence are more highly correlated than those from surfaces at different fluences, due to the presence of a limited number of surface features in the 703 Å 2 surface. The presence of a single, weakly bound −C 2 D 5 surface moiety could substantially enhance the ethyl radical sputtering yield for one particular surface, for example. The true 'steady state' surface is the ensemble of all such surfaces generated after the initial, prepared surface has stopped evolving systematically. Based on the evolution of the species-specific sputtering yields, it appears that the surface begins to sample configurations from this 'steady state' surface for fluences of >10 20 m −2 . In the MD simulation, different surfaces from this ensemble are sampled temporally, as the surface is progressively damaged at higher fluences. In the experiment, of course, the coverage of the ion beam is many orders of magnitude larger, of order of mm 2 , and there is also a large degree of spatial averaging over these different microstructural features. Thus, in order to make an effective comparison between experimental and simulated sputtering yields, it is necessary to perform an average over the simulated sputter yields, and to eliminate the contribution from the initial, transient 'wearing in' period.
For this reason, the sputtering yields, discussed below, were averaged over six surfaces obtained for fluences between 400-1400 per cell area (between 400 and 1400 ps of cumulative D 2 bombardment). In the following we focus on these averages and do not specifically address fluence-dependent properties. Six different surfaces are sampled at each set of conditions for the impinging particle (energy, type and internal state). Thus, the sputtering yield for each set of impact conditions was obtained from the impact of N = 24 000 deuterium atoms. By sampling over different surfaces, as well as over different impact sites on the surface, sputtering yields can be obtained with low statistical error. Although the sputter yields from each of the N = 24 000 D impacts are not rigorously independent (due to the correlation of ejecta with structure for each individual surface) we calculate the statistical error for the yield n/N, where n particles of a particular type are sputtered, as ε = √ n/N, which is correct in the limit where each impact is truly independent.
The variation of the calculated sputtering yields with energy is shown in figure 8 for various species, for impact of D ( * ) 2 , where ( * ) refers to a vibrationally excited state with bond length of 1.3 Å, i.e. 2 eV below the vibrational continuum (as compared to 4.5 eV dissociation energy from the ground vibrational state, of bond length 0.734 Å). At all energies, the sputtered species are consistent with the arguments presented above for impacts at 20 eV D −1 . Noteworthy trends include the increasing prevalence of acetylene and CD ejecta with increasing energy (note that CD may be a byproduct of the decomposition of energetic acetylene). There is a corresponding decrease of C 2 D x species that does not quite rise to the level of statistical significance. Remarkably, the behaviour of the yields of CD, CD 2 , and CD 3 appear to be correlated with the behaviour of the R-CD, R-CD 2 , and R-CD 3 moieties in figure 6 (as predicted in the structural analysis of the bombarded surface), with the methyl radical sputtering yields exhibiting a distinct peak at 15 eV. Thus, the terminal groups (R-CD x ) bond-breaking (TGBB) mechanism emerges as the leading subgroup in the well-known swift C-C bond-breaking [9,25], associated to the chemical erosion in hydrogenated carbon.
The importance of preparing and sputtering the surface by with mutually consistent types and energies of particles in both stages, as described in section 2, is shown in figure 9 (solid lines), for the case of 'methane' i.e. CD 3 + CD 4 sputtering. The dashed lines show the sputtering results for surfaces prepared with mutually consistent impact energies, but with one type of impact particles, D ( * ) 2 , that were not consistent with the particles used to calculate the sputter yield. The differences are biggest at lowest energies, reaching a factor of 2-4, while each pair of curves converges (within the statistical uncertainty) at higher energies (>15 eV). The 'self-consistent' results with D ( * ) 2 sputtering of the surfaces prepared by D ( * ) 2 impacts (black solid line) will be used as an example in a detailed analysis of sputtering spectra. In order to compare our results with the beam experiments measuring methane yield upon impact of D + , D + 2 and D + 3 , we perform also self-consistent calculations with impacting D 3 (dissociating, with initial bond length of 2 Å; see [22] for details), as shown in figure 10(a). Interestingly, the experimental results for the three different impact particles obtained at the same energy per D diverge at low energies (with D + giving the lowest, and D + 3 the highest yield). The MD simulation results in figure 10(a) show the Figure 10. (a) Comparison of the experimental data for methane sputtering by D + [7], D + 2 [4,7] and D + 3 [7] with calculated sputtering yields of CD 3 + CD 4 by impact of D atoms and by dissociating D 2 and D 3 molecules, for 'self-consistently' prepared surfaces. (b) As in (a) for acetylene (C 2 D 2 ) yields, experimental data from [5]. same trend and order. We note that these qualitative results cannot be obtained if the surfaces are not prepared in a fully self-consistent way. Figure 10(b) shows a comparison of experimental and MD simulation results for the sputtering of acetylene. A divergence of the results of simulation by various projectiles is present here, as in the case of methane, though to a smaller extent. Again, as with methane sputtering, the slope of the MD sputtering yield curve which corresponds to dissociating D 2 impact is closest to the experimental one, indicating vibrationally highly excited molecular ions in the beam, together with the neutralization mechanism discussed in the section 2.
Because of the significant difference in results when the sputter yields are not calculated on surfaces prepared self-consistently with the impacts under the same conditions, we will consider only such self-consistently prepared surfaces in what follows. A total of more than 150 such surfaces were prepared (i.e. six surfaces sampled from the steady state at different fluences for each of five different impact particles at each of five different energies).
Our calculations for (CD 3 + CD 4 ) are compared with recent measurements of the sputtering yield of methane by Vergara et al [5]- [7], who were able to achieve very low impact energies DEUTSCHE PHYSIKALISCHE GESELLSCHAFT Figure 11. Total carbon yields, from experimental impact by D + 2 [4], and D + 3 [1], compared with simulations of sputtering from 'self-consistent' surfaces. The dashed line shows the total calculated hydrocarbon yield by D ( * ) 2 impact. The MD simulation results from [7] by D 2 impact of a 'virgin' surface are also shown.
for molecular ions D + 2 and D + 3 , and also for D + ions. For comparison with this experiment, the methane (CD 4 ) and methyl radical (CD 3 ) yields from the MD simulations have been combined. The reasons are some uncertainties in both experiment and simulations. On one hand, the sputtered and reflected particles coming from the target surface are reflected at least once from the stainless steel walls or baffle in the experimental geometry [5]. These walls are likely coated with deuterium and deuterated hydrocarbons, remaining from earlier experimental runs [5]. As a consequence, radicals such as CD 3 could be passivated by wall collisions between the sputtering event and the mass spectrometric detection. On the other hand, even though we use the 2002 version of REBO potential [12] which includes terms in the bond order that are designed to give better description of radical energetics, inaccuracies in the relative energetics of these species could affect the passivation of the free CD 3 radicals-resulting in a poor CD 4 /CD 3 ratio. Thus, the discrepancy between the apparent dominance of methane sputtering products in the experiments versus the dominance of methyl in the simulations remains an open question.
In spite of these uncertainties, the agreement between the simulated and experimental methane yields is remarkable, and comparable to the statistical error bars of the respective measurements, particularly for the impacts by vibrationally excited D 2 . In order to present the dispersion of the results over the six surfaces used at each energy, for each particle and at each internal (i.e. vibrational) state of the molecule (when applicable), we show in figures 9 and 10, the error bars corresponding to one standard error (standard deviation σ normalized to the sample size n s = 6, i.e. σ/ √ n s ). The standard error is shown also for D (g) 2 in figure 11. The total yield of all carbon-containing species is shown in figure 11 and is directly related to the erosion of the surface and thus is of significant interest in fusion applications. The total yield rises smoothly and nearly monotonically with impact energy. While variation of the yield with energy over the energy range of 7.5 to 30 eV D −1 is small with dissociating D 2 or atomic D (less than a factor of 7), the variation with energy for D ( * ) 2 impact is a factor of 30 and still ten times larger with D (g) 2 . This can be attributed to the energy needed to dissociate the projectile, since 'geometric' size of the non-dissociated particle is not an advantage for its penetration into the surface. Figure 11 also shows experimental results from [1] and [4] which are in reasonable agreement with the calculations. We have performed a detailed analysis of the methane and acetylene yield as functions of the vibrational state of the impacting molecule in [22], using only surfaces prepared by D ( * ) 2 . Figures 9, 10 and 11 show that calculations with 'self-consistent' surfaces further emphasize the dependence of the sputtering yields on the initial vibrational state.
We have also studied the energy distributions of sputtered particles. This information, which is much needed for modelling of fusion plasma-surface interactions, has not been available for sputtering from hydrogenated carbon, primarily due to poor statistics in the energetic and angular decomposition of a small number of ejecta, being shown for the first time in [26]. By using a statistical sample of 24 000 impacts at each set of conditions, our sputtering simulations are able to provide sufficient statistics to evaluate meaningful product distributions for some of the more prevalent sputtered particles, including protomethane (methyl + methane) and acetylene and certainly for the total population of sputtered hydrocarbons. Figure 12 shows the distribution of ejection energies, E s , averaged across the six surfaces prepared by 20 eV D −1 impacts with vibrationally excited D 2 . These are translational KEs, obtained from centre-of-mass velocities, and do not include contributions from rotational or internal vibration energy. To a reasonable approximation, the ejection KEs can be fitted to a Boltzmann distribution, ρ(E s ) =Ē −1 s exp(−E s /Ē s ), withĒ s = E s = 0.6 eV [26]. Similar Boltzmann-like behaviour can be observed for all impact energies, as well as for other individual species with sufficient statistics. Our choice of a near-Boltzmann distribution of energies is an attempt to describe sputtering processes in which the energy of the incoming particle is largely dispersed and randomized by the initial stages of the collision cascade before the sputtering products are ejected. The average ejection energy of 0.6 eV represents only 1.5% of the incoming KE for a 20 eV D −1 impact by molecular D 2 . However, 0.6 eV is also still much higher than thermal energies, clearly indicating that the incoming energy has not been fully thermalized. A molecular kinetic energy of 0.6 eV corresponds to a thermal ejection from a collision cascade that is locally still at a temperature of over 6500 K, far above the 300 K equilibrium substrate temperature. Our energy distribution can also be fitted to a Thompson distribution ∼E/(E + E B ) 3 like the one used by Marian et al [17] provided that a very different binding energy is used E B ∼ 0.2 eV. It is not clear whether this difference is due to the increased physical sputtering component of the 100 eV deuterium impacts in [17], or to the fact that they have calculated sputter yields for a non-self-consistently prepared surface. We note, however, that a Boltzmann distribution may fit their data as well as the Thompson distribution.
Our interpretation is further supported by figure 13, which shows that the translational KE of a sputtered particle does not vary strongly with the bombardment energy. If the sputtered particles represent long-delayed desorption from a fully thermalized surface, ejection energies would be at thermal energies, and independent of incident energy. That this kinematic ejection energy is weakly dependent on incident energy, much lower than incident energies, but much higher than thermal energies, confirms the view that ejection takes place after a largely, but incompletely thermalized collision cascade, locally hot due to the particle impact.
The dependence of ejection energy on particle mass is also very weak, providing further evidence for substantial, but incomplete, thermalization. Figure 14 shows the averaged translational KE for ejecta of all the different species observed, at individual bombardment energies from 7.5 to 30 eV D −1 . A number of very interesting features can be observed in this plot. First is the relative invariance of average ejection energy with particle mass; it varies somewhat, but typically with just statistical variation. Two species, however, have anomalously low ejection energies: C atoms and C 2 molecules. This persists for all bombardment energies, indicating a real physical effect rather than a statistical artifact. (The low energy at mass 50 is statistical in nature, resulting from only a single ejection.) It is not clear from the current analysis whether the anomalously low ejection energies for C and C 2 result from characteristic ejection at later times when the impact region has been more fully thermalized, or from a different mechanism of ejection of bare carbon species, as proposed by Marian et al [17]. The average ejection energies at each mass (for the statistically significant species) are also weakly dependent on impact energy, consistent with the data in figure 13.  The decomposition by mass reveals one additional trend, which is made more obvious by the mass-dependent yields presented in figure 15: higher impact energies generate more ejecta of larger mass. At the lowest incident energy of 7.5 eV D −1 , no species heavier than C 2 D 5 (mass 34) were observed. At 10 eV D −1 , a small number of heavier species are observed, and the C 2 D y species are comparable in total yield with the C 1 D y ejecta. Between the energies of 20 and 30 eV D −1 , the C 1 D y and C 2 D y yields have stabilized, and much of the growth in the sputter yield is due to increases from more complex species. This indicates that not just the yield but the distribution of sputtered species changes quite noticeably with impact energy.
The distribution of the x-and y-components of the ejected particle velocities is uniform across all azimuthal angles, with approximately zero average vectorial velocity in the x-y-plane. However, the distribution of polar emission angles is not uniform, with preferential ejection in the z-direction, for all impact projectile energies.A quantitative picture is provided in figure 16, which shows the distributions of polar angles θ for total hydrocarbons (figure 16(a)) and methyl radicals ( figure 16(b)), for various energies of the vibrationally excited D 2 projectiles. These are expressed as differential counts per solid angle, dn/d = dn/2π sin θ dθ. The behaviour is similar in all cases, with a pronounced preference for ejection at small angles from the surface normal, decaying approximately exponentially with a width of approx. 40 • . These angular distributions are well described by cotangent distributions, which yield cosine distributions for dn/dθ.
Most of the impacting D and D 2 are reflected. The most complex case is that for D 2 projectile, which can dissociate, in addition to reflecting or implanting. The variation of this branching with impact energy is illustrated in figure 17 for the case of the D ( * ) 2 and D reflection probability, P r , is here defined as the total number of 'emitted' D 2 per an impinging D 2 .
In figure 17(a) we do not distinguish between true reflections and sputtered D 2 , as they would not be distinguished experimentally. The dissociation probability, P d , is defined as the number of emitted D per D 2 projectile, and is shown in figure 17(b). Finally, we define the implantation probability, P i , as the number of D atoms per impacting D that were not counted as either reflected or dissociated, i.e. P i = 1-P r -P d /2. This is shown in figure 17(c). The implanted particles may be emitted back to the plasma as hydrocarbons. Not surprisingly, because the ground state D 2 molecules are more compact and more strongly bound, they have larger reflection probability, smaller dissociation probability and lower implantation probability at low energies. The dissociation probability is not too sensitive to the impact energy, while reflection is quite probable at low energies (∼0. 8), decreasing to about 0.2 at the highest energies considered here. It is interesting to consider the energy distributions of the reflected particles (D 2 ) and dissociated products (D) upon impact of D ( * ) 2 , as illustrated in figure 18, for the lowest impact energy considered, 15 eV. Considering that the energy spectrum of emitted particles contains both (inelastically and elastically) reflected and sputtered particles we expect two superimposed Boltzmann-like distributions: one corresponding to low sputtering energies (0.6 eV) characteristic of the chemically sputtered particles from the partially thermalized collision cascade, and another one, with a higher average energy, representing inelastically and elastically scattered particles. The data support this intepretation, as illustrated in figure 18, with both high-and low-energy components to the energy distribution. The inelastically reflected particles have an average energy of 1.75 eV for D 2 and ∼3 eV for D atoms. We note that the distributions of the emitted D and D 2 cannot be fitted well to a Thompson type of distribution.

Conclusions
We have shown that when proper attention is given to the preparation of the deuterated carbon target surfaces as a function of impact fluence, energy and type of projectile species and to the internal (rovibrational) state of the molecular projectile, mimicking as nearly as possible the particular experimental conditions, one can expect good quantitative agreement of the sputtering yields between simulation and experiment. This allows MD simulations to be used with some predictive power. Furthermore, by analysing the microstructure of the prepared bombarded 24 DEUTSCHE PHYSIKALISCHE GESELLSCHAFT surfaces one can obtain good qualitative estimates of the behaviour of the sputtering yields and their variation with impact energy and type of sputtered species. We have also shown that one has to account carefully for the vibrational state of the impacting molecule in low energy sputtering, since the sputtering yields are correlated with the degree of excitation of the molecule. The vibrationally excited D 2 molecules are expected to be produced by neutralization of the impinging D + 2 molecular ions above the surface by resonant charge exchange with deuterium that has accumulated at the surface. Although not a subject of the present study, it is quite possible that a similar dependence can be found on the initial rotational state of the molecular projectile; this will be the subject of a forthcoming study.
We have presented the distribution of sputtered species, as well as their KE and angular spectra for impacts in the energy range of 7.5-30 eV D −1 . The energy distributions are well described by Boltzmann distributions. The average KEs of the ejecta are, surprisingly, distributed about means of 0.5-1 eV that are insensitive to the mass of the emitted particles and only weakly dependent on the impact energy. The velocity vectors of the sputtered particles are predominantly oriented in the z-direction, perpendicular to the surface of the simulation cell. The D 2 and D particles, emitted by reflection, sputtering and dissociation of the impact D 2 follow the sum of two Boltzmann distributions: a low-energy one, corresponding to sputtering processes, and a higher energy one, corresponding to the inelastic reflection of the particles.
Here, we are mimicking beam experiments. We have stressed the importance of the surface preparation under specific conditions (type of impinging particles, energy, angle, internal state and surface history) that replicate the experiment. Consequently, it is to be expected that sputtering in beam experiments will differ somewhat from that due to a plasma, due to the spread in angle, energy, density of impacting particles and surface history. To properly model sputtering in a plasma environment, one would need to prepare the surface with the appropriate energy, angle and particle distributions along with their time-resolved profiles. This will be a subject of our future studies.