Reactive and electron force field molecular dynamics simulations of electric field assisted ethanol oxidation reactions

In this research, a combination of reactive force field (ReaxFF) and electron force field (eFF) molecular dynamics (MD) simulations is constructed to reveal the fundamental mechanisms for the influence of the electric field on ethanol oxidation reactions at atomic and subatomic scales. In total, 21 ReaxFF MD sim- ulations and 35 eFF MD simulations have been conducted. ReaxFF MD results indicate that the ethanol oxidation reaction is a two-stage process where the electric field plays varied roles in each stage. The first stage features the decomposition of ethanol molecules, in which the electric field influences the decomposition reaction rate by changing the kinetic energy of carbon-containing molecules/radicals on the order of 100–1000 kJ/mol and altering the molecular conformation and thereby the bond dissociation energy. At the second stage where oxygen molecules participate in the reaction, the electric field affects reactions by mod- ifying the reaction pathways. The application of the eFF MD simulations, for the first time, extends our understanding of the electric field effects on ethanol oxidation reaction to subatomic scales. The results in- dicate that the electric field modifies the electron energy on the order of 10–100 kJ/mol. The present study also offers interpretation of previous findings on electric field effects on reaction pathways and fluorescence experimental observations, and provides support for both “ionic wind” and chemistry-driven hypotheses. This research provides unprecedented insight into reactions aided by the electric field, which potentially can facilitate the design of realistic field-assisted combustion systems.

There are two plausible hypotheses interpreting the influence of an electric field on flame behavior. One is the ionic wind hypothesis where the alteration in flame behavior is attributed to the momentum exchange between ions and neutral atoms via collisions, causing macroscopic movement. The ionic wind acts as either a sole cause or an important contributor to the flame changes [4 , 10-12] . The "ionic wind" has also been visualized in laminar jet flames [13] . The other hypothesis states that flame behavior alteration by the electric field is chemistry-driven rather than momentumdriven, as the pressure exerted by the ionic wind (0.0004 atm [14] ) is insufficient to account for some critical changes in the flame front [15] . The electric field facilitates the movement of ions to (or from) the reaction zone where ions undergo dissociative recombination with electrons and generate combustion radicals (like OH [15] and H 3 O [16] ) that activate chain reactions [17 , 18] .
In our previous paper [19] , we have demonstrated that the electric field can modify reaction pathways of ethanol oxidation reactions and the alteration in reaction pathways accounts for the discrepancy in species between situations with and without an electric field. However, the mechanisms behind modification of reaction pathways by the electric field were not revealed in that study, due to insufficient information about subatomic events. A few questions remain unanswered to date. For example, will the electric field change subatomic properties as well as "observable" phenomena? Will the electric field modify the structure or conformation of the molecules? And will the electric field affect the properties of electrons, which are not calculated in previous reactive force field molecular dynamics (ReaxFF MD) research? To answer these questions, we have greatly extended our previous research by incorporating additional in silico experiments and employing new techniques such as the electron force field molecular dynamics (eFF MD).
In this research, we conduct a series of molecular dynamics (MD) studies to quantify the influence of the electric field with varying strength on ethanol oxidation reactions. The factors which might be responsible for the influence of the electric field are individually examined, including atomic velocity and kinetic energy distributions, molecular conformational changes, collision frequency and electron energy differences. The aim is to elucidate the detailed mechanisms for reactions assisted by the electric field from the atomic/molecular and even subatomic (i.e., electronic) perspectives.

Computational methods
In the present study, MD simulations with a reactive force field are conducted to investigate the atomic/molecular behavior of ethanol oxidation reactions, and MD simulations with an electron force field are introduced to elucidate the dynamics of electrons. Both methods are a step up compared with the classic MD that is incapable of simulating chemical reactions.

