Commentary: A heterocyst antisense RNA controlling its sense mRNA target to develop a cell-specific regulation D. Kirilovsky, Institute of Integrative Biology of the Cell, I2BC, CNRS, CEA, Univ Paris-Saclay, 91198 Gif sur Yvette M. Crespi, Institute of Plant Sciences Paris-Saclay, IPS2, 91192 Gif sur Yvette, France Cyanobacteria, prokaryotes doing oxygenic photosynthesis, are important fixers of CO 2 in oceans and lakes. While most of cyanobacteria assimilate only combined nitrogen compounds (such as nitrate, ammonium, urea), certain strains are able to fix atmospheric nitrogen when other sources of N are lacking. However, the enzyme responsible for N 2 fixation is inactivated by oxygen and then certain filamentous strains develop specialized cells, called heterocyst, in which there is no production of O 2 via photosynthesis (reviews about heterocysts: [START_REF] Flores | Compartmentalized function through cell differentiation in filamentous cyanobacteria[END_REF][START_REF] Herrero | The multicellular nature of filamentous heterocyst-forming cyanobacteria[END_REF][START_REF] Magnuson | Heterocyst Thylakoid Bioenergetics[END_REF][START_REF] Meeks | Cellular differentiation in the cyanobacterium Nostoc punctiforme[END_REF]. Also key enzymes involved in the photosynthetic fixation of CO 2 like Rubisco and phosphoribulokinase are absent whereas others are largely diminished. Furthermore, the heterocysts envelope is different and decreases the entry of gases inside these cells together with activation of terminal oxidases to cope with any residual oxygen interfering with its nitrogen-fixing function. When vegetative cells evolve into heterocysts, the differentiation process presents multiple structural and functional changes in these cells that require an extensive modification of gene expression under the control of the global nitrogen regulator NtcA and HetR, an specific regulator of its differentiation.

Non-coding RNAs (sRNAs) have been sown to regulate complex transcriptional networks in all life kingdoms. Although the majority of non-coding RNA is present in eukaryotes, many non-coding RNAs have been shown to regulate virtually any adaptive response of bacteria to environmental changes (Jimenez-Zurdo and Robledo, 2017). The canonical activity of prokaryotic non-coding RNAs involves base pairing to coding mRNAs through which their stability or translation is regulated and, in many cases, involves extensive perfect complementarity as it is the case for cis-acting antisense sRNAs (asRNAs). Natural antisense RNAs (asRNA) are a diverse group of bacterial RNAs first discovered many decades ago however their biological relevance has been appreciated more recently. Antisense transcript regulatory mechanisms affect different levels of gene expression including: transcription interference, transcription attenuation, translation stimulation or inhibition, and RNA stability [START_REF] Saberi | Natural antisense RNAs as mRNA regulatory elements in bacteria: a review on function and applications[END_REF]. Several advantages can be referred to asRNA-based regulation such as shorter regulatory times (direct targeting of the sense mRNA), specificity of gene regulation based on base pairing (including potentially simultaneous trans effects at different homologous loci) and the capacity of asRNA to modulate conformation interactions depending on temperature and other environmental changes. Abundant nitrogen-regulated antisense transcripts have been detected in cyanobacteria strains including the nitrogen fixing Nostoc PCC7120 [START_REF] Flaherty | Directional RNA deep sequencing sheds new light on the transcriptional response of Anabaena sp. strain PCC 7120 to combined-nitrogen deprivation[END_REF][START_REF] Mitschke | An experimentally anchored map of transcriptional start sites in the model cyanobacterium Synechocystis sp[END_REF]. Only few (6) of them have been characterized: IsiR, regulates the degradation of the isiA mRNA; As1_flv4, downregulates the flv2-4 operon; PsbA2 and 3, upregulate the expression of psbA2 and psbA3 by stabilizing their coding mRNAs; RblR, a positive regulator of the RubisCO; pilR, a negative regulator of pilA11 and finally furA which modulates the concentration of FurA, a regulator of iron assimilation is the only characterized in Nostoc PCC 7120.

In this issue, Olmedo-Verd et al. characterized in Nostoc PCC 7120, a first heterocystspecific asRNA, as_glpX, antisense to the glpX gene encoding the Calvin cycle enzyme sedoheptulose-1,7-biphosphatase/fructose-1,6-biphosphatase (SBPase) involved in CO 2 assimilation. First, the authors demonstrated its existence and transcriptional induction upon removal of combined nitrogen. This 458 nucleotide long asRNA requires the presence of HetR to be transcribed. By introducing a gfp gene under the control of the as_glpX promoter, the authors demonstrated that the asRNA is expressed only in heterocysts (and not in vegetative cells) since the very early stages of differentiation. In the as_glpX promoter a DIF1 motif was found, a cis-sequence presents also in other promoters of early differentiation genes [START_REF] Muro-Pastor | The heterocyst-specific NsiR1 small RNA is an early marker of cell differentiation in cyanobacterial filaments[END_REF][START_REF] Muro-Pastor | A combinatorial strategy of alternative promoter use during differentiation of a heterocystous cyanobacterium[END_REF].

By constructing a strain overexpressing the as_glpX in all cells and not only in heterocyts (only 10-15% of cells), the authors demonstrated its repressive effect on the glpX mRNA. Smaller RNAs were detected with the glpX5'-probe suggesting the triggering of mRNA degradation by the antisense transcript. Also, the concentration of the SBPase protein was lower in the overexpressing asRNA mutant confirming that as_glpX induces the degradation of glpX mRNA and repression of SBPase protein accumulation in heterocyst cells. This perturbation of antisense expression diminished the concentration of SBPase protein and largely affects culture growth and photosynthesis under high light conditions (400 µE m -2 s -1 ). The cells grew slower than control cells and their photosynthesis saturated at lower light intensities and the maximal rate was about 50% of that of control cells.

In conclusion, as_glpX is a first heterocyst specific asRNA with a known function in early stages of heterocyst differentiation to control photosynthetic O 2 production and permit nitrogenase function. Other HetR and NtcA regulated asRNAs may also contribute to fine-tune gene expression regulation during cyanobacterial differentiation. These results open wide perspectives for the research of regulatory mechanisms involved in bacterial differentiation.


Teemu T T Myllymäki

Aynur Guliyeva

Antti Korpi

Mauri A Kostiainen

Ville Hynninen

Nonappa

Patrice Rannou

Olli Ikkala
email: olli.ikkala@aalto.fi

Sami Halila
email: sami.halila@cermav.cnrs.fr

Lyotropic liquid crystals and linear supramolecular polymers of end-functionalized oligosaccharides

We synthesized permethylated maltoheptaose oligosaccharides, whose both ends have been functionalized with supramolecular binders 2-ureido-4[1H]-pyrimidinones (UPy) after single ringopening of -cyclodextrin counterpart. In 1,1,1,3,3,3-hexafluoro-2propanol (HFIP), they show lyotropic liquid crystallinity. In the dried state they allow linear saccharide-based supramolecular polymers by UPy-dimerization.

Polysaccharides consist of a wide range of naturally occurring biopolymers, including starch and cellulose. They and their modifications have attracted a major revived interest due to their abundance, safety, diversity, and sustainability as they are not extracted from oil-based resources.

1 Their classic applications are in pharmaceutics, cosmetics, foods, beverages, detergents, textiles, paper, and paint, but recently new high-tech functional application has been suggested by nanocelluloses.

By contrast to conventional polymers with covalently connected repeat units, supramolecular polymers are based on low molecular weight repeat units, which are connected to form one dimensional arrays via supramolecular interactions. [3][4][5][6] Therein, 2-ureido-4[1H]-pyrimidinone (UPy) is particularly feasible, undergoing dimerization by four cooperative hydrogen bonds. 7 The dimerization can be controlled by solvent exchanges. 8,9 More generally, UPy is also feasible in supramolecular constructions. 10-12 Surprisingly, linear supramolecular polymers based on oligosaccharide building blocks have not yet been reported. Still they could allow relevant design options due to their non-covalently connected chains, e.g., to tackle the microplastic pollution problem arising from the slow degradation of oil-based covalent polymers or to develop bioactive and adaptive glycomaterials with remarkable low-temperature processability and recyclability. 13 Such arguments encouraged us to explore routes for oligosaccharides where both ends are decorated with the UPy supramolecular binding motifs enabling the fabrication of supramolecular glycomaterials with promising applications. Note that this is synthetically untrivial as it is well-known that selectively reacting the nonreducing end of oligo-or polysaccharides remains very challenging. This urges us to identify key-enabling and synthetically accessible protection/deprotection paths.

On the other hand, lyotropic liquid crystal systems have been reported for various polysaccharides including cellulose, chitin, amylose or their derivatives. [15][16][17][18][19][20] This has allowed specific benefits, such as improved alignment in fiber spinning. Therefore, it can be asked whether also the repeating units of polysaccharide-based supramolecular polymers could form lyotropic liquid crystals with appropriately selected processing solvents.

In this work we present a concise synthetic scheme for monodisperse oligosaccharides symmetrically end-capped with UPy quadruple hydrogen bonding (H-bonding) motifs, aiming at lyotropic liquid crystallinity and supramolecular polymers. The synthesis is illustrated in Fig. 1 and fully characterized (ESI †). The synthesis was achieved from the commercially available trimethyl-β-cyclodextrin (TM-CD) (1), which is a derivative of amylose, or more generally, of starch. Therein the lateral OH-groups are completely protected by permethylation. It was treated with a 30% HClO 4 aqueous solution at 0 °C and warmed up to room temperature for 6 h to trigger a ring opening reaction yielding a mixture of linear malto-oligomeric chains. It consisted of a majority (79%) of maltoheptaose derivative (2) with a degree of polymerization of 7 (DP 7 ), accompanied by lower homologues DP 6 (8%), DP 5 (6%), DP 4 (5%), and DP 3 (2%), The lateral methyl protection of the OH-groups allows to perform topochemical ring-opening reactions of the two OH end-groups in the permethylated maltoheptaose (MeMH), still involving the challenge to accomplish both of the end-group reactions.

The reaction conditions were first tuned using n-hexyl isocyanate C 6 -NCO (4) (Fig. 1c) to modify the two hydroxyl groups at the reducing end C-1 I of the maltoheptaose derivative ( 2), and at the nonreducing end C-4 VII . Dibutyltin dilaurate (DBTL) and Et 3 N were needed to catalyze the reaction and to achieve high -stereoselectivity 22,23 for C 6 -MeMH-C 6 (7) (Fig. 3). The expected product was confirmed by Such reaction conditions were next used to react UPy-C 6 -NCO (3) with MeMH (2) to afford the targeted product (5). Note importantly, that the lower than DP 7 oligosaccharides present in the crude mixture resulting from ring-opening of TM-CD (vide supra: Fig. 1a) were removed during purification by flash chromatography on silica gel (ESI †). Finally, the synthesis yielded monodisperse permethylated maltoheptaose derivative (5), where both ends have been equipped with "sticky" UPy groups with a hexyl spacer, abbreviated as UPy-C 6 -MeMH-C 6 -UPy (5). The hexyl linker can create hydrophobic pockets to promote UPy dimerization even in polar solvents like water. [24][25][26][27] The carbamate linker between the hexyl and sugar in (5) can also mediate lateral hydrogen bonds between the neighboring molecules. Therefore, the building block ( 5) is reminiscent of a pentablock architecture, each block having a specific functional role. Note also, that the central MeMH block is not totally rigid, as the internal hydrogen bonds have been blocked by methylation.

The dissolution properties of UPy-C 6 -MeMH-C 6 -UPy (5) are subtle. The molecule is poorly water-soluble. It dissolves in polar solvents such as dimethyl sulfoxide, dimethyl formamide, methanol, ethanol, and acetone. Surprisingly, UPy-C 6 -MeMH-C 6 -UPy (5) dissolves also in chloroform, toluene, and dichloromethane, despite the dimerizing UPy-groups.

Considering the liquid crystallinity, none of above solvents or their mixtures led to birefringent textures characteristic for lyotropic LCs, even if a wide range of conditions (mixtures of solvents and concentrations) were screened. Upon searching potential solvents, the previously reported lyotropic liquid crystallinity of related polysaccharides turned instructive: Halogenated acids, such as trifluoroacetic acid, allow lyotropic liquid crystallinity of methylcellulose and methylated amylose, whose chemical structures are closely related to the present MeMH. [28][29][30][31] However, the strong acidity of trifluoroacetic acid (pKa = 0.3) and other halogenated acids can cause degradation of the sugars. 1,1,1,3,3,3-Hexafluoro-2-propanol (HFIP) turned out, so far, the only solvent that allowed lyotropic LC of UPy-C 6 -MeMH-C 6 -UPy (5). It is a strong hydrogen bonding donor, due to the strong electron-withdrawing fluorine atoms attached to the α-carbon of the OH group. 32 This has previously allowed dissolution of several polymeric ethers, including methylated sugars and hydrocarbons. 33 However, there are two adverse properties, the safety and especially the high volatility. The high vapor pressure (16 kPa at 20 °C) and the low boiling point (58.2 °C) makes it exceedingly difficult to prevent the solvent from evaporating and maintaining the desired concentrations. This limited severely quantitative characterizations.

The lyotropic liquid crystallinity of UPy-C 6 -MeMH-C 6 -UPy (5) in HFIP was first studied using polarized optical microscopy (POM). The sample was first dissolved in HFIP at a relatively high concentration (500 mg/mL). It was placed between two glass slides, sealed with thermosetting acrylic glue with inherent pinholes to reduce the solvent evaporation rate (Fig. 2a). Special efforts were made to keep the sample and the sealant separate, i.e., to prevent their mixing. The HFIP solution was initially nonbirefringent (Fig. 2b). Upon a slow evaporation of HFIP through the polyacrylate sealant between the glass plates, a birefringent texture started to evolve after ca. 18h (Fig. 2c). The fact that the birefringent texture is really due to liquid crystallinity, is shown in Movie S1 ESI †, which shows flow of the texture upon poking with a spatula. After ca. 5 days, the birefringent LC texture covered the whole sample area (Fig. 2d). After removing the cover slip and evaporating all the remaining solvent, the sample showed some residual birefringence (Fig. 2e).

The POM patterns of Fig. 2d were instructive to suggest the LC structure of UPy-C 6 -MeMH-C 6 -UPy (5) in HFIP. It suggests a spherulitic texture with smooth domains of a fluid columnar mesophase, nucleating from the lower concentration isotropic areas which appear black under crossed polarizers. Worth noticing is the uniaxial symmetry of these defects with bright (orange vs. blue) colored domains separated by sharp (black) isogyres defining large size (up to ca. 100m in diameter) of Maltese crosses which follow the extinction of the analyzer/polarizer (A/P) directions. This reveals a high-quality orthogonal alignment of columnar aggregates radiating from the melatope of the Maltese crosses. Interestingly, the long-Please do not adjust margins Please do not adjust margins range order of the columnar self-assembly remained in the higher concentration spherulitic textures (Fig. 2d). A wave Fig. 2 The sample cell and POM images of UPy-C6-MeMH-C6-UPy 5 in HFIP. a) The sample cell allowing slow evaporation of HFIP. b) Isotropic solution (concentration of ca. 500 mg/mL) immediately after mixing, imaged between crossed polarizers with a POM using a full λ-wave plate. c) Emergence of the first birefringent texture typical of a lyotropic columnar mesophase after ca. 18 hours. d) Birefringent spherulitic texture covering large areas after 5 days, also illustrating the effect of the 90° rotation of a full λ-wave plate. Complete drying leads to amorphous glassy material, with some residual birefringence (e), (f) DSC indicating glass transition at Tg of ca. 74 o C, and (g) broad XRD-peaks suggesting that crystallization does not take place. plate discriminated the columnar aggregates aligned perpendicular and along the wave plate axis by retarding (by one full wavelength) the green component of polarized white light. Figs 2d illustrate the orthogonal orientations of the axes of columnar aggregates which manifest through the inversion of orange vs. blue domains of the Maltese crosses when rotating the full wave plate by 90 °.

