# Molecular dynamics simulations of reactive neutral chemistry in an argon-methane plasma

Glenn Christopher Otakandza Kandjani, Pascal Brault, Maxime Mikikian, Gautier Tetard, Armelle Michau, Khaled Hassouni

► To cite this version:

Glenn Christopher Otakandza Kandjani, Pascal Brault, Maxime Mikikian, Gautier Tetard, Armelle Michau, et al.. Molecular dynamics simulations of reactive neutral chemistry in an argon-methane plasma. Plasma Processes and Polymers, 2023, 20, pp.2200192. 10.1002/ppap.202200192. hal-03917195

HAL Id: hal-03917195

https://hal.science/hal-03917195

Submitted on 1 Jan 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
RESEARCH ARTICLE

PLASMA PROCESSES
AND POLYMERS

# Molecular dynamics simulations of reactive neutral chemistry in an argon-methane plasma

Glenn O. Kandjani<fnref n="1" />

Pascal Brault<fnref n="1" /> ID

Maxime Mikikian<fnref n="1" />

Gautier Tetard<fnref n="2" />

Armelle Michau<fnref n="2" /> ID

Khaled Hassouni<fnref n="2" />

Correspondence

Pascal Brault, GREMI, CNRS - Université d'Orléans, 45067 Orléans, France.

Email: pascal.brault@univ-orleans.fr

Funding information

Agence Nationale de la Recherche, Grant/Award Number: ANR-18-CE30-0016; Fédération CaSciModOT, Grant/Award Number: User account to HPC LETO; National Scholarship Agency of Gabon (ANBG), Grant/Award Number: PhD Fellowship grant

Abstract

Molecular dynamics simulations are performed to study the reactions in the volume of a low-pressure methane plasma diluted in argon. In the first step, a 1D fluid model is used to determine the initial molar fractions of initial species. The obtained composition thus becomes the input for the reactive molecular dynamics simulations. The study is carried out at 300, 400, 500, and 1000 K. Increasing the temperature increases the range of different molecules formed. The time evolution

of C₂H and CH₃ and the differ-
ent reaction pathways leading
to larger molecules show that
C₂H is the main precursor and
CH₃ is the main intermediate
precursor for the formation of
large molecules.

<figure><img src="image_2.png" /></figure>

KEYWORDS

carbon cluster, hydrocarbon, methane, molecular dynamics, plasma chemistry

1 | INTRODUCTION

Reactive low-temperature plasmas containing hydro-
carbon gases such as methane or acetylene are weakly
ionized gases containing atoms, molecules, ions, elec-
trons, and sometimes, solid particles of nanometer or
micrometer size.[1–4] These plasmas are of great scientific
and industrial interest because they are used in many
applications. Deposition of thin films of amorphous
hydrogenated carbon used for tribological materials,[5,6]
passivation layers,[7] cold-field emission cathodes for flat-
screen displays,[8,9] or for the deposition of other carbon-
based structures such as diamond-like carbon (DLC)
films,[7,10] carbon nanotubes, nanowalls, and so on[11,12]

have been addressed. Research in this field is mainly
devoted to the understanding of the mechanisms
involved in film growth and formation of particles in
the volume homogeneous phase, from the injection of
the hydrocarbon monomer to the formation of large
molecular clusters, to determine an optimal operating
mode. Essentially, gaseous hydrocarbon molecules
react with each other to create larger molecules, then
molecular clusters, being starting species for ultimate
soot or dust nucleation and growth.[13–15] The exact
mechanism of molecular species formation in the gas
phase and soot/dust nucleation remains a challenging
task owing to the complexity of the phenomena.[16] Since
these processes are both atomic and molecular in nature,

This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

© 2022 The Authors. *Plasma Processes and Polymers* published by Wiley-VCH GmbH.

<footnote n="1">¹GREMI, UMR7344, CNRS - Université d'Orléans, Orléans, France</footnote>

<footnote n="2">²LSPM, UPR3407, CNRS - Université Sorbonne Paris Nord, Villetaneuse, France</footnote>
reactive molecular dynamics (RMD) simulations are relevant to gain insight into some neutral fundamental phenomena of plasma chemistry. RMD is a simulation method for analyzing the physical and chemical changes of atoms and molecules interacting with each other through interatomic chemical potentials by solving Newton's equations of classical dynamics. With recent advances in powerful and efficient computer architectures as well as the availability of relevant reactive force fields,[17,18] RMD is becoming increasingly popular for studying the fundamental processes of reacting hydrocarbon molecules.

In this work, the objective is to use RMD simulations to study the volume reactivity of the main neutral species in a low-pressure Ar/CH₄ plasma. Present RMD simulations are targeted to account for the conditions of parallel plate capacitively coupled radiofrequency discharges. Temperature dependence is studied at 300, 400, 500, and 1000 K to identify the precursors initiating the formation of new species and large molecular clusters that are the first steps before nucleation.

Increasing the temperature is known to activate some chemical reactions. The present temperature range is considered typical for capacitively coupled radio-frequency discharges, usually operating at 300–400 K[19] Temperatures up to 1000 K are also considered as methane is a molecule of interest for deposition of diamond requiring larger gas temperature.[20]

This simulation technique requires the knowledge of all the initial species (obtained here thanks to a 1D fluid model) with their positions and velocities, and the corresponding appropriate interatomic potentials to obtain robust results.[17,18,21] Reactive interatomic potentials are necessary to study reactions in the gas phase, as these various processes require bond breaking/formation. Fortunately, several reactive potentials are available for addressing molecular and radical hydrocarbons.[22-28] The simplest and most robust reactive potential for hydrocarbons is the reactive empirical bond order (REBO) potential,[29] which has been extended to adaptive intermolecular reactive empirical bond order (AIREBO) to include long-range and torsional interactions.[27,30] There is another class of reactive interatomic potentials that includes partial charge dynamics in addition to bond order: Reactive variable charge force fields charge optimized many-body (COMB) and reactive force field (ReaxFF) being the most commonly used.[29,30] A recent study compared density-functional theory (DFT) with REBO, AIREBO, and ReaxFF potentials for the study of hydrocarbon reactions leading to hydrocarbon polymer growth on a silver substrate.[31] In this study, the REBO potential gave very similar results to DFT calculations. Especially because the ReaxFF parametrization-induced

energy barriers are not present in DFT and REBO.[30] However, an Aryanpour parameterization of the ReaxFF potential,[32] gave equivalent reactions and structures to the REBO potential after 6 ns. Since the ReaxFF potential is still very computationally expensive compared to the REBO potential without substantial gain, the REBO force field was used for all simulations. However, although the RMD offers a good level of reaction precision, the direct simulation involving the interaction of all particles present in the plasma, that is, ions, radicals, gaseous monomers, and electrons, is beyond the simulation capabilities, due to the fact that the REBO potential does not account any electronic processes.[31] So the present study is limited to volume reaction processes involving only neutral molecules of an Ar/CH₄ plasma. A 1D fluid model calculates the initial neutral composition of the plasma.[3,33]

The following section briefly describes the 1D fluid model giving the initial neutral species and their mole fractions. Then, the molecular dynamics simulation using the REBO potential is described, and the procedure for setting up the simulation box based on the fluid model is detailed. Section 3 presents the simulation results, such as the production rate of the newly formed species as a function of time and temperature. The time evolution of the concentration ratio of each of the initial neutral species is obtained to provide insight into their role in the formation of new species. Mass spectra and size distributions of the carbon clusters are calculated and analyzed. The large molecules formed, and the reaction pathways are then studied for different temperatures. Bond order ratios and structural characterization in terms of aliphatic, alicyclic, and aromatic compounds are then determined for each temperature.

## 2 | METHODOLOGY AND COMPUTATIONAL DETAILS

### 2.1 | 1D model description

The capacitively coupled Ar/CH₄ radiofrequency plasmas are generated in a 2.54 cm discharge gap between 20 cm diameter electrodes. These plasmas are, therefore almost uniform along the electrode radius and may be described using a one-dimensional model along the gap direction. The details of the model used in this work are given in reference,[3,33] and we will limit ourselves here to describe its main characteristics. The plasma is generated under a pressure above 10 Pa. It is, therefore, collisional enough to assume quasi-stationary electron momentum, neglect the inertial terms in the momentum balance, and use a drift-diffusion model for electron transport.[34] Even though the ion momentum is likely
nonstationary and the ions do not follow the RF field, we also made use of a drift-diffusion model with an effective electric field as described in[34] for these species. The electric field that governs the charged species' transport is governed by the plasma space charge and determined by Poisson's equation.

