Dynamics and kinetics of reversible homo-molecular dimerization of polycyclic aromatic hydrocarbons

Physical dimerization of polycyclic aromatic hydrocarbons (PAHs) has been investigated via molecular dynamics (MD) simulation with the ReaxFF reactive force field that is developed to bridge the gap between the quantum mechanism and classical MD. Dynamics and kinetics of homo-molecular PAH collision under different temperatures, impact parameters, and orientations are studied at an atomic level, which is of great value to understand and model the PAH dimerization. In the collision process, the enhancement factors of homo-molecular dimerizations are quantified and found to be larger at lower temperatures or with smaller PAH instead of size independent. Within the capture radius, the lifetime of the formed PAH dimer decreases as the impact parameter increases. Temperature and PAH characteristic dependent forward and reverse rate constants of homo-molecular PAH dimerization are derived from MD simulations, on the basis of which a reversible model is developed. This model can predict the tendency of PAH dimerization as validated by pyrene dimerization experiments [H. Sabbah et al., J. Phys. Chem. Lett. 1(19), 2962 (2010)]. Results from this study indicate that the physical dimerization cannot be an important source under the typical flame temperatures and PAH concentrations, which implies a more significant role played by the chemical route.


I. INTRODUCTION
3][4][5][6] Among these processes, it is generally believed that the particle nucleation starts from polycyclic aromatic hydrocarbon (PAH) dimerization as supported by experimental evidences, such as the bimodality in the soot particle size distributions (PSD) of premixed flames, stacked PAH conformations observed by transmission electron microscopy (TEM). 7,8As a consequence, PAH dimerization is regarded as the source term of the soot particle dynamics in the population balance models. 8,9n the irreversible PAH dimerization model, sticking coefficient (γ) is introduced as a tunable constant given by scaling the collision rate to the PSD in experimental observations. 10,11hough this sticking coefficient can be easily utilized in a) Author to whom correspondence should be addressed: K.Luo@ucl.ac.uk.
population balance models, it does not fully reflect the underlying physics in the dimerization process due to its independence on temperature.1][12][13][14] Recently, studies from Thomson et al. considered the reversibility in modeling soot nucleation and condensation, and the newly developed model was reported to give the best agreement with experimental data. 9,15In their models, the forward collision rate is calculated based on the collision theory, with a sticking coefficient in the range of 0.75-1.For the reverse process, some PAH dimers formed by van der Waals (vdW) interactions are assumed to instantaneously convert to chemically bonded dimers, while the rest dimers decompose to PAH monomers, which is lack of theoretical and experimental validations.
The van der Waals (vdW) interaction among particles can enhance their collision rate, which is corrected by an enhancement factor multiplied by the original collision rate based on the hard-sphere model. 16Based on the assumption of the spherical balls, previous studies reported that the enhancement factor for homo-molecular PAH collision is 2.2 being independent of PAH characteristics and temperature, [17][18][19] which is extensively adopted in soot modeling.As the enhancement factor is important to deduce the collision rate, 10 it is unclear whether the spherical ball assumption is appropriate for planar PAHs.Therefore, quantification of the enhancement factors for the planar PAHs at different temperatures requires more investigation.
To facilitate the study of soot modeling, more accurate understandings of the dynamics and kinetics of PAH dimerization are demanding.Studies of the PAH dimerization and nucleation have been aided by advanced modeling based on quantum mechanical (QM) methods 4 and molecular dynamics (MD) simulations. 20,21Schuetz and Frenklach utilized nonequilibrium QM/MD with semi-empirical quantum forces, named parametric method 3 (PM3), to study the energy transfer during binary collisions and the survival lifetime of PAH dimers at combustion temperatures. 21Limited by computational cost, the simulation time was no more than 20 ps.Thus, statistic quantities such as sticking coefficient, decomposition constant could not be derived.
In this study, detailed analyses on both the dynamics and kinetics of the reversible PAH dimerization at a range of temperatures with varied PAH characteristics have been conducted by MD simulations with the ReaxFF force field.The objective of this investigation is threefold: (i) to demonstrate the dynamic processes of the PAH collision and dimer decomposition under different temperatures and impact parameters, (ii) to develop a temperature and PAH characteristic dependent reversible model by assessing the survival lifetime of PAH dimers, and (iii) to test the model performance by comparing with previous irreversible model and experimental observations.