Incidentally these POM observations are reminiscent to those observed for the lyotropic LC mesophases of "short" 4-20 base pair DNA and RNA duplexes with "sticky" ends.

34-37

Therein, the observed phases are the chiral nematic, uniaxial columnar, and higher-order columnar. A key feature for this behavior is the end-to-end adhesion into columnar aggregates by the sticky ends. This can only develop if a subtle balance between competing interactions exists. In the present case, the self-assembling UPy dimerize in less polar environments which the C 6 -linkers can provide. In the initial isotropic state with UPy-C 6 -MeMH-C 6 -UPy (5) in HFIP, the UPy end-groups form hydrogen bonds preferentially with HFIP, instead of UPy dimerization. Upon HFIP evaporation, the balance gradually shifts towards the UPy dimerization.

After these qualitative findings, considerable efforts were made to assess the critical concentration for the appearance of the lyotropic LC mesophases. Due to the extremely volatility of HFIP, reliable estimates could not be achieved. Also, X-ray measurements were tried, however, efforts to allow controlled evaporation in time-resolved X-ray scattering measurements turned unsuccessful.

To shed light on the role of the end-group UPy-moieties, the C 6 -linkers, and the central MeMH-moieties on the observed mesomorphic behaviors, reference compounds were studied. Fig. 3a shows maltoheptaose MH with free OH groups (6). It was not soluble in HFIP, but was soluble in water. Permethylated maltoheptaose with OH-end groups, MeMH (2), and permethylated maltoheptaose with hexyl end-groups without self-associating UPy moieties, C 6 -MeMH-C 6 (7), were soluble in the same solvents as UPy-C 6 -MeMH-C 6 -UPy (5). Note that none of these reference samples formed lyotropic LCs in any tested solvent or in any tried concentration and they all crystallized upon drying.