The electron average energy is determined by a governing equation that expresses the balance between the time-variation of the electron energy, the power received from the RF field, the rate of electron energy dissipated in collisional processes, and the energy transport fluxes by conduction and enthalpy diffusion. The rates of electron energy dissipation processes required to estimate the power dissipation are determined as a function of the electron average energy by solving the two-term expansion of the electron Boltzmann equation. In fact, this equation yields a two-Maxwellian electron energy distribution functions, a first one that describes the low energy electron population, that is, typically below 10 eV, which represents the major fraction of the electron population and a second one for the high energy population that governs the high energy threshold processes such as electronic excitation, dissociation, and ionization. The electron temperature inferred from the solution of the average electron energy corresponds to the temperature of the low-energy electron population.[35]

The EEDFs obtained from the solution of the Boltzmann equation were also used to determine the electron mobility as a function of the electron temperature. The electron diffusion coefficient was determined using the Einstein equation. For ions, we used Langevin mobility values, and the diffusion coefficients were determined from the Einstein equation assuming that ions are in thermal equilibrium with the neutral gas species, that is, $T_i = T_g = 300$ K.

The model described above results in a discharge module that solves for the electron and ion continuity equations, electron energy equation, Poisson's equation, and ion effective field equations. Besides this discharge module, the one-dimensional plasma model includes a reactive transport module that describes the coupled chemistry and transport for the neutral species that show characteristic times much longer than the RF period. Due to the low pressure, diffusion transport is relatively fast, and convection fluxes were neglected. The chemical model of Ar/CH₄ plasmas used in this work is a simplified version of the one published by De Bie et al.[36] Compared to this later, we only took into account up to two carbon-containing species, which results in a 22 species/56 reaction model. The diffusion coefficients were estimated from Lennard-Jones collision

integrals. The reaction rate constants were estimated using Arrhenius expression and gas temperature for reactions between heavy species. The rate constants for electron-heavy species collisions depend on the electron temperature and were determined from the solution of the Boltzmann equation as described previously.

The chemical model involves 12 neutral species and 10 charged species. Therefore, the discharge module includes 10 continuity equations, while the reactive transport module involves 12 continuity equations. The boundary conditions required for the continuity equations, the electron energy equation, and Poisson's equations are thoroughly discussed in Tetard et al.[3,33] Basically, we assume that the ions recombine totally on the electrodes; that the stable neutral species do not react on the electrode surfaces, while the radical species recombine on the electrode with probability values given by De Bie et al. and reference therein.[36]

Since the discharge characteristic time is much smaller than the neutral species characteristic times, the plasma composition may evolve over durations corresponding to a huge number of RF periods, that is, much larger than the number of RF periods required to reach a permanent discharge regime. Therefore, simulation of both discharge dynamics and plasma composition in terms of neutral species is performed using an iterative procedure that involves a two-step cycles. The first step consists of solving the coupled set of the charged species continuity equations, electron energy equations, Poisson's equations, and ion effective field equations while assuming a constant neutral species composition. This is performed by a time-integration of the discharge module equation over a few hundred RF periods until a permanent regime is reached. This first step yields time-averaged charged species density, electron temperature, and electric field, as well as electron-heavy species reaction frequencies. These are used in the reactive transport module that solves the stationary continuity equations of the long characteristic time-neutral species. The newly obtained neutral species composition is then used to solve the discharge dynamic equations for the next iteration. The procedure is iterated until a permanent discharge regime, and stationary neutral composition are achieved.

## 2.2 | Molecular dynamics simulations

Molecular dynamics simulations (MD) is a technique aimed at solving Newton's equations of motion (Equation 1) to numerically describe the displacement and interactions over time of all particles present in a
system (either solid, liquid, or gas), to interpret and predict possible changes in the system.[17,18,37]

$$ \frac{\partial^2}{\partial t^2} \vec{r}_i = \frac{1}{m_i} \vec{f}_i, \vec{f}_i = -\frac{\partial}{\partial r} V(\vec{r}_i(t), \vec{r}_i(t), \dots, \vec{r}_N(t)), \quad (1) $$

where $\vec{r}_i(t)$ is the position of atom $i$ with mass $m_i$ at time $t$, and $V$ is the interaction potential between all involved N species and $\vec{f}_i$ is the total force on the particle $i$.

The interactions between particles are described by appropriate interatomic potentials. The choice and quality of this potential are crucial to obtain accurate results. There are several potentials to describe the interactions between atoms of hydrocarbon molecules. In this work, the C–C, C–H, and H–H interactions are described by the REBO potential due to its agreement with density functional theory simulations (DFT).[31] For the Ar–C, Ar–H, and Ar–Ar interactions, we used the Lennard-Jones potential.[38]

$$ V_{LJ}(r_{ij}) = 4\epsilon_{ij} \left[ \left( \frac{\sigma_{ij}}{r_{ij}} \right)^{12} - \left( \frac{\sigma_{ij}}{r_{ij}} \right)^{6} \right], \quad (2) $$

where $r_{ij}$ is the interatomic distance between atoms $i$ and $j$. $\sigma_{ij}$ and $\epsilon_{ij}$ are the Lennard-Jones parameters corresponding to the potential range and well depth between atoms $i$ and $j$, respectively.

The $\sigma_{ii}$ and $\epsilon_{ii}$ parameters of the Ar and C atoms are presented in Graves and Brault[17] and those of the H atom in Delcorte and Garrison.[39] Knowledge of these parameters allows us to apply the Lorenz-Berthelot mixing rules below to determine approximately $\sigma_{ij}$ and $\epsilon_{ij}$[17] These values are summarized in Table 1.

$$ \sigma_{ij} = \frac{\sigma_{ii} + \sigma_{jj}}{2}, \quad (3) $$

$$ \epsilon_{ij} = \sqrt{\epsilon_{ii} \cdot \epsilon_{jj}}. \quad (4) $$

The REBO potential developed by Brenner is an extension of the Tersoff potential characterized by the addition of a many-body bond order functional to the two-body interatomic potential.[23] The following

<table><caption>TABLE 1 Lennard-Jones parameters of interactions between Ar, C, and H</caption><thead><tr><th>Atom pair</th><th>σ (Å)</th><th>ε (eV)</th></tr></thead><tbody><tr><td>Ar–Ar</td><td>3.4</td><td>0.01</td></tr><tr><td>Ar–H</td><td>2.42</td><td>0.0049</td></tr><tr><td>Ar–C</td><td>3.4</td><td>0.0014</td></tr></tbody></table>

analytical form, derived from Abell's pseudopotential theory, defines the total chemical binding energy $E$ between nearest neighbors.[40]

$$ E = \sum_i \sum_{j>i} [V_R(r_{ij}) - b_{ij} V_A(r_{ij})], \quad (5) $$

where $V_R$ and $V_A$ stand for the repulsive and attractive parts of the pair potential, respectively, and $b_{ij}$ is an empirical term of the bond order. The forms of the functions $V_R(r)$ and $V_A(r)$ are given by the following relations:

$$ V_R(r) = f_c \left(1 + \frac{Q}{r}\right) A e^{-\alpha r}, \quad (6) $$

$$ V_A(r) = f_c \sum_{n=1}^{3} B_n e^{-\beta_n r}, \quad (7) $$

where $f_c$ is the cutoff function that restricts the pair potential to the nearest neighbors. A, $B_n$, $\alpha$, and $\beta_n$ represent the parameters of the Morse function, and the term $Q$ is the parameter of the screened Coulomb function.[29] The empirical bond order term $b_{ij}$ between atoms $i$ and $j$ considers the local coordination, bond angles, radical character, conjugated bonds, and dihedral angle for C–C double bonds. It can therefore reproduce the reaction with the correct bond breaking/formation. It is composed of several subparameters whose details are available in the reference.[29] The open-source code large-scale atomic/molecular massively parallel simulator (LAMMPS), which implements the velocity Verlet algorithm for integrating Newton's equations, was used to perform the MD simulations.[37] The visualization of simulation results is carried out using the software visual molecular dynamics (VMD)[41] and OVITO,[42] postprocessing of the results make use of python scripts.[43]

## 2.3 | Simulation method

The first step of this study is to properly define the MD simulation box and the number of individual molecules consistently with the previously described 1D fluid model (Section 2.1). Therefore, to obtain a description of the plasma core consistent with the fluid simulation, we need a simulation box leading to the same number of collisions. This procedure will allow downsizing of the RMD simulation box size from the cm range (of the fluid model simulation box) to the nanometer range. This reduced size should not be less than the largest interaction cutoff length.[21]

The number of collisions $n_{coll}$ along a path from a simulation box/reactor side to the opposite one, is
proportional to the pressure-distance product by the following relation[21]:

$$n_{\text{coll}} = \frac{P d \sigma}{k_B T_g}, \qquad (8)$$

where $T_g$ is gas temperature, $d$ is the reactor or simulation box size, $P$ is the gas pressure and $\sigma$ is the collision cross-section.

Then, the connection between the simulation boxes of RMD and the fluid model is established by assuming that the number of collisions in both cases is equal,[21] which leads to