A. Simulation method
The ReaxFF force field utilized in MD simulation is developed to bridge the gap between the QM methods and classical MD, which allows MD to simulate both the physical and chemical processes. 22,23Energy calculation in the ReaxFF force field is composed of bonded and non-bonded interactions.The components of the bonded energy are bond, valence angle, torsion angle, over-coordination penalty, undercoordination stability, lone pair energy, etc.The van der Waals and Coulombic energy belong to the non-bonded part, and partial charge on every atom is updated every step following charge equilibration (QEq).All the parameters in the ReaxFF force field are trained against QM-based training sets.A more detailed description of ReaxFF energy terms can be found in the literature. 24The free energy profile of the coronene dimer calculated by the ReaxFF force field with C/H description is compared with that calculated by semi-empirical quantum force of PM3 and the density functional (U) M06-2X 25 with basis set def2SVP as presented in Fig. S1 in the supplementary material.The ReaxFF force field possesses the similar trend with M06-2X/def2SVP and has better accuracy than the PM3 method.Moreover, ReaxFF enables much longer MD simulations as it is about 700 times faster than PM3 and 150 000 times faster than M06-2X/def2SVP (see Table S1 in the supplementary material).Though the binding energy calculated by the ReaxFF force field has some deviation from the high accuracy DFT method, it is similar to the force field developed by Herdman and Miller, 26 which has been widely accepted and utilized to predict the thermodynamics of PAHs. 4,27In our previous study, this ReaxFF force field has been applied to study the soot nucleation mechanism and the boundary dividing the physical nucleation and no nucleation agrees well with the boiling/sublimation temperatures. 20etailed parameter sets for the ReaxFF C/H description can be found in the supplementary material.

B. Computational setup
According to the latest kinetic model of KAUST PAH mechanisms 2 (KM2), formation and growth of PAH have been extended to coronene. 28Besides the six-membered rings, five-membered rings are also commonly reported during the growth of PAH. 29,30Therefore, PAHs constituted of both sixmembered and five-membered rings were selected for this study, and their structures are displayed in Fig. 1.Before studying the bimolecular collisions, PAH monomers are firstly geometrically minimized under the ReaxFF force field via a conjugate gradient algorithm.Then, each PAH monomer is vibrational equilibrated to a desired temperature by 1 × 10 6 iterations with a time step of 0.25 fs through using the Nose-Hoover thermostat. 31he initial configuration for binary homo-molecular collision is shown on Fig. 2. PAH on the left is termed as PAH1 which is selected from the vibrational equilibrated state, and the other PAH, termed as PAH2, is replicated from PAH1.The size of the simulation box in the Z direction is 2 L d .PAH1 and PAH2 locate at 0.5 L d and 1.5 L d with an initial distance of L d .Normally, PAH collisions are assumed in the free molecular regime, which can be approximately estimated by where P and T are the partial pressure of the PAH monomer and temperature, respectively, R is the gas constant, and d p is the collision diameter expressed as where d A denotes the size of a single aromatic ring and equals to 1.395 coronene at 400 K when assuming the PAH mole fraction is of 10 3 at atmospheric pressure. 4According to Fig. S1 of the supplementary material, the binding energy profile of the coronene dimer diminishes to zero as the mass center distance increases to 8 Å.As such, an initial distance (L d ) of 30 Å is selected for all cases in this work.Perpendicular to the initial translational velocity, different impact parameters (L 0 ), i.e., offset distances of the two PAHs as shown in Fig. 2, continuously increase in the MD simulations to study the dynamics and determine the critical value at which collision no longer occurs.This critical value is considered to be the capture radius (r f ) for PAH dimerization.The black straight arrow and curve arrow indicate the initial relative translational and rotational velocities of PAH1, while these two components of PAH2 are set Based on the principle of equipartition of energy, translational (v) and rotational (w) velocities of PAHs can be described as follows: where K tran and K rot are the translational and rotational kinetic energies, k B is the Boltzmann constant, J and m denote the moment of inertia and mass of the PAH, and the angle brackets mean the ensemble average.Since PAH2 is set at rest, the relative velocity between PAH1 and PAH2 is √ 2 times of the original velocity of PAH1.The PAH dimerization process is significantly influenced by the relative collision velocity and collision orientation.Therefore, both translational and rotational stochastic velocities are generated satisfying the Gaussian distributions of each temperature for different PAHs at different impact parameters.Figure 3 shows the distributions of the translational and rotational velocities of pyrene at 400 K. Then these values are assigned to PAH1 in different collision events at a certain impact parameter.For each PAH at each temperature, a total of 10 000 simulation cases are performed, which is enough to get statistics significance.The simulation of the homo-molecular binary collisions is conducted adiabatically using microcanonical ensemble (NVE) under periodic boundary conditions with a time step of 0.25 fs and a total simulation time of 100 ps.By solving Newton's equations of motion, trajectories of PAHs are derived.MD simulations are performed using Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) package implemented with reax/c. 32Snapshots and movies in this study are prepared by Ovito. 33

