Ab initio molecular dynamics of atomic-scale surface reactions: insights into metal organic chemical vapor deposition of AlN on graphene

Metal organic chemical vapor deposition (MOCVD) of group III nitrides on graphene heterostructures oﬀers new opportunities for the development of flexible optoelectronic devices and for the stabilization of conceptually-new two-dimensional materials. However, the MOCVD of group III nitrides is regulated by an intricate interplay of gas-phase and surface reactions that are beyond the resolution of experimental techniques. We use density-functional ab initio molecular dynamics (AIMD) with van der Waals corrections to identify atomistic pathways and associated electronic mechanisms driving precursor/surface reactions during metal organic vapor phase epitaxy at elevated temperatures of aluminum nitride on graphene, considered here as model case study. The results presented provide plausible interpretations of atomistic and electronic processes responsible for delivery of Al, C adatoms, and C–Al, CH x , AlNH 2 admolecules on pristine graphene via precursor/surface reactions. In addition, the simulations reveal C adatom permeation across defect-free graphene, as well as exchange of C monomers with graphene carbon atoms, for which we obtain rates of B 0.3 THz at typical experimental temperatures (1500 K), and extract activation energies E exca = 0.28 (cid:2) 0.13 eV and attempt frequencies A exc = 2.1 ( (cid:3) 1.7 (cid:2) 1 ) THz via Arrhenius linear regression. The results demonstrate that AIMD simulations enable understanding complex precursor/surface reaction mechanisms, and thus propose AIMD to become an indispensable routine prediction-tool toward more effective exploitation of chemical precursors and better control of MOCVD processes during synthesis of functional materials.


Introduction
Metal organic chemical vapor deposition (MOCVD) is acknowledged as the prime rational and scalable deposition process for obtaining (ultra)thin films and heterostructures of semiconductor materials applied in main areas of established and emerging new technologies. MOCVD of heterostructures of group III nitrides (AlN, GaN, InN) with graphene offers new opportunities for the development of flexible optoelectronics (e.g. InGaN/GaN/graphene, ref. 1) and for the stabilization of conceptually-novel group III nitrides (see, e.g., synthesis of two-dimensional GaN via graphene encapsulation 2 ). Recent ab initio investigations, motivated by the perspectives for highimpact applications such as room-temperature spintronics at the nanoscale, 3 reported on the structural, vibrational, and electronic properties of single and few-layer graphitic-like AlN, 4 including in heterostructures with graphene. 5 The MOCVD of AlN on epitaxial graphene has also been attempted. 6 MOCVD processes are regulated by a complex interplay between precursor/precursor and precursor/surface reaction thermodynamics and kinetics. 7 In this respect, the MOCVD of AlN represents the most challenging case among group III nitrides (including BN) 8 and group III-V semiconductor materials (e.g., AlGaInAsP, see ref. 9) in general. MOCVD of highquality AlN films requires elevated process temperatures, B1500-1700 K, 10 considering the high Al-N bond energy of 66 kcal mol À1 (compared to, e.g., the Ga-N bond strength, 45 kcal mol À1 ) 11 and the need for sufficiently rapid surface diffusion. 12,13 The MOCVD of AlN is further complicated by the fact that typical precursors -trimethylaluminum, (CH 3 ) 3 Al, and ammonia, NH 3 -form highly-reactive (CH 3 ) 3 Al:NH 3 adducts, subsequently triggering intricate precursor/precursor and precursor/ surface reactions. 14 Atomistic reaction pathways are not readily determinable via MOCVD experiments. 7 However, ab initio calculations have substantially contributed to the understanding of MOCVD processes 8,9,14 by predicting the energetics of gas and surface reactions at 0 Kelvin. [15][16][17][18] The effects of pressure and temperature on reaction rates can be estimated via thermodynamic and kinetic models 19 or, whenever possible, combining ab initio calculations with classical molecular dynamics (CMD) based on empirical potentials. 20,21 Nevertheless, static 0 K ab initio techniques such as density functional theory (DFT), or even more accurate first-principles electronic-structure-based methods, 22,23 are inherently unable to resolve the intricate molecule/surface dynamics, especially when several degrees of freedom are involved. On another hand, efficient deterministic CMD simulations present disadvantages such as questionable reliability/ transferability and require fitting a large set of potentialparameters for describing the large variety of atomic interactions.
Density-functional ab initio molecular dynamics (AIMD) 24 yields reliable descriptions of the time-evolution of solidstate systems directly embedded in realistic environmental conditions. 25,26 AIMD simulations often reveal non-intuitive reaction pathways 27,28 and system configurations 29 while providing corresponding kinetic rates at finite temperatures. Despite its known limitations for treating electron transfer during chemical reactions and/or the energetics of transition states for molecule splitting, 30 AIMD is an excellent tradeoff between accuracy + reliability vs. computational cost for identifying chemical reaction pathways and estimating the relative occurrence of competing reactions. To this point, AIMD has been proven useful computational tool for investigating reactions of graphene with different functional groups including N-containing amino radicals. 26 The purpose of this work is to identify reaction pathways relevant for MOCVD processes of AlN on graphene and interpret the reactions on the base of associated variations in electron densities. AlN on graphene serves as a model to rationalize atomistic and electronic mechanisms of importance to deposition of group III nitrides on graphene, in general. Our results demonstrate the great potential of AIMD simulations for revealing precursor/surface dynamics and for guiding experiments toward the realization of prototype group III nitrides heterostructures with graphene, already attempted by MOCVD. 1,2