$$P_{\text{fluid}} d_{\text{fluid}} = P_{\text{MD}} d_{\text{MD}}. \qquad (9)$$

Thus, the number of molecules in the RMD simulation can be deduced from $P_{\text{MD}}$, assuming a $d_{\text{MD}}$ value, typically a few nm. This is chosen for allowing enough large molecule numbers necessary for statistical significance.

Preliminary simulations were performed at a temperature of 300 K using two simulation boxes, referred to as box 1 and box 2. The goal of these preliminary simulations is to choose the best simulation box in terms of computation time for our study. The total number of species and the size of the different simulation boxes were deduced from relation (9) and considering the molar fractions of each species (see Section 2.1). The ratios between the different hydrocarbon species are similar in both simulation boxes.

Box 1, with a size of $13.7 \times 13.7 \times 13.7 \text{nm}^3$ contained a total of 1519 molecules distributed approximately as follows: 66% Ar, 18% H$_2$, 8% CH$_4$, 3% C$_2$H$_4$, 2% C$_2$H$_2$, 2% C$_2$H, and 1% CH$_3$. Unreactive Ar atoms gain and carry out to the walls the excess energy of hydrocarbon bond formation: it means that Ar gas is operating as a thermostat maintaining the temperature of the neutral species. For this system, we chose the Langevin thermostat, which does not affect the time integration,[44] but only modifies the forces to perform the thermostating.[45] It is used with a separate time integration corrector as a microcanonical ensemble to update the velocities and positions of the atoms.[45] The Langevin thermostat was applied to argon at a temperature of 300 K, with a damping step of 10 fs which determines the rate at which the temperature is relaxed in the microcanonical ensemble. For other molecules, the simulations are carried out without any constraint, that is, in the NVE ensemble with $N$, $V$, and total energy $E$ kept constant. Box 2 has the same size as box 1 and does not consider argon atoms. It contains a total of 519 species distributed approximately as follows: 54% H$_2$, 24%CH$_4$, 9%C$_2$H$_4$, 5% C$_2$H, 4% C$_2$H$_2$, and 4% CH$_3$. For this box, the simulations

were performed in the NVT ensemble at the temperature of 300 K with a damping time of 100 fs. In these simulations, the average volume and temperature of the molecular system remain constant over time, but the average energy can vary. For our MD simulations at a constant temperature, the Nose-Hoover thermostat is the one implemented in the LAMMPS code for using the canonical NVT ensemble.[37,46–48] Periodic boundary conditions in the directions {$x, y, z$} are applied to both simulation boxes to mimic infinite volume. The initial velocities of the molecules are randomly chosen for each temperature according to the corresponding Maxwell-Boltzmann distribution, and the total simulation time is limited to 4 ns with a time step of 0.1 fs. The calculations are performed to determine the evolution of the mass spectrum during the simulation. At the end of the simulations, the species formed/destroyed are identified. Cluster growth was identified by assigning a cluster identifier to each atom so that atoms that have the same identifier at the end of the calculation belonging to the same cluster. If an atom has no neighbors within the prescribed cutoff distance, it remains a single atom.[37] The cutoff distance is generally chosen as the minimum between the first neighbor peak (r1) and the second neighbor peak (r2) of the radial distribution function.[49,50] In the present calculations, the cut-off radius of the carbon atom was set at 1.7 Å, since typical values for r1 and r2 for the C–C bond are at 1.45 and 2.55 Å, respectively.[49] Each atom has an initial identifier at the beginning and during all the simulations. Therefore, it is possible to track the reaction that initiated the formation of a molecule.

Five simulations were repeated in both cases to generate statistics. Figure 1 shows the concentration ratio (CR) (i.e., the ratio between the final number and the initial number of individual molecules) of the common molecules in the two boxes. This quantity provides information about the reactivity of each initial molecule during the simulation (i.e., if CR is strongly less than 100,

<figure><img src="image_3.png" /><figcaption>FIGURE 1 Concentration ratio of the initial neutral molecules in the two simulation boxes (blue bar box 1 and orange bar box 2).</figcaption></figure>
it means that the molecule was strongly involved in the reactivity of the system). It can be noted that the CR values of H₂ in both boxes are about 98% and that those of the molecules CH₄ and C₂H₄ are higher than or equal to 100% in both cases. The values of these molecules are also clustered around their respective average values for the five repeated simulations in both simulation boxes. These results indicate that these molecules behave very similarly in the two simulation boxes and that they react very weakly during the simulation. The CR values of the CH₃, C₂H, and C₂H₂ molecules show a rather scattered behavior compared to the other initial molecules. For the CH₃ molecule, a CR value around 77% in box 1 and 89% in box 2 is obtained, with a larger standard deviation in box 1, which means that the CR values in box 1 are more scattered with the five repeated simulations. For the C₂H₂ molecule, we have a CR value of around 114% in box 1 and 103% in box 2, which means that there is a slight production of C₂H₂ in both cases. In addition, the standard deviation of CR for this molecule is quite large and approximately the same in both cases. For the C₂H molecule, a very low CR value is obtained in both cases, around 7% in box 1 and 11% in box 2, with a slightly larger standard deviation in box 2. These results indicate that the reactivity in both boxes is conducted by the C₂H radical and that there is no statistically significant difference in the behavior of the CH₃, C₂H₂, and C₂H molecules in these boxes.

Figure 2 shows the superposition of the five repeated simulations of the simulated mass spectra of the molecules in box 1 (Figure 2a and box 2, Figure 2b) at 4 ns. These simulated spectra are obtained by accumulating the number of the different species in all the simulations, and it is sorted by the masses instead of the species names. In these spectra, only the hydrocarbon molecules are shown for clarity, since argon and molecular hydrogen are the main species in box 1, and molecular hydrogen is the main species in box 2. We

have previously observed from the value of CR that the majority of the initial molecules, in particular CH₃, CH₄, C₂H₂, and C₂H₄, have not changed significantly compared to their initial amount, so these are the main molecules in the two spectra.

Both simulations give the same types of formed molecules. Some discrepancies between the two boxes concern isolated molecules appearing in only one copy. These differences are not relevant and are considered to be due to the limited statistics. Among the formed molecules, the most significant (number greater than or equal to 5) are those with masses 40, 50, 51, 75, and 76, corresponding to C₃H₄, C₄H₂, C₄H₃, C₆H₃, and C₆H₄ molecules, respectively. These molecules are present in both boxes with roughly the same amount statistically. This allows us to conclude that the obtained mass spectra can be considered similar in both systems. Through this comparison between the two boxes, it can be concluded that the reactivity of the two systems is very similar. It shows that applying the Langevin thermostat to argon in an NVE ensemble yields essentially similar results than without argon and with a Noose-Hover thermostat in an NVT ensemble. This means that NVT simulations capture the essential part of the argon role, that is, dissipating bond energy formation and cluster cooling. Thus, performing simulations without argon has the advantage of reducing the simulation time, while keeping similar results for the formation of new species. Based on these results, argon is excluded from the simulation box, and simulations are run with the same parameters as in box 2 at temperatures of 300, 400, 500, and 1000 K. A total of 2683 molecules with the same composition ratios as in box 2 were introduced into a simulation box of size 10 × 10 × 10nm³. This box size and the number of molecules were determined by the relation (9) considering a pressure of 4% of the total gas pressure since 4% of methane is dissolved in 96% of Ar (fluid model). This procedure is valid since the mean

<figure><img src="image_4.png" /><figcaption>FIGURE 2 Superposition of mass spectra of hydrocarbon species obtained in (a) box 1 and (b) box 2 from five repeated simulations</figcaption></figure>
distance between molecules remains greater than the largest potential cutoff radius. The total simulation time for each temperature was set to 40 ns with a time step of 0.1 fs, which corresponds to $4 \times 10^8$ timesteps.

## 3 | RESULTS AND DISCUSSION

In this section, we first determine the production rate of our system at different temperatures. This rate represents the fraction of new molecules formed during the simulation, excluding the initial molecules. The following relation defines the production rate:

$$ \text{Production rate (\%)} = \frac{N}{N_t} \times 100, \quad (10) $$

where N is the number of new molecules, that is, excluding the molecules of the initial type (i.e., H₂, CH₄, C₂H₄, C₂H₂, C₂H, and CH₃), and Nt is the total number of starting molecules. Figure 3 shows the temporal production rate at temperatures of 300, 400, 500, and 1000 K. This rate is around 3% at the four temperatures. Between 0 and 5 ns, a rapid growth of the production rate is observed, with a slightly larger value at 1000 K. After 5 ns, the production rate starts to decrease slowly at all temperatures, indicating that several created molecules react for the formation of larger molecules. At 40 ns, the production rate at the respective temperatures of 300, 400, 500, and 1000 K is approximately equal to 2.84%, 2.74%, 2.81%, and 2.84%, indicating that the total amount of new molecules produced is approximately the same at all temperatures. These results indicate that the molecules that are strongly involved in the formation of new molecules are the same at all temperatures. These