III. RESULTS AND DISCUSSION
Based on the above simulation method, dynamics and kinetics of reversible PAH dimerization are discussed separately.Section III A mainly discusses dynamics of PAH collision and dimer decomposition at different impact parameters and temperatures.Then, from the perspective of kinetics, a reversible dimerization scheme is developed in Sec.III B. Finally, Sec.III C tests the model performance by comparing the reversible dimerization scheme with the previous irreversible model and experimental results.
FIG. 3. Distributions of the initial translational and rotational velocities of pyrene monomers at 400 K with 500 collision events at L 0 = 0.

A. Dynamics of PAH dimerization
Dynamics of PAH collision and dimer decomposition at a range of temperatures and impact parameters are studied and those of pyrene are displayed in Fig. 4. First, center-tocenter (L 0 = 0) collisions at 400 and 1200 K are scrutinized with snapshots of key stages recorded in Figs.4(a) and 4(b), respectively.Before collision, PAH2 is stationary, and PAH1 moves towards PAH2 at its initial translational direction and velocity as indicated by the straight mass-center trajectory highlighted in blue.For case (a), the formed PAH dimer continues to travel in its original translational direction accompanied by rotation as indicated by the mass-center trajectory after collision.The dimer formed in this condition is quite stable and survives the rest of the simulation.While for case (b), PAH1 has relatively higher translational velocity, which leads to a faster collision.After collision, the PAH dimer, however, is quite unstable and decomposes into PAH monomers after about 20 ps, with PAH1 and PAH2 traveling in different velocities.
In general, centrifugal collisions are more common compared with the center-to-center collisions.Two distinctive PAH collision events of pyrene with different impact parameters, that is, 1.5d p and 1.6d p , are examined, and key stages are recorded in Figs.4(c) and 4(d).Before collision (t < 11.25 ps), dynamics and trajectories of PAHs are similar.Upon approaching, PAH1 and PAH2 start to bend towards each other and subsequently coalesce into a dimer for case (c).Afterwards, the PAH dimer moves together in a straight line for a while.
After 60 ps, the PAH dimer takes apart into PAH monomers, which travel at different translational velocities.Nevertheless, for case (d), upon approaching, though PAH2 is attracted by PAH1 with a small rotation, they still do not collide.Then, PAH1 almost moves at its original velocity and direction while PAH2 keeps stationary.Thus, the capture radius for pyrene at 400 K is between 1.5d p and 1.6d p .
To quantitatively evaluate the exact value of critical capture radius and enhancement factor, collisions at different impact parameters are examined for PAHs at each temperature.The enhancement factor is evaluated by determining the critical capture radius of two approaching PAHs, where r o is the sum of the radius of two colliding monomers and r f is the critical capture radius.First, collisions at mean collision velocity (v p = 8RT πM ) of each temperature is performed to clarify the temperature effect.Here, we take pyrene, for example.From 400 to 1600 K, the enhancement factor decreases from 2.465 to 2.22 as shown in Fig. 5(a).This agrees with the theoretical study by Tikhomirov et al. 16 as the enhancement factor is proportional to 1/k B T. Then the enhancement factors for PAHs at 400 K are illustrated in Fig. 5(b).It is noteworthy that the enhancement factor decreases significantly with increasing of PAH size.For naphthalene, the enhancement factor for homo-molecular collision is about 3.6.As for coronene, it reduces to about 2.2.This trend  is similar to TiO 2 nanoparticles as the enhancement factor for smaller particles is larger than that of the larger ones. 34

