A reactive molecular dynamics study of the effects of an electric ﬁeld on n-dodecane combustion

A reactive Molecular Dynamics (MD) study of n-dodecane combustion at high temperatures under externally applied electrostatic ﬁelds is performed to investigate their effect on chemical kinetics. A local charge equilibration method is used to enable charge transfers up to the overlap of the atomic orbitals and introduce molecular polarization induced by an electric ﬁeld. The atomic charges of an isolated n-dodecane molecule with and without external electrostatic ﬁelds are ﬁrst compared with Density Functional Theory (DFT) computations, to assess the accuracy of the charge equilibration method and its ability to capture polarization. Then, the impact of external electrostatic ﬁelds on the reaction kinetics of fuel, oxidizer and products is studied for a range of ambient temperatures and densities. The activation energy and pre-exponential factor of Arrhenius-type reactions under various electrostatic ﬁelds are also investigated by performing Nudged Elastic Band (NEB) computations on selected reactions’ Minimum Energy Path (MEP) and by analysing the collision frequency, respectively. Results show that the atomic charge transfers due to close interactions and molecular polarisation are relatively weak in all investigated conditions, leading to the necessity of strong external electric ﬁelds to induce changes to chemical kinetics. The consumption rate of n-dodecane decreases for strong electrostatic ﬁelds, whereas for low values of the electrostatic ﬁeld strength no clear trend is observed. In addition, at high temperature and density conditions, oxygen consumption increases under strong electrostatic ﬁelds, whereas the opposite trend is observed as the temperature and density decrease. NEB analysis shows alterations of the activation energy up to 2.3 kcal/mol for oxygen compound reactions with varying strength of the external electrostatic ﬁeld. Furthermore, analysis of the translational, rotational and vibrational kinetic energy modes and collision frequency reveals the inﬂuence of translational motion and molecular stabilization on the reaction rates. The kinetics of oxygen molecules was found to be of primary importance to determine the reaction behaviour under external electrostatic ﬁelds, as oxygen molecules have a direct effect on the oxidation reactions and also indirectly affect n-dodecane pyrolysis when an electrostatic ﬁeld is present. This study provides fundamental understanding of the interactions between heavy hydrocarbons and electrostatic ﬁelds for the development of future hybrid thermal-electrical technologies.