<figure><img src="image_5.png" /><figcaption>FIGURE 3 Time evolution of the production rate (%) of new molecules formed during the simulation at temperatures of 300, 400, 500, and 1000 K</figcaption></figure>

precursors are present in small amounts, and some of them react very quickly in the first times of the simulation while others react all along the duration of the simulation.

### 3.1 | Evolution of the initial species

In this section, the temporal concentration ratio (TCR) of the initial species is determined, that is, the ratio between the concentration at time t of each species and the initial concentration. Figure 4a–f shows TCR of CH₄, C₂H₄, H₂, C₂H₂, CH₃, and C₂H, respectively. The TCR of CH₄ stays very large, with values over time of about 98%, 97.8%, 97.8%, and 96.8% at temperatures of 300, 400, 500, and 1000 K, respectively. It shows that the amount of CH₄ is not significantly different from the initial amount, implying that CH₄ has reacted very weakly with the other species over time. It can be noted that its consumption slightly increases with temperature. In Figure 4b, the TCR of the C₂H₄ molecule over time is about 99% at temperatures of 300, 400, and 500 K and about 96% at 1000 K. It shows that, as for CH₄, the C₂H₄ molecule is almost not consumed at temperatures from 300 to 500 K and that it reacts only weakly with the other molecules at 1000 K. In Figure 4c, TCR of the hydrogen molecule is around 98.8%, 98%, 97.5%, and 96% at 300, 400, 500, and 1000 K, respectively. These values are approximately constant over time, indicating that molecular hydrogen is only very weakly involved in the reactivity of the system and produces more atomic hydrogen in the time interval of [0, 1 ns] (the illustration of the evolution of the atomic hydrogen is not shown here) as the temperature increases. Figure 4d shows the concentration ratio of C₂H₂. At 300 K, this ratio first decreases and then becomes constant with time after reaching a minimum of 95%. In contrast, at temperatures of 400, 500, and 1000 K, we observe an initial increase in the concentration ratio around the respective maximum values of 108%, 115%, and 136%, followed by a decreasing trend before stabilizing around the respective values of 102%, 108%, and 112%. The formation of C₂H₂ during the simulation originates from the reaction:

$$ C_2H + H \rightarrow C_2H_2. \quad (11) $$

The reaction (11) is observed at all temperatures, and the amount of C₂H₂ created at the very beginning increases with increasing temperature. These results suggest that C₂H₂ serves as an intermediate molecule for the formation of other hydrocarbon molecules and that its contribution is more significant with increasing temperature. A similar observation of concentration
<figure><img src="image_6.png" /><figcaption>FIGURE 4 Temporal concentration ratio of initial neutral species CH₄ (a), C₂H₄ (b), H₂ (c), C₂H₂ (d), CH₃ (e), and C₂H (f), at temperatures of 300, 400, 500, and 1000 K</figcaption></figure>

increase and a decrease was also made by Majumdar et al.[51] in the case of a CH₄/Ar plasma in a high-voltage dielectric barrier discharge where they observe that the sum of the concentration of C₂Hₘ species (m = 2, 4, 6) over time initially increases and shows a decreasing trend after reaching a maximum concentration. The authors mentioned that it indicates that C₂Hₘ molecules (m = 2, 4, 6) serve as intermediate steps for the formation of larger hydrocarbons. This behavior has also been observed in low-pressure radiofrequency discharges in pure CH₄.[52] Although not under the same conditions, in the case of our study, we can observe (Figure 4d) that

C₂H₂ molecules play an important role as intermediates for the formation of new molecules at 1000 K. Figure 4e shows that the TCR of the CH₃ radical decreases with time and with the increase of temperature. The temporal evolution is consistent with the evolution of the production rate observed in Figure 3, implying that CH₃ is the initial molecule continuously involved in the formation of larger hydrocarbon molecules. Figure 4f shows that the TCR of the C₂H radical decreases with increasing temperature between times 0 and 1 ns meaning that the reactivity of C₂H becomes faster with increasing temperature. From 1 ns, the TCR of C₂H
becomes null at all temperatures, showing that the C₂H molecules were completely and very quickly consumed during the first times of the simulation, recombining with each other or reacting with the other initial molecules to form new molecules. C₂H is indeed a very well-known radical in the chemistry of acetylene-based hydrocarbon plasmas for its reactive character.[⁴,⁵³,⁵⁴] In a series of papers,[⁵⁵⁻⁵⁷] the temporal evolution of the different neutral species in a He/Ar/C₂H₂ plasma was studied, and the authors identified the main mechanism of polymerization of neutrals by reactions with the C₂H. In Benedikt,[⁴] a reaction table between the most common hydrocarbon radicals in C₂H₂ plasmas and stable molecules is shown at a temperature of 400 K. Large differences between the reactivities of the hydrocarbon radicals were observed, and radicals with a low number of hydrogen atoms, such as C₂H, appeared to be the more reactive and lead to extremely fast polymerization rates. This is consistent with our study, although it is not exactly the same plasma conditions. By comparing the TCR of each initial molecule (Figure 4a–f), it can be suggested that the reactions of the stable molecules CH₄, H₂, C₂H₄, and C₂H₂ with the CH₃ and C₂H radicals are probably very weak and that the reactions between the radicals are most favored at all temperatures. This suggests that the C₂H radical is the main precursor for the formation of new molecules at all temperatures.

## 3.2 | Newly formed molecular species

Mass spectrometry (MS) is a very versatile diagnostic method that plays a prominent role in the characterization of reactive plasmas such as hydrocarbon plasmas.[¹³,¹⁹,⁵⁸⁻⁶⁰] Nevertheless, before any comparison between simulated and experimental spectra, care should be taken to possible molecule fragmentation during mass spectrometry measurements. Present MD simulations do not account for such a phenomenon. Thus, corrected experimental mass spectra should be provided before comparison with MD simulations is planned.

At the end of each simulation, a mass spectrum is calculated to identify the different hydrocarbon molecules that were created during RMD. As observed in Section 3.1 from the evolution of the concentration ratio, the amounts of CH₄, C₂H₄, C₂H₂, and H₂ have only slightly changed compared to their initial amounts, while the amounts of C₂H and CH₃ changed drastically. Therefore, only the mass spectra of the newly formed molecules are presented in Figure 5a–d after 40 ns for the five repeated simulations at temperatures of 300, 400, 500, and 1000 K, respectively. For clarity, the spectra end

<figure><img src="image_7.png" /><figcaption>FIGURE 5 Superposition of the obtained mass spectra of hydrocarbon species for temperatures of (a) 300, (b) 400, (c) 500, and (d) 1000 K, excluding initial molecules</figcaption></figure>

at mass 180, since only isolated masses are present above this value.

