Dimerization of Polycyclic Aromatic Hydrocarbon Molecules and Radicals under Flame Conditions

This work presents a dynamic and kinetic study on the dimerization of polycyclic aromatic hydrocarbon (PAH) molecules and radicals under flame conditions using reactive force field (ReaxFF) molecular dynamics (MD) simulations. The accuracy of the ReaxFF force field is evaluated through comparing with quantum chemistry (QC) calculations of the barrier heights and species concentrations of PAHs reacting with H and OH radicals. A series of homobinary collisions between PAH molecules/radicals are performed to reveal the influence of temperature, molecular size, PAH composition, and the number of radical sites on the dynamics and kinetics of PAH dimerization. Instead of directly forming the strong covalent bonds, the majority of the binary collisions between PAH radicals are bound with weak intermolecular interactions. Effects of oxygen on PAH radical dimerization are also investigated, which indicates that the oxygenated PAH radicals are less likely to contribute to soot nucleation. In addition, the temperature, PAH characteristic, and radical site dependent collision efficiency for PAH radical−radical combinations is extracted from this study.


INTRODUCTION
Incomplete combustion of hydrocarbon fuels under fuel-rich conditions leads to the generation of the carbonaceous particulates known as soot. 1 In addition to its adverse effect on engine fuel efficiency, soot is regarded as one of the main environment contaminants, which poses a significant risk to human lives. 2 Studies of soot formation have been performed for many practical and laboratory-scale combustion systems, and it is commonly accepted that soot formation experiences the following stages: aromatic formation, soot nucleation, coagulation, and surface growth. 3 Among them, soot nucleation is a critical step yet still poorly understood. According to the previous experimental study, 4 the soot nucleation pathway is through physical PAH dimerization to stacks. However, both theoretical 5−7 and experimental 8 studies report that the intermolecular interaction between PAH molecules is too weak to stabilize the PAH dimers at flame temperatures. Recent studies on soot nucleation have shifted from physical dimerization to chemical dimerization as PAHs are connected by covalent bonds. 9, 10 Elvati et al. 11 proposed the role of radical−radical combinations in nascent soot formation and evaluated the sensitivity of soot nucleation to collision efficiency. And they 11 also reported the oxygen chemistry has a significant impact on the chemical growth of soot. Moreover, Johansson et al. 12 proposed the clustering of hydrocarbons by a radical-chain reaction (CHRCR) mechanism.
In the model of homodimerization of PAH molecules/ radicals, the reaction rates are estimated from the collision theory by taking into account the specific flame temperature and collision efficiency. 10,13,14 The collision efficiency for PAH radical−radical combinations is a parameter introduced to account for the reversibility, which is usually given arbitrarily and independent of temperature and PAH characteristic. 10,11 Therefore, an accurate collision efficiency for PAH molecules/ radicals is essential.
Quantum chemistry (QC) calculations of large PAH molecules/radicals are still challenging due to the prohibitively high computational cost despite impressive progress in computing hardware and software in recent decades. To simulate such large systems, some have tried to use a force field approach instead of the QC method. 6,15,16 ReaxFF is a bond order based force field and parametrized against data from QC calculations, which is capable of describing chemically reacting systems without a priori knowledge of predefinition of reactive sites. 17 ReaxFF-based reactive molecular dynamics (MD) has matured over the past decade and become a powerful yet affordable tool for the development of kinetic mechanisms for systems of large molecules and complex reactions. 18−20 The objectives of this study are threefold: first, to assess the reliability of the ReaxFF MD method in describing reactive systems composed of PAH and H, OH radicals by comparing against QC calculation based kinetic studies; second, to understand the dimerization mechanism of PAH molecules/ radicals and extract the temperature and PAH characteristic dependent collision efficiency for radical−radical combina-tions; and, finally, to evaluate the importance of the radical− radical combination to PAH growth and soot nucleation.