We discuss next the dried UPy-C 6 -MeMH-C 6 -UPy (5). Based on thermogravimetric analysis (ESI †), the DSC measurements were limited to 200 °C. The DSC curve of dried UPy-C 6 -MeMH-C 6 -UPy (5) shows a glass transition at ca. 74 o C (Fig. 2f) and no crystallization/melting transitions. WAXS of the dried sample shows only broad reflections instead sharp crystalline peaks, also indicating amorphous disordered solid with only correlation hole peaks. In conclusion, by contrast to the crystal-forming reference compounds of Fig. 3, (5) showed a tendency for glass formation and residual birefringence upon drying. It showed suppressed crystallization, suggesting kinetic trapping into anisotropic glassy state, indirectly hinting formation of polymer-like chains with free volume. This is not surprising, as the UPy dimerization is expected, when the hydrogen bond donating HFIP has evaporated. Note that the binding constant of the parallel UPy-groups is substantially smaller than that of their end-to-end arrangements, whereupon the UPy dimerization preferably leads to supramolecular polysaccharide instead of their lateral bindings. The FTIR spectra shows the N-H stretching of carbamate at 3730-3750 cm -1 (Fig. S11 in ESI †), which can indicate lateral hydrogen bonds between the carbamates. Importantly, UPy-C 6 -MeMH-C 6 -UPy (5) allowed to draw fibers (Figs 4d) before completely dried, which one would not expect if the supramolecular polymer does not form. However, the fibers turned brittle due to weak lateral connections between the molecules in the case of permethylated maltoheptaose. In this sense, it would be ultimately desirable to allow for a final removal of the lateral methyl-protecting group that was needed to allow the topochemical reactions at