We first observe a wider diversity of newly formed molecules at 1000 K compared to the other temperatures. Among the main identified masses, masses 30 and 44 correspond to the alkanes C₂H₆ and C₃H₈, respectively. The other main identified masses, namely 40, 42, 56, 65, 80, 82, 106, and 130, correspond respectively to the
<table><caption>TABLE 2 Summary of the formed molecules still present at 40 ns and corresponding reaction pathways for each temperature (first and last reactions, three dots replace intermediate reactions)</caption><thead><tr><th rowspan="2">Product</th><th rowspan="2">Initial molecules leading to the product</th><th colspan="4">Number of molecules formed versus temperature</th></tr><tr><th>300 K</th><th>400 K</th><th>500 K</th><th>1000 K</th></tr></thead><tbody><tr><td>C<sub>2</sub>H<sub>6</sub></td><td>CH<sub>3</sub><sup>*</sup> + CH<sub>3</sub><sup>*</sup> → C<sub>2</sub>H<sub>6</sub></td><td>95</td><td>100</td><td>106</td><td>89</td></tr><tr><td>C<sub>3</sub>H<sub>4</sub></td><td>C<sub>2</sub>H<sup>*</sup> + CH<sub>3</sub><sup>*</sup> → C<sub>3</sub>H<sub>4</sub><sup>*</sup></td><td>79</td><td>59</td><td>64</td><td>41</td></tr><tr><td>C<sub>3</sub>H<sub>6</sub></td><td>C<sub>2</sub>H<sub>2</sub> + CH<sub>3</sub><sup>*</sup> + H<sup>*</sup> → ... → C<sub>2</sub>H<sub>3</sub><sup>*</sup> + CH<sub>3</sub><sup>*</sup> → C<sub>3</sub>H<sub>6</sub></td><td>4</td><td>8</td><td>13</td><td>20</td></tr><tr><td>C<sub>3</sub>H<sub>8</sub></td><td>C<sub>2</sub>H<sub>4</sub> + CH<sub>4</sub> → C<sub>3</sub>H<sub>8</sub></td><td>6</td><td>9</td><td>6</td><td>14</td></tr><tr><td>C<sub>4</sub>H<sub>8</sub></td><td>C<sub>2</sub>H<sub>2</sub> + 2CH<sub>3</sub><sup>*</sup> → ... → C<sub>3</sub>H<sub>5</sub><sup>*</sup> + CH<sub>3</sub><sup>*</sup> → C<sub>4</sub>H<sub>8</sub></td><td>4</td><td>6</td><td>17</td><td>18</td></tr><tr><td>C<sub>5</sub>H<sub>5</sub></td><td>2C<sub>2</sub>H<sup>*</sup> + CH<sub>5</sub><sup>*</sup> → ... → C<sub>4</sub>H<sub>2</sub> + CH<sub>3</sub><sup>*</sup> → C<sub>5</sub>H<sub>5</sub><sup>*</sup></td><td>26</td><td>13</td><td>13</td><td>12</td></tr><tr><td>C<sub>6</sub>H<sub>8</sub></td><td>2C<sub>2</sub>H<sup>*</sup> + 2CH<sub>3</sub><sup>*</sup> → ... → C<sub>5</sub>H<sub>5</sub><sup>*</sup> + CH<sub>3</sub><sup>*</sup> → C<sub>6</sub>H<sub>8</sub></td><td>24</td><td>29</td><td>16</td><td>14</td></tr><tr><td>C<sub>6</sub>H<sub>10</sub></td><td>C<sub>2</sub>H<sup>*</sup> + 2CH<sub>3</sub><sup>*</sup> + C<sub>2</sub>H<sub>2</sub> + H<sup>*</sup> → ... → C<sub>5</sub>H<sub>7</sub><sup>*</sup> + CH<sub>3</sub><sup>*</sup> → C<sub>6</sub>H<sub>10</sub></td><td>—</td><td>2</td><td>—</td><td>11</td></tr><tr><td>C<sub>8</sub>H<sub>10</sub></td><td>2C<sub>2</sub>H<sup>*</sup> + 2CH<sub>3</sub><sup>*</sup> + C<sub>2</sub>H<sub>2</sub> → ... → C<sub>7</sub>H<sub>7</sub> + CH<sub>3</sub><sup>*</sup> → C<sub>8</sub>H<sub>10</sub></td><td>12</td><td>18</td><td>7</td><td>7</td></tr><tr><td>C<sub>10</sub>H<sub>10</sub></td><td>4C<sub>2</sub>H<sup>*</sup> + 2CH<sub>3</sub><sup>*</sup> → ... → C<sub>9</sub>H<sub>7</sub><sup>*</sup> + CH<sub>3</sub><sup>*</sup> → C<sub>10</sub>H<sub>10</sub></td><td>12</td><td>9</td><td>9</td><td>6</td></tr></tbody></table>

Note: Reactions with more than 2 interacting molecules summarize successive reactions occurring during a short period of time.

molecules C<sub>3</sub>H<sub>4</sub>, C<sub>3</sub>H<sub>6</sub>, C<sub>4</sub>H<sub>8</sub>, C<sub>5</sub>H<sub>5</sub>, C<sub>6</sub>H<sub>8</sub>, C<sub>6</sub>H<sub>10</sub>, C<sub>8</sub>H<sub>10</sub>, and C<sub>10</sub>H<sub>10</sub>, are mainly unsaturated hydrocarbons. The amount of each identified molecule can vary from one temperature to another.

A summary of the main newly formed molecules (for C<sub>n</sub> [n ≤ 10]) identified after 40 ns is given in Table 2 for each temperature. The reaction by which each molecule was formed is also given, as is its amount as a function of temperature. The three dots indicate the intermediate reactions before the last reaction that formed the molecule, and the horizontal lines mean that the molecule is not present at the indicated temperature. Only the C<sub>6</sub>H<sub>10</sub> molecule is not present at all temperatures, its formation appears to be favored at 1000 K. By limiting the calculation to the molecules presented in Table 2, the total number of molecules per temperature is 262, 253, 251, and 232 molecules at 300, 400, 500, and 1000 K, respectively. Thus, considering that the production rate at 40 ns is nearly the same at all temperatures (Figure 3), it means that increasing the temperature increases the diversity of formed molecules, as shown in Figure 5d. From this table, we can see that the main molecules involved in the reactions are the C<sub>2</sub>H and CH<sub>3</sub> radicals and that CH<sub>3</sub> contributes to the intermediate and final reactions. This is consistent with the temporal evolution of the concentration rate of these radicals shown in Figure 4e–f.

The reactions leading to C<sub>n</sub>H<sub>m</sub> with n ≤ 5 have been included in kinetic simulations.[61] Our results are consistent with the reactions obtained by these simulations. It should be pointed out that the MD simulations are sensitive to initial conditions, and thus, direct comparisons between MD simulations and experiments

(or kinetic simulations) need a very high accuracy on these conditions. Nevertheless, an agreement between kinetic and MD simulations can be considered as correct when the produced (small) n ≤ 5 molecules are the same. In this case, with such an initial reasonable agreement, MD simulations will predict the formation of larger molecules with confidence. In combustion studies, this is exactly what has been done, leading to hierarchical models of the molecular, cluster, and, finally, soot growth.[62,63] C<sub>2</sub>H<sub>6</sub> and C<sub>3</sub>H<sub>4</sub> hydrocarbons are the predominant new molecules formed at all temperatures. These two molecules were formed by the reactions identified in our RMD simulations listed below:

$$CH_3 + CH_3 \rightarrow C_2H_6, \qquad (12)$$

$$C_2H + CH_3 \rightarrow C_3H_4. \qquad (13)$$

Reactions (12) and (13) are gas-phase recombination reactions in hydrocarbon plasmas based on pure methane or argon-methane mixtures. They have been addressed in various experimental and simulations works, for different kinds of plasmas.[34,51,61,64,65]

The production of C<sub>2</sub>H<sub>6</sub> and C<sub>3</sub>H<sub>4</sub> molecules is lower at 1000 K compared to other temperatures, suggesting either that these molecules contribute as fast intermediates to the “temperature-accelerated” formation of other hydrocarbon molecules, or that the molecules that initiated their formation react less with each other at 1000 K. The temporal evolution of C<sub>2</sub>H<sub>6</sub> and C<sub>3</sub>H<sub>4</sub> molecules formed on average for the five repeated simulations are shown in Figures 6a,b, respectively. At temperatures from 300 to 500 K, a gradual increase of C<sub>2</sub>H<sub>6</sub> is observed, with typical
<figure><img src="image_8.png" /><figcaption>FIGURE 6 Temporal evolution of C₂H₆ (a), C₃H₄ (b), and of the C₂nH₃ (n = 1, 2, 3) family and C₄H₂ (c–f) at temperatures of 300, 400, 500, and 1000 K</figcaption></figure>

time scales consistent with the exponential decrease of CH₃ (Figure 4e). At 1000 K, the amount of C₂H₆ increases very quickly between 0 and 5 ns before decreasing very slowly with time, indicating that ethane served only very weakly as an intermediate for the formation of larger molecules at this temperature. Figure 6b shows that the amount of C₃H₄ increases very quickly at all temperatures during the first times of simulations, and then its production no longer occurs because C₂H molecules are no longer available (Figure 4f).

In Figure 2, where the simulations were run for only 4 ns, masses 50 and 51, corresponding to C₄H₂ and C₄H₃

molecules, respectively, were detected while they are no more significantly present after 40 ns as seen in Figure 5. These molecules are formed during our simulations by the following reactions:

$$C_2H + C_2H \rightarrow C_4H_2, \qquad (14)$$

$$C_2H_2 + C_2H \rightarrow C_4H_3. \qquad (15)$$

These molecules are very well known in hydrocarbon plasmas.[4,53,54,66] The C₄H₂ molecule and the C₂nH₃ (n = 1, 2, 3) molecule family are important intermediates
for the formation of larger hydrocarbon molecules. It explains why they are present at the early stages of the simulation and disappear very quickly due to their strong involvement in the growth of larger molecules. This behavior is also confirmed in our simulations, as shown in Figure 6c–f, where the time evolution of C₄H₂ and C₂nH₃ (n = 1, 2, 3) is plotted as a function of the temperature. C₄H₂ and C₄H₃ molecules are the most formed, and their amount decreases when temperature increases. The time scale of their initial formation coincides with the decrease of C₂H (Figure 4f). The amount of C₆H₃ is smaller than C₄H₂ and C₄H₃ but behaves in the same way versus temperature, while C₂H₃ slightly increases with temperature. These results indicate that C₄H₂ and C₄H₃ molecules strongly contribute to the growth of larger molecules. This contribution is reduced as temperature increases.

At 40 ns, the size distribution of the carbon clusters is calculated and summed up for each set of five repeated simulations, that is, the number of carbon atoms in each C_nH_m (n ≥ 3) molecule (C₂H₆ being the only identified carbon molecule with n = 2 at all temperatures). The corresponding cumulative distributions are shown in Figure 7a–d at temperatures of 300, 400, 500, and 1000 K,