Reactive force field MD simulations
The potential of an MD system is calculated by the sum of bonded and non-bonded energies, as shown in Eq. (1) .
Specifically, in the ReaxFF MD methods, E bonded includes bond energy, overcoordination energy penalty, undercoordination stability, lone pair energy, energy of valence angle strain and energy of four-body torsional angle; E non-bonded usually comprises energy from Coulombic and van der Waals interactions between atoms [20] . For detailed explanation of a ReaxFF, the readers are referred to Refs. [21][22][23] . For ethanol oxidation reactions, the chemical elements in the system are C, H and O. Thus, a specific parameter set for CHO [20] is selected for ReaxFF MD simulations in the present study.
To study the effects of the electric field on ethanol oxidation, six electric strength values (i.e. 2 × 10 4 V/m, 5 × 10 4 V/m, 2 × 10 5 V/m, 5 × 10 5 V/m, 2 × 10 6 V/m, and 5 × 10 6 V/m) are applied in the x -direction of the system. These chosen values are based on previous experimental and computational studies [3 , 24 , 25] . For comparison, a base situation without any external electric field is included. For each situation, three replicates with different initial configurations are established. Therefore, 21 ReaxFF MD simulations have been conducted and analyzed.
The simulation configuration is a 100 × 100 × 100 Å 3 cubic box, where the reactions of fiv e ethanol (C 2 H 6 O) and 100 oxygen (O 2 ) molecules occur. A model initial configuration for the oxidation system is shown in Fig. 1 a. In each simulation, a canonical (NVT) ensemble under the Nosé-Hoover thermostat algorithm with a damping constant of 100 fs is performed. The periodic boundary conditions are applied to all three directions. Before heating the reactants, the system is energy minimized for 40 ps at 300 K to eliminate simulation artifacts that could cause simulation collapses [19] . After the energy minimization, the system is gradually heated from 300 K to 3000 K with a rate of 5.4 K/ps without any electric field. A temperature of 3000 K is maintained after 500 ps, and simultaneously an electric field with one of the aforementioned strengths is applied to the + x  direction. The time step for the ReaxFF MD simulations is 0.1 fs, and data are sampled every 0.1 ps.

Electron force field MD simulations
The electron force field (eFF) is an ab initio molecular dynamics model that includes calculations of electrons by solving a Schrödinger equation in a simplified way. In eFF, both the electrons and nuclei are mobile, and the interactions between electrons and nuclei are calculated via an effective potential that includes a repulsion between pairs of electrons due to the Pauli exclusion principle, a kinetic energy for individual electrons arising from the Heisenberg uncertainty principle and electrostatics [26] . In eFF, electrons are represented by electron gas which can grow and shrink in size over time. Detailed introduction about the development, application and validation of the eFF MD method can be found in Ref. [27] . In the present research, effective core potential models were used to simulate O and C atoms.
The aim of conducting an eFF MD simulation is to investigate the dynamics of electrons under varying electric field strength. Thus, structures of ethanol molecules at 500 ps of the ReaxFF MD simulations are extracted as input for eFF MD calculations. The structures of ethanol molecules are then studied under varying electric field strength. For each electric field, fiv e eFF MD simulations are constructed. Thus, 35 eFF MD simulations are conducted and post-processed. The temperature is maintained at 3000 K throughout the eFF MD simulations in a Nosé-Hoover NVT ensemble with a damping constant of 0.1 fs. The time step of the eFF simulation is 0.005 fs, and the simulated physical time is 5 ps.

Parallel computing and post-processing
All the simulations are carried out in the LAMMPS platform [28] . Specifically, a REAXC package is used for ReaxFF MD simulations, and an EFF package for eFF MD simulations. ARCHER, the UK national supercomputing service.
In each ReaxFF MD situation, results are averaged over three replicates. A 0.3 bond order cutoff [29] is selected for molecular recognition for every ReaxFF MD output. In eFF MD simulations, results are averaged over five replicates. All the post-processing of results is accomplished using in-house scripts. The VMD [30] package is used for the visualization of the molecular structures.

Validation of methods
The ReaxFF MD method for ethanol oxidation reactions has been validated in our previous publication [19] , and the validity of the eFF MD method for calculating the electron dynamics of C, H, and O atoms is demonstrated in Ref. [27] .

Results and discussion
This section starts with a phenomenological description of the influence of the electric field on ethanol oxidation reactions. To explore atomic scale events whereby the electric field affects the reactions, the velocity and kinetic energy distributions, collision frequency, and molecular conformations of ethanol molecules under various electric field strengths are scrutinized based on the ReaxFF MD results. To shed light on how an electric field influences subatomic-scale events, the electron energy from eFF MD simulations is analyzed.

