Molecular dynamics analysis of incipient soot morphology

Understanding mechanism of the soot formation process is important for reduction of harmful emissions from combustion and also for synthesis of carbon nanostructures. However, at the moment, both the mechanisms of carbon cluster formation and its morphology are poorly understood. In this paper, we present the results of molecular dynamics simulation of the formation and growth of carbon clusters formed at high temperatures from polyaromatic hydrocarbons.


Introduction
Due to enormous number of possible chemical ways for polycyclic aromatic hydrocarbon (PAH) formation and aggregation into nascent soot during hydrocarbon combustion, straightforward exhaustive analysis of the corresponding kinetic pathways is an unreachable goal. General mechanisms of these processes are formulated only on the phenomenological level and typically include such steps as PAH dimerization, growth of the resulting hydrocarbon complexes and their subsequent condensation into solid-state particles.
Recently, it was suggested that the initial stages of PAH dimerization and the formation of their covalently bound complexes are regulated by radical reactions [1,2]. From a computational point of view, these early stages of hydrocarbon growth can be studied using kinetic Monte Carlo approaches, which allow choosing the most energetically favorable pathways and products [3]. Another possibility is a direct atomistic modeling of combustion processes via classical molecular dynamics (MD) [4].
The condensation of gas-phase nanoparticles precursors into soot results from the interaction between physical and chemical aggregation: the first one implies clusterization of molecules via weak non-covalent interactions such as van der Waals forces, while the second includes formation of much more energetically stable covalent bonds. One can expect that at different temperature conditions these two mechanisms will have different contribution into the nucleation process. It can affect morphology and thermodynamic properties of the particles including their thermal stability [5]. Recently high resolution transmission electron microscopy (HRTEM) combined with image-processing algorithms gave important insights on the internal structure of primary soot particles [6,7]. Nevertheless scientific community is still far from developing a detailed 2 model of soot nanoparticles morphology and its evolution at particular temperatures. In this article on the basis of classical MD we study the atomistic structure of nascent soot formed from hydrocarbon precursors at T = 2500, 3000, 3500 K, analyze the possible products of its thermal decomposition and estimate the corresponding heat of sublimation.

Methods and model
Our MD calculations are performed using the LAMMPS program package [8] with the reactive interatomic potential ReaxFF. Computations are conducted in the periodic boundary conditions at density 0.078 g/cm 3 (the simulation cell volume is 100 × 100 × 100Å 3 ). N V T ensemble was chosen for simulations. Canonical ensemble is a typical choice for MD modeling of chemically active carbon systems at high temperatures [4,9,10], where drift of the temperature due to exothermicity of some reaction can significantly alter the initial conditions.
ReaxFF model [11] is an interatomic potential depending on the general bond order. Unlike non-reactive force fields, which are suitable for the modeling of hydrocarbon systems near the room temperature [12][13][14], reactive potentials take into account breaking and formation of covalent bonds and change of atom hybridization. It makes them one of the most widely used tools for atomistic simulations of combustion and related processes. Apart from the modeling of carbonaceous systems at extreme conditions, reactive potentials also give precise data on the elastic properties of carbon materials at moderate temperatures and pressures [15,16]. The ReaxFF potential divides the system potential energy into various contributions that depend on the bond order (except for the van der Waals and Coulomb terms): where E bond describes the bond energy between atoms, E val describes the interaction of atoms through valence angles, E tors is the energy of torsion (four-particle) interactions, E vdW and E Coulomb are the energies which represent the non-bonded interactions.

Incipient soot nucleation from PAH precursor
The modeling of cluster formation was divided into two stages. At the first stage, PAH molecules with mass in the range 600-10 4 Da (smallest molecules are product of coronene dimerization) were generated from coronene gas at 2000 and 2250 K during 5 ns simulations in N V T ensemble. On the second stage these molecules were placed into the main computational cell where simulations of clusterization process were performed at temperatures 2500, 3000 and 3500 K. Total amount of atoms in system was 5141, of which 3777 were carbon atoms and 1364 were hydrogen atoms. Total number of atoms remains constant during simulation. Simulations of nucleation process were performed in three cells at temperatures of 2500, 3000 and 3500 K.
Total simulation time is about 4 ns with a timestep of 0.1 fs. Figure 1 schematically shows transition from the initial particles [in figure 1(a), a small fragment of the initial system is presented] to the large particle [see figure 1(b)] at 3000 K (blue spheres indicate carbon atoms, white spheres indicate hydrogen). It was expected that this molecules will unite and form large cluster. Main mechanisms of PAH clusterization is physical (without breaking and formation of chemical bonds) and chemical nucleation [4]. We calculate number of carbon atoms in biggest cluster to analyze clusterization kinetics. It is assumed that atom belongs to the cluster if the distance to at least one of the particles in cluster is less than r cut = 1.85Å. Figure 2 which of the mechanisms prevails, we analyzed number of bonds in largest cluster between two atoms that have belonged to different molecules on initial step. This parameter can show how quickly molecules rebuild bonds and show total amount of new covalent bonds. One can see in figure 2(b) that carbon compounds at 2500 K almost do not rebuild bonds which corresponds to more flaked and amorphous structure of cluster compared to systems at 3000 and 3500 K.

