Atomic nitrogen is formed in the high-temperature shock layer of hypersonic vehicles and contributes to the ablation of their thermal protection systems (TPSs). To gain atomic-level understanding of the ablation of carbon-based TPS, collisions of hyperthermal atomic nitrogen on representative carbon surfaces have recently be investigated using molecular beams. In this work, we report direct dynamics simulations of atomic-nitrogen [N(4S)] collisions with pristine, defected, and oxidized graphene. Apart from non-reactive scattering of nitrogen atoms, various forms of nitridation of graphene were observed in our simulations. Furthermore, a number of gaseous molecules, including the experimentally observed CN molecule, have been found to desorb as a result of N-atom bombardment. These results provide a foundation for understanding the molecular beam experiment and for modeling the ablation of carbon-based TPSs and for future improvement of their properties.
I. INTRODUCTION
Entry of space vehicles into Earth’s atmosphere at high velocities generates a thin, high-temperature shock layer above the vehicle’s leading surface, typically a thermal protection system (TPS) material. Within this shock layer, temperatures as high as 10 000 K may be reached, leading to dissociation of O2 and N2.1,2 The resulting O and N atoms diffuse through the boundary layer and react with the TPS surface. If the TPS is ablative, reaction products are transported back into the flow, which, in turn, affects the chemical state of the boundary layer.3 The degradation, or ablation, of the TPS depends on the coupled gas-surface and gas-phase processes that add heat to the TPS and result in gas-surface reactions at high temperatures. To address the challenges of atmospheric entry, predictive models are used to guide TPS and vehicle design. The fidelity of such models is improved by increasing the understanding of the molecular-level chemistry under relevant extreme conditions.
Many TPS materials are based on carbon (e.g., carbon–carbon composite) or become carbonized through charring (e.g., pyrolysis of a phenolic ablator), so a clear understanding of oxidation and nitridation reactions on carbon is important for many TPS systems.4,5 As a result, there has been strong recent interest in experimental and theoretical studies of collisions of the relevant gaseous species with carbon-based materials such as highly oriented pyrolytic graphite (HOPG) and vitreous carbon.6–11 While the reactions between O2 and O with carbon surfaces have been extensively studied,6,8,9,11,12 fewer studies have been conducted on nitrogen reactions on carbon.
Nevertheless, there have been some investigations of the effect of nitrogen atoms on carbon surfaces, many of which were conducted in the context of TPS ablation. N atoms were reported by McCarroll and McKee to etch graphite surfaces anisotropically above 1000 K, leaving characteristic hexagonal pits.13 These pits are similar to those found with O-atom etching, suggesting that the two species could react through similar reaction mechanisms at surface defects.6 However, no chemistry was observed with impinging molecular nitrogen. The ablation of a graphite surface by N atoms was also reported by Suzuki et al., who found that for the range of surface temperatures (TS), 1351 K–2184 K, the probability of reaction increased from 1.4 × 10−3 to 3.2 × 10−3.14 Another study of mass-loss measurements from samples of purified graphite (TS = 873 K–1373 K) exposed to N atoms showed reaction probabilities of 0.2–9.8 × 10−3.15 To obtain information on the mechanism of the mass loss, products resulting from such reactions have also been detected using various techniques.16–18 It was found that the dominant product was C2N2, which might have been produced via condensation of the nascent CN radicals produced by the reaction of atomic N with the carbon surface. However, there has been a lack of quantitative data on the CN formation.19
In addition to its relevance to TPS chemistry, nitridation of graphene is also of considerable interest in materials science.20 The motivation here is to create N-doped graphene materials that have desired semiconducting properties.21 One such approach is to bombard graphene with neutral or charged atomic nitrogen in plasmas, which causes reactions of these species with the carbon network.22,23 However, the detailed mechanisms of these processes are still poorly understood.24
In order to gain mechanistic insights into the nitridation of carbon, detailed gas-surface scattering investigations have recently been performed. Molecular beam experiments by Murray and Minton studied the scattering of hyperthermal N(4S) atoms and N2 molecules with vitreous carbon at high temperatures between 1023 K and 1923 K.10 The molecular beam contained N and N2 with translational energies of 110 kcal/mol and 193 kcal/mol, respectively. They found that N2 scattered non-reactively, while some of the N(4S) reacted with the surface to form CN with an Arrhenius activation energy of 49.5 kcal/mol. It was concluded that the CN is formed by a Langmuir–Hinshelwood mechanism as CN was found to desorb in thermal equilibrium with the surface. Subsequent work by Murray et al. with a much higher flux of lower-energy N atoms (∼8 kcal/mol) is more relevant to the conditions in the shock layer on the leading edge of hypersonic vehicles and is thus more relevant for the modeling of TPS chemistry. This study also observed the CN reaction product and further refined the CN pathway activation energy to 41.1 kcal/mol. These authors also reported evidence for the formation of N2 from recombination of N atoms. Again, it was believed that N2 was produced via a Langmuir–Hinshelwood mechanism as the translational energy and angular distributions of the molecular product were consistent with thermal desorption. Overall, N atoms were found to be less than 5% as reactive with the surface as O atoms when the surface temperature was below 1900 K.11
To understand the microscopic details of the aforementioned experimental observations, it is necessary to carry out theoretical investigations of the reactive scattering. Although non-reactive scattering of N2 from carbon surfaces has been investigated,25,26 there have been few chemical dynamics investigations that allow the breaking and formation of chemical bonds. The only exception was a recent direct dynamics simulation, in which up to 100 N atoms were directed to a pristine graphene sheet in order to understand nitrogen doping resulting from nitrogen plasma exposure.27 With incident energies of 23 kcal/mol, 46 kcal/mol, and 92 kcal/mol, atomic nitrogen was found to react with the graphene sheet and release various molecular species including N2, CN, CN2, and C2N2. These incident energies are higher than desired for modeling TPS chemistry.
In this work, we report direct dynamics studies of the reactive scattering of hyperthermal atomic nitrogen in the ground electronic state [N(4S)] on a variety of pristine, defected, and oxidized carbon surfaces, simulated by using the corresponding graphene sheets. The modeling of N-atom reactive scattering from oxidized graphene is of particular significance as it is relevant to air-carbon ablation chemistry, and no experiment has yet been performed that probed N-atom reactions on an oxidized carbon surface. These calculations, which use an efficient semi-empirical tight-binding density functional theory (TB-DFT), allow not only the simulation of non-reactive events but also breaking and formation of chemical bonds. The choice of this approximate DFT method is primarily due to its efficiency, which is important for large systems and for long-time trajectories. The aim here is to explore various reaction channels under conditions similar to the aforementioned molecular beam experiments10,11 and to quantify the reactivity of each channel. To this end, the chemical dynamics calculations were performed for atomic N scattering from model graphene surfaces at a relevant experimental temperature (1375 K) including a pristine, non-defected structure [model 1, Fig. 1(a)], a single carbon vacancy (SV) [model 2, Fig. 1(b)], and two oxygenated structures based on the SV surface [models 3 and 4, Figs. 1(c) and 1(d), respectively]. The graphene surfaces are modeled with periacene, and its modified analogs are modeled with artificial edge constraints. Various gas-phase species as well as surface species have been found as a result of the reactive scattering, and the reaction mechanisms are discussed.
Periacene (5a, 6z) graphene models: (a) pristine (model 1), (b) single carbon vacancy (SV) defect (model 2), (c) SV defect with one oxygen atom (model 3), and (d) SV defect with two oxygen atoms (model 4).
Periacene (5a, 6z) graphene models: (a) pristine (model 1), (b) single carbon vacancy (SV) defect (model 2), (c) SV defect with one oxygen atom (model 3), and (d) SV defect with two oxygen atoms (model 4).
II. METHODS
A. Electronic structure calculations
Direct dynamics simulations were performed in which the energy and forces are computed on the fly using the semi-empirical self-consistent charge density functional tight-binding (SCC-DFTB) method.28 The SCC-DFTB (DFTB for short) method is an approximate density functional theory (DFT) method based on a second-order expansion of the Kohn–Sham total energy in terms of the charge density fluctuation relative to a reference density. It is much more efficient than the conventional DFT method as the Hamiltonian, which contains only terms for nearest neighbors (tight-binding), is parameterized semi-empirically. This method has been successfully employed in many applications, ranging from materials29 to biological systems.30 Indeed, it was also used in the aforementioned work by Moon et al. in a direct dynamics simulation of N-atom scattering from a graphene surface.27 In our calculations, the third-order correction31,32 was utilized with the 3ob-3-1 Slater–Koster parameter set33 and the Lennard-Jones (LJ) dispersion model.34
Because of the semi-empirical nature of the SCC-DFTB method, it is important to assess its accuracy and reliability in describing the systems of our interest. To this end, a detailed comparison has been performed here against the conventional DFT method with the three-parameter Becke–Lee–Yang–Parr (B3LYP) functional,35,36 with the 6-31G* basis set.37 In addition, the D3 dispersion correction38 with Becke–Johnson damping (D3BJ)39 was used to handle the dispersion corrections. Fermi smearing was utilized to allow for quick convergence of the self-consistent processes. The quartet spin state was selected in all calculations. In such comparisons, potential energy curves were generated first by optimizing geometry with fixed C–N or C–CN bond distances at several intervals with the B3LYP-D3BJ method. The DFTB3-LJ/3ob-3-1 energies were then determined at these geometries for comparison. Hydrogen atoms were kept fixed to simulate the geometry restrictions of a larger graphene sheet. All DFT calculations were performed with the Gaussian 16 program,40 and the DFTB calculations performed with the DFTB+ program.41
B. Direct dynamics simulations
Although there have recently been several theoretical studies of atomic and molecular scattering from carbon surfaces, many were restricted to non-reactive events.25,26,42–45 These studies used empirical force fields that may or may not have an accurate descriptions of the reactive events. To explore the chemical channels in atomic-nitrogen collisions with graphene, we chose to use a direct dynamics approach,46 in which the energy and forces for the nuclear motion are computed on the fly. Such approaches have been successfully used by several groups recently on scattering dynamics involving HOPG.9,47–49 Here, because of the relatively large number of atoms (and electrons), we have chosen to take advantage of the efficient DFTB Hamiltonian.
Direct chemical dynamics simulations were carried out using the VENUS chemical dynamics program50 externally interfaced to the DFTB+ program. Namely, the atomic positions of the whole system provided by VENUS were read by the DFTB+ program and used to calculate energies, which were then used by VENUS to generate the next trajectory steps. The atomic forces were determined numerically using finite differencing (step size of 0.01 atomic units). In this way, it is possible to use all the quasi-classical trajectory features present in VENUS (for initial condition generation and for the propagation of an ensemble of trajectories) with the DFTB method. The VENUS/DFTB+ interface can be obtained upon request.
In the work reported here, the DFTB3-LJ/3ob-3-1 method was used for all the on-the-fly direct dynamics simulations. The calculations were performed with graphene models based on periacene (pristine and modified with representative defects) and the ground state N(4S) atom. Initial conditions utilized the options in the VENUS program package for gas-surface interactions. The mass of every hydrogen atom at the edge of the periacene molecule was set to 1000 amu in order to simulate the restrictions present in a larger graphene sheet. The surface normal modes were populated from a Boltzmann distribution with TS = 1375 K. The algorithm for selecting the initial conditions for a gas-surface collision has been described in an earlier work,51 and the coordinates are displayed in Fig. 2(a). The initial position and momentum of the N(4S) atom were chosen with respect to a surface plane and defect site. The surface plane was defined by three carbon atoms, NN1, NN2, and NN3, at the edge of the periacene sheet [Fig. 2(b)]. The reaction site was defined by displacing the origin of the surface plane, determined by NN1, by distances RX, RY, and RZ. The initial atomic coordinates and momenta were chosen randomly from the Boltzmann distribution at the surface temperature. Five coordinates were used to characterize the initial position and orientation of the N(4S) velocity vector with respect to the surface: b, θ, φ1, φ2, and s. The coordinates, b and θ, designate an aiming point on the surface. The impact parameter, b, sets the distance from the defect site for the graphene sheet (models 1–4) and varies from 0 Å to 3 Å in intervals of 0.5 Å. For each value of b, 50 trajectories were calculated. The angle, φ1, was chosen randomly from a uniform distribution from 0° to 360°. The angle, θ, which sets the angle of the velocity vector of the N(4S) atom projectile relative to the surface normal, was fixed at 0°. The initial collision energy of the projectile used to impact the surface was 14.9 kcal/mol and 110 kcal/mol, designed to simulate collision velocities for hypersonic and low-Earth-orbital vehicles, respectively. (We note in passing that the latter could induce electronic excitations, which are not included in our theoretical model.) The angle, φ2, of the N(4S) atom projectile was chosen randomly between 0° and 360° but was entirely redundant in normal collisions, having no effect on the reaction dynamics. The separation between the graphene sheet and the nitrogen atom, s, was chosen to be 10 Å, sufficiently large to ensure no interaction.
Schematic for (a) the coordinate definitions used for initial conditions, and (b) carbon atoms, NN1, NN2, NN3, (red circles) define the surface plane in model 1 (these atom choices were made in the same way for all models).
Schematic for (a) the coordinate definitions used for initial conditions, and (b) carbon atoms, NN1, NN2, NN3, (red circles) define the surface plane in model 1 (these atom choices were made in the same way for all models).
All trajectories were integrated in time with the sixth order symplectic integration algorithm,52,53 utilizing a time step of 0.5 fs and total simulation time of 3 ps. A trajectory is considered non-reactive if the nitrogen atom impacts the surface and leaves without forming a bond or altering the surface in a permanent way within the simulation time, which is confirmed with visualization of the trajectory. Occurrence of novel surface group functionalization, alteration, and gaseous product formation was also ascertained from trajectory visualization.
III. RESULTS
A. Assessment of the SCC-DFTB model
A (5a, 6z) periacene was used to model a larger graphene sheet such as that utilized in a previous dynamics study with O2.9 Full optimization of the periacene with DFT and DFTB methods leads to root mean square difference (RMSD) values between identical atomic centers that never exceed 0.058 Å for any structure (Fig. S1), suggesting that the DFTB method provides an excellent geometric description of the model graphenes.
For the interaction of an N(4S) atom with pristine graphene (model 1), two routes were found from the direct dynamics simulations: N bonding to a single carbon center or N bonding to two carbon centers at a C–C bond. In the first case of an N atom bonding to a single carbon atom, the optimized B3LYP-D3BJ/6-31G* potential energy curve (PEC) features a barrier and a minimum for the bonded structure, which are 11.0 kcal/mol and −14.2 kcal/mol, respectively, with respect to the dissociation limit (Fig. S2 of the supplementary material), respectively. The minimum features an sp3-like carbon connected with the nitrogen. Along this PEC, the bonding carbon atom can be seen to rise smoothly out of the graphene plane as the N atom approaches the surface and makes the C–N bond (Fig. S3). The corresponding DFTB3-LJ/3ob-3-1 energy along the DFT PEC follows the same trend, but the barrier and minimum are both higher in energy, suggesting underbinding.
In the second N + graphene route of N approaching a C–C bond, the nitrogen atom bonds with both carbon atoms and the corresponding PEC is shown in Fig. S4. The potential minimum at the DFT optimized geometry is below the dissociation limit by about 20.4 kcal/mol. The corresponding DFTB binding energy is at 8.4 kcal/mol above the dissociation limit. The barrier energy of about 15.5 kcal/mol found with DFT at a C–N separation of 2.1 Å is larger in DFTB at 24.0 kcal/mol. The bonding carbon atoms again smoothly rise out of the plane as the nitrogen atom approaches the sheet (Fig. S5). This comparison further suggests that DFTB underestimates the binding of nitrogen with graphene.
The PEC of a product molecule, CN, leaving an oxygenated surface in model 3 was also examined. In this case, a nitrogen atom would have bonded to a carbon atom at the inserted oxygen defect, pulled out of the surface, breaking the C–CN bond and leaving the surface as CN. The potential minimum for CN bonded to a carbon at the defect site is −104.7 kcal/mol with DFT and is higher in energy, −90.2 kcal/mol, with DFTB (Fig. S6). No potential barrier is observed along the PEC. Selective structures from the PEC in Fig. S7 show the dissociation of CN. As with the other cases, the bonding carbon attaching CN to the sheet falls back into the plane as the CN leaves the sheet (Fig. S8). Analysis of the CN bond distance shows that it does not vary greatly with respect to the C–CN fixed bond distance, decreasing by only about 0.01 Å compared to the free CN at the dissociation limit, 1.178 Å (Fig. S9). This compares well to the literature value for a CN triple bond of 1.157 Å.54
Finally, the PEC for NO leaving model 4 after an impinging N atom reacts with the carbonyl oxygen is found in Fig. S10. A 17.9 kcal/mol barrier is found at a C–ON distance of 1.9 Å, while the potential minimum is located above the complete dissociation limit by 7.28 kcal/mol and a C–ON distance of 1.39 Å. DFTB under binds by about 19.3 kcal/mol at the DFT minimum energy structure. The fully optimized DFTB structure finds the minimum to be 19.0 kcal/mol with a C–ON distance of 1.36 Å. The bonding carbon atom (Fig. S11) goes back into the plane as NO leaves forming the model 3, single oxygen insertion defect.
Based on these comparisons, we conclude that DFTB provides qualitatively correct interaction patterns for nitrogen interaction with graphene but underestimates the interaction energy.
B. Reactive scattering
1. Model 1—Pristine graphene sheet
Most trajectories in model 1 with a translational energy of 14.9 kcal/mol (Fig. 3, red) result in the nitrogen atom rebounding from the surface without reacting. The velocity distribution for the scattered nitrogen atoms is shown in Fig. 4. It is clear that the scattered N atoms exchanged much energy with the surface, mostly losing energy to the surface but with a small fraction gaining some energy. The average kinetic energy of scattered N is 9.42 kcal/mol, and the corresponding velocity distribution compares reasonably well with a Maxwell–Boltzmann distribution at 4750 K, much higher than the surface temperature in the experiment.10
Probability of reaction, Pr(b), with respect to the impact parameter, b, from DFTB3-LJ/3ob-3-1 direct dynamics simulations of model 1.
Probability of reaction, Pr(b), with respect to the impact parameter, b, from DFTB3-LJ/3ob-3-1 direct dynamics simulations of model 1.
Velocity distribution of scattered N(4S) atoms from DFTB3-LJ/3ob-3-1 direct dynamics simulations of model 1 with N(4S) translational energy set to 14.9 kcal/mol. Maxwell–Boltzmann distributions shown for 1923 K and 4750 K.
Velocity distribution of scattered N(4S) atoms from DFTB3-LJ/3ob-3-1 direct dynamics simulations of model 1 with N(4S) translational energy set to 14.9 kcal/mol. Maxwell–Boltzmann distributions shown for 1923 K and 4750 K.
The reactivity is quite small (only five trajectories react in total) and is mostly found with large impact parameters, b, nearing the edge of the periacene model, but never increases beyond a probability [Pr(b)] of 0.06 (3/50 trajectories) for a given b value. The only reactive event observed for model 1 at this incidence energy is nitrogen insertion into the graphene sheet breaking a C–C bond, as shown in Fig. 3. In addition to five such reactive trajectories, there are 11 additional trajectories that react with the surface, with N inserting into the sheet, but finally ejecting from the surface within 3 ps simulation time. In light of this, it is possible that these reactive trajectories contributing to the calculation of Pr(b) could eventually be ejected from the sheet with a longer simulation time if energy dissipation does not remove the energy in the vicinity of the impact site.
Increasing the translational energy of the N(4S) atom to 110 kcal/mol is found to enhance the reactivity (Fig. 3, blue). The probability of reaction for a given b value increases slightly to a maximum of 0.12 (six trajectories). The character of the reactive trajectories does not change though, producing nitrogen insertions into the graphene sheet. The number of trajectories where the nitrogen atom projectile reacts with the surface by inserting, but then is ejected before the total 3 ps passes, increases greatly to 110 of the total 350 trajectories. This accounts for the greater reactivity for the 14.9 kcal/mol trajectories at b = 3.0 Å: more of those are likely to stay bonded to the surface during the simulation time, whereas the N atom with 110.0 kcal/mol of incidence energy is much more likely to leave the surface after reacting. This is the only other change in reactivity associated with the higher translational energy N(4S) atom impacting the pristine graphene sheet. At both incidence energies, the nitrogen atom is not observed to undergo significant diffusion on the graphene surface.
2. Model 2—SV graphene
The SV defect of model 2 displays increased reactivity compared to the pristine sheet of model 1. It has been well established that the SV defect structure is more reactive compared to pristine graphene, owing to the presence of dangling bonds.55,56 The probability of reaction at low incidence energy (14.9 kcal/mol, Fig. 5, black) for values of b from 0.0 Å to 1.0 Å is either 1.0 or very close. The number of reactive trajectories then decreases with increasing b as the number of trajectories resulting in N-atom scattering increases.
(Top left) Probability of reaction, Pr(b), with respect to the impact parameter, b, from direct dynamics simulations of model 2 with N(4S) translational energy set to 14.9 kcal/mol (red) and 110.0 kcal/mol (black). Probabilities for the two primary N-functionalization products of the 14.9 kcal/mol (middle left) and 110.0 kcal/mol (bottom left) simulations. Typical configurations of these major reaction scenarios are shown in (a) N insertion in SV defect, (b) singly bonded CN group, (c) N insertion in C–C bond, (d) SW insertion, (e) SQ insertion, and (f) CN product formation.
(Top left) Probability of reaction, Pr(b), with respect to the impact parameter, b, from direct dynamics simulations of model 2 with N(4S) translational energy set to 14.9 kcal/mol (red) and 110.0 kcal/mol (black). Probabilities for the two primary N-functionalization products of the 14.9 kcal/mol (middle left) and 110.0 kcal/mol (bottom left) simulations. Typical configurations of these major reaction scenarios are shown in (a) N insertion in SV defect, (b) singly bonded CN group, (c) N insertion in C–C bond, (d) SW insertion, (e) SQ insertion, and (f) CN product formation.
The trajectories at 14.9 kcal/mol incidence energy display a variety of different nitrogen functionalized configurations in addition to N insertion between into a C–C bond [Fig. 5(c)]. For values of b = 0.0 Å, 0.5 Å, and 1.0 Å, almost all of the trajectories result in the N atom inserting into the SV defect [Fig. 5(a)]. As expected, the probability of this process decreases with b thereafter. Beginning with values of b ≥ 1.5 Å, a singly bonded CN group becomes another major configuration when the N atom bonds to a carbon at the defect and pulls it out of plane [Fig. 5(b)]. In addition, N-atom insertions at the defect site can also result in Stone–Wales (SW) [Fig. 5(d)] or square (SQ) configurations [Fig. 5(e)]. Gaseous CN species are also observed coming off the surface, leading to a double carbon vacancy after forming in a similar way to the singly bonded CN group [Fig. 5(f)].
Increasing the incidence translational energy to 110.0 kcal/mol decreases the probability of reaction for b = 0.0 Å and 0.5 Å because the incoming N atom can pass through the SV defect without reacting or scattering. For larger values of b, the probability of all other reactive configurations increases. A difference is seen for the singly bonded CN and CN product trajectories for small b values. The N atom bonds to a carbon atom at the SV defect site and pulls it out of plane as it passes through the vacancy. These trajectories result in functionalization of the opposite side of the sheet.
There is a limitation in the current model’s ability to fully characterize this reaction with 110.0 kcal/mol of incidence energy as the N atoms are found to pass through the SV defect into vacuum. In an extended system with additional layers or in a real bulk sample, the N atoms would encounter another layer with which to react or from which to rebound. In the case of rebounding, the N atom could then react with the top layer from underneath. Reactions of rebounding N atoms could play a role in the pristine case described earlier in this work and would certainly impact the results of the defect structures discussed below.
3. Model 3—SV graphene with one oxygen atom
The presence of an oxygen atom at the SV defect of graphene has a large effect on the reactivity of the N(4S) atom, as illustrated in Fig. 6 (red and black columns for incidence energies of 14.9 kcal/mol and 110 kcal/mol, respectively). At the lower incidence energy, reactivity of the sheet toward the nitrogen atom decreases with increasing b (increasing distance from the defect). For trajectories aiming directly at the defect (b = 0.0 Å), Pr(b) is the largest at 0.8 and smoothly decreases to 0.18 for trajectories aiming furthest from the defect (b = 3.0 Å).
(Top left) probability of reaction, Pr(b), with respect to the impact parameter, b, from direct dynamics simulations of model 3 with N(4S) translational energy set to 14.9 kcal/mol (red) and 110.0 kcal/mol (black). Probabilities for the two primary N-functionalization products of the 14.9 kcal/mol (middle left) and 110.0 kcal/mol (bottom left) simulations. Typical configurations of these major reaction scenarios are shown in (a) for the singly bounded CN group, (b) the doubly bonded CN group at the defect site, (c) N insertion into the C–C bond, (d) N sticking on surface, and (e) replacement of O by N.
(Top left) probability of reaction, Pr(b), with respect to the impact parameter, b, from direct dynamics simulations of model 3 with N(4S) translational energy set to 14.9 kcal/mol (red) and 110.0 kcal/mol (black). Probabilities for the two primary N-functionalization products of the 14.9 kcal/mol (middle left) and 110.0 kcal/mol (bottom left) simulations. Typical configurations of these major reaction scenarios are shown in (a) for the singly bounded CN group, (b) the doubly bonded CN group at the defect site, (c) N insertion into the C–C bond, (d) N sticking on surface, and (e) replacement of O by N.
No gaseous species were found to leave the graphene sheet though the nitrogen impact functionalized the surface resulting in two primary products at the defect site: a singly bonded CN group [Fig. 6(a)] and a doubly bonded CN group [Fig. 6(b)]. The former is found to be the major product and forms after the nitrogen atom attaches to a carbon atom at the oxygen defect, which sticks out of the plane, though it still remains bonded to the graphene sheet in all cases. Conversion between the two forms has been observed. Comparing the black and green columns in Fig. 6 (top left and middle left plots), this reaction is found to account for almost all of the reactions at each impact parameter, b. The formation of a CN bond has a Pr(b) smaller than that of the attached CN group but still has a probability of greater than 0.1 for b = 0.5 Å, 1.0 Å, and 1.5 Å. This reaction occurs in a similar way to the attached CN group, but the nitrogen atom only inserts in the CO bond and does not pull a carbon atom out of plane. The reaction probability for this product decreases to zero by b = 2.5 Å.
Increasing the incidence energy of the impinging N(4S) atom to 110 kcal/mol increases Pr(b) significantly, especially for trajectories with b > 0.0 Å (Fig. 6, red). As was found for the 14.9 kcal/mol simulations, most trajectories result in the singly bonded CN group, but this group stays attached to the graphene for the duration of the 3 ps integration time. In contrast to the 14.9 kcal/mol case, at large values of b (1.5 Å–3.0 Å), a small number of trajectories (eight) is found to have a gaseous CN molecule ejected from the graphene sheet. These gaseous CN molecules are always formed from the attached CN groups. It is possible that the simulation time is insufficient to capture more of these CN groups detaching from the surface.
A greater variety of nitrogen-surface functionalization is found with incident atomic N with 110.0 kcal/mol of translational energy (Fig. 6). The major configurations are presented in Figs. 6(a)–6(c). The primary surface functionalization remains the attached CN group as mentioned before, but many minor reaction pathways were observed and at higher probabilities than those at the lower translational energy. Additional surface events such as N-atom insertions into a C–C bond and N-atom bonding to a single carbon center on the surface (N stick) were found. At small values of b, in a small number of trajectories, the incoming nitrogen atom may displace the oxygen atom defect forming a three-fold coordinated nitrogen dopant.
4. Model 4—SV graphene with two oxygen atoms
Sampling model 4 with N(4S) atoms with 14.9 kcal/mol of incidence energy shows a rather different reactivity profile compared to the reactions with model 3 (Fig. 7, top left, black). Pr(b) values are much lower especially for small values of b aiming at or near the defect site. The mechanism governing this involves the carbonyl oxygen atom playing a strong directing role, knocking away the incoming nitrogen atom as it approaches the defect site. Away from the defect site, this carbonyl oxygen plays less of a role in reducing the Pr(b) values, and the reactivity becomes similar to that of model 3. The reactions that do occur (Fig. 7 middle left) are mostly N insertions into a C–C bond [Fig. 7(b)] or pulling a carbon atom out of the plane as a singly bonded CN group [Fig. 7(a)] as seen in the previous models. The latter mechanism serves to increase the size of the defect as well by damaging the surface. No gaseous product molecules are ejected from the surface.
(Top left) probability of reaction, Pr(b), with respect to the impact parameter, b, from direct dynamics simulations of model 4 with N(4S) translational energy set to 14.9 kcal/mol (red) and 110.0 kcal/mol (black). Probabilities for the two primary N-functionalization products of the 14.9 kcal/mol (middle left) and 110.0 kcal/mol (bottom left) simulations. Typical configurations of these major reaction scenarios are shown in (a) for the singly bounded CN group, (b) N insertion, (c) gaseous NO product formation, (d) gaseous CO product formation, (e) gaseous NCO product formation, and (f) N sticking on surface.
(Top left) probability of reaction, Pr(b), with respect to the impact parameter, b, from direct dynamics simulations of model 4 with N(4S) translational energy set to 14.9 kcal/mol (red) and 110.0 kcal/mol (black). Probabilities for the two primary N-functionalization products of the 14.9 kcal/mol (middle left) and 110.0 kcal/mol (bottom left) simulations. Typical configurations of these major reaction scenarios are shown in (a) for the singly bounded CN group, (b) N insertion, (c) gaseous NO product formation, (d) gaseous CO product formation, (e) gaseous NCO product formation, and (f) N sticking on surface.
Increasing the incidence energy to 110.0 kcal/mol greatly increases the reaction probability at all values of b (Fig. 7, red). While the carbonyl oxygen still directs the incoming N(4S) projectile as discussed above, the increase in energy results in greater probability of reaction closer to the center of the defect (small values of b). The probability of reaction increased with increasing b in contrast to model 3. An important contribution to the increased values of Pr(b) comes from the variety of gaseous product molecules, which decreases with b values.
Several gaseous products were observed to be ejected from the graphene sheet after the impact of the nitrogen atom projectile (Fig. 7, bottom left). Especially at values of b = 0.0 Å, 0.5 Å, and 1.0 Å, the increase in reaction probability comes primarily from the production of NO [Fig. 7(c)]. The reaction mechanism is Eley–Rideal (ER) recombination involving the direct reaction of gas phase species (N) and a surface species (O). The nitrogen atom strikes the carbonyl oxygen atom, and rather than being directed away from the sheet as was observed for 14.9 kcal/mol incidence energy, NO forms and is immediately ejected from the graphene surface. CO also forms [Fig. 7(d)] at larger b values after the nitrogen atom inserts into a C–C bond involving the carbonyl. The C–CO bond is then cleaved to form a C–N–C insertion at the SV defect in the sheet, while the CO is ejected [Fig. 7(e)]. Another major product is NCO, which is formed in much the same way as CO, but instead of CO leaving the surface, the nitrogen atom bonds to the CO, before the C–NCO bond breaks and NCO is ejected from the surface. Two minor products, oxygen atom produced by N–O replacement and CN, are also formed but at very low probabilities.
IV. DISCUSSION
The lack of reactivity in model 1 underscores the chemical inertness of pristine graphene. This point has been recognized in several previous studies.11,15 As discussed above, the N–C bond is minimally stable in DFTB relative to the asymptote, which is consistent with the simulation results. The underbinding of N in DFTB relative to the DFT results suggests that the nitrogen atom should perhaps be significantly more sticky and reactive than revealed in the current simulations. More quantitative simulations might have to await future investigations within the more reliable (and more costly) DFT framework. In addition, we note that these chemisorbed nitrogen species could eventually encounter another surface N species and react to form N2 via a Langmuir–Hinshelwood mechanism, which will then desorb. This could account for the observed N2 products in the recent experiment.11 Such events might take a long time, however, and are thus beyond the limit of the current direct dynamics simulations. However, the experimentally observed N2 product is very likely to be both internally and translationally hot due to the large exothermicity of the reaction. The fact that the N2 was found to have a Maxwell–Boltzmann distribution at the surface temperature11 supports the notion that the scattered/desorbed species undergo further collisions with the rough carbon surface that was used in the experiment.
Despite significant energy loss to the surface, interestingly, the scattered atomic nitrogen from defect-free graphene (model 1) possesses a kinetic energy distribution that is markedly hotter than the experimental surface temperature of 1923 K.11 As discussed above, the majority of the N-atom scattering trajectories are direct, without chemisorption on the graphene surface, and as such, N is unlikely to equilibrate with the surface in a single collision. This could be an artifact of the underbinding in DFTB or the small periacene model used to simulate graphene. On the other hand, we note that inelastic scattering of N atoms on vitreous carbon with 110 kcal/mol of incidence energy also showed high translational energy characteristic of impulsive, or direct, scattering.10 The experiment with 8.2 kcal/mol N atoms incident on a vitreous carbon surface found that non-reactive N atoms desorbed in a Maxwell–Boltzmann distribution at the surface temperature, but the surface in this experiment was very rough, and it is thus likely that the scattered N atoms suffered many collisions before escaping to the gas phase, leading to scattering angle randomization and significant energy transfer, ultimately yielding scattering dynamics that are characteristic of thermal desorption. This speculation could be verified experimentally by using pristine HOPG instead of vitreous carbon.
Extensive nitridation of the graphene surface was found in all models, but the precise functionalization differs. Our simulations suggest chemisorption at a bridge site, insertion into the C–C bond, the formation of CN moiety, and the insertion into the defect. Some of the motifs have been reported in a recent simulation of N collisions with graphene.27 These surface species underscore the wide range of bonding the nitrogen can have with the carbon surface. The ultimate fate of these species was not explored in the current study because of the long simulation time required. We plan to examine the diffusion and reaction barriers of these species in future work, which will shed light on the degradation mechanisms of carbon surfaces.
Our simulations found gaseous CN formed in models 2, 3, and 4 with graphene defects. The CN moiety is strongly bound with graphene, as shown in Fig. S6, with a binding energy of 104.7 kcal/mol and 90.2 kcal/mol at the DFT and DFTB levels of theory, respectively. Since the potential energy curve did not exhibit any barrier, the activation energy of 41 kcal/mol determined in the recent experiment11 could be due to thermal desorption of CN. However, the current calculated binding energy is significantly larger than the experimental value, and this disagreement could conceivably stem from the highly defected carbon surface. This speculation of course requires further investigation, particularly by calculating the corresponding free-energy change for the CN bond in different defect sites of graphene.
Our simulations further predicted several other gas phase molecules, such as NO, CO, NCO, and O, that have not been observed in experimental studies. NO, for example, represents a major product channel in model 4. However, NO is only formed at the higher incidence energy, so its relevance to the gas-surface chemistry on the TPS of a hypersonic vehicle in dense air might not be a significant concern. Unfortunately, this species was not investigated in the recent experiments, so the theoretical prediction has yet to be confirmed.
V. CONCLUSIONS
This work reports an exploratory study of the physical and chemical processes initiated by the collisions of energetic nitrogen atoms in their ground electronic state with pristine graphene and its modified forms, aiming to understand the underlying reaction dynamics related to recent experiments. The collision dynamics was computationally simulated with a model graphene sheet, using a direct dynamics method based on a semi-empirical electronic Hamiltonian. Apart from non-reactive scattering of the impinging N atom, several reactive events have been discovered. The nitrogen atom was found to stick to the graphene surface, insert into the C–C bond and defect, functionalize the graphene to form various surface species, and form products that desorb from the surface. In particular, the experimentally observed scattered N atoms and CN molecules were confirmed. Furthermore, several other gaseous species were predicted. The detailed information provides insight into the chemistry in the shock layer at the gas-surface interface above the TPS, which is invaluable for the design of TPSs for re-entry vehicles.
Despite the wealth of information provided by the simulations, it is also recognized that a complete and quantitative understanding of the chemical processes is still far from being realized. The semi-empirical method used here is shown to be reasonable in describing bond forming and bond breaking, but it is far from quantitative. Improvements are sorely needed to provide a more reliable characterization of the gas-surface interactions, presumably by using the more expensive and more accurate density functional theory. The modeling of the graphene sheet using periacene is another potential shortcoming of the current model, which might provide a poor description of the energy transfer process between the projectile and the surface, and totally ignores energy exchange between the graphene layers. These energy dissipation processes are expected to impact reactivity and dynamics in a significant way. A more realistic model is needed to simulate both scattering and energy transfer. Experimentally, a better characterization of the carbon surface is needed to compare with theoretical simulations. Work in these directions is underway in our groups.
SUPPLEMENTARY MATERIAL
See the supplementary material for comparison of the results calculated from DFT and DFTB methods.
ACKNOWLEDGMENTS
This work was supported by the Air Force Office of Scientific Research (Grant No. FA9550-18-1-0413 to H.G.), the Robert A. Welch Foundation (Grant No. D-0005 to W.L.H.), Texas Tech University (to W.L.H.), University of Colorado (to T.K.M.), and ANR DynBioReact (Grant No. ANR-14-CE06-0029-01 to R.S.). The computer time from the High Performance Computational Center at Texas Tech University and the Center for Advanced Research Computing at the University of New Mexico is gratefully acknowledged.
DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request.