B. Kinetics of PAH dimerization
PAH dimerization is a reversible process.However, most previous models only consider the forward collision, namely PAH + PAH → 2PAH, and neglect the dimer decomposition, i.e., 2PAH → PAH + PAH.Kinetics of PAH collision and dimer decomposition are related to temperatures, collision impact parameters, and PAH types.
To shed light on the temperature effect, 500 independent center-to-center collisions (L 0 = 0) at temperatures from 400 to 2000 K are investigated.The lifetime of a PAH dimer is recorded starting from the dimer formation and ending with the PAH dimer decomposition.Distributions of the number of pyrene dimers with different lifetimes at each temperature are demonstrated in Fig. 6.It is apparent that at low temperatures, especially at 400 K, the majority of the collisions lead to the formation of pyrene dimers with survival lifetime exceeding FIG. 6. Distributions of the number of surviving pyrene dimers formed from center-to-center (L 0 = 0) collisions at temperatures from 400 to 2000 K according to different lifetimes.80 ps.As the temperature increases to 600 K, the overall survival lifetime of the pyrene dimer decreases compared to that at 400 K. Specifically, the number of pyrene dimers with lifetime exceeding 80 ps decreases significantly.While that of short survival lifetime (less than 10 ps) increases.When the temperature rises to above 1200 K, none of the 500 dimers can survive over 80 ps.As the temperature increases to 1600 K, nearly all of the pyrene dimers have survival lifetimes shorter than 10 ps.Distributions of the number of PAH dimers at different lifetimes for naphthalene, acenaphthylene, phenanthrene, pyrene, cyclopean[cd]pyrene, and coronene from center-tocenter collisions are added to the supplementary material in Fig. S2.Overall, the temperature has a negative correlation with the survival lifetime of PAH dimers, and the number of PAH dimers with long survival lifetime increases with larger PAH size.
Apart from the temperature, impact parameter also significantly affects the lifetime.The number of surviving pyrene dimers formed from collisions at different impact parameters FIG. 7. Changes of the number of surviving pyrene dimers with respect to lifetime, temperature, and impact parameter.Dashed lines are the mean number of surviving pyrene dimers at different lifetimes.and temperatures is illustrated in Fig. 7.It is noteworthy that the number of surviving PAH dimers is non-uniform within the critical capture radius.At each survival time, the number of PAH dimers decreases at larger impact parameters, and it decreases to zero when the impact parameter is larger than the critical capture radius.As a result, the mean number of PAH dimers, ndimer (T , t), at each temperature and lifetime is proposed as where n dimer (T , t, L 0 ) is the number of surviving PAH dimers at temperature T, impact parameter L 0 , and lifetime t.The mean numbers of PAH dimers of different survival lifetimes at 400 K are added to Fig. 7 as shown in the dashed lines.Afterwards, mean numbers of surviving PAH dimers for each type of PAH at different temperatures and lifetimes are summarized in Fig. 8.The rate of the dimer decomposition can be expressed by where k reverse is the rate constant for the dimer decomposition.
Relation between ndimer and t obeys an exponential expression and follows where n 0 is the mean number of surviving PAH dimers at t = 0. Therefore, explicit values of n 0 and k reverse are quantified for a certain type of PAH at a certain temperature according to the fitting lines on Fig. 8.It is found that n 0 is not equal to the total number of initial collisions, which indicates the existence of a collision efficiency ( β) in the expression of the forward collision rate.Based on Figs. 5 and 8, values of the enhancement factor (ε), collision efficiency ( β), and decomposition constant (k reverse ) are obtained.The relation between the above parameters with respect to temperatures and PAH characteristics is summarized in Fig. 9. Specifically, the left panel shows the products of the collision efficiency and enhancement factor of PAHs at each temperature, which can be fitted by It is reported that the decomposition rates are sensitive to the ratio of the intermolecular energy to the thermal motion of the atoms. 15The intermolecular energy is proportional to the second power of the monomer size, 26 while the thermal energy

C. Model performance
To validate our dimerization scheme, we establish a simple kinetic model that contains the forward and reverse reactions between monomers and dimers The model developed in this study is compared to the irreversible model 10 and experimental results. 35In experiments, pyrene monomer concentrations in low-temperature supersonic flows have been detected by photoionization mass spectrometry.Dimerization rates in irreversible models only contain the forward collisions by utilizing a temperatureindependent sticking coefficient to evaluate the collision efficiency, which equals to 0.025 for pyrene.In our model, forward collision efficiency and decomposition rate at 60 and 235 K are directly derived from Eqs. ( 9) and (10).The reciprocal pyrene monomer concentrations (1/[n py ]) from our reversible model, previous irreversible model, 10  simulation cases, if there is one pyrene dimer taking decomposition within the 60 ps, the decomposition constant (k reverse ) is found to be 1.6667 × 10 6 ps 1 .Therefore, the decomposition constant for pyrene at 60 K is definitely smaller than 1.6667 × 10 6 ps 1 .Thus, we use the decomposition rate of 1 × 10 7 ps 1 , and the reciprocal pyrene monomer concentration is shown in solid brown line in Fig. 10, which matches well with the experimental results.The decomposition constant at 235 K is found to be about 0.000 408 7 ps 1 according to the MD simulation results.The reciprocal pyrene monomer concentration in green dashed line is found to be consistent with the result of the proposed model in this study as well as the experimental result.
Finally, we utilize the reversible PAH dimerization model to predict the PAH dimerization at a typical flame temperature of 1500 K. Initial concentrations are 3.5 × 10 14 , 3.5 × 10 19 , and 3.5 × 10 20 molecules cm 3 , which correspond to PAH partial pressure of ∼1 × 10 4 , 10, and 100 bars, respectively.Traces of reciprocal pyrene monomer concentration fractions (1/x py ) for the above three conditions are shown in Fig. 11.Apparent physical dimerization between pyrene monomers only happen when the partial pressure of the pyrene monomer is above about 10 bars.On the basis of these observations, we affirm that physical dimerization of PAH only occurs at low temperatures and cannot contribute to the dimerization source in soot formation under typical combustion conditions.A more stable or strong binding interaction between PAHs should exist during PAH dimerization at flame temperatures to suppress the decomposition of PAH dimers.[38]