Fig. 1

 1 Fig. 1 The synthesis scheme. a) The reaction of permethylated β-cyclodextrin (TM-CD) (1) with perchloric acid HClO4 to achieve ring-opening and liberate two free OH groups, one in the both ends of the permethylated maltoheptaose (MeMH) (2). b) 2-ureido-4[1H]-pyrimidinone isocyanate with hexyl spacer (UPy-C6-NCO) (3). c) Hexyl isocyanate C6-NCO (4), which was used as a model compound to tune the synthesis conditions, and to synthesize a reference compound for the LC formation. d) The reaction between (2) and (3) to achieve the permethylated maltoheptaose, where both ends have been equipped with UPy 4-H bonding motifs (UPy-C6-MeMH-C6-UPy) (5).

21

 21

1 H

 1 NMR where -stereoselectivity of the carbamate linker was observed as a doublet with a coupling constant of J 1,2 = 8.1 Hz at 5.43 ppm. MALDI-TOF-MS further confirmed the bisfunctionalization where a peak at 1723.89 m/z was detected and assignable to the molecular ion (m/z, [M+Na] + ) (ESI †).

Fig. 3

 3 Fig. 3 Reference compounds, not leading to lyotropic liquid crystallinity (ESI †).

both of the end-group reactions. This issue is currently under investigation and various tracks are explored for alternative protection/deprotection paths.