METHODOLOGY
2.1. Quantum Chemistry Calculation. According to the transition state theory, the barrier height of a reaction is the key factor to determine the reaction rate constant. In this study, geometries of reactants and transition states of hydrogen abstraction from PAH and PAH radicals by H and OH radicals are optimized by the density functional B3LYP method with the basis set 6-311G(d,p), 21 which has been widely used in previous studies involving chemical reactions of PAHs. 22−24 Transition states are confirmed to connect the reactants and the products by inspection of the normal modes of the corresponding imaginary frequencies. On the basis of the optimized geometries, single point energies are then calculated using the density functional theory (DFT) methods of M06-2X 25 and B3LYP with the 6-311G(d,p) basis set. All of the QC calculations are carried out using the Gaussian 09 software package. 26 2.2. Chemical Kinetic Modeling. A chemical reaction mechanism for reactions between PAHs and H, OH radicals is built including 19 species and 27 reactions in total. Details can be found in Table 1. Reaction rate parameters in the form of the Arrhenius equation are derived on the basis of QC calculations and the transition state theory. Coronene, one of the common PAHs, is selected for study. Specifically, R1 in Table 1 represents the reaction of an H atom on a coronene molecule (A7) abstracted by a H radical. R2 means the successive hydrogen abstraction takes place on R71 to form an A7 diradical named R72. Here, the different structures of R72 are not differentiated due to the fact that the per-site rate constants of H abstraction from A7 and A7 radicals are similar according to our recent study 24 at the level of M06-2X/6-311G(d,p). Moreover, simplified reacting systems of PAH and H and PAH and OH mixtures are then investigated in a 0-D homogeneous reactor model at constant volume and temperature. Chemical kinetic modelings are performed using the Chemkin-Pro software 27 to evaluate the kinetics of the formation of PAH radicals through hydrogen abstraction by H and OH radicals.
2.3. ReaxFF MD Simulations. The ReaxFF force field utilized in this study is the most up-to-date version involving interactions among C, H, and O, which are derived through training against DFT calculations based on the B3LYP/6-311G(d,p) method. 32−34 The van der Waals interaction for carbon is additionally described by the PBE exchangecorrelation functional with the DFT-D2 parameters. 35,36 This force field has been successfully applied to investigate thermal fragmentation of fullerene, 36 defect formation of graphene, 37 PAH dimerization and nucleation, 6,7 carbon nanotube nucleation from hydrocarbon precursors, 38 etc. Moreover, prediction of the potential energy surface (PES) of PAH dimers by the current ReaxFF force field is only 2−3 kcal/mol lower than that by the DFT method of M06-2X, while the predictions of the classical force field in the form of Buckingham 1 and Lennard-Jones formulations 39 are about 8  29 and the equilibrium constants calculated at the B3LYP/6-311G(d,p) level. A7 represents the coronene molecule, and R7X represents the coronene radical with X equal to the total number of radical sites.
The Journal of Physical Chemistry A Article kcal/mol lower. Details of the parameters of the force field are given in the Supporting Information in Table S1.
Reaction barriers of H abstraction from PAH molecules/ radicals by H and OH radicals are calculated using the ReaxFF force field on the basis of the optimized geometries obtained in section 2.1. Afterward, the kinetics of the chemical reactions between PAH and H/OH radicals are studied by ReaxFF MD simulations under the same initial conditions with the chemical kinetic modeling demonstrated in section 2.2.
Pyrene (A4) and coronene (A7) are simple but common precursors for soot inception. 2 Thus, investigation of the PAH dimerization is performed using the above two PAH molecules and their radicals, namely, R41, R42, R71, and R72 shown in Figure 1. Besides, the R71_O radical is also considered in our study, as the oxygen has a significant impact in the chemical growth of soot. 11 Before studying the binary collisions, structural equilibrium of PAH molecules and radicals is performed in the canonical ensemble (NVT) at temperatures of 1500, 2000, and 2500 K, respectively. Constant temperature is maintained by the Nose−Hoover thermostat with a damping constant of 100 fs. 40 Afterward, binary collisions between PAH molecules/radicals are performed adiabatically in the microcanonical ensemble (NVE) with a time step of 0.1 fs to ensure the energy conservation during the collisions, which is similar to our previous studies. 7,41 To achieve statistical significance, a total of 120,000 binary collisions are performed in the present study. Specifically, the initial relative mass center distance between the two colliding PAH molecules/radicals ranges from 50 to 80 Å, which is much larger than the effective intermolecular interaction distance between PAHs. 7 The relative translational and rotational stochastic velocities are sampled from a Gaussian distribution. The impact parameter of the colliding molecules/radicals is sampled from 0 to ε d p , where ε is the enhancement factor for PAH dimerization 7,42 and d p is the collision diameter of a colliding molecule. 13 The orientation of a molecule is rotated from 0 to 360°in three axes. In this study, MD simulations are performed using the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) package 43 implemented with reax/c. Snapshots and movies in this study are prepared by Visual Molecular Dynamics (VMD). 44