Introduction
To reduce the environmental impact of aviation and meet the increasingly stringent regulations on emissions, new propulsion technologies must be developed.Aircraft and engine manufacturers are moving towards electrification of propulsion systems, which offers the possibility of reducing polluting emissions by optimising aircraft architectures and by using clean electrical energy.New design concepts have been introduced in the past few years [1] , including a more synergetic integration of the propulsion system with the aircraft and distribution of propulsion [2] .Although fully electrical solutions could be implemented for shortdistance flights, hybrid thermal-electrical propulsion appears to be the most viable solution for long-distance transportation [1,3] .Therefore, the development of low-emission combustion systems integrated with the electrical propulsion units is an important task to be addressed.The on-board availability of a large amount of electrical energy could open up new possibilities for the control of combustion through electromagnetic fields, possibly improving the combustion performance in terms of flame stability [4,5] and emissions [6] .
Electrostatic fields can influence the combustion process through interactions at the molecular level, resulting in measurable effects at the macro-scales.Experimental studies have shown effects on flame liftoff and stabilization [7][8][9] , blowoff limits [10][11][12] and possibility of altering the flame shape [13][14][15][16] .The impact of electrostatic fields on the reduction of pollutant emissions and soot formation has also been demonstrated [6,[17][18][19][20] .In general, the effect of external electric fields on the flame characteristics can be attributed to two main contributions: (i) the bulk movement of ionic species and electrons induced by the externally applied electrostatic forces and subsequent movement of neutral species due to collisions with charged species (known as ionic wind, a non-reacting effect) [21] and (ii) the influence of the electric field on the kinetics of chemical reactions involving polarized molecules [22] , ions and electrons [23,24] .It has been shown that strong electric fields can affect the catalysis or inhibition even of non-redox reactions based on the orientation of the electric field [22,25] .This is known as electrostatic catalysis and can influence both polar and apolar reactants.Numerical studies conducted so far have mainly focused on the investigation of ionic wind and charge density distribution effects on the flame behaviour, showing some success in predicting the flame response to external electric fields [26][27][28][29][30] .However, such studies usually neglect the direct effect of electrostatic fields on reaction kinetics, not least because of a lack of understanding of charge transfers and the impact of electrostatic forces on molecule dynamics.
To increase the understanding of the effects of electric fields on the chemical reaction pathways, species dynamics under externally applied electrostatic fields have been investigated at the atomistic scale using reactive Molecular Dynamics (MD) simulations.Previous studies using reactive MD [31][32][33][34] showed effects of external electrostatic fields on reaction rates, species consumption and initiation time of the combustion process.Results discussed in these studies were mainly obtained for high pressure and high temperature conditions with relatively strong electrostatic fields, which are hardly found in practical applications.Furthermore, the investigated fuels are not representative of typical aviation kerosene.Therefore, further investigations with systems much closer to aviation applications are needed.It should also be noted that the aforementioned studies employed the Charge Equilibration (QEq) method [35] to compute the charge distribution in the MD domain.Such method could be inadequate to capture the charge transfers between molecules [36] , with a potential impact on the prediction of the electric field effects on the combustion kinetics [37] .Moreover, no assessment of the capability of the method to predict the charge distribution, e.g. through comparisons between MD and reference Quantum Mechanics (QM) computations, is typically reported.
In this study reactive MD simulations of n-dodecane combustion under externally applied electrostatic fields have been performed.The Charge Transfer with Polarization Current Equalization (QTPIE) method [38] is used to compute the atomic charge distribution.This method enables molecular polarization induced by external electric fields and limits long-distance charge transfers [37] .n-Dodecane was selected as a single-component surrogate for aviation fuels, such as Jet A-1 and JP-8.Previous studies in absence of external electrostatic fields have shown that combustion of ndodecane could follow different reaction pathways depending on the conditions of the reacting environment [39,40] .At high temperatures n-dodecane primarily follows a pyrolytic reaction process.Initially, n-dodecane undergoes pyrolysis by C-C bond fission followed by β−scissions, and subsequently the decomposed products are oxidized.At lower temperatures a different reaction mechanism is typically observed.n-Dodecane decomposition is initiated by H abstraction by an oxygen molecule.This leads to the forma-tion of an alkyl radical which subsequently combines with an oxygen molecule leading to the formation of a peroxy radical [39] .In this study we focus on conditions at relatively high temperatures, where pyrolysis is expected to be the most frequent reaction path.Note that low-temperature combustion simulations with the low densities investigated in this study would lead to an unfeasible computational cost.Alternative methods, such as hyperdynamics [41] , must be considered in the future to reduce the simulation time and make simulations at low temperature feasible.
The first objective of this work is (i) to assess the capability of the QTPIE method [38] to compute the atomic charges and to predict molecular polarization induced by an external electrostatic field.To achieve this objective, the charge distribution of an isolated n-dodecane molecule predicted by MD simulations is compared with DFT computations.In addition, results obtained with the QTPIE method are compared with results obtained with the widely used QEq method to assess the effect of charge shielding on the charge distribution.The second objective is (ii) to analyse the electrostatic field effects on the reaction kinetics and reaction pathways of n-dodecane combustion at relatively high temperatures.This objective is achieved by performing large-system simulations of n-dodecane combustion for a range of densities, temperatures and external electrostatic fields, and by evaluating the time evolution of species concentration and reaction kinetics.A further objective is (iii) to provide insight into the effect of the electrostatic fields on the activation energy and the pre-exponential factor of Arrhenius-type reactions.To do that, the activation energy is examined by performing NEB computations for selected reactions under oriented electrostatic fields.Additionally, the translational, rotational and vibrational kinetic energy modes of the system, as well as the collision frequency, are investigated to unveil the electrostatic field effects on the system's kinetics.The structure of the paper is as follows.First, the numerical methods are described, followed by results and discussion.A summary of the work and main conclusions close the paper.

Reactive molecular dynamics
Reactive MD is based on the use of reactive force fields able to accurately describe bond breaking and formation in MD simulations.Force field parameters are tuned using QM computations.The use of reactive force fields enables the simulation of reactive systems containing a large number of atoms with reasonable accuracy at a low computational cost.ReaxFF [42] is a bond order potential that can predict reactions during simulation, since it computes the bond order based on the interatomic distance.The computation of atomic forces comes from the system's total energy [42,43] : E bond corresponds to the bond formation/dissociation energy.E over , E under are the overcoordination and undercoordination energy penalties, respectively; E val and E tor are the valence and torsion angle energies.The non-bonded energy terms E vdWaals and E Coulomb correspond to the van der Waals and Coulomb long-range interactions.Additional correction terms can be included depending on the investigated system to account for accurate description of lone pairs, C 2 molecules and H bonds to name a few.These additional energy terms are referred here to as E additional .Further details on the various energy terms can be found in Ref. [44] .ReaxFF has been widely used in combustion simulations of various hydrocar-bons, showing high accuracy in the prediction of reaction rates and large system dynamics [45][46][47][48][49] .
In this study the CHO-2016 ReaxFF force field developed by Ashraf et al. [50] is used, with a bond order cutoff of 0.3 for the recognition of molecules.This force field has an improved description of small hydrocarbon oxidation and oxidation initiation temperature compared to the previous counterpart by Chenoweth et al. [44] .

Computation of atomic charges
At the atomistic level, the existence of an external electric field directly affects the charge distribution of the atom or molecule, since the positively charged core and the negatively charged electron cloud are pushed in opposite directions.This perturbation of the atomic charge distribution leads to polarization of nonpolar molecules or alteration of the dipole moment of polar molecules.The dipole moment induced by the external electrostatic field follows the direction of the electric field lines (for isotropically polarizable molecules).As a result, the electrostatic forces applied on atomic charges can influence the bond energies of molecules, as well as the rotational and vibrational kinetic energies of neutral species, and the translational kinetic energy of charged species.Strong electric fields could also lead to excitation or even ionization of atoms and molecules.Therefore, redox reactions can become significant and their effect on the overall chemistry evolution of the system should be taken into consideration.
The focus of this study is on the effects of polarisation and local charge transfers between inter-and intra-molecular pairs of atoms, which occurs when there is a spatial overlap of the electronic orbitals [36,51,52] .It is noted that the excitation and ionization states of molecules are not considered in this study, since they cannot be described by standard reactive MD methods due to the absence of an explicit electron description.A few cases of infrequent reactions involving electron collisions can be neglected for low-energy electric fields.However, when electrons acquire sufficient acceleration and thus energy, many avalanche-like reactions can lead to ionization of the neutral species.In the context of gaseous reactive flows the main processes that lead to ionization are ionization by collision, electron transfer between atoms, ionization by transfer of excitation energy and chemi-ionization reactions [21] .An analysis of the electric field strengths for which ionization by collision cannot be neglected, for the conditions investigated in this work, can be found in Section S1.3 of the Supporting Information (SI).Note that possible alternatives to ReaxFF to explicitly include free electrons in MD simulations, such as eReaxFF [53] , are part of current research and development.Their use to simulate ionic species is left for future study.
To study the effect of electric fields on chemical reactions, an accurate prediction of atomic charge distribution is of primary importance.In reactive MD the atomic charge distribution of the domain is calculated at each timestep based on the atomic position and calibrated element properties by applying Sanderson's Electronegativity Equalization (EE) principle [54] with Parr's definition of electronegativities [55] .The electrostatic energy of the system is minimized for fixed atomic positions under the charge neutrality constrain ( i q i = 0 ).Two of the pioneering and most widely used methods for the computation of charge distribution based on this approach are the QEq method [35] and the Electronegativity Equalization Method (EEM) [56] .Compared to the EEM method, the QEq method uses an additional iterative scheme for the computation of the charge-dependent self-Coulomb interaction.Also, the QEq method employs a single normalized n s Slater orbital to describe the atomic density instead of a shielded Coulomb potential, which is used by EEM [35] .The electrostatic energy of a sys-tem of N atoms is given by [57] : where q is the atomic charge, k c = 1 / (4 π 0 ) is the Coulomb constant and 0 the vacuum permittivity.χ and η represent the Mulliken electronegativity and the atomic Parr-Pearson hardness, respectively.J is the screened Coulomb interaction matrix which is modelled with a shielded Coulomb potential.Tap (r i j ) is a sevenorder taper function which allows for the suppression of any nonbonded energy discontinuities above a cutoff distance and γ i j is a shielding parameter between the i -th and j-th atoms.The cutoff distance of long range interactions is usually set to 10 Å.This may affect the accurate description of electrostatic interactions in large systems.When large systems are considered, other methods such as Ewald or particle-particle particle-mesh (PPPM) solvers may be more appropriate.
During the charge equilibration process both the EEM and QEq methods perform a long-range charge transfer despite the interatomic distance.Charge transfers take place even between atoms or molecules at infinite separation distance, therefore treating the system as an ideal conductor [36,38] .As a consequence, neutral molecules with a zero overlap integral with nearby atoms can acquire a net charge and diffuse in the direction of an electric field.Although in some applications this treatment of the charges may be appropriate, in the study of electric field effects in gaseous combustion these charge equilibration methods could lead to an unphysical charge distribution with a significant impact on the results [37] .
In order to shield charge transfers and improve the prediction of the charge distribution, a number of methods have been proposed in the literature, such as the Split Charge Equilibration (SCE) method [58] , the QTPIE method [38] and the Atom-condensed Kohn-Sham density functional theory to Second order approximation (ACKS2) method [59] .These methods locally restrict the interatomic interactions to atoms in close proximity, therefore shielding the charge transfers.If a perfect charge transfer shielding is achieved, isolated neutral molecules do not acquire a net charge, thus the external electric force influences the rotational and vibrational energy components of these molecules but not their translational energy.In this study the QTPIE method is used, which has been shown to be successful in shielding charge transfers in hydrocarbon combustion [37] .The QTPIE method was initially derived using partial charge variables p i j ( q i = j p i j , where q i is the atomic charge) in order to compute diatomic charge transfer, similarly to the SCE method.However, a transformation was later applied in order to shift from the bond space to the atom space description allowing for a mathematical formulation similar to QEq method [60] .In contrast, the bond-space formulation of the SCE method increases the computational cost and requires additional bond-dependent parameters to be calibrated for the system, which are similar to the atomic hardness and atomic electronegativity of the QEq method.The ACKS2 method is based on a more sophisticated derivation from Kohn-Sham density functional theory (KS-DFT) and, to the best of the authors' knowledge, represents the state-of-the-art charge equilibration method.However, its use for the investigation of the effect of electrostatic fields on hydrocarbon combustion is left for future study as the use of such a method requires a complete re-calibration of the CHO-2016 ReaxFF, which was developed using the EEM method.
In the QTPIE method the computation of the electrostatic energy is performed using effective electronegativities in place of atomic electronegativities.The effective electronegativity depends on the difference of atomic electronegativities between pair of atoms.The overlap integral is used as a weighting factor to limit the charge transfers to the atomic orbitals interactions [60] .Specifically, the electrostatic energy is given by Eq. ( 2) , but the electronegativity is replaced with the effective electronegativity given by: where S i j is the overlap integral of the atomic orbitals [61] .In the present work, the overlap integrals are computed with the primitive Gaussian Type Orbitals (GTO).The exponential constants are taken from Ref. [62] , as they were calibrated to reproduce the twoelectron Coulomb integral of Slater Type Orbitals (STO) with a negligible absolute error.
The formulation of the QTPIE method is similar to the QEq method, therefore the same electronegativity and hardness parameters can be used without significant discrepancies [38] .Nevertheless, re-calibration of the ReaxFF with the QTPIE method could further improve the accuracy of the computed charges.A shortcoming of all aforementioned methods, except ACKS2, is that they suffer from the superlinear polarizability issue [63,64] .In particular, as the molecular size increases, the dipole polarizability tensor, α i j = ∂ μ j /∂ E i , does not follow a linear behaviour.However, this is not expected to be an issue for the relatively short molecules investigated in this study (including n-dodecane, which is unlikely to be fully elongated during the simulation).Another limitation of all the methods discussed above is the inability to describe out-ofplane polarization for planar molecules.

Modelling of the electric field interactions
In addition to the global charge transfer, another shortcoming of EEM and QEq methods is the lack of description of polarization induced by external electrostatic fields (however, note that polarization due to close interactions is included in all charge equilibration methods).To introduce the influence of external electric fields on the charge distribution, some studies [65][66][67] utilized the method proposed by Chen and Martinez [68] .In this method, the atomic electronegativity is modified by an additional term r • E , where r is the vector of nuclei position and E is the electric field vector.However, this approach is insufficient to describe polarization for systems of more than one molecule when charge equilibration methods unable to shield the charge transfer are used, due to an unrealistic spatially dependent charge distribution.In the QTPIE method, selected in this study to compute the atomic charges, the effect of external electric fields on molecular polarization is introduced by modifying the effective electronegativity as [62] : where χ is the electronegativity in eV and E is the electric field in V/ Å.The resulting formulation allows the electric field to affect the effective instead of the intrinsic atomic electronegativities, satisfying the nonlocality constraint [69] .As previously mentioned, the QTPIE method does not include a description of out-of-plane polarization for planar molecules.Therefore, considering for example an oxygen molecule, no polarization is induced when the molecule's bond is placed at 90-degree angle with respect to the electric field lines.
The atoms that are immersed in an external electric field experience a Coulomb force equal to: The work done by the electrostatic force to move a charged particle i from a position in space r 1 to a position r 2 corresponds to the negative of the variation of the electric potential energy, which is given by [70] : Therefore, a movement of positively charged atoms in the direction of the electric field leads to a decrease in electrical potential energy.

ReaxFF molecular dynamics simulation setup
Simulations of n-dodecane/oxygen mixtures are performed using LAMMPS software [71] integrated with the USER-REAXC package, and using the QTPIE method for the charge equilibration [37] .A timestep of 0.1 fs is used, which is recommended for reactive MD simulations at high temperatures [44] .The equations of motion are advanced in time using the Verlet algorithm [72] .The canonical NVT ensemble, i.e. constant number (N) of atoms at a constant volume (V) and temperature (T), was chosen.The temperature is kept constant by applying the Nosé-Hoover thermostat [73,74] on the translational degrees of freedom of atoms, with a temperature damping parameter of 100 fs.For post-processing, the atomic connectivities and trajectories are outputted every 20 fs and 1 ps, respectively.The sampling period for the atomic connectivities is kept sufficiently low to avoid merging consecutive reactions [46] .
Cubic domains are filled with 30 n-C 12 H 26 and 555 O 2 molecules (stoichiometric conditions) at densities of 250, 150, 50 mol/m 3 .A large number of molecules is used in order to have statistically representative data to perform averages, but also to observe the influence of the products on the system's kinetics.Periodic boundary conditions are applied in all directions.Simulations are performed at temperatures of 180 0, 210 0, and 240 0 K for the 250 mol/m 3 domain and 2100 K for the 150 and 50 mol/m 3 domains.Simulations were also performed at T = 1500 K for densities equal to 250 and 150 mol/m 3 showing little to no reactivity during a 10 ns simulated time.Thus, such results are not presented herein. .Relatively low values of density and temperature compared to previous studies [46,75] were chosen to be more representative of realistic combustor engine conditions (however, note that combustion with pure oxygen instead of air is investigated).Furthermore, relatively low values of density and temperature allow us to better investigate the effect of external electrostatic fields as the combustion process is not significantly accelerated by the underlying ambient conditions.Cases with strength of the external electrostatic field equal to 0 (base case), 10 −7 , 10 −6 , 10 −5 , 10 −4 , 10 −3 , 10 −2 , 10 −1 , 1 .5 × 10 −1 , 2 × 10 −1 V/ Å are investigated.The external electrostatic field is assumed to act uniformly along a Cartesian coordinate direction.
To generate the initial configurations, energy minimized ndodecane and oxygen molecules were randomly placed inside the domains using Packmol software [76] .Afterwards, energy minimization was applied to the whole configuration, in order to remove any residual artificial potential energy.The conjugate gradient method [77] was used for energy minimization.For each energy minimized configuration, 10 different velocity distributions were generated, sampled from a Maxwell-Boltzmann distribution at the temperature used for equilibration.All systems, prior to simulations with the externally applied electric field, were equilibrated at 10 0 0 K for 4 ns for the densities of 250 and 150 mol/m 3 and 10 ns for the density of 50 mol/m 3 .During the equilibration process, the C-O and H-O bonds were switched-off to avoid oxidation reactions.

Density functional theory computation setup
To assess the accuracy of the charge equilibration methods, reference DFT computations of the charge distribution of an isolated n-dodecane molecule were performed.All geometries are calculated using a hybrid density functional method employing Beckes three-parameters approach, B3LYP-D3 [78] , and the 6-311G * * ++ basis sets [79] .Note that the computation of partial atomic charges is an active area of research [80] .Multiple classes of methods can be used to calculate partial atomic charges [81] , such as population analysis, partitioning the electron density, and electrostatic fitting.Hence, DFT charge calculations are performed using three methods from different classes for comparison.We present the charges calculated using Mulliken population analysis [82] , Stockholder charges from Hirshfeld partitioning [83] and the ESP method [84,85] , corresponding to the classes of methods in the same order as stated above.Please note that the Stockholder charges are calculated using the 6-31G basis set.The n-dodecane molecule was first geometry optimized using the 6-311G * * ++ basis set and the resulting geometry was used to calculate partial atomic charges.The Jaguar software package [86] is used for all DFT computations without external electrostatic fields.
DFT computations are also performed under external electrostatic fields, using the ADF software suite [87] .The B3LYP-D3 method is employed again with the QZ4P basis set.The charge calculation is performed using the Hirshfeld partitioning.Note that DFT methods may lead to significant errors in charge calculations.Calculations using more accurate methods such as Coupled Cluster Theory (CCT) are required to reduce these errors [88] .These are left for future study.

Computational details of Nudge elastic band analysis with climbing image method
In order to assess the influence of the external electric force and inter-and intra-molecular polarization on the transition state of reactions, NEB [89] computations are performed with the Climbing Image (CI) [90] method.Given an energy minimized initial and final state, the MEP of a reaction can be obtained via the simultaneous minimization of intermediate replicas.A spring force is applied between intermediate replicas to fix their position on the transition path.The image with the highest energy further climbs the potential energy surface to find the saddle point that corresponds to the activation energy of the reaction.
NEB computations were performed on LAMMPS [71] .A total of 6 replicas are used in the analysis using a spring force of 1 kcal/(mol • Å) between adjacent images.The additional potential energy due to the electric forces is included in the system's total potential energy, which is the quantity that is minimized.Both pyrolytic and oxidation reactions were studied.

Reaction kinetics analysis
ChemTraYzer code [46] is used to compute the reaction rates from MD simulations.ChemTraYzer is based on first-order kinetics to derive the rate constants of all reactions occurring in a single simulation.The reaction rate constants derived from each simulation are Poisson-averaged over all simulations.Poisson average is used as rare reaction events follow the Poisson distribution rather than the normal distribution [91] .Therefore the rate constant of a reaction is computed according to the following equation: where N is the number of reaction events, V is the domain volume and C is the product of the reactants' concentrations.t is the time interval over which the reactants concentrations, and thus the reaction rate, are constant.S is the total number of simulations and M is the number of pieces (divisions) of the i -th simulation where the reaction rate is constant.The minimum and maximum values of the error bars shown in the rate constant graphs represent the confidence intervals with a 90% confidence.
For the recognition of reactions a recrossing filter is applied.Backward and forward reactions that take place within the filtering period are assumed to be recrossing events and are excluded from the analysis.The filtering period was chosen based on the collision frequency of the mixture from the kinetic theory of gases.For a mixture of two molecules A and B , the collision frequency based on the Hard Sphere (HS) model is given by [92] : where N A and N B are the number of A and B molecules, respectively.d is the kinetic diameter.The kinetic diameters of the fuel and oxygen are d f = 7 .12 Å [93] and is the reduced mass.The filtering period is computed based on the collision frequency of the fuel and oxidizer (initial concentration) and is chosen to be 0.2 ps, 0.3 ps and 0.9 ps for the density of 250 mol/m 3 , 150 mol/m 3 and 50 mol/m 3 , respectively, keeping the collision probability around 60%.The filtering period and collision probability are summarised in Table S6 of the SI.
To examine the consumption rate of fuel and oxygen, first-order kinetics is used.Therefore, the rate constant, k , of species consumption is computed via the first-order integrated rate law given by: ln where N 0 and N t are the number of molecules at the initial state and at time t, respectively.The final value of the consumption rate constant is the slope resulting from the linear curve fitting of the ln (N 0 / N t ) term over time.

Collision frequency analysis
Collision theory links the pre-exponential factor of the Arrhenius expression for the reaction rate with the number of collisions.Specifically, the reaction rate of a bimolecular reaction between species A and B can be expressed as [61] : where r is the reaction rate, Z = f /V is the collision frequency per unit volume , ρ s is a steric factor and A and E a are the preexponential factor and the activation energy, respectively.Therefore, the reaction rate is directly proportional to the number of collisions in the simulation.The steric factor is a correcting factor between experimental and theoretical results, which takes care of the fact that the reaction probability depends on the certain mutual orientations of the reactant molecules and that all the collisions do not have the same probability of leading to a reaction.To assess whether the external electrostatic field affects the reaction rate through a change in collision frequency, the collisions between species are directly counted in MD simulations.The counting of collisions between molecules was directly implemented into LAMMPS and performed at runtime in each simulation.Molecular collisions, i.e. collisions between nearby molecules, are measured based on interatomic distances.Specifically, a binary collision between molecules is counted if the distance between an atom of the first molecule and an atom of the second molecule is below a threshold value.To define that threshold value, the van der Waals atomic radii are used [95] .Therefore, the atomic radius of collision of C, H and O atoms are taken as 1.7, 1.2 and 1.52 Å, respectively.This approach of interatomic distance was chosen instead of measuring the momentum change in order to be able to track binary, triple, quadruple and higher order collisions, which are expected for the relatively long n-dodecane molecule initially surrounded by many oxygen molecules.Taking into account the analysis based on the kinetic theory of gases presented in the SI, we chose a time interval of 10 fs for the detection of collisions, sufficiently small to account even for the fast moving particles.To accurately analyse the collision events, selected cases were rerun with a data collection dumping period of 10 fs.In particular, results were collected for the base case ( E = 0 V/ Å) and a strong ( E = 0 . 2 V/ Å) electrostatic field at T = 2100 K and ρ = 250 mol/m 3 .
In addition, the kinetic energy components, namely the translational, rotational and vibrational kinetic energies, were directly computed during simulation.The total kinetic energy of the system is computed from the atomic masses and velocities.The translational energy is computed based on the molecular masses and velocities.The molecular velocity is found by projecting the atomic velocities on the molecule's center of mass.The rotational energy is calculated from the angular velocity of each molecule, and the vibrational energy is computed from the residual atomic motion.During bond formation and dissociation abrupt changes in the kinetic energy components occur.For example, consider the case where the bond order of an n-dodecane molecule's atom and a nearby oxygen atom surpasses the 0.3 threshold value set in the current study.In that case, a bond is created and the oxygen atom is included in the molecule's analysis.Given that the oxygen atoms are more likely to have a higher velocity, an increase in the kinetic energy components of the newly formed molecule will be observed.These peaks were filtered out in the analysis presented in this work.The filtering method is based on the time evolution of one of the kinetic energy components, and then applied simultaneously to the other components, in order to conserve the total kinetic energy.Specifically, the time evolution of the rotational energy was used to guide the filtering operation.When peaks in rotational energy surpassed a difference of 10 0 0 kcal/mol between the current dumping timestep and the average energy value of the previous 100 dumping timesteps, peaks of all energy components were removed.The very small dumping time interval chosen for the recognition of collisions was also used for the kinetic energy modes.This guarantees that the extreme values, which are filtered, are not due to actual system's kinetic energy changes, but due to the aforementioned variations arising from the bond recognition.

Zero-dimensional reactor analysis
To assess the MD results, the rate constants obtained with MD simulations are compared with results from homogeneous zerodimensional reactor simulations.Zero-dimensional reactor simulations, performed using the Cantera software [96] , were advanced for the same time period as the MD simulations, i.e. 10 ns simulated time.Computations are initialised with a stoichiometeric mixture of n-dodecane and oxygen.The temperature and volume of the reactor are held constant throughout the simulation to match the conditions of the MD simulations.The chemical mechanisms used in the computations are the JetSurF2.0[97] and the CaltechMech [98] for high temperature n-C 12 H 26 combustion.
The JetSurF2.0 has 44 reactions involving n-dodecane dissociation, 7 of which are unimolecular pyrolysis reactions.The rest of these reactions correspond to reactions of n-dodecane and oxygen compounds.By neglecting the position of the hydrogen atom extracted from the n-dodecane molecule, reactions that lead to the same products can be grouped together leading to 5 unique reactions between n-dodecane and oxygen compounds.The Caltech-Mech mechanism has 6 reactions of n-dodecane where C-C β- Lines are added for ease of following data points.The atomic id's are shown in Fig. 1a.scission occurs, whereas reactions involving H abstraction by an oxygen molecule are not included in the mechanism.Moreover, the methyl group thermal cracking reaction is also not part of the Cal-techMech mechanism.

Charge distribution
We begin the analysis by comparing the charge distribution of isolated molecules computed by MD with the QEq and QTPIE charge equilibration methods (based on the CHO-2016 ReaxFF parameters) and by different DFT simulations (Stockholder, Mulliken and ESP charge partitioning methods).The atomic charges of an n-dodecane molecule are presented in Fig. 1 b, where each atom is identified by an integer number as labeled in Fig. 1 a.The same molecular conformation was used in all computations.The atomic charges predicted by MD with the two charge equilibration methods are in good agreement with the Stockholder charge distribution.On the contrary, the Mulliken and ESP methods predict a very different trend of charge distribution.The Mulliken charges are approximately 8 times larger than the charges obtained with MD.However, it is known that the Mulliken charges may not reproduce the dipole and higher moments as accurately as ESP [44] .On the other hand, ESP is typically very sensitive to the conformation of large molecules, which may explain the relatively large fluctuations of atomic charges observed in this study.It should be noted that the differences between the charges obtained with QTPIE and QEq methods appear negligible when MD results are compared with DFT results.Only minor differences between the charges predicted by the two MD methods are observed, with the QTPIE method predicting lower atomic charges as a consequence of the shielding of charge transfers.This suggests that using the same ReaxFF pa- rameters with both QEq and QTPIE methods leads to reasonably accurate results.Even though the charges predicted by the two charge equilibration methods are very similar for most molecules, larger differences could be found for highly polar molecules.This is for example the case of H 2 O .As shown in Table S4, the QEq method provides a better estimation of the reference value of the H 2 O dipole moment compared to the QTPIE method.Nevertheless, the concentration of H 2 O and other highly polar molecules is sufficiently low in all investigated cases (e.g., at ρ = 250 mol/m 3 and T = 2100 K there are less than 4 H 2 O molecules, averaged over the 10 simulations).Thus, the lack of accuracy to predict the dipole of highly polar molecules is not expected to significantly affect the results.Note that this lack of accuracy could also affect, to some extent, the charges of polar transition states of reactions.The related effects on the kinetics of the system should be carefully studied in future work.Re-calibration of ReaxFF using the QTPIE method instead of the EEM could improve the prediction of dipole moments in MD computations with the QTPIE method.
Figure 2 shows the atomic charges and the related (internal) electrostatic potential of an isolated n-dodecane molecule under different external electrostatic fields.Results obtained with both the QTPIE method and DFT computations are presented.The elec-trostatic potential is shown on the accessible surface area, and it is derived from the atomic charges using the Multiwfn software [99] .All results were obtained using the same molecular conformation.Charges predicted with the QEq method are not presented here since no significant effects of external electrostatic fields on the charge distribution were found.As already discussed, the QEq method does not include polarization induced by an external electrostatic field and thus the external electric field effects are limited to the application of Coulomb forces.However, Coulomb forces do not considerably alter the bond length to result in significant differences of the charge distribution, even under E = 0 . 2 V/ Å. Conversely, polarization is evident with the QTPIE method under strong external electrostatic fields (note that polarization is negligible for electrostatic fields below 10 −3 V/ Å).Polarization follows the direction of the electric field lines and is proportional to the field's strength.The dipole moments of various neutral molecules are reported in Table S4 with and without external electrostatic fields.The bond lengths of isolated molecules under 0 and 0.2 V/ Å external electrostatic fields were also investigated.The Coulomb force acting on a bonded atom may lead to elongation or shrinking of the bond, thus causing variations in bond energies [44] .As shown in Tables S2 and S3, variations in bond lengths are negligible (of the order of 0.1% for n-C 12 H 26 and close to 0% for O 2 ).Therefore, no significant variation of the bond energy is expected.
MD-QTPIE and DFT results show a very similar variation of partial charges due to induced polarization compared to the charges computed with no external electric fields (see Fig. 2 ).Specifically, the QTPIE method predicts an average 84% increase or decrease of the end H atoms' charge (atoms H14 and H37 in Fig. 1 a) for an electric field of 0.1 V/ Å along the x-coordinate.The corresponding variation with the Hirshfield DFT method is 89% on average.Differences in the polarization behaviour predicted by DFT and MD-QTPIE computations with electric field oriented along the xcoordinate are limited to the charges of the two methyl groups.As shown in Fig. 2 (and in more detail in Figure S4 of the SI), MD-QTPIE computations predict a weak local polarization in some of the atoms that are located close to the end H atoms of the ndodecane molecule.The same behaviour is not evident in the DFT computations, which give a more uniform variation of the charges at the two ends of the n-C 12 H 26 molecule when an electric field is applied.The difference between the charge distribution predicted by MD and DFT can be attributed to the shielding technique of the QTPIE method, i.e. to the overlap of the atomic orbitals which for the studied atoms is around 3-4 Å, as a similar behaviour is not evident for small molecules or for n-dodecane with an electric field in the y-and z-directions.
Further investigations on isolated molecules show that the effects of polarization on the internal electrostatic field (and consequently on long range interactions) of an n-dodecane molecule are confined to within a 10 Å distance from the molecule (see Section S1.1).It is also shown that the QTPIE method successfully predicts neutral oxygen atoms during O 2 dissociation (Section S1.2).
We conclude the analysis of charge distributions by comparing the charges predicted by MD with the QEq and QTPIE methods for the full system, i.e. a domain comprising n-dodecane and oxygen molecules in stoichiometric conditions, at ρ = 250 mol/m 3 .
Figure 3 shows the charge distribution of the domain under various external electrostatic fields after performing a charge equilibration on the energy equilibrated atomic positions.Figure 3 demonstrates that the charge transfers have been successfully shielded with the QTPIE method.The QEq method produces significantly high net molecular charges, whereas the QTPIE method preserves (with good accuracy) the charges of the isolated molecules when these do not interact with nearby atoms.When an external electric field is applied, the QTPIE method predicts polarization of the molecules, which depends on both the strength and the direction of the electric field.This is the reason of the dispersed distribution of charges observed in Fig. 3 for O 2 molecules, as they acquire their maximum and zero polarization when they are placed at an angle of 0 and 90 degrees, respectively, with respect to the electric field lines.It should be noted here that the net charge of formally neutral molecules predicted by the QTPIE method, although close to zero, is not exactly zero.A net charge of the order of 10 −2 e for the n-dodecane molecule and 10 −3 e for the oxygen molecule is found for the zero electric field case of Fig. 3 a, even when the molecules are placed far away from all the others (distance where the overlap integral tends to zero).This weak charge transfer occurs due to the use of a global charge equilibration process while forcing a locally shielded charge transfer.Additionally, the charge neutrality constrain is imposed for the system as a whole instead of individual neutral atoms.Although substantially small compared to the charge transfers due to close interactions (see Figure S3), this residual net charge on isolated molecules may contribute to their diffusion under the external electrostatic forces.

Species consumption/production
The time evolution of n-dodecane and oxygen concentration is shown in Fig. 4 for cases at ρ = 250 mol/m 3 and T = 2100 K. n-Dodecane is generally consumed within 4 ns for all investigated cases.For weak external electrostatic fields, below a threshold value of about 10 −3 V/ Å, no significant effect of the electric field on the consumption rate of n-dodecane is observed.It is interesting to note that some of the results for weak electric fields suggest a slight increase of the average consumption rate of n-dodecane (see Fig. 5 ).However, results are characterised by a relatively high dispersion (high standard deviation) making it difficult to draw definitive conclusions on the trend.Above the threshold value of 10 −3 V/ Å the consumption rate of n-dodecane drastically decreases.Simulations of systems with only n-dodecane molecules (see Section S2.3) show that the electric field effects on n-dodecane pyrolysis are less significant when oxygen molecules are not included in the simulation domain, even under E = 0 . 2 V/ Å.Therefore, the  kinetics of the oxygen indirectly affect the pyrolysis process of ndodecane.
Results for various intermediate hydrocarbons (propene, methyl radical, ethylene, methane and acetylene) are presented and discussed in Section S2.1.In general, from the simulations performed here, we can conclude that for strong electrostatic fields ( E ≥ 0 . 1 V/ Å) the production rate of hydrocarbons (consumption rate in the case of n-dodecane) consistently decreases compared to the value found for lower strengths of the electric field.Conversely, for weak electrostatic fields ( E < 0 . 1 V/ Å) the behaviour of the production/consumption rate depends on the given species, although no clear trend can be identified due to the high dispersion of the data.
Contrary to hydrocarbons, oxygen molecules (see Fig. 4 ) experience a systematic increase in reactivity at high electric field strengths.The consumption rate of oxygen, shown in Fig. 5 , remains relatively unaffected, within the confidence intervals, for low external electrostatic fields, whereas for E > 0 . 1 V/ Å a drastic increase of the reaction rate is observed.Results for oxygen compounds are shown and discussed in Section S2.1.Similar to oxygen molecules, the production of CO and HO 2 molecules seems to increase under strong electrostatic fields.An increase of the CO 2 to CO ratio and HO 2 to HO ratio under strong electrostatic fields is also observed, suggesting that the presence of an external electric field facilitates complete combustion.The high reactivity of oxygen under strong external electric fields can be attributed to the strong polarization of the highly electronegative oxygen molecule and oxygen compounds.This strong induced polarization leads to variations in the activation energy of reactions involving oxygen, as will be further discussed in the following sections.The different effect of the electric field on the consumption rates of n-dodecane and oxygen could affect the time evolution of intermediate species (such as of H 2 , see Section S2.1).It should be noted that the strong external electric fields studied in this work are above the breakdown threshold (as described in Section S1.3).Therefore, at such conditions, a complete description of the physics requires to consider additional phenomena, such as excitation and ionization of atoms, as well as free electrons.However, it is also important to note that the computation of atomic charges is characterised by a relatively high uncertainty, as demonstrated by DFT computations.Since the Coulomb force on the atoms is directly proportional to the charge values, the results found here under strong electrostatic fields could be indicative of the system's behaviour at lower electric fields in case of higher atomic charges.
The response to external electric fields is also influenced by the temperature of the system.Figure 5 , top, shows the n-dodecane and oxygen first-order consumption rates for different temperatures at ρ = 250 mol/m 3 .Under strong electric fields n-dodecane experiences a decrease in consumption rate for all temperatures investigated here.This decrease becomes larger as the temperature of the system reduces.Conversely, no significant changes of n-dodecane consumption rate are observed for weak electric fields for all investigated temperatures.Note that a higher temperature implies an accelerated kinetics, which reduces the effect of the molecule acceleration due to the external electrostatic forces.As far the oxygen molecules are concerned, at weak electric fields no clear trends for the consumption rates are observed.Under strong electric fields, a systematic increase of oxygen consumption with increasing strength of the electric field is found for T = 2100 K and 2400 K, whereas a decrease is noted for T = 1800 K.This decrease in oxygen consumption at T = 1800 K is caused by the significantly reduced availability of n-dodecane pyrolysis products for the oxygen molecules to react with, which also has an effect on oxygen compound reactions.
The effect of ambient density on the response of the system to different strengths of the external electrostatic field has also been evaluated.Figure 5 , bottom, shows the first-order consumption rates of n-dodecane and oxygen molecules under different densities for a temperature of 2100 K.For all studied densities, the n-dodecane consumption rate shows a significant reduction only for strong electrostatic fields, whereas for weak electrostatic fields no clear trend is observed compared to the case with no electric field.This is shown in more detail in Figure S7.The effect of strong electric fields becomes more evident as the density of the system reduces.Oxygen consumption rate at 250 mol/m 3 and 150 mol/m 3 seems to be unaffected when weak electric fields are applied, as any changes of the reaction rate are within the error bars.Conversely, for a density of 50 mol/m 3 weak electric fields lead to an increase of the oxygen consumption rate compared to the value found for no electric fields.At high electric fields there is a significant increase in the number of consumed oxygen molecules for densities of 250 mol/m 3 and 150 mol/m 3 .On the contrary, for 50 mol/m 3 the oxygen consumption is decelerated.This can be attributed to the lack of pyrolysis products leading to less reactions with oxygen molecules, thus decreasing the reaction rate of oxidation reactions.
A comparison between the fuel and oxidizer consumption rates predicted with QTPIE and QEq methods is presented and discussed in Section S2. 4. Results suggest that the different charge distribution predicted by the two methods can greatly influence the consumption rate of oxygen.In particular, polarization induced by the external electrostatic field, neglected by QEq method, could affect the activation energy of oxygen compound reactions as will be discussed in Section 3.4 .Future work will aim at improving the charge description in MD simulations under externally applied electrostatic fields with a focus on oxygen molecules, which seems to be crucial for the progress of the system's reactivity.
The influence of external electrostatic fields on the initiation time of reactions involving n-dodecane, oxygen and selected products is analysed next through Fig. 6 .The initiation time denotes the time when the first fuel or oxygen dissociation occurs.For products it is the time till the first molecule is formed.To extract the initiation time, the output from ChemTraYzer used.Any recrossing events are excluded from the analysis.The external electrostatic field seems not to affect the initiation time of both n-dodecane and oxygen, suggesting that the external electrostatic forces do not considerably affect the initial stages of reactions.This result is in agreement with the study of Ref. [37] .Conversely, the effects of electrostatic fields on the initiation time are more evident for smaller hydrocarbons where clear variations on the timing of their first appearance are observed.These variations are caused by the change of reaction dynamics of fuel, oxidizer and other products participating in the combustion process.

Reaction pathways
Figure 7 shows the mechanism of pyrolysis down to 4-carbon hydrocarbons.Branches of at least 4% branching ratio are shown for the initial stage of pyrolysis and above 30% for the lateral stages.The branching ratio of branch 7 was up to 4%, thus any subsequent reactions are excluded from the figure.The n-dodecane pyrolysis undergoes mainly a β-scission and secondary a C-H bond dissociation.By examining the branching ratios (see Table S5), no influence of the electric field on reaction pathways was found.The reason of this behaviour may be attributed to the fact that the electric field, even at high strengths, does not affect the thermal cracking energies of n-dodecane as discussed in the following.
A comparison between the reaction rates of the main pyrolytic reactions for different strengths of the external electrostatic field is presented in Fig. 8 a.The results are also compared with the widely used JetSurF2.0[97] and CaltechMech [98] mechanisms.Very good agreement between MD results and the rate constants from such mechanisms is found under no external electric fields.The external electrostatic field affects the rate constants of ndodecane pyrolytic reactions with a trend similar to the one observed for the first order consumption rate of n-dodecane.Under strong electric fields the thermal cracking of n-dodecane slows down, whereas for weak electric fields no clear trend can be identified.
The reaction rate constants of selected reactions involving oxygen molecules are shown in Fig. 8 b.Although good agreement with the reference reaction mechanisms is found for HCO and n-CH 3 O dissociation reactions, and n-CH 2 O + H bimolecular reaction, the rest of the reactions are significantly accelerated (higher rate constants) in the MD results.However, it should be noted that differences of an order of magnitude can also be observed between the reaction rate predicted by the two different chemical mechanisms (e.g., rate constants of the n-CH 3 O → n-CH 2 O + H reaction).For weak electric fields the selected reactions show different trends.For example, the formation of water and consumption of hydroperoxyl radical is slightly accelerated by an increase of the strength of the electrostatic field, whereas the formaldehyde and hydrogen reaction is decelerated.At strong electric fields almost all reactions of oxygen compounds exhibit an increase of the reaction rate compared to the values found for weak electric fields.Results for the other temperature and density conditions investigated in this study can be found in Figure S10.It should be noted that very few oxidation events of n-dodecane take place in the present MD simulations (Figure S5).In addition, reactions between n-dodecane and oxygen compounds are not observed in the MD simulations since no oxygen compounds are formed before n-dodecane molecules undergo pyrolysis.This is consistent with the fact that at high temperatures the decomposition of large hydrocarbons is more likely to follow a pyrolysis path rather than direct oxidation [39,40] .

theory interpretation
Considering the collision theory, the electric field may have two distinct effects in reactive MD simulations.On the one hand, the electric field may alter the activation energy of a reaction, due to changes in polarized molecules, polar species, polar transition state conformations, bond energy and electrostatic forces.Previous studies showed that strong oriented electric fields can affect the electrostatic catalysis or inhibition of non-redox reactions depending on the direction of the electric field for both polar and apolar reactants [22,25] .On the other hand, the electric field may lead to variations of the molecular translational, rotational and vibrational kinetic energy, therefore affecting the collision frequency.These two possible interaction mechanisms are further examined below.As discussed in Section 3.1 , the external electric field perturbs the molecular charge distribution leading to polarization of molecules.To examine whether the variation in inter-and intramolecular charge distribution and applied forces has an effect on the activation energy of selected reactions, NEB computations were performed.The accuracy of the method is first evaluated by comparing the present results with the literature in the case of no external electric fields.Activation energies of n-dodecane pyrolysis reactions are given in Table S8.Due to symmetry of the ndodecane molecule, it is sufficient to examine the C-H bond dissociation from 7 instead of 26 H atoms. NEB results seem to be in relatively good agreement with the literature.Specifically, energy barriers are in good agreement with data from the first method of Ref. [100] , whereas approximately a 10 kcal/mol difference is found when results are compared with the other cited data.The influence of the external electrostatic field on the MEP of n-dodecane pyrolysis reactions was also analyzed with NEB computations.Computations were performed for E = 0 , 0 .01 , 0 .05 , 0 . 1 , 0 . .No significant influence of the external electrostatic field on the MEP of n-dodecane pyrolysis reactions is observed.Therefore, the activation energy of n-dodecane pyrolysis reactions seem not to be affected by neither the magnitude nor the orientation of the external electrostatic field.As a result, according to the collision theory, the changes that are observed on n-dodecane reactivity under different external electrostatic fields should be attributed to changes in the collision frequency.
Differently from pyrolysis reactions, the MEP of reactions involving oxygen can be affected by external electrostatic fields.Figure 9 shows the MEP of the methane and oxygen reaction forming a methyl and hydroperoxyl radical.The atomic charges of the O 2 molecule and the H atom that result in the formation of HO 2 are also shown.Under no external electrostatic fields, the activation energy of the reaction is 61.2 kcal/mol, which is close the theoretical value of 52-56 kcal/mol of Ref. [101] .Under strong external electrostatic fields, the activation energy of this reaction increases or decreases depending on the orientation of the electric field.Variations up to 0.9 kcal/mol are observed.When the external electrostatic field is oriented in the y-direction larger effects on the activation energy are observed compared to the case with the external electrostatic field in the x-direction.Therefore, the effects of the electrostatic field are more evident in the configuration where the oxygen molecule has a very weak polarization (placed at 11 degrees angle with respect to the x-direction).Figure S12 also shows the analysis of the energy variation of the same reaction when the oxygen molecule has fixed charges, corresponding to the maximum O 2 polarization under a 0.2 V/ Å electrostatic field.Results show no significant energy variations.Therefore, we could conclude that the variations of the MEP with the external electrostatic field can be attributed to the atomic charges that the O 2 molecule and H atom acquire during the formation of the HO 2 molecule.Specifically, as observed in Fig. 9 , the atomic charges increase significantly during the formation of the H-O bond.The consequent higher electric forces result in variations of the activation energy.The relatively small charges of O 2 due to polarization also explain why the external electrostatic field seems not to in- fluence the initiation time even at high strengths of the electric field.On the contrary, the relatively high polarization of the products results in significant variation of their initiation time.Results for other selected reactions, typically investigated in many combustion studies, are discussed in detail in Section S4.2.In general, under no external electric field the activation energy of the reactions is in good agreement with the literature.Additionally, a clear effect of the external electric field on the activation energy of selected oxygen compound reactions is found, with variations up to 2.3 kcal/mol.Also, it is evident that changes in the MEP are not only affected by the strength of the external electric field but also by the orientation of the oxygen compound molecules and therefore by the orientation of the dissociating and forming bonds.
To further examine the influence of the external electrostatic field on the reactions kinetics, the translational, rotational and vibrational kinetic energies of the system are analysed.Results are presented in Fig. 10 for the case at T = 2100 K and ρ = 250 mol/m 3 and for zero and a strong external electrostatic field.The total kinetic energy of the system is approximately 14080 kcal/mol.Under no external electric field, the vibrational energy represents the largest contribution to the total kinetic energy, whereas the rotational energy contributes to the smallest percentage.The translational energy of the system exhibits an initial decrease due to the perturbation of atomic velocities after the energy equilibration, which leads to an increase in vibrational energy.Afterwards, during the n-dodecane pyrolysis and the production of smaller molecules, the translational energy gradually increases until asymptotically reaches a maximum value of about 4650 kcal/mol.The increase in number of products also leads to a rapid increase in rotational energy, up to 3600 kcal/mol.The vibrational energy follows the opposite trend, decreasing from 7300 to 5800 kcal/mol.
The application of a strong external electrostatic field ( E = 0 . 2 V/ Å) significantly changes the time evolution of the kinetic energy components compared to the case with no external electric field.An overall increase of translational kinetic energy and decrease of the other kinetic energy components are observed (as shown in Fig. 10 ).The comparison between cases with and without an external electrostatic field is further conducted by analysing the time-averaged kinetic components in three consecutive time periods (each corresponding to 1/3 of the total simulated time, as also indicated for reference in Fig. 10 ).The results are presented in Table 1 .The translational kinetic energy is on average 10 0 0 kcal/mol higher when the external electrostatic field is applied, whereas a decrease of rotational and vibrational energies of the order of 300 and 700 kcal/mol, respectively, is observed.The increase in translational kinetic energy occurs mainly due to frequent charge transfers between molecules that closely interact.As products with a smaller size are formed, charge transfers are more frequent, which leads to larger effects of the electrostatic forces on the translational movement of smaller molecules.The acquired net charges, which can be significant, introduce a net electrostatic force that diffuses the molecule in the direction of the electric field, if the net charge is positive, or to the opposite direction, if the net charge is negative.This movement of charges leads to a decrease of the electric potential energy.The decrease of the rotational and vibrational energies suggests a stabilization of the molecules under external electrostatic fields.Although this result could be partly affected by the use of the Nosé-Hoover thermostat to impose a constant temperature, external electrostatic fields could also provide a stabilisation mechanism.In more detail, molecules with an inherent dipole moment or molecules with a polarisation induced by the external electric field have the tendency to orient their dipole moment parallel to the electric field lines due to the applied external electrostatic forces.In addition, the continuously applied electric force leads to damping of the bond oscillations, similarly to the decay curve of a common spring-  mass system when an external force is applied.Electric field stabilization effects and their isolation from the thermostat technique should be further investigated in future studies.It is noteworthy that for most of the simulated time (from 0 to 7 ns), fewer reactions take place under strong electric fields compared to the case with no external electric fields.This is also evident in Fig. 11 where with strong external electrostatic fields the increase of the system's potential energy is low compared to the zero electric field case.Fig. 12 shows the time evolution of the collision frequency.In absence of electric fields, the collision frequency is approximately 10 12 Hz.This value is in good agreement with the collision frequency of 5 ×10 12 Hz predicted by the kinetic theory of gases for a mixture of two molecules using the HS model, which assumes elastic collisions.For an external electrostatic field of 0.2 V/ Å less collisions take place in the initial 7 ns of simulated time compared to the zero electric field case.This decrease in collision frequency leads to the decrease in n-dodecane consumption rate discussed in the previous sections.On the contrary, above 7 ns an increase in the number of collisions is observed.This result can be attributed to the increase of oxygen consumption and oxygen compounds formation, leading to more product species in the domain and thus more collisions.On average there is a 4.1% reduction of the total collision count between the 0.2 V/ Å electric field case and the case with no external electric field.The change in number of collisions affects the pre-exponential factor of the Arrhenius-type equations.Therefore, the effect of the electrostatic field on the reaction kinetics is primarily due to changes in the number of collisions, in addition to changes in the activation energies of the oxygen compounds reactions.

Summary and conclusions
A reactive MD study of n-dodecane combustion under external electrostatic fields has been performed.The QTPIE charge equilibration method was used, which shields the charge transfers up to the overlap of the atomic orbitals.The charge distribution of ndodecane for various strengths of the external electrostatic field was compared with DFT computations and with results obtained with the widely used QEq charge equilibration method to verify the accuracy of the MD method used in this study.Comparison with DFT computations showed that the QTPIE method can predict shielding of charge transfers and molecular polarization due to an external electrostatic field with sufficient accuracy.
Stoichiometric simulations of large systems of molecules for a wide range of electrostatic fields were performed to analyse the effect of external electric fields on the species and reaction kinetics.The electric field effects were studied under different ambient temperatures and densities.Results suggest that charge transfers due to close interactions and molecular polarisation induced by the external electrostatic field lead to relatively low atomic charges.Therefore, low-energy electric fields (below 10 −3 V/ Å) do not significantly affect the system's kinetics or chemical reactions.For high strengths of the external electrostatic field (above 10 −3 V/ Å) the consumption of n-dodecane is considerably hindered.The same trend is observed for all temperatures and densities.The oxygen consumption rate under strong external electric fields depends on the ambient conditions.Under high temperature and density conditions, oxygen reactivity increases with increasing strength of the electrostatic field.For lower temperature and density conditions, the inverse behaviour is observed, which is related to the lower availability of pyrolysis products in the system.It is noted that the initiation time of the fuel and oxidizer reactions, and therefore the initial stage of combustion, is not affected by the external electrostatic field.On the contrary, a clear effect of the electric field on the initiation time of products is observed.
To give further insight into the effect of electrostatic fields on the system's kinetics, the activation energy of key reactions and kinetic energy components were investigated.Analysis of the reaction transition paths using NEB computations shows that the electric field does not affect the reaction pathways of n-dodecane pyrolysis, whereas reactions involving oxygen show a difference in the activation energy up to 2.3 kcal/mol.The analysis of the translational, rotational and vibrational kinetic energy modes reveals an increase in translational movement in the direction of the electric field, together with a decrease in rotational and vibrational energies.These are linked to changes in the collision frequency.
These findings provide a fundamental understanding on the effects that charge transfers and molecular polarization have on the combustion process of heavy hydrocarbons when external electrostatic fields are applied.Results suggest that relatively high atomic charges are necessary to enable manipulation of the system's reactivity with relatively low energy electric fields.For the atomic charges predicted by the present MD study of n-dodecane combustion, strong electric fields are necessary to affect the kinetics and chemical reactions of the system.

Fig. 1 .
Fig. 1.(a) Atomic labels and orientation of n-dodecane molecule used as reference for the study of the isolated molecule; (b) atomic charges of n-dodecane molecule computed with QEq and QTPIE charge equilibration methods and various DFT charge partitioning schemes (Mulliken population analysis, Stockholder charges from Hirshfeld partitioning and the ESP method).Lines are added for ease of following data points.The atomic id's are shown in Fig. 1a.

Fig. 2 .
Fig. 2. Atomic charges in multiples of electron charge (1.0 is a proton) of an n-dodecane molecule under external electrostatic fields of different strengths and directions.Results are presented for MD-QTPIE (left column) and DFT-Hirshfield (right column) computations.The electrostatic potential derived from the atomic charges is coloured on the accessible surface area.ˆ e x , ˆ e y , ˆ e y and ˆ e z are the unit vectors in the x, y and z direction, respectively.

Fig. 3 .
Fig. 3. Atomic charges at timestep 0 after the energy equilibration process for ρ = 250 mol/m 3 ; all O atoms form O 2 molecules and all C and H atoms belong to the n-C 12 H 26 molecule.Note the y-axis scaling difference between the results from QEq and QTPIE methods.

Fig. 4 .
Fig. 4. Time evolution of the n-dodecane and oxygen number of molecules.The system's conditions are ρ = 250 mol/m 3 and T = 2100 K.

Fig. 5 .
Fig. 5. Fuel (left) and oxidizer (right) first-order consumption rates under various temperature and density conditions.

Fig. 6 .
Fig. 6.Initiation time for fuel/oxidizer (left) and selected products (right) for ρ = 250 mol/m 3 and T = 2100 K; the dash denotes the average value over the 10 simulations per external electric field case.

Fig. 7 .
Fig. 7. Initiation pyrolysis mechanism of n-dodecane.The number in brackets correspond to the branch label.

Fig. 8 .
Fig. 8. Reaction rate constants of (a) n-dodecane main pyrolysis reactions and (b) key oxygen compound reactions.Solid lines with error-bars are used for MD results; dashed lines for JetSurf2.0[97] and dotted lines for CaltechMech [98] combustion mechanisms.Note that the lines indicating the rate constants of HCO → CO + H and O 2 + H → HO 2 reactions computed with JetSurf2.0and CaltechMech mechanisms are overlapping.

Fig. 9 .
Fig. 9. MEP of CH 4 + O 2 → CH 3 + HO 2 reaction (top) and corresponding atomic charges of oxygen molecule and bonded hydrogen atom (bottom) under different strengths and orientations of external electrostatic fields.The oxygen with label O(2) bonds with the H atom to form HO 2 .The electric field lines of figures (a) and (c) are parallel to the x-direction, whereas of those of figures (b) and (d) are parallel to the y-direction.
2 V/ Å, with the electrostatic field aligned along the Cartesian directions ˆ e x , − ˆ e x , ˆ e y , − ˆ e y , ˆ e z , − ˆ e z for the molecule configuration presented in Fig. 1 a. Results are shown in Figures S13 to S21

Fig. 10 .
Fig. 10.Kinetic energy components under various electric fields at 2100 K and 250 mol/m 3 .Results were initially averaged over 10 simulations per electric field and subsequently with a time interval of 50 ps.

Fig. 11 .
Fig. 11.Potential energy (solid line) and electric potential energy (dotted line) of the system under different electric fields at 2100 K and 250 mol/m 3 .Results were initially averaged over 10 simulations per electric field and subsequently with a time interval of 50 ps.

Fig. 12 .
Fig. 12. Collision frequency under different electric fields at 2100 K and 250 mol/m 3 .Results were initially averaged over 10 simulations per electric field and subsequently with a time interval of 100 ps.

Table 1
Translational, rotational and vibrational kinetic energies difference between the case with an external electrostatic field and the case without any electrostatic field.The values are averaged over the first τa , second τ b and third τc third of the simulated time as indicated in Fig.10.The energies are measured in kcal/mol.