In conclusion, we synthesized monodisperse permethylated maltoheptaose, where both ends have been functionalized using quadruple H-bonding UPy motifs allowing supramolecular polymerization of this saccharide building block. Both ends were simultaneously reacted with UPy-C 6 -NCO (3) with high -stereoselectivity and efficiency. UPy allows supramolecular polymerization by H-bonding upon solvent removal. We foresee potential for linear supramolecular polysaccharides owing to the renewable and bioactive carbohydrate building-blocks and the ability to undergo reversible switching of structure, morphology, and properties under exposure to certain external stimuli.

Funding from Academy of Finland, Centre of Excellence in Molecular Engineering of Biosynthetic Hybrid Materials, and ERC Advanced Grant are acknowledged. This work made use of Aalto University Nanomicroscopy Center premises. Maria Godinho is acknowledged for useful discussions. OI, TTTM, PR, and SH thank Univ. Grenoble Alpes (UGA) for the partial support under the Project BioTriCo. AG, SH and PR acknowledge Polynat Carnot Institute (No. 16-CARN-025-01), Glyco@Alps project (ANR-15-IDEX-02), ICMG FR 2607 platform, UPR5301-CERMAV (CNRS/UGA), and UMR5819-SyMMES (CNRS/CEA/UGA).

Conflicts of interest