RESULTS AND DISCUSSION
On the basis of the above simulation methods, reliability of the ReaxFF MD method in describing reacting systems containing PAH molecules/radicals and H, OH radicals is evaluated by QC calculation based kinetic studies in section 3.1. Section 3.2 discusses the dynamics of homodimerization of PAH molecules/radicals and extracts the collision efficiency for radical−radical combinations through forming covalent bonds. Finally, section 3.3 assesses the contribution of the radical− radical combination to the PAH growth and soot nucleation.
3.1. Chemical Kinetics. Table 2 illustrates the structure of the A7 molecule and the barrier heights of hydrogen abstraction from the A7 molecule and A7 radicals by H and OH radicals. The barrier heights at different C−H sites are calculated by two DFT methods, namely, B3LYP/6-311G(d,p) and M06-2X/6-311G(d,p), and the ReaxFF force field. It is found that the barrier heights predicted by the ReaxFF force field situate between the results from B3LYP/6-311G(d,p) and M06-2X/6-311G(d,p), with deviations around 1 and −4.5 kcal mol −1 , respectively. Our previous study indicates that the M06-2X/6-311G(d,p) is found to obtain reaction barriers for H abstraction from PAHs closest to the high-level ab initial method CCSD(T)/CBS among 11 tested DFT methods. 24 For hydrogen abstraction from coronene molecule by OH radicals, the van der Waals complex (A7···OH) is typically formed first as discussed extensively by Shiroudi and Deleuze. 45 For simplicity, we assume the barrier height for hydrogen abstraction from A7 by OH radical is equal to the energy difference between its complex and the sum of A7 and   311G(d,p)) and the ReaxFF Force Field (Units: kcal mol −1 ) a a Radicals with one abstracted site and two abstracted sites are named as R71 and R72(xy), respectively [R72(xy) represents a coronene radical of two radical sites with a sequence number of x and y, as shown in the inserted snapshot].
The Journal of Physical Chemistry A Article OH radical, which is negative, as listed in Table 2. Generally, the barrier heights of H abstraction reactions by H and OH radicals at different sites predicted by the ReaxFF force field have only slight deviation from DFT calculations, which indicates the reliability of the ReaxFF force field to describe a system with PAHs and radicals. Furthermore, it is noteworthy that the computational cost for a single point energy by the ReaxFF force field is at least 100,000 times faster than the DFT calculations, which enables some complex systems to be studied.
Afterward, we compare the kinetics of reactions between PAH molecules and H/OH radicals. Figure 2 shows the mole fractions of species in a 0-D homogeneous reactor at the constant volume and temperature with initial molecule fractions of x H(OH) = 0.33 and x A7 = 0.66, respectively. The mole fractions of reactants and key products from the ReaxFF MD simulation are shown in dots, and those from kinetic modeling are shown in solid lines. Though there is some deviation especially at the beginning due to the discrete nature of MD simulations, results from MD simulations globally agree well with those from the chemical kinetic modeling. According to Figure 2, it is found that hydrogen abstraction from PAHs by OH radicals is faster than that by H radicals. Formation of species R71 is the most predominant, and then, species R72 follows with two H abstractions. It is noted that the mole fraction of (R71) 2 (R71 dimer connected by covalent bonds) is not negligible in the system composed of PAH and OH radicals according to the chemical kinetic modeling, which differs greatly from the ReaxFF MD simulation. Similar results happen for systems with initial molecule concentrations ranging from 0.125 to 0.913 mol/cm 3 . In fact, the collision efficiency of PAH radical−radical combinations is given arbitrarily in the range 0.1−0.85 in the chemical kinetic models, 10,11,31 and it is desirable to determine a more precise value for radical−radical combinations.
3.2. Dimerization of PAHs. To obtain the collision efficiency for PAH radical−radical combinations through forming covalent bonds, a series of binary collisions need to be studied. Initially, binary collisions between homo PAH molecules/radicals are studied at temperatures ranging from 1500 to 2500 K. In this study, the products from both the molecule−molecule and radical−radical dimerization are uniformly named as dimers. Firstly, representative dynamics including collision, dimerization, combination, and decomposition of PAH molecules, radicals at 2000 K are displayed in Figure 3. For the A7 molecule, the formed dimer is quite unstable, as it falls apart into molecules with a lifetime of about 10 ps according to Figure 3a. For PAH radicals, such as R71 radicals, they travel and rotate around each other after collision, as seen in Figure 3b. Subsequently, a covalent bond is formed as the two radical sites come closer, which represents the combination of the PAH radicals. However, the majority of the dimers composed of R71 radicals cannot take combinations through forming covalent bonds before the dimer falls apart, as illustrated in Figure 3c. Furthermore, as demonstrated in Figure 3d, the edge-to-edge collision is more conducive to producing a dimer connected by the covalent bond.
Through examining the behaviors of PAH molecules and PAH radicals after binary collisions, quantitative studies, such as the lifetime distribution of PAH dimer and the collision efficiency of the PAH radical combination, are scrutinized. The lifetime of a dimer is evaluated from the moment of dimer formation to that of the dimer decomposition. For PAH dimers formed through intermolecular interactions, the interlayer displacement ranges from 3.5 to 3.9 Å, 46 and our previous study indicated that the criterion of 4 Å is acceptable to distinguish PAH dimers from molecules, 6 which is also used in this study. Parts a−c of Figure 4 summarize the percentage of dimerization from identical PAH molecules/radicals versus the lifetime. First, behaviors of dimerization of A4, A7, and R71_O are quite similar, as the percentage of dimerization has an exponential correlation with the corresponding dimer lifetime shown in Figure 4a. Globally, the lifetime of PAH dimer decreases with the increasing temperature, and is sensitive to both the molecular size and composition. Because of the stronger intermolecular interactions between A7 molecules than that of the A4 molecules, 2,47 the A7 dimer has a longer lifetime than the A4 dimer at the same temperature. Though R71_O and A7 have a high similarity in size, the lifetime of the R71_O dimer is slightly shorter than the A7 dimer. This observation agrees with the result from Violi and co-workers, 11 as the oxygenated PAH prohibits the dimerization to some extent. According to Figure 4b and c, the relationship between the percentage of PAH radical dimerization and the lifetime of PAH dimer is markedly different from that of the PAH molecules, such as A4 and A7 in Figure 4a. That is, two remarkably distinctive relations are observed. Initially, the lifetime of PAH radical dimer has a negative correlation with temperature in the first 20−30 ps shown in Figure 4b and c, which is similar to the behaviors of A4, A7, and R71_O in Figure 4a. This indicates that those PAH dimers are bound with intermolecular interactions. Afterward, the percentage of PAH radical dimerization remains mostly constant versus dimer lifetime at 1500 and 2000 K due to the formation of covalent bonds through radical−radical combinations. Contrary to the negative temperature dependent relation for the physical dimerization, the percentage of PAH dimerization through forming covalent bonds at 2000 K exceeds that at 1500 K. Even though the R71 radicals have a higher possibility of forming a covalent bond than R41 radicals, The Journal of Physical Chemistry A Article it is still much lower compared to the suggested collision efficiency in previous studies. 10,11 Figure 4d illustrates the percentages of fragmentation of PAHs after binary collisions at different temperatures. It is found that, when the temperature rises to 2500 K, nearly all of the R72 and R42 radicals fragment, which leads to the percentage of dimerization decreasing instead of remaining constant, as shown in Figure  4b and c at 2500 K. This agrees with the previous experimental and theoretical studies of thermal fragmentation of PAHs at about 2200 K. 6,48 It is noteworthy that the thermal stability decreases with the increasing radical site and decreasing PAH radical size. That is, the percentage of the decomposition of R42 is higher than that of R41 and R72 at the same temperature.
Additionally, heterodimerization between PAH radicals is investigated. Percentages of the dimerization versus lifetime of PAH dimer from binary collision between an R71_O radical and a R71 radical at temperatures of 1500, 2000, and 2500 K are presented in Figure 5. It is found that the relation between the dimerization percentage and temperature is similar to that of the homodimerization of R71_O radicals, as illustrated in Figure 4a. Instead of forming covalence bonds, the dimers formed by the above two radicals are associated with intermolecular interactions. This once more confirms the  where β is the collision efficiency extracted from radical− radical combinations from Figure 4. For example, β for R41 is roughly equal to 0.001 and R71 is 0.003 at 2000 K. The collision efficiency (β = 0.1) for PAH radical combinations in the study of Violi and co-workers 11 is selected and adopted for comparison. The typical concentration of A4 is about 10 13 − 10 14 cm −3 , and that for A7 is 10 11 −10 12 cm −3 as detected in flames, 2 while the concentration of the PAH radicals is 1−2 orders of magnitude smaller than their PAH molecules. 10 On the basis of the relation between the PAH radical and dimer expressed above, the normalized concentrations (n/n 0 ) of the dimer from radical−radical combinations at varied initial concentrations are deduced and plotted as a function of time in Figure 6. Through applying the collision efficiency of 0.1, 11 the concentration of (R41) 2 to R41 equals about 0.3 at an initial R41 concentration of 6.5 × 10 12 cm −3 . However, it is about 1− 2 orders of magnitude smaller using the collision efficiency extracted from the current MD simulations under the same initial condition. The (R41) 2 to R41 is about 0.1 when the initial R41 concentration increases to 6.5 × 10 13 cm −3 , which is acceptable to explain the soot nucleus concentration detected in experiments. 10 The concentrations of (R71) 2 from R71 radicals at an initial number concentration of 6.5 × 10 10 and 6.5 × 10 11 cm −3 roughly equal 10 7 and 10 9 cm −3 . This is too small to account for the observed soot nuclei at the concentration of 10 11 −10 13 cm −3 . 22 Therefore, soot nucleation from radical−radical combinations is accessible to the situation where the concentration of PAH radical is high.