Reaction rate and start time
The consumptions of ethanol and oxygen molecules in situations with and without an electric field have been recorded from t = 500 ps when an electric field is first applied. Figure 1 b and c illustrates the time-evolutions of ethanol and oxygen molecules, respectively. As shown in Fig. 1 b and 1 c, the ethanol oxidation reactions comprise two stages: at the initial stage (stage I, between 500 and 1300 ps), the electric field accelerates the reaction of ethanol molecules, as faster consumptions of ethanol and oxygen molecules are observed in the cases with an electric field; at the subsequent stage (stage II, 1300 ps and afterwards), the consumption rate of ethanol molecules in the non-electric-field case exceeds its electric field counterparts, and the reaction of oxygen molecules in the non-electricfield case occurs. The reaction rate coefficients for ethanol molecules at both stages are provided in the Supplementary Information. Figure 1 b and c suggests that the principal reaction at stage I is the decomposition of ethanol molecules, while oxygen molecules mainly participate in the reactions at stage II. This reaction sequence is consistent with our previous findings about reaction pathways [19] , indicating the ethanol oxidation reaction starts with the decomposition of ethanol molecules, followed by the oxidation of radicals from ethanol decomposition. Reaction pathways in pyrolysis include C 2 H 6 O ↔ C 2 H 5 + OH, C 2 H 6 O ↔ C 2 H 4 + H and C 2 H 5 ↔ C 2 H 4 + H, with H + O 2 ↔ HO + O being a typical oxidation pathway.
The starting instants when the decomposition of the first ethanol molecules occurs in each case are further examined. As compared in Fig. 1 d, an electric field slightly postpones the starting time of the ethanol decomposition.

Atomic velocity and kinetic energy distributions
The initial stage (stage I) witnesses the acceleration of ethanol decomposition by the electric field. To check whether the alterations in reaction rates is attributed to the atomic velocity changes when the electric field is applied, the x -component velocity (in the same direction as the electric field) and the kinetic energy (KE) distributions of molecules/radicals containing carbon atoms are scrutinized. For simplicity, translational movements of these molecules/radicals are considered. The velocities are represented by x -component velocities of carbon atoms ( v x,C ), and the KE distributions are simplified by the translational KE dis- , m C is the atomic mass and v y,C and v z,C are the y -and z -component velocities, respectively). In Fig. 2 a and b, the velocities/KEs of carbon atoms are averaged over three replicates in each time frame, and a distribution of velocities/KEs is obtained based on the array formed by the average velocities/KEs from every single time frame. The mean values of velocity/KE distributions of the electric field cases are individually compared with that of E = 0 case, and T -test results suggest the differences in mean values between cases with and without electric field are statistically significant with p -value < 0.05. As shown in Fig. 2 a, the electric field does alter the x -component velocities of C atoms in terms of the mean value of velocity distribution (bar in the middle of a violin plot) and velocity variations, and the alteration is non-linear. For example, an electric field of 5 × 10 6 V/m accelerates C atoms in the x direction, whereas a strength of 2 × 10 6 V/m decelerates the x -direction motion. However, when the velocities of the other two components ( y and z ) are included, the electric field, indeed, energizes the molecules/radicals containing C atoms, as Fig. 2 b shows that higher mean values of KEs are observed in the cases with an electric field. The imposition of an electric field induces flow motion via collisions between gas molecules and electric field-driven radicals, which is termed the "ionic wind". The ionic wind increases mass transport and mixing [31] . In the present study, the higher KEs suggest that stronger ionic winds from Please cite this article as: X.Z. Jiang   C-containing radicals could exist in the electricfield cases, and thus higher mass consumption of ethanol molecules are expected at stage I as shown in Fig. 1 b. The KEs are also compared between the starting (500-550 ps) and concluding (1250-1300 ps) periods of the decomposition-dominated stage (stage  Fig. 2 c, in the cases with the electric field, the KEs are attenuated as the decomposition continues, whilst the KE in the nonelectric-field case slightly increases. A similar analysis has also been performed for oxygen atoms and the details can be found in the Supplementary Information.

Collision frequency
The slightly increased KE at the concluding period of the decomposition-dominated stage (stage I) in the E = 0 case can contribute to the faster ethanol consumption at stage II. In the meantime, the higher quantities of C 2 H 6 O and O 2 molecules in the E = 0 case in the transient period from stage I to stage II enable a higher collision frequency, thereby facilitating the faster consumption of ethanol at stage II. The collision frequencies in the transient period from stage I to stage II (1250 to 1350 ps) are compared among the cases as shown in Fig. 3 . (The calculation of collision frequency is shown in the Supplementary Information.) The E = 0 case outnumbers the electric-field cases in the collision frequency, which is mainly due to the slightly higher number of ethanol and oxygen molecules at the end of Stage I. According to the simple collision theory [32] , the reaction rate is in proportion to the product of the collision frequency and a probability factor, q , which is a function of temperature ( q = exp(-E c / RT ), R is the ideal gas constant, and E c = E a − 0.5 RT with E a being the activation energy of the reaction). In the present study, the temperatures of all the cases are set as a constant of 3000 K, and all the cases share the same q . Thus, more effective collisions occur in the E = 0 case, accounting for the higher consumption rate of ethanol at stage II.

Molecular conformation and bond energy
The dynamics of molecules/radicals introduced in the previous two sections explain the discrepancies in reaction rates among cases with and without  an electric field at various stages. In this section, the effect of the electric field on the molecular conformation will be investigated.
The conformations of ethanol molecules are further scrutinized under different electric field situations. The C -C bond length ( d C -C in Fig. 4 a), Fig. 4 Fig. 4 c) are calculated before the earliest ethanol decomposition reaction occurs (at t = 505 ps). The electric field attenuates d C -O , whereas the changes in d C -C and α C -C -O by the electric field appear to be irregular. The distinct effects can be attributed to the polarity of the ethanol molecule: the ending -OH group forms the main polarity of the ethanol molecule, resulting in sensitivity of C -O bond to the imposition of the electric field; compared with the C -O bond, the remaining part of the ethanol molecule is relatively insensitive to the electric field due to the weak polarity.