respectively. First, even and odd clusters between C₃ and C₂₀ are observed at all temperatures, which is not very common in hydrocarbon plasmas, especially acetylene plasmas, where the even clusters are the most identified compared to the odd clusters.[4,53] C₃ clusters are the most formed at all temperatures. Beyond C₂₀, a few larger clusters appear and are limited at 1000 K. It could be due to the presence of only a very small amount of intermediate molecules, such as C₄H₂ and C₄H₃ at 1000 K (Figure 6f), where most of the C₄ peak is due to C₄H₈ molecules as shown in Table 2.

To get more insights into the reaction pathway leading to cluster growth, the detailed reactions leading to the largest clusters, indicated by arrows in Figure 7, are analyzed in the following section.

## 3.3 | Evolution of large clusters and their reaction pathway

The growth pathways of the largest molecules C₃₂H₂₄, C₅₀H₃₉, C₄₂H₃₇, and C₃₄H₃₂, identified in Figure 7 at each temperature are presented in Figure 8a–d. The black boxes represent the initial precursors that are at the

<figure><img src="image_9.png" /><figcaption>FIGURE 7 Size distribution of hydrocarbon clusters for C_nH_m (n ≥ 3) as a function of the numbers of carbon, at temperatures of 300 K (a), 400 K (b), 500 K (c), and 1000 K (d)</figcaption></figure>
**FIGURE 8** Reaction pathways of large molecules at 300 K (a), 400 K (b), 500 K (c), and 1000 K (d). The black, blue, and green boxes represent the initial, intermediate, and final molecules, respectively. The numbers on the left of the squares indicate the number of molecules. The red, yellow, and green arrows represent the initial, intermediate, and final reactions, respectively. The numbers on the arrow lines represent the number of molecules that contributed to the reaction.

<figure><img src="image_10.png" /></figure>

origin of the pathway leading to the final large molecules (green boxes). The blue boxes represent the intermediate molecules formed during the growth of the large molecule. The C₂H, C₂H₂, CH₃, and H species are the initial precursors present at all temperatures, with C₂H molecules being the more numerous. These molecules

were identified as the initial precursors of all large molecules formed at all temperatures with a crucial role in the C₂H radical. At 1000 K, other initial precursors for the growth of large molecules, such as C, C₂, and CH, were identified in small amounts. The red arrows represent the initial reactions (i.e., the reactions at
0.1 ns), and the yellow and green arrows represent the intermediate and final reactions, respectively. The numbers on the arrow indicate the number of molecules of the species that contributed to the reaction. As the temperature increases, there are more first reactions because the number of collisions between particles also increases. The first reactions are mainly initiated by the C₂H radical, which is consistent with the temporal evolution of its concentration rate (Figure 4f) and results in the literature.[4,53,54,66] A large fraction of the C₂H recombines to form C₄H₂ by the reaction (14). The C₄H₂, C₄H₃ (Figure 6), and C₆H₄ molecules are intermediate species formed in large amounts and have been shown to contribute to the formation of larger molecules. Indeed the C₆H₄ molecule was formed in the first times of simulation and then very quickly consumed. It is notably identified in the mass spectrum of Figure 2b for the simulations made at 4 ns. The C₆H₃ radical (Figure 6), formed by the interaction of C₄H₂ with C₂H, also contributes to the growth of large molecules from 300 to 500 K, while C₂H₃ is only formed at 1000 K. These results show that the C₂H radical is indeed the most important precursor involved not only in the formation of small molecules but also in the formation of large clusters. Moreover, C₄H₂ appears as a good intermediate for the growth of large molecules, and at least one molecule of the C₂nH₃ (n = 1, 2, 3) family was involved in the growth of large molecules. These results also suggest that the low production rate observed in Figure 3 is mainly due to the low amount of highly reactive molecules, such as the C₂H radical, compared to the stable molecules in the initial simulation box.

These results, summarized in Figure 8, are consistent with pathways reproduced in Heijkers et al.,[61] dealing with CH₄ conversion in a DBD reactor. The comparison can only be done by observing that the same molecules and radicals are produced up to C₅ and that generating reactions are similar.

It seems reasonable to consider the not so much addressed formation of large molecules with a C number greater than 5, bringing new information on neutral chemistry. Moreover, consistency between MD simulations and kinetic modeling up to C₅ species, suggests that MD simulations are relevant for exploring the formation for Cₙ, with n > 5.

## 3.4 | Structural properties of the formed molecules

The radial distribution function (RDF), g(r), is ubiquitous in the analysis of the internal structure of particle systems.[67] It provides important information on molecule structure, especially, the nature of the carbon bonds involved: single, double or triple bonds. The radial distribution function g(r) describes the probability of finding an atom j at a distance r from another atom i.[68] In the context of molecular dynamics simulations, g(r) is calculated using a finite number of trajectories that sample the thermodynamic ensemble of interest.

The RDFs at 40 ns of C–C bonds for C_n (n ≥ 3) are shown in Figure 9a at temperatures of 300, 400, 500, and 1000 K. The initial simulation box contains only double and triple bonds between carbon atoms, with a majority of C≡C triple bonds (C₂H, C₂H₂). Peaks 1, 2, and 3 correspond to intervals [1.16 Å; 1.3 Å], [1.3 Å; 1.5 Å], and [1.5 Å; 1.7 Å], respectively attributed to triple C≡C, double C=C, and single C–C bond lengths, respectively. The peaks are sharp at 300, 400, and 500 K, which means that the bonds of carbon atoms with their nearest neighbors are well-defined. Peak 2 is the highest one meaning that C=C bonds are the majority, but their occurrence decreases as the temperature increases. The height of the triple bond peak is very small at all temperatures and tends to flatten out at 1000 K, implying that these bonds

<figure><img src="image_11.png" /><figcaption>**FIGURE 9** Radial distribution function between pairs of carbon atoms C_n (n ≥ 3) (a) and Bond order ratio (b) at temperatures of 300, 400, 500, and 1000 K.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>FIGURE 10 Structural properties of the formed $C_nH_m$ ($n \ge 3$) molecules (a) and 5-, 6-, and 7-membered aromatic rings (b) at temperatures of 300, 400, 500, and 1000 K</figcaption></figure>

are little encountered. Peak 3 appears to broaden as temperature increases and tends to overlap with peak 2 at 1000 K, suggesting that increasing thermal motion favors identification as single bonds. The bond order is calculated as the area of each peak (Figure 9b). From these results, it appears that the majority of the newly formed molecules are mainly unsaturated. The appearance of single bonds, which were initially absent, corresponds to the growth of alternating –C=C–C= chains, which is strongly favored as temperature increases.

The post-processing utility MAFIA-MD[69] is used to characterize the spatial structures of the formed $C_n$ ($n \ge 3$) hydrocarbons. It identifies existing structures and generates statistics of alicyclic and aromatic hydrocarbons based on the number of carbon atoms at 40 ns elapsed time. It should be noted that the MAFIA-MD utility does not strictly verify the four classical conditions of aromaticity according to Huckel's rule,[70] because it is impossible to extract the exact electronic structure from the trajectory files, so the conditions of $(4n + 2)\pi$ electrons and perpendicularity of p-orbitals cannot be checked.[69] Moreover, in molecular dynamics, the atoms of an aromatic ring are often not in the same plane at a given time. Therefore, the utility focuses only on the bond distances and the closed nature of the ring structures.[69]

The idea is to extract chemical and structural information and to identify the presence of 5-, 6-, and 7-membered alicyclic and aromatic ring structures in an atomic cluster. The identification of the ring structures is important because the stability of these structures, especially the aromatic rings, is considered a key to the formation and growth of dust or soot particles, as identified during combustion processes.[16,62,63,71,72] Since the utility cannot distinguish between alicyclic and aromatic rings, we mention both alicyclic and aromatic compounds in our ring characterization. Figure 10a shows the chemical

characterizations as a function of temperature, with the green and blue bars representing the ratio of aliphatic compounds and the ratio of alicyclic and aromatic compounds, respectively. Aliphatic compounds are most abundant at all temperatures, but their percentage decreases as temperature increases leading to the presence of more alicyclic and aromatic compounds. Since double bonds are the most numerous (Figure 9), it suggests that aliphatic compounds are mainly unsaturated, especially alkenes and alkynes, and that high temperatures favor the formation of closed rings. Figure 10b shows the number of 5-, 6-, and 7-membered rings of alicyclic and aromatic compounds formed at each temperature. The number of 5-, 6-, and 7-membered rings increases with the temperature. Rings with 6 carbon atoms are the most numerous. A large increase in the number of 5 and 6-membered rings is observed between 500 and 1000 K, while this increase is lesser for the 7-membered rings. These results indicate that the formation of alicyclic and aromatic compounds is highly favored at high temperatures (above 500 K) and that 6-membered rings are formed most frequently.