Results and discussion
The MOCVD of AlN implements trimethylaluminum, (CH 3 ) 3 Al, and ammonia, NH 3 , as typical group III and group V precursors, respectively. It is known that a chain of gas-phase chemical reactions is opened by an easy dimerization of the adductderived species (CH 3 ) 2 AlNH 2 following the rapid methane (CH 4 ) elimination from the adduct (CH 3 ) 3 Al:NH 3 . 14 AIMD simulations were carried out to reveal how each of the individual precursorsammonia NH 3 , and trimethylaluminum (CH 3 ) 3 Al -and the adduct-derived species (CH 3 ) 2 AlNH 2 with direct Al-N bonding, evolve/react with a graphene sheet. The interaction of each precursor with graphene is treated in separate sets of simulations. In order to observe different relevant and possibly competing reaction mechanisms, we also carry out two, or more, separate and independent AIMD runs for each of the investigated precursor/graphene reactions. The description of AIMD results is focused on reactions which we retain of primary relevance for understanding MOCVD of AlN on graphene. The equations describing chemical reactions are based on our empirical interpretation of events observed in AIMD videos and schematically represented in the way we believe these to be more likely to occur. Duration of precursor/graphene physisorption events and frequency of precursor/graphene collisions prior to molecule transformations allow for qualitatively comparing the reactivity of different precursors with the graphene sheet.
Electron transfer maps are used to identify the variations in charge distribution associated with molecule/molecule and molecule/surface interactions, or any reaction of interest. These are calculated by subtracting the DFT self-consistent electron densities of individual molecules and/or the graphene sheet (by using the atomic positions generated during AIMD) from the DFT self-consistent charge density of the entire system.