Structural decomposition of soot nanoparticles using clustering algorithm
According to the experimental data [5,7,[17][18][19], carbon soot sublimation temperature varies in wide range from approximately 2000 to almost 4500 K. In the previous section, examples of soot were obtained using MD methods at different temperatures. Thus, the temperature of sublimation can be directly obtained through quasi-stationary heating of carbon cluster. However, such approach has high computational cost. In spite of using direct MD heating, it is better to transfer the logic of the experiment into calculations. In the experiment, the particles are heated with laser impulses until they fall apart into pieces with sizes on the order of several nm or less, the critical resolution of particle detection [5]. Provided that splitting into subclusters during the sublimation process corresponds to the minimum energy supplied, the energy could be estimated as energy of a minimum C-C covalent bonds cut in cluster. Therefore, energies of thermal decomposition were calculated based on the energy of a single C-C bond. This is a rough estimation; however, it allows us to understand the order of sublimation energy. Moreover, it allows us to distinguish structures by thermal stability factor among each other, which is an important aspect of our study. Hence, the major step of  To provide heterogeneity in atomic distribution, the uniform graphene list was perforated. Then agglomerative clustering algorithm was applied to allocated subclusters with the minimum number of connections between each other. the technique that we use is to split the biggest particle from MD calculations into subclusters with minimum number of connections to cut. For that need, we use agglomerative clustering algorithm [20]. The purpose was to split the particle into pieces, containing less than 100 carbon atoms (which approximately corresponds to the 1 nm particle size), providing minimum  cut connections. In figure 3 the demonstration of agglomerative clustering is shown. For that illustration we use randomly perforated graphene followed by clustering. Such approach allows us to estimate the number of connections to be torn during sublimation, taking into account that the energy of such process should be minimal.
From a physical point of view, the more connections in the cluster have to be cut, the more stable the cluster is. Following that logic, we calculated the number of connections to be torn to decompose the cluster per atom for three temperatures, and evolution of that parameter during the clusterization. The results are shown in figure 4.

Discussion
It should be noted that the temperature regimes discussed in this paper (T = 2500-3500 K) slightly overestimate the typical values for combustion processes. At lower temperatures small rates of clusterization and graphitization do not allow us to obtain sufficiently large and wellequilibrated soot nanoparticles on the scale of several nanoseconds, which is the upper limit for MD simulation. Nevertheless, it is known that at T < 4000 K carbon does not overcome any solid-liquid phase transitions (according to both experimental [21,22] and computational data [23][24][25][26][27]). Therefore in our calculations we expect to observe the same qualitative picture for the dependence between temperature and carbon nanoparticle morphology evolution as at 1500 < T < 2500 K. Recent MD calculations suggest that at T = 2500 [4], 3000 [9] and 3500 K [10] small fullerene-like particles do not demonstrate any signs of thermal decomposition, at least on the timescale of several nanoseconds. The same is also true for our system [see, e.g., figure 2(a)].
Visual analysis of soot nanoparticles formed at 2500 and 3000 K shows that they consist of chaotically interconnected graphene flakes with the sizes up to several nm (much larger than initial PAHs). This structure reminds of the incipient soot outer shell structure, for which experimental HRTEM analysis gives the size of the flakes (typically called "fringes") also on the scale of several nm [28,29]. If we assume that for such loosely interconnected structures sublimation process starts not from C 2 /C 3 ejection (as it was observed for giant fullerenes [30,31]) but from decomposition onto individual graphene flakes, we should expect that the value of the sublimation energy is significantly affected by the manner these flakes are covalently bonded to each other. Figure 4 gives basic estimations on the contribution of this covalent bond energy: we can see that decomposition of nanoparticle obtained at T = 3000 K requires almost twice as much bond breaking events as it is for T = 2500 K. With rough assumption for the single C-C bond energy E C-C = 350 kJ/mol it gives quite sufficient ∆E ≈ 30 kJ/mol (right axes in figure 4). This result could be interpreted in favor of the idea that higher temperatures on the stage of incipient soot formation results in more covalently interconnected and therefore more stable structures. This in turn should lead to higher temperatures required for their thermal decomposition in pulse-heating experiments [5].

Conclusion
Isothermal MD simulations were performed on ns timescale to investigate the influence of the formation temperature on the structure of incipient carbon soot nanoparticle. We observe that at higher temperatures covalent bond rebuilding during clusterization is more intense which results in more covalently interconnected and therefore more stable structures. Rough estimations suggest that contribution of this effect into sublimation energy could be on the order of several tens of kJ/mol, but additional MD calculations in wider range of temperatures are required for more precise statistics.