4 | CONCLUSION

In this study, reactive molecular dynamics is used to identify the key precursors that initiate the formation of new hydrocarbon molecules and the growth of large clusters in the volume phase of an Ar/CH₄ plasma. The initial input data for the RMD simulations were obtained using a 1D fluid model of a capacitively coupled radiofrequency discharge operating at low pressure. The main neutral species H₂, CH₄, C₂H₄, C₂H, C₂H₂, and CH₃ were included in the RMD simulation, and the study was performed at temperatures of 300, 400, 500, and 1000 K. The results showed that the increase in
temperature increased the diversity of newly formed
molecules. The very low values of the production rate,
observed at all temperatures, are because some of the
initial major molecules, such as CH₄, C₂H₄, C₂H₂, and
H₂, did not react significantly. However, the formation of
small species and large hydrocarbon molecules is
observed at all temperatures. The molecules C₂H₆ and
C₃H₄, formed by the recombination of CH₃ and the
interaction of CH₃ with C₂H, respectively, are the main
molecules at all temperatures. Their amount is slightly
smaller at 1000 K because, at this temperature, the
molecules that initiate their formation are involved in
several other reactions. The newly formed species
CₙHₘ (n ≥ 3) are mainly unsaturated hydrocarbons, as
deduced from the predominance of carbon double bonds
at all temperatures. The structural characterization of
these species was carried out, and it was found that a
small amount of alicyclic and aromatic compounds is
present. This amount becomes more important as the
temperature increases, with a preferential formation of
rings with 6 carbon atoms. The C₂H radical was
identified as the main precursor for the reactions that
led to the formation of these different molecules at all
temperatures. The C₂H₂ molecule and the CH₃ radical
also contributed to the reactions for the formation of
larger molecules, and CH₃ was identified during the
simulation as the main intermediate for the formation of
small molecules and the growth of large molecules. The
C₄H₂ molecule and the C₂nH₃ (n = 1, 2, 3) family were
also identified in the early stages of the simulations and
contributed as intermediates to the formation of large
molecules.

ACKNOWLEDGMENTS

We are grateful to the Fédération Calcul Scientifique et
Modélisation Orléans Tours (CaSciModOT, LETO HPC
facility) for the numerical resources made available to us.
This work has been supported by the National Scholar-
ship Agency of Gabon (ANBG) and the French National
Research Agency (ANR) through the MONA project
(ANR-18-CE30-0016)

CONFLICT OF INTEREST

The authors declare no conflict of interest.

DATA AVAILABILITY STATEMENT

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

ORCID

Pascal Brault http://orcid.org/0000-0002-8380-480X

Armelle Michau http://orcid.org/0000-0002-7669-0277

REFERENCES

[1] M. Mikikian, L. Couëdel, M. Cavarroc, Y. Tessier, L. Boufendi, *Eur. Phys. J. Appl. Phys.* **2010**, *49*(1), 13106. https://doi.org/10.1051/epjap/2009191

[2] J. Winter, J. Berndt, S.-H. Hong, E. Kovačević, I. Stefanović, *Plasma Sources Sci. Technol.* **2009**, *18*(3), 034010. https://doi.org/10.1088/0963-0252/18/3/034010

[3] G. Tetard, A. Michau, S. Prasanna, J. Mougenot, P. Brault, K. Hassouni, *Plasma Process. Polym. 2022*, **19**, 2100204. https://doi.org/10.1002/ppap.202100204

[4] J. Benedikt, *J. Phys. D Appl. Phys.* **2010**, *43*(4), 043001. https://doi.org/10.1088/0022-3727/43/4/043001

[5] A. Grill, *Wear* **1993**, *168*(1-2), 143. https://doi.org/10.1016/0043-1648(93)90210-D

[6] B. K. Daniels, D. W. Brown, F. M. Kimock, *J. Mater. Res.* **1997**, *12*(9), 2485. https://doi.org/10.1557/JMR.1997.0328

[7] J. Robertson, *Mater. Sci. Eng. R Rep.* **2002**, *37*(4-6), 129. https://doi.org/10.1016/S0927-796X(02)00005-0

[8] A. N. Obraztsov, A. P. Volkov, K. S. Nagovitsyn, K. Nishimura, K. Morisawa, Y. Nakano, A. Hiraki, *J. Phys. D Appl. Phys.* **2002**, *35*(4), 357. https://doi.org/10.1088/0022-3727/35/4/311

[9] J. Robertson, *Thin Solid Films* **1997**, *296*(1-2), 61. https://doi.org/10.1016/S0040-6090(96)09381-9

[10] M. Frenklach H. Wang, *Phys. Rev. B Condens. Matter Mater. Phys.* **1991**, *43*(2), 1520. https://doi.org/10.1103/PhysRevB.43.1520

[11] M. Tanemura, K. Iwata, K. Takahashi, Y. Fujimoto, F. Okuyama, H. Sugie, V. Filip, *J. Appl. Phys.* **2001**, *90*(3), 1529. https://doi.org/10.1063/1.1382848

[12] K. Ostrikov, *Rev. Mod. Phys.* **2005**, *77*(2), 489. https://doi.org/10.1103/RevModPhys.77.489

[13] I. B. Denysenko, E. Wahl, S. Labidi, M. Mikikian, H. Kersten, T. Gibert, *Plasma Process. Polym.* **2019**, *16*(6), 1800209. https://doi.org/10.1002/ppap.201800209

[14] H. A. Michelsen, M. B. Colket, P. E. Bengtsson, A. D'Anna, P. Desgroux, B. S. Haynes, J. H. Miller, G. J. Nathan, H. Pitsch, H. Wang, *ACS Nano* **2020**, *14*(10), 12470. https://doi.org/10.1021/acsnano.0c06226

[15] A. Sharma, K. M. Mukut, S. P. Roy, E. Goudeli, *Carbon* **2021**, *180*, 215. https://doi.org/10.1016/j.carbon.2021.04.065

[16] H. Wang, *Proc. Combust. Inst.* **2011**, *33*(1), 41. https://doi.org/10.1016/j.proci.2010.09.009

[17] D. B. Graves P. Brault, *J. Phys. D Appl. Phys.* **2009**, *42*(19), 194011. https://doi.org/10.1088/0022-3727/42/19/194011

[18] E. C. Neyts P. Brault, *Plasma Process. Polym.* **2017**, *14*(1-2), 1600145. https://doi.org/10.1002/ppap.201600145

[19] C. Deschenaux, A. Affolter, D. Magni, C. Hollenstein, P. Fayet, *J. Phys. D Appl. Phys.* **1999**, *32*(15), 1876. https://doi.org/10.1088/0022-3727/32/15/316

[20] F. Silva, K. Hassouni, X. Bonnin, A. Gicquel, *J. Phys. Condens. Matter* **2009**, *21*(36), 364202. https://doi.org/10.1088/0953-8984/21/36/364202

[21] P. Brault, *Front. Phys.* **2018**, *6*, 59. https://doi.org/10.3389/fphy.2018.00059

[22] J. Tersoff, *Phys. Rev. B Condens. Matter Mater. Phys.* **1988**, *37*(12), 6991. https://doi.org/10.1103/PhysRevB.37.6991

[23] D. W. Brenner, *Phys. Rev. B Condens. Matter Mater. Phys.* **1990**, *42*(15), 9458. https://doi.org/10.1103/PhysRevB.42.9458
[24] A. C. T. van Duin, S. Dasgupta, F. Lorant, W. A. Goddard, *J. Phys. Chem. A* **2001**, *105*(41), 9396. https://doi.org/10.1021/jp004368u

[25] P. Erhart K. Albe, *Phys. Rev. B Condens. Matter Mater. Phys.* **2005**, *71*(3), 035211. https://doi.org/10.1103/PhysRevB.71.035211

[26] N. A. Marks, *Diamond Relat. Mater.* **2005**, *14*(8), 1223. https://doi.org/10.1016/j.diamond.2004.10.047

[27] T. C. O'Connor, J. Andzelm, M. O. Robbins, *J. Chem. Phys.* **2015**, *142*(2), 024903. https://doi.org/10.1063/1.4905549

[28] D. Zhang, M. R. Dutzer, T. Liang, A. F. Fonseca, Y. Wu, K. S. Walton, D. S. Sholl, A. H. Farmahini, S. K. Bhatia, S. B. Sinnott, *Carbon* **2017**, *111*, 741. https://doi.org/10.1016/j.carbon.2016.10.037

[29] D. W. Brenner, O. A. Shenderova, J. A. Harrison, S. J. Stuart, B. Ni, S. B. Sinnott, *J. Phys. Condens. Matter* **2002**, *14*(4), 783. https://doi.org/10.1088/0953-8984/14/4/312