Simulations of trimethylaluminum (CH 3 ) 3 Al with graphene
Run 1: Al adatom formation via CH 3 -AlQ Q QCH 2 dissociation on graphene. In this case, the system comprised by the (CH 3 ) 3 Al molecule and the graphene sheet evolves by sequential collisions of the molecule with the carbon layer. In most cases, a collision leads to temporary physisorption timeframes (of duration ranging between B2 and B6 ps) with the (CH 3 ) 3 Al molecule rolling and translating on the graphene surface. The fifth molecule/surface collision, at a simulation-time (t sim ) of approximately 29 ps, excites scissoring/bending vibrational modes of Al-methyl intramolecular bonds. This assists proton transfer from one methyl group to the other [ Fig. 1(a and b)] and leads to methane elimination from the trimethylaluminum (CH 3 ) 3 Al molecule (eqn (1) and Fig. 1c). The products of this reaction is a methane molecule and a CH 3 -AlQCH 2 molecule: (1) The formation of a CH 3 -AlQCH 2 molecule is followed by a rapid sequence of reactions. At t sim = 37 ps, the CH 3 -AlQCH 2 molecule collides for the sixth time with, and adsorbs onto, the graphene sheet by its Al-atom side (Fig. 2a). This can be explained by the fact that Al has an empty orbital of sp 2character (eqn (1)), which acts as electron acceptor from the p-cloud overlying the graphene sheet. The effect is evidenced by the electron accumulation generated between the Al atom with the two closest graphene C atoms (see electron-transfer density cut on the C 1 G -Al-C 2 G plane in Fig. 2b). Remarkably, attachment to the graphene sheet rapidly leads to CH 3 -AlQCH 2 dissociation (at t sim = 39 ps) caused by an electrostatic repulsion between the graphene p-cloud and the electrons of intramolecular Al-methyl bonds (Fig. 2c). This repulsion induces CH 3 -AlQCH 2 bending, which brings the C atoms of the CH 3 and CH 2 groups closer to each other (see Fig. 3a and electron accumulation between C A and C B atoms in Fig. 3b). The subsequent rapid formation of a dehydrogenated ethane molecule CH 3 -: CH 2 leaves an Al adatom (Al ad ) on graphene (Fig. 3c). 31 Finally, at t sim = 42 ps, Al ad /CH 3 -: CH 2 recombination is observed. The delivering of Al ad onto graphene is an insightful It is noteworthy that this simulation run indicates a straightforward mechanism for supply of Al atoms to the graphene surface. The Al adatom remains physisorbed for as long as it does not interact with other gas species resulting from the (CH 3 ) 3 Al/graphene reactions. Additional AIMD runs that we have carried out at temperatures ranging from 300 to 4300 K, each with duration of 0.1 ns, confirm that individual Al adatoms are stable species, unlikely to spontaneously desorb from graphene (see Appendix).
For Al adatom diffusivities on graphene, our simulations yield diffusion coefficient D(T) results that follow two distinct Arrhenius-like trends: D(T) = [5.6 Â 10 À4 cm 2 s À1 ] Â exp[À0.03 eV/(k B T)] (Â3 AE1 ) for 300 K r T r 1200 K; D(T) = [1.1 Â 10 À2 cm 2 s À1 ] Â exp[À0.34 eV/(k B T)] (Â3 AE1 ) for 1200 K r T r 4300 K. The high lifetime and diffusivity of Al adatoms on graphene might be important for promoting AlN nucleation, with the Al adatom itself acting as initial seed, whereby high surface mobility is considered as a factor that improves growth efficiency. 32 Alternatively, adsorbed Al atoms may intercalate at graphene/SiC interfaces; a mechanism analogous to that  leading to formation of two-dimensional GaN via graphene encapsulation. 2 Dissociative patterns for metal organic precursors at MOCVD conditions, as well as delivery of metal atoms to the growing surface -similarly to the reactions observed in this work -have been predicted via DFT and discussed for, e.g., trimethylgallium and trimethylindium on GaN and InN surfaces, respectively. 18,19 However, although 0 K ab initio calculations provide reliable evaluation of reaction pathways and Gibbs free energies of reactants and products, they are unable to determine the intricate atomistic dynamics underlying MOCVD processes.
Run 2: Al adatom formation catalyzed by C adatoms and carbon mass transport on, and across, graphene. Similar to the first simulation run, during the initial 35 ps of the simulation (prior to the first (CH 3 ) 3 Al/graphene reaction), the trimethylaluminum (CH 3 ) 3 Al molecule collides 8 times with the graphene surface, where exhibits periods of physisorption which amount to a total of B10 ps. We observe methane elimination (see eqn (1) as described in the previous section (Run 1)), which leads to formation of a CH 3 -AlQCH 2 molecule directly after the eighth (CH 3 ) 3 Al/graphene collision.
During a timeframe of 15 ps following the methane elimination, the CH 3 -Al QCH 2 molecule collides 6 times with the graphene sheet, with interposed physisorption periods lasting for approximately 2 ps each. Then, the CH 3 -Al QCH 2 molecule reacts with graphene as illustrated in Fig. 4. First, the CH 3 -Al QCH 2 molecule bonds with the graphene carbon atom C 1 G via its Al-atom side, while the C atom of the CH 2 group (C A ) attaches onto a graphene carbon atom (C 2 G ) adjacent to C 1 G (Fig. 4a). After 0.2 ps, the C A /Al bond breaks, and the CH 2 binds even more strongly to the underlying C 2 G atom (Fig. 4b). The average C A -C 2 G bond length of 1.45 Å is comparable with the nearest neighbor C G -C G distance in graphene. Approximately 0.1 ps later, the Al-CH 3 molecule desorbs, while the CH 2 group remains as admolecule on graphene (Fig. 4c). The CH 2 admolecule adsorbs preferentially atop C G atoms and migrates among these sites, through bridge positions, above C G -C G bond centers, with a migration frequency of 4.3 (Â1.4 AE1 ) THz at 4300 K. After B5.5 ps, we observe the formation of a C adatom consequent to dehydrogenation of the CH 2 admolecule, as illustrated in the snapshot sequence of Fig. 5. The formation of C adatoms has been reported to significantly enhance the reactivity of graphene by forming localized states in the vicinity of the adspecies, similar to the effects induced by defects in the graphene sheet. 33 In our simulations, the carbon adatom (C ad ) diffuses among atop-C G sites via the center of C G -C G bonds with jump frequency of 10.0 (Â1.5 AE1 ) THz. Carbon adatoms located above C G -C G bonds cause the temporary formation of C 7 rings, with the adatom embedded between two C G atoms, and ripple the graphene sheet with a cusp in correspondence of the adatom position. Approximately 1 ps later, the C ad and the underlying C G atom undergo an exchange reaction, as shown in Fig. 6a and b. The C adatom (cyan sphere) moves inward the graphene layer pushing the underlying C G atom (gray sphere) to the other side of the sheet.
The exchange reaction occurs while the AlCH 3 molecule (which had previously formed, see Fig. 4c) is approaching the graphene layer from the gas phase (Fig. 6c). The simulation    molecule due to C ad /Al attraction which leads to Al-C ad bond formation and elimination of the methyl group from the gas molecule (Fig. 6d). The reaction leaves a C-Al admolecule anchored to graphene by its C atom. However, the C-Al admolecule rapidly splits, with the Al atom adsorbing on graphene and migrating away from the carbon adatom. Thus, in the second (CH 3 ) 3 Al/graphene simulation, the formation of a free Al adatom is catalyzed by the presence of a C adatom on the graphene surface. This is an alternative atomistic mechanism for supplying graphene with Al adatoms. However, this pathway is considerably less likely to occur (entropically unfavored) than the direct attachment of CH 3 -AlQCH 2 molecules on graphene (see Run 1).
During the subsequent B13 ps of this simulation, frequent C ad /C G exchange reactions (eight events for 13 ps) as well as permeation of a C adatom through graphene (two events) are observed. This yields an overall rate of 0.75 (Â1.5 AE1 ) THz for C atom transport across pristine graphene at 4300 K. The presence of a carbon adatom locally weakens C G -C G interatomic bonds and promotes bond stretching at elevated temperatures (note C 2 G -C 4 G bond length 42 Å in Fig. 7a and b), thus favoring C permeation across an in-plane threefold-coordinated position ( Fig. 7c and d). The reaction ends with the C becoming an adatom on the opposite site of graphene sheet, Fig. 7(e-h). DFT modeling of C atom penetration though large aromatic C x H y planar molecules yields extremely high energy barriers (B7 eV). 34 AIMD simulations of this work describe C exchange (Fig. 6) and permeation (Fig. 7) phenomena in a more realistic way, providing global activation energies for C transport across pristine graphene of only 0.28 eV (see below).