The authors declare no conflict of interests.

Notes and references


Reed Nieman

Riccardo Spezia

Bhumika Jayee

Timothy K Minton
email: tminton@colorado.edu

William L Hase

Hua Guo
email: hguo@unm.edu

H J Smead

Exploring Reactivity and Product Formation in N( 4 S) Collisions with Pristine and Defected Graphene with Direct Dynamics Simulations

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( 4 S)) collisions with pristine, defected, and oxidized graphene. Apart from non-reactive scattering of nitrogen atoms, various forms of nitridation of the 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, hightemperature 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 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][7][8][9][10][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 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 -1373 K) exposed to N atoms showed reaction probabilities of 0.2 -9.8 × 10 -3 . 15 To obtain This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. PLEASE CITE THIS ARTICLE AS DOI:10.1063/5.0028253 information on the mechanism of the mass loss, products resulting from such reactions have also been detected using various techniques. [16][17][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 et al. studied the scattering of hyperthermal N( 4 S) atoms and N2 molecules with vitreous carbon at high temperatures between 1023 and 1923 K. 10 The molecular beam contained N and N2 with translational energies of 110 and 193 kcal/mol, respectively. They found that N2 scattered nonreactively while some of the N( 4 S) 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 the 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 This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. 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, 46, 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.

PLEASE CITE THIS ARTICLE AS

In this work, we report direct dynamics studies of the reactive scattering of hyperthermal atomic nitrogen in the ground electronic state (N( 4 S)) 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 the simulation of not only 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 experiments 10,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) include a pristine, nondefected structure (Model 1, Figure 1a), a single carbon vacancy (SV) (Model 2, Figure 1b), and two oxygenated structures based on the SV surface (Models 3 and 4, Figure 1c andd, respectively).

The graphene surfaces are modeled with periacene and its modified analogs 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. This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS

II. Methods

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 materials 29 to biological systems. 30 a graphene surface. 27 In our calculations, the third-order correction 31,32 was utilized with the 3ob-3-1 Slater-Koster parameter set 33 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 threeparameter Becke-Lee-Yang-Parr (B3LYP) functional, 35,36 with the 6-31G* basis set. 37 In addition, the D3 dispersion correction 38 with Becke-Johnson damping (D3BJ) 39 was used to handle the dispersion corrections. Fermi smearing was utilized to allow for quick convergence of the selfconsistent 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

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][43][44][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][48][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 program 50 externally interfaced to the DFTB+ program. Namely, the atomic positions of the whole system provided by VENUS were read by 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( 4 S) 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 Figure 2a. The initial position and momentum of the N( 4 S) 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 (Figure 2b). 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( 4 S) 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( 4 S) 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 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( 4 S) atom projectile was chosen randomly between 0° to 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.

This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. All trajectories were integrated in time with the 6 th 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 nonreactive 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.

PLEASE CITE THIS ARTICLE AS

III. Results

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 (Figure S1), suggesting that the DFTB method provides an excellent geometric description of the model graphenes.

This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI:10.1063/5.0028253

For the interaction of an N( 4 S) 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 and -14.2 kcal/mol, respectively, with respect to the dissociation limit (Figure S2 of the Supplementary Information (SI)), respectively. The minimum features an sp 3 -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 (Figure 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 Figure 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 (Figure S5). This comparison further suggests that DFTB underestimates the binding of the 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 (Figure S6). No potential barrier is observed along the PEC. Selective structures from the PEC in Figure 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 (Figure 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 Å (Figure S9). This compares well to the literature value for a CN triple bond of 1.157 Å. [START_REF] Herzberg | Electronic Spectra of Polyatomic Molecules[END_REF] This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI:10.1063/5.0028253

Finally, the PEC for NO leaving Model 4 after an impinging N atom reacts with the carbonyl oxygen is found in Figure S10. A 17.9 kcal/mol barrier is found at C-ON distance of 1.9 Å, while the potential minimum is located above the complete dissociation limit by 7. Based on these comparisons, we conclude that DFTB provides qualitatively correct interaction patterns for nitrogen interaction with graphene but underestimates the interaction energy.

Reactive Scattering Model 1 -Pristine Graphene Sheet

Most trajectories in Model 1 with a translational energy of 14.9 kcal/mol (Figure 3, red) result in the nitrogen atom rebounding from the surface without reacting. The velocity distribution for the scattered nitrogen atoms is shown in Figure 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 the 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 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 Figure 3. In addition to five such reactive trajectories, there are eleven 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.

This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. Increasing the translational energy of the N( 4 S) atom to 110 kcal/mol is found to enhance the reactivity (Figure 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 are 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( 4 S) atom impacting the pristine graphene sheet. At both incidence energies, the nitrogen atom is not observed to undergo significant diffusion on the graphene surface.

PLEASE CITE THIS ARTICLE AS

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. [START_REF] Nieman | [END_REF]56 The probability of reaction at low incidence energy (14.9 kcal/mol, Figure 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.

This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. 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 (Figure 5c). 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 (Figure 5a). As expected, the probability of this process decreases with b thereafter.

PLEASE CITE THIS ARTICLE AS

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 (Figure 5b). In addition, N-atom insertions at the defect site can also result in Stone-Wales (SW) (Figure 5d), or square (SQ) configurations (Figure 5e). 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 (Figure 5f). This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. 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.

PLEASE CITE THIS ARTICLE AS

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( 4 S) atom, as illustrated in Figure 6 (red and black columns for incidence energies of 14.9 and 110 kcal/mol, respectively). At the lower incidence energy, reactivity of the sheet towards 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 Å).

This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. 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 (Figure 6a) and a doubly bonded CN group (Figure 6b). 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 Figure 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 This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. 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 Å.

PLEASE CITE THIS ARTICLE AS

PLEASE CITE THIS ARTICLE AS

Increasing the incidence energy of the impinging N( 4 S) atom to 110 kcal/mol increases Pr(b) significantly, especially for trajectories with b > 0.0 Å (Figure 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 (Figure 6). The major configurations are presented in Figure 6a-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.

Model 4 -SV Graphene with Two Oxygen Atoms

Sampling Model 4 with N( 4 S) atoms with 14.9 kcal/mol of incidence energy, shows a rather Several gaseous products were observed to be ejected from the graphene sheet after the impact of the nitrogen atom projectile (Figure 7, bottom left). Especially at values of b = 0.0, 0.5, and 1.0 Å, the increase in reaction probability primarily from the production of NO (Figure 7c). 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 (Figure 7d) 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 (Figure 7e). 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 the 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 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 temperature 11 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 Natom 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 energies characteristic of impulsive, or direct, scattering. 10 The experiment with 8.2 kcal/mol N atoms incident on a vitreous carbon surface found that nonreactive 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 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 the graphene, as shown in Figure S6, with a binding energy of 104.7 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 experiment 11 could be due to thermal desorption of CN. However, the current calculated This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. Our simulations further predicted several other gas phase molecules, such as NO, CO, NCO and O, that have not been observed 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.

PLEASE CITE THIS ARTICLE AS

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

  DOI:10.1063/5.0028253

  DOI:10.1063/5.0028253

Figure 1 :

 1 Figure 1: 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.

  Indeed, it was also used in the aforementioned work by Moon et al. in a direct dynamics simulation of N-atom scattering from This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. PLEASE CITE THIS ARTICLE AS DOI:10.1063/5.0028253

Figure 2 :

 2 Figure 2: 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).

  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 (Figure S11) goes back into the plane as NO leaves forming the Model 3, single oxygen insertion defect.

Figure 3 :

 3 Figure 3: Probability of reaction, Pr(b), with respect to impact parameter, b, from DFTB3-LJ/3ob-3-1 direct dynamics simulations of Model 1.

Figure 4 :

 4 Figure 4: Velocity distribution of scattered N( 4 S) atoms from DFTB3-LJ/3ob-3-1 direct dynamics simulations of Model 1 with N( 4 S) translational energy set to 14.9 kcal/mol. Maxwell-Boltzmann distributions shown for 1923 and 4750 K.

Figure 5 :

 5 Figure 5: Top left, probability of reaction, Pr(b), with respect to impact parameter, b, from direct dynamics simulations of Model 2 with N( 4 S) translational energy set to 14.9 (red) and 110.0 kcal/mol (black). Probabilities for the two primary N-functionalization products of the 14.9 (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.

  DOI:10.1063/5.0028253 Increasing the incidence translational energy to 110.0 kcal/mol decreases the probability of reaction for b values = 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.

Figure 6 :

 6 Figure 6: Top left, probability of reaction, Pr(b), with respect to impact parameter, b, from direct dynamics simulations of Model 3 with N( 4 S) translational energy set to 14.9 (red) and 110.0 kcal/mol (black). Probabilities for the two primary N-functionalization products of the 14.9 (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.

  DOI:10.1063/5.0028253

different reactivity profile compared to the reactions with Model 3 (

 3 Figure 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 (Figure 7 middle left) are mostly N insertions into a C-C bond (Figure 7b) or pulling a carbon atom out of the plane as a singlybonded CN group (Figure 7a) 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.

Figure 7 :

 7 Figure 7: Top left, probability of reaction, Pr(b), with respect to impact parameter, b, from direct dynamics simulations of Model 4 with N( 4 S) translational energy set to 14.9 (red) and 110.0 kcal/mol (black). Probabilities for the two primary N-functionalization products of the 14.9 (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 (Figure 7, red). While the carbonyl oxygen still directs the incoming N( 4 S) 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

  DOI:10.1063/5.0028253 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 CN bond in different defect sites of the graphene.

This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.PLEASE CITE THIS ARTICLE ASDOI:10.1063/5.0028253

Acknowledgments:

This work was supported by the Air Force Office of Scientific Research (Grant No. FA9550-15-1-0305 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 statement:

The data that support the findings of this study are available from the corresponding author upon reasonable request.

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 Materials: Comparison of results calculated from DFT and DFTB methods.