[30] S. J. Stuart, A. B. Tutein, J. A. Harrison, *J. Chem. Phys.* **2000**, *112*(14), 6472. https://doi.org/10.1063/1.481208

[31] M. Zarshenas, K. Moshkunov, B. Czerwinski, T. Leyssens, A. Delcorte, *J. Phys. Chem. C* **2018**, *122*(27), 15252. https://doi.org/10.1021/acs.jpcc.8b01334

[32] M. Aryanpour, A. C. T. van Duin, J. D. Kubicki, *J. Phys. Chem. A* **2010**, *114*(21), 6298. https://doi.org/10.1021/jp101332k

[33] G. Tetard, A. Michau, S. Prasanna, J. Mougenot, P. Brault, K. Hassouni, *Plasma Sources Sci. Technol.* **2021**, *30*(10), 105015. https://doi.org/10.1088/1361-6595/ac2a17

[34] E. Gogolides, H. H. Sawin, *J. Appl. Phys.* **1992**, *72*(9), 3971. https://doi.org/10.1063/1.352250

[35] K. Hassouni, C. D. Scott, S. Farhat, A. Gicquel, M. Capitelli, *Surf. Coat. Technol.* **1997**, *97*(1-3), 391. https://doi.org/10.1016/S0257-8972(97)00201-6

[36] C. De Bie, B. Verheyde, T. Martens, J. van Dijk, S. Paulussen, A. Bogaerts, *Plasma Process. Polym.* **2011**, *8*(11), 1033. https://doi.org/10.1002/ppap.201100027

[37] A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in 't Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, R. Shan, M. J. Stevens, J. Tranchida, C. Trott, S. J. Plimpton, *Comput. Phys. Commun.* **2022**, *271*, 108171. https://doi.org/10.1016/j.cpc.2021.108171

[38] P. Brault, M. Ji, D. Sciacqua, F. Poncin-Epaillard, J. Berndt, E. Kovacevic, *Plasma Process. Polym.* **2021**, *19*, 2100103. https://doi.org/10.1002/ppap.202100103

[39] A. Delcorte, B. J. Garrison, *J. Phys. Chem. B* **2000**, *104*(29), 6785. https://doi.org/10.1021/jp001374h

[40] G. C. Abell, *Phys. Rev. B Condens. Matter Mater. Phys.* **1985**, *31*(10), 6184. https://doi.org/10.1103/PhysRevB.31.6184

[41] W. Humphrey, A. Dalke, K. Schulten, *J. Mol. Graphics* **1996**, *14*(1), 33. https://doi.org/10.1016/0263-7855(96)00018-5

[42] A. Stukowski, *Modell. Simul. Mater. Sci. Eng.* **2010**, *18*(1), 015012. https://doi.org/10.1088/0965-0393/18/1/015012:

[43] J. Sundnes, *Introduction to Scientific Programming with Python*, Springer International Publishing, **2020**. https://doi.org/10.1007/978-3-030-50356-7

[44] O. Farago, *Phys. A Statist. Mech. Appl.* **2019**, *534*, 122210. https://doi.org/10.1016/j.physa.2019.122210

[45] N. Grønbech-Jensen, *Mol. Phys.* **2020**, *118*(8), e1662506. https://doi.org/10.1080/00268976.2019.1662506

[46] D. J. Evans, B. L. Holian, *Nose-Hoover Thermostat* **1985**, *83*, 7.

[47] W. G. Hoover, *Phys. Rev. A.* **1985**, *31*, 1695.

[48] S. Nosé, *J. Chem. Phys.* **1984**, *81*(1), 511. https://doi.org/10.1063/1.447334

[49] S. Wang, K. Komvopoulos, *Sci. Rep.* **2020**, *10*(1), 8089. https://doi.org/10.1038/s41598-020-64625-w

[50] K. Kupka, A. A. Leino, W. Ren, H. Vázquez, E. H. Åhlgren, K. Nordlund, M. Tomut, C. Trautmann, P. Kluth, M. Toulemonde, F. Djurabekova, *Diamond Relat. Mater.* **2018**, *83*, 134. https://doi.org/10.1016/j.diamond.2018.01.015

[51] A. Majumdar, J. F. Behnke, R. Hippler, K. Matyash, R. Schneider, *J. Phys. Chem. A* **2005**, *109*(41), 9371. https://doi.org/10.1021/jp053588a

[52] I. Géraud-Grenier, M. Mikikian, F. Faubert, V. Massereau-Guilbaud, *J. Appl. Phys.* **2019**, *126*(6), 063301. https://doi.org/10.1063/1.5099326

[53] K. De Bleecker, A. Bogaerts, W. Goedheer, *Phys. Rev. E* **2006**, *73*(2), 026405. https://doi.org/10.1103/PhysRevE.73.026405

[54] J. R. Doyle, *J. Appl. Phys.* **1997**, *82*(10), 4763. https://doi.org/10.1063/1.3663333

[55] J. Benedikt, A. Consoli, M. Schulze, A. von Keudell, *J. Phys. Chem. A* **2007**, *111*(42), 10453. https://doi.org/10.1021/jp072892w

[56] A. Consoli, J. Benedikt, A. von Keudell, *J. Phys. Chem. A* **2008**, *112*(45), 11319. https://doi.org/10.1021/jp8042413

[57] A. Consoli, J. Benedikt, A. von Keudell, *Plasma Sources Sci. Technol.* **2009**, *18*(3), 034004. https://doi.org/10.1088/0963-0252/18/3/034004

[58] J. Benedikt, A. Hecimovic, D. Ellerweg, A. von Keudell, *J. Phys. D Appl. Phys.* **2012**, *45*(40), 403001. https://doi.org/10.1088/0022-3727/45/40/403001

[59] X. Gillon, L. Houssiau, *Plasma Sources Sci. Technol.* **2014**, *23*(4), 045010. https://doi.org/10.1088/0963-0252/23/4/045010

[60] A. Baby, C. M. O. Mahony, P. D. Maguire, *Plasma Sources Sci. Technol.* **2011**, *20*(1), 015003. https://doi.org/10.1088/0963-0252/20/1/015003

[61] S. Heijkers, M. Aghaei, A. Bogaerts, *J. Phys. Chem. C* **2020**, *124*(13), 7016. https://doi.org/10.1021/acs.jpcc.0c00082

[62] A. Violi, A. Kubota, T. N. Truong, W. J. Pitz, C. K. Westbrook, A. F. Sarofim, *Proc. Combust. Inst.* **2002**, *29*(2), 2343. https://doi.org/10.1016/S1540-7489(02)80285-1

[63] C. A. Schuetz, M. Frenklach, *Proc. Combust. Inst.* **2002**, *29*(2), 2307. https://doi.org/10.1016/S1540-7489(02)80281-4

[64] H. Kojima, H. Toyoda, H. Sugai, *Appl. Phys. Lett.* **1989**, *55*(13), 1292. https://doi.org/10.1063/1.1016363

[65] K. Hiraoka, K. Aoyama, K. Morise, *Can. J. Chem.* **1985**, *63*(11), 2899. https://doi.org/10.1139/v85-481

[66] V. De Vriendt, S. M. Miladinovic, J. L. Colaux, F. Maseri, C. L. Wilkins, S. Lucas, *Langmuir* **2011**, *27*(14), 8913. https://doi.org/10.1021/la2003035

[67] B. A. F. Kopera, M. Retsch, *Anal. Chem.* **2018**, *90*(23), 13909. https://doi.org/10.1021/acs.analchem.8b03157

[68] B. G. Levine, J. E. Stone, A. Kohlmeyer, *J. Comput. Phys.* **2011**, *230*(9), 3556. https://doi.org/10.1016/j.jcp.2011.01.048

[69] K. M. Mukut, S. Roy, E. Goudeli, *Comput. Phys. Commun.* **2022**, *276*, 108325. https://doi.org/10.1016/j.cpc.2022.108325

[70] Y. B. Apriliyanto, S. Battaglia, S. Evangelisti, N. Faginas-Lago, T. Leininger, A. Lombardi, *J. Phys. Chem. A* **2021**, *125*(45), 9819. https://doi.org/10.1021/acs.jpca.1c06402
[71] K. O. Johansson, M. P. Head-Gordon, P. E. Schrader, K. R. Wilson, H. A. Michelsen, *Science* **2018**, *361*(6406), 997. https://doi.org/10.1126/science.aat3417

[72] N. L. P. Andrews, J. Z. Fan, R. L. Forward, M. C. Chen, H.-P. Loock, *Phys. Chem. Chem. Phys.* **2017**, *19*(1), 73. https://doi.org/10.1039/C6CP06015C

**How to cite this article:** G. O. Kandjani, P. Brault, M. Mikikian, G. Tetard, A. Michau, K. Hassouni, *Plasma. Process. Polym. 2022*, e2200192. https://doi.org/10.1002/ppap.202200192