Carbon adatom dynamics on graphene: surface migration and exchange reactions
The observation of the phenomenon of carbon transport through defect-free graphene motivated additional AIMD simulations (270 ps in total) to investigate more in detail the dynamics of carbon adatoms on graphene at different temperatures. Previously, C adatom migration on graphene and on carbon nanotubes (both internally and externally) has been studied via ab initio and via Monte Carlo simulations. [35][36][37] Our AIMD results show that both C ad /C G exchange and C ad migration events are extremely frequent (with respect to the timescales feasibly accessible to AIMD) at all investigated temperatures (between 1000 and 4100 K) and follow a fairly Arrhenius-like behavior. Adsorption of isolated carbon atoms on graphene can be imaged via conventional transmission electron microscopy. 38 To our knowledge, however, theoretical investigations of C ad /C G exchange reactions had not been carried out.
Simulations at relatively low temperatures (1000 and 1500 K), allow for easily identifying the transition state configuration for C ad /C G exchange (Fig. 8), with the two carbon atoms occupying specular positions with respect to the sheet. The dumbbell-like C ad /C G configuration entails that both the carbon adatom and its twin C G atom are bonded to each other (bond length of E1.7 Å) as well as to three neighboring graphene atoms with bond lengths of E1.5 Å; slightly larger than the calculated average C G -C G nearest neighbor distance (E1.45 Å) at 1000 K.
Linear regression of C ad /graphene AIMD migration rates [ln(rate)] vs. 1/T yields attempt frequencies A mig = 24.4 (Â1.3 AE1 ) THz and activation energies E mig a = 0.40 AE 0.09 eV, consistent with C adatom energy barriers (E0.55 eV) computed via DFT at 0 K. 35 For C transport across graphene, we obtain E exc a = 0.28 AE 0.13 eV and A exc = 2.1 (Â1.7 AE1 ) THz. It should be noted that, although carbon transport across graphene requires lower activation energies, the rate of these reactions is smaller than that for C ad surface migration at all investigated temperatures ( Fig. 9) due to A exc attempt frequencies being much smaller than A mig . C ad /C G exchange rates at typical experimental temperatures (1500 K) 10 are B0.3 THz. This result indicates that formation of C adatoms during MOCVD may significantly affect reactivity and stability of defect-free graphene areas and consequently influence the deposition of AlN on graphene, here considered as case study. We note that exchange of C adatoms with graphene atoms is never observed to disrupt the honeycomb graphene structure during our AIMD simulations. The presence of carbon adatoms and the occurrence of C ad /C G exchange reactions may contribute to weaken graphene C-C bonds vicinal to existing defects, e.g., Stone-Wales, pentagonal, and double pentagon defects, and assist formation of holes in the bonding network, as previously shown for graphene-like carbon-based systems. 39 Consequently, the increased reactivity of unsaturated bonds in the carbon sheet would lead, favored by attachment of gas radical molecules, to partial or full disintegration of graphene.