b), and C -O bond length ( d C -O in
The attenuation of the C -O bond length further leads to the variations in C -O bond energy which can be calculated by the equation provided in Ref. [22] . As Fig. 4

Energy of electrons
Sections 3.2 to 3.4 have revealed the atomicscale behavior of ethanol molecules when an electric field is imposed. How the electric field could influence subatomic-scale dynamics of ethanol molecules is the focus of this section. The average total energy of electrons of individual ethanol molecules is calculated by eFF MD simulations. As Fig. 5 suggests, the electric field generally energizes ethanol molecules by charging the electrons, except for the E = 2 × 10 4 V/m case (will be mentioned in the next paragraph). The ethanol oxidation reaction is a redox reaction, where electron transfer occurs and the transfer rate is modulated by energy levels of dipolar states [33] . The modified electron  energy by the electric field revealed in this study can result in changes in the energy level of dipolar states and consequently modify the electron transfer rate. It is noteworthy that the changes in electron energy by the electric field is of an order of 10-100 kJ/mol, while the changes in KE is of an order of 100-1000 kJ/mol ( Fig. 2 ). Thus, the principal contribution to the elevated reaction rate at stage I in the presence of an electric field stems from the increased KE. Compared with the non-electricfield case, the KE increase in the E = 2 × 10 4 V/m case is far greater than the electron energy decrease, and thus a faster ethanol decomposition rate is observed in the E = 2 × 10 4 V/m case.

Discussion
In the present study, atomic and subatomic events are unveiled to shed light on the effect of the electric field on ethanol oxidation reactions. The ethanol oxidation reaction is a two-stage process: stage I features the decomposition of ethanol molecules with the electric field accelerating the ethanol consumption; at stage II, oxygen molecules participate in the reactions and the case with no electric field has a higher ethanol consumption rate. The electric field plays varied roles in the two stages of the ethanol oxidation reaction. At stage I, the electric field alters the ethanol molecular conformations, thereby changing the bond dissociation energy and slightly retarding the reactions. The electric field also energizes the molecules/radicals by increasing their kinetic energy. The changes in the electron energy of the ethanol molecule by the electric field also contributes to the reaction rate differences at stage I. The kinetic energy increase is about 10-fold higher than the electron energy changes and is thus the principal contributor to the elevated reaction rates in the electric field cases. At stage II, with the participation of O 2 molecules, the effects of the electric field reside in modifying the reaction pathways as demonstrated in our previous study [19] . (We have checked the KE changes of C-containing molecules/radicals at stage II, but the changes are not in a similar trend to changes in ethanol consumption rate at this stage. Details are shown in the Supplementary Information.) The results from the present research offer interpretations for previous findings. Our previous study [19] has identified C 2 H 6 O → C 2 H 5 + OH as the principal first step for ethanol oxidation reactions, with C 2 H 6 O → CH 3 + CH 3 O as secondary. The prevalence of the reaction C 2 H 6 O → C 2 H 5 + OH can be attributed to the polarity of the ethanol molecule. As Fig. 4 suggests, the separation of electric charge by hydroxide (OH) enables C -O bond to be the most sensitive point of an ethanol molecule to respond to the external electric field, whereas the impact of the electric field on other parts of the ethanol molecule (e.g. C -C bond length and C -C -O angle) can be neglected in terms of energy variations. Besides, the novelty of the present research also includes the elaborate analysis of dynamics of atoms and structural changes of Please cite this article as: X.Z. Jiang  molecules. The eFF MD calculation of electron energy should also be highlighted, as it enhances our understanding of the influence of the electric field on the ethanol oxidation reactions at the subatomic scales.
A previous experimental study has reported an enhanced fluorescence of flames in the presence of an electric field of ∼10 5 V/m [4] , indicating that the molecules/radicals are at a higher energy level. The present eFF results ( Fig. 5 ) reveal that an electric field of ∼10 5 V/m can energize molecules by elevating electron energies, and thus contribute to a stronger fluorescence. In the introduction section, two plausible hypotheses for the influence of an electric field on flame behavior are mentioned. Indeed, the present study suggests that both hypotheses are valid in the electric-field-assisted ethanol oxidation reactions: stage I satisfies the momentum-driven hypothesis and stage II verifies the chemistry-driven idea.
It is noteworthy that both ReaxFF and eFF consider bonding and bond geometry preferences; they differ from each other in the expressions of energy terms: in ReaxFF, an atom is treated as an unity and the energy is calculated as a function of atom (mainly nucleus) coordinates; whereas in eFF, the factors from electrons are considered and the energy is a function of three variables, i.e. nucleus positions, average electron positions, and average electron sizes. Both force fields aim to retain the accuracy of quantum mechanics (QM) calculations of energies and barriers for reactions while obtaining the results at much faster speeds. For example, the bond dissociation energy of O -H calculated from ReaxFF and eFF are ∼100 kcal/mol [20] and ∼110 kcal/mol [27] , respectively, which are in good agreement with the QM result of 100-120 kcal/mol.

Conclusions
The present research has pioneered the use of eFF MD simulations in combination with ReaxFF MD simulations to provide insights into the influence of the electric field on ethanol oxidation reactions at the atomic and subatomic scales. The ReaxFF MD results have identified the ethanol oxidation reaction as a two-stage procedure where the effects of an electric field on reactions are varied. At stage I, ethanol decomposition is the principal reaction and the electric field affects the reaction by altering the kinetic energy on the order of 100-1000 kJ/mol and influencing molecular conformation and bond dissociation energy. At stage II, oxygen molecules participate in the reaction, and the influence of the electric field is to modify the reaction pathways. The application of the eFF MD simulations, for the first time, extends our understanding of the electric field effects on ethanol oxidation reaction to subatomic scales. The results indicate that the electric field modifies the electron energy on the order of 10-100 kJ/mol. Moreover, the study provides a new interpretation of electric field effects on reaction pathways and experimental observations of enhanced flame fluorescence through eFF MD simulations. It also gives support for the "ionic wind" and chemistry-driven hypotheses, which are confirmed in Stage I and in Stage II of the ethanol oxidation reaction, respectively, in the ReaxFF/eFF MD simulations. This fundamental study adds unique atomic and subatomic insights into the electric field effects and can potentially facilitate the design of realistic electric-fieldassisted combustion systems. In the meantime, it calls for more experimental and QM studies in the field that may provide comparable information at the relevant scales.

Declaration of Competing Interest
None.