CONCLUSIONS
This study utilizes ReaxFF MD simulations to investigate the PAH molecule−molecule dimerization and the PAH radical− radical combination. Prediction of the energy barriers and species concentrations of PAHs reacting with H and OH radicals by the ReaxFF MD simulations is quantitatively consistent with those by quantum chemistry calculations based kinetic studies. To reveal the underlying relation among the lifetime of the PAH dimer, temperature, and PAH characteristic, a total of 120,000 binary collisions between two identical/ different PAH molecules/radicals are performed. Behaviors of homodimerization of PAH molecules, such as A4 and A7, and oxygenated PAH radicals, such as R71_O, are similar, as the formed dimers are attracted by intermolecular interactions. The lifetime of these PAH dimers is found to decrease with increasing temperature and decreasing PAH size. Moreover, the molecular composition affects the lifetime of the PAH dimer as the existence of oxygen decreases its lifetime. For PAH radicals, such as R41 and R71, the majority of the binary collisions lead to the formation of PAH dimers bound with intermolecular interactions, which is similar to the dimers formed from PAH molecules. And few of them are connected through covalent bonds from radical−radical combinations after collision. Below the thermal fragmentation temperature thresholds of the PAH radicals, collision efficiencies of radical− radical combinations of pyrene and coronene are extracted from ReaxFF MD simulations. The collision efficiency becomes larger with increasing PAH size, radical site, and temperature. This study indicates that collision efficiency prescribed in previous studies was generally overestimated and the contribution from PAH radical combinations to PAH growth and soot nucleation is highly dependent on the initial radical concentration.
■ ASSOCIATED CONTENT

* S Supporting Information
The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpca.8b07102. ReaxFF force field parameters for the C/H/O system (PDF) Figure 5. Percentage of the dimerization from heteromolecular binary collisions between a R71 radical and a R71_O radical with respect to dimer lifetime. Figure 6. Normalized concentrations (n is the concentration of PAH dimer, and n 0 is the initial concentration of PAH radical) of PAH dimer formed from R41 and R71 at different initial concentrations and collision efficiencies (from our current study and Violi and coworkers 11 ) as a function of time at 2000 K.