Simulations of (CH 3 ) 2 AlNH 2 with graphene
As indicated, the general case of MOCVD of AlN could be considerably complicated and dominated by rapid ({1 s) irreversible formation of (CH 3 ) 3 Al:NH 3 adducts, which nearly proceeds at gas-kinetic collision rates, with no activation energy upon mixing of the precursors trimethylaluminum (CH 3 ) 3 Al and ammonia NH 3 in the gas phase. 40 This is followed by a rapid methane (CH 4 ) elimination from the adduct (CH 3 ) 3 Al:NH 3 , which in turn initiates the formation of adduct-derived species such as (CH 3 ) 2 AlNH 2 with direct Al-N bonding 14 (3) For sake of completeness, we present below (CH 3 ) 2 AlNH 2 / graphene simulations (alternatively written as H 2 N:Al(CH 3 ) 2 / graphene).
Run 1: Al-NH 2 and CH 3 admolecule formation. Analogous to the behavior observed for trimethylaluminum (CH 3 ) 3 Al precursors, the adduct-derived species H 2 N:Al(CH 3 ) 2 exhibits relatively high affinity to graphene. The initial two collisions of H 2 N:Al(CH 3 ) 2 with graphene result in relatively long physisorption periods of 3 and 6 ps, respectively, with H 2 N:Al(CH 3 ) 2 rototranslating on the carbon sheet. The third collision (t sim = 16 ps) rapidly leads to detachment of a : CH 3 radical: which flies away from the surface. The reaction is enabled by the high vibrational energy at 4300 K. 41 The remnant H 2 N: : AlCH 3 radical physisorbes at the graphene sheet, on which it remains for approximately 3 ps (up to t sim = 19 ps). Thus, H 2 N: : AlCH 3 loses the other methyl radical group : CH 3 , which attaches onto graphene for the subsequent 2 ps while the AlNH 2 species desorbs from the surface (Fig. 10).
At a simulation time t sim = 24 ps, the gaseous AlNH 2 species absorbs at the graphene surface by its Al side for 3 ps (figure not shown). This reaction supplies the graphene surface with an (already formed) Al-N bond. It should be noted that, also in the case of AlNH 2 physisorption, graphene exhibits strong affinity for the Al atom of the adspecies.
Run 2: delivery of Al adatoms to the graphene sheet. The second simulation of H 2 N:Al(CH 3 ) 2 with graphene reveals a significantly different sequence of events in comparison to the first run. Ten H 2 N:Al(CH 3 ) 2 /graphene collisions occur before the H 2 N:Al(CH 3 ) 2 species initiates a rapid sequence of chemical reactions. At the 11th collision (t sim = 38 ps), the H 2 N:Al(CH 3 ) 2 species undergoes methane elimination activated by the two methyl groups coming close enough to assist proton transfer from one methyl to the other: This differs from the observation of the first H 2 N:Al(CH 3 ) 2 / graphene simulation, where a reactive CH 3 methyl radical detaches from the adduct at the third H 2 N:Al(CH 3 ) 2 /graphene collision. The H 2 N:AlCH 2 species attaches to the graphene surface by its Al-atom side. Very few fs later, the NH 2 and CH 2 groups bond to each other, while breaking their bond with the Al atom, and the H 2 NCH 2 species desorbs from the surface (the mechanism is analogous to the one illustrated in Fig. 3 for delivery of Al adatoms by (CH 3 ) 3 Al/graphene reactions), while the Al atom remains physisorbed at the surface for the remaining simulation time (from t sim = 38 to 41 ps). As in the case of (CH 3 ) 3 Al/graphene reactions, also this simulation suggests a different reaction of delivery of isolated Al adatoms on graphene: H 2 N:AlCH 2 -H 2 NQCH 2 + Al ad (6) During the final 3 ps of the simulation, the H 2 NCH 2 gas species undergoes several reactions of dehydrogenation and recombination with hydrogen. However, AIMD results indicate that H x N-CH y molecules are not likely (at least at elevated temperatures) to adsorb on graphene.
Simulations of ammonia NH 3 with graphene NH 3 is a strongly-bonded molecule: the N-H bond energy calculated via DFT is B420 kcal mol À1 , while the Al-N and Al-C bond energies in (CH 3 ) 3 Al, the adduct, and their derivatives vary between 40 and 120 kcal mol À1 . 42 Thus, the decomposition of NH 3 is expected to occur in the gas phase at very high temperatures (typically above 1800 K), or on the graphene surface in presence of other species or at defect sites. 43 This is consistent with our AIMD results, as explained below.
AIMD simulations show 64 NH 3 /graphene collisions during 100 ps (Fig. 11). Although ammonia exhibits physisorption periods on graphene (B3 ps in total), no reaction is recorded in the gas phase or on the surface, despite the elevated temperature. This is in direct contrast to what observed (see above) for (CH 3 ) 3 Al and H 2 N:Al(CH 3 ) 2 reactions, which occur after very few collisions with the graphene sheet (initial precursor reactions take place within 40 ps).
Overall, AIMD results show a considerable greater affinity of Al-containing species as well as of methyl groups toward graphene in comparison with NH 3 . While Al has an empty orbital, prone to accept graphene p electrons, the N atom has a completely filled electronic shell that repels the graphene p-cloud and thus prevents ammonia attachment on graphene. Formation of Al and C adatoms consequent to (CH 3 ) 3 Al/graphene and H 2 N:Al(CH 3 ) 2 /graphene collisions provides surface sites of higher reactivity which can assist subsequent gas-molecule/ graphene reactions.
Ammonia is known to functionalize graphene at defective sites at temperatures between 1000-1300 K. 43 Molecular dynamics simulations performed using quantum-chemistry based forcefields showed that excess of nitrogen impurities and the presence of N-N bonds within the graphene layer undermine graphene stability. 44 However, in the case of MOCVD of AlN on defect-free (high-quality) graphene, treatment of graphene with ammonia prior to initiating MOCVD is not expected to provide significant differences in surface reactivity. Nevertheless, Al-N bond splitting in H 2 N:Al(CH 3 ) 2 or ammonia/atomic-H reactions -H 2 dissociation is relevant at the typical synthesis temperatures, B1700 K, used for growth of high-quality AlN thin films 10 -would  promote formation of : NH 2 radicals which are expected to exhibit considerably higher affinity toward graphene. The presence of an ammine group on defect-free regions of the graphene sheet was indicated to assist various reactions, including the possibility of nitrogen-atom delivery to the graphene surface. 26

Conclusions
In summary, density-functional molecular dynamics simulations with van der Waals corrections have been employed for achieving deeper understanding at atomistic level of the rich chemistry of reactions occurring during MOCVD of AlN on graphene at elevated temperatures. Trimethylaluminum (CH 3 ) 3 Al, ammonia NH 3 , and an adduct-derived species (CH 3 ) 2 AlNH 2 have been employed as chemical precursors. The simulations reveal atomistic pathways, with associated changes in electron densities, responsible for reactions in the gas phase (e.g. methane elimination, dehydrogenation) as well as for delivery of Al, C adatoms, and C-Al, CH x , Al-NH 2 admolecules on defect-free graphene consequent to precursor/surface collisions. DFT charge-density results suggest that attachment and reactivity of Al-containing molecules on graphene is relatively rapidly activated owing to facile graphene p-electron transfer to empty Al orbitals. In contrast, NH 3 precursors are never observed to react with the carbon layer, thus indicating that pre-treatment of defect-free graphene via exposure to ammonia does not significantly contribute to N-atom delivery or creation of reactive sites on the carbon sheet. The simulations show that individual Al adatoms exhibit considerable stability and high diffusivity on graphene, and might therefore contribute to increase surface reactivity, promote AlN nucleation, or may rapidly intercalate at grain boundaries. Irrespective of temperature, carbon adatoms undergo frequent exchange reactions with graphene carbon atoms, thus leading to C transfer across the sheet with rates G = 2.1 (Â1.7 AE1 ) THzÁexp[(À0.28 AE 0.13 eV)/(k B T)]. On the base of the results presented, we argue that density-functional molecular dynamics is an efficient method to unveil precursor/surface chemical reactions and thus to guide MOCVD experiments toward a more effective exploitation of gas precursors.

Computational methods
Ab initio molecular dynamics simulations are carried out with the VASP code 45 using the local density approximation (LDA) within the scheme of Ceperley-Alder, 46 as parameterized by Perdew and Zunger, 47 and the projector augmented wave 48 method. The approximation proposed by Grimme 49 is used to model the non-locality of electron correlation. At each AIMD time step (0.1 fs for modeling dynamics in all simulations involving hydrogen atoms and 1 fs for studying diffusion of isolated Al and C adatoms on graphene), the total energy is evaluated to an accuracy of 10 À5 eV per atom employing a plane-wave energy cutoff of 300 eV and G-point sampling of the Brillouin zone. Electronic thermal excitations are accounted for by employing Gaussian smearing energies equal to k B T.
The present AIMD simulations employ boxes with periodic boundary conditions in all directions. The supercells contain 27 Å-thick vacuum regions separating graphene-sheet replicas along the vertical direction. The graphene layer consists of 6 Â 6 unit cells, for a total of 72 carbon atoms. The AIMD simulations are carried out at a temperature of T = 4300 K, which is below the melting point of graphene (B4500 K). 50 The use of a supercell with relatively small graphene surface area, (B13.8 Å 2 ), prevents populating long-wavelength (B80 Å) phonons which cause formation of defects from T Z 3500 K. 51 The choice of a simulation temperature (4300 K) well above the typical temperature implemented in the MOCVD of AlN (T E 1500-1700 K), 10 dictated by the fact that defect-free graphene exhibits extremely low reactivity at the timescales accessible to AIMD, 43 significantly reduces the simulation times required to obtain qualitative understanding of surface reaction pathways. Chemical reactions occur considerably more rapidly at defective sites than on defect-free graphene regions. 43 Nevertheless, the investigation of the inherent reactivity of pristine graphene is a necessary initial step for understanding MOCVD growth on its surface.
In all our simulations, the graphene lattice parameter a is maintained fixed to 2.465 Å. This value corresponds to that determined at 2400 K via atomistic Monte Carlo simulations 52 based on reliable bond-order interatomic potentials for carbon. 53,54 Considering that, for T Z 900 K, graphene has a positive thermal expansion coefficient, 52 our choice of a causes a slight in-plane compression 55 which, reducing the degree of vibrational anharmonicity, ensures structural stability of the graphene sheet at 4300 K.
The geometries of chemical precursors are optimized via 0 Kelvin DFT energy minimization. Then, random initial velocities (corresponding to a translational kinetic energy of B300 K) are assigned to each atom. Thus, microcanonical NVE dynamics is followed during 2 ps to ensure stability of the molecules. Coupling directly the internal degrees of freedom of the molecules with a thermostat would alter their dynamics, thus possibly causing artifacts such as sudden molecule splitting. The precursor internal atomic positions and velocities obtained via AIMD NVE sampling are used as input for subsequent precursor/graphene modeling.
The precursors are placed at random initial positions, approximately 13 Å above the graphene sheet, previously equilibrated in a separate AIMD simulation at a temperature T = 4300 K. Then, the molecule center of mass is gently pushed (during B1 fs) toward the surface using weak (o0.05 eV Å À1 ) bias forces. 56 This is necessary to induce precursor/graphene collision since the molecule center of mass has no drift after NVE sampling. The temperature of graphene is controlled via coupling with an isokinetic thermostat (velocity rescaling at each time step for all carbon atoms), which mimics canonical NVT phase-space sampling, while the precursors rotational and vibrational degrees of freedom are equilibrated at T = 4300 K upon frequent collisions with the graphene layer. The combination of NVE (for precursors) and NVT (for graphene) sampling within the same AIMD simulation was recently implemented in an own version of VASP. 56 AIMD snapshots are generated using the software Visual Molecular Dynamics. 57