IV. CONCLUSIONS
This study investigates both the dynamics and kinetics of the reversible homo-molecular PAH dimerization by molecular dynamics (MD) simulations with the ReaxFF force field.By performing both center-to-center and centrifugal collisions, the quantified enhancement factors for different types of PAHs are found to decrease with increasing PAH size and temperature.The survival lifetime of PAH dimers is highly dependent on temperature.At low temperatures, the majority of the PAH collisions lead to the formation of PAH dimers which can survive over 80 ps.When the temperature reaches above 1600 K, nearly all of the collision events have survival lifetime shorter than 10 ps.Besides the temperature effect, the number of surviving PAH dimers also decreases at larger impact parameters.To account for this effect, the number of surviving PAH dimers at different impact parameters is integrated into a mean number of surviving PAH dimers at a certain temperature.
Based on the relationship between the mean number and lifetime of surviving PAH dimers, we have developed a reversible PAH dimerization model composed of both forward and reverse processes.Expressions of the dimerization model have considered the effects from temperature and PAH characteristic.The reversible PAH dimerization model agrees well with experimental data in terms of the PAH dimerization trend at different temperatures, giving correct physics and better predictions compared to previous irreversible dimerization models.On the basis of these analyses, we find that the physical dimerization of PAHs can hardly contribute to soot nucleation at typical flame temperatures.More stable or strong binding interactions, such as Coulombic interactions between PAHs or PAHs connected by covalent bonds, should exist to suppress the rapid decomposition from PAH dimers to monomers, which needs further investigation.

SUPPLEMENTARY MATERIAL
See supplementary material for the free energy profiles of coronene dimers (Fig. S1), computational cost for single point energy (Table S1), distributions of the number of PAH dimers at different lifetimes (Fig. S2), the number of surviving pyrene dimers formed from different collision impact parameters and lifetimes at 60 and 235 K (Fig. S3), mean numbers and corresponding exponential fitting lines for the surviving pyrene dimers at 60 and 235 K (Fig. S4), and parameter sets for C/H descriptions of ReaxFF force field.Animations of the PAH collision processes at different temperatures and impact parameters are also attached.
FIG. 1. Structures of PAHs investigated in this study.

FIG. 4 .
FIG. 4. Trajectories and dynamics of the pyrene monomers at varied temperatures and impact parameters.(a) T = 400 K, L 0 = 0; (b) T = 1200 K, L 0 = 0; (c) T = 400 K, L 0 = 1.5d p ; (d) T = 400 K, L 0 = 1.6d p .Initial settings for cases (a), (c), and (d) are the same except for the impact parameter.The blue and red lines represent the mass center trajectories of PAH1 and PAH2, respectively.Snapshots of PAH1 and PAH2 are recorded with the evolution of time.Grey and white spheres represent carbon and hydrogen atoms.Movies of the above processes are attached in the supplementary material.

FIG. 8 .
FIG. 8. Mean numbers and corresponding exponential fitting lines for surviving PAH dimers at temperatures from 400 to 2000 K for different PAHs with respect to time.
FIG.10.1/[n py ] as a function of time at 60 and 235 K with initial concentrations of [n py ] 0 = 6.5 × 10 13 and 3.5 × 10 14 molecules cm 3 .The red solid line is the result from the reversible PAH dimerization model developed in the current study.The green solid line is the result directly extracted from the MD simulations as the decomposition constant of zero at 60 K.The brown solid line is the result as decomposition constant equals to 1 × 10 7 ps 1 at 60 K.The green dashed line is the result directly extracted from the MD simulations at 235 K.The blue solid and dashed lines are results based on the irreversible PAH dimerization model.10The black filled circles, triangles, and grey regions are the experimental data and error bars from Klippenstein et al.'s experiments.35