Shear-driven reactions of organosulfur compounds on ferrous surfaces: A molecular dynamics study

Reactive molecular dynamics simulations were used to investigate reactions between di-tert -butyl disulfide and ferrous surfaces. Tribochemical reactions were characterized with and without a model base oil, n -dodecane, and on either Fe(100) or H-passivated Fe 2 O 3 surfaces. Reaction yield increased with both temperature and pressure for all model systems. The presence of the base oil did not significantly affect the reaction yield or pathway, but the base oil impeded transmission of shear stress to the reactants. Replacing the ideal Fe(100) with H-passivated Fe 2 O 3 surfaces enabled reaction pathways involving oxygen from the surface and decreased the reaction yield. The rate-limiting step of the reactions was analyzed in the context of the Bell model as a step towards better understanding tribochemical reactions.


Introduction
Mechanical stress can couple with chemical reactions at the molecular scale to lower the energy barrier that must be overcome for reactions to occur.Such mechanochemical reactions have several potential advantages, including higher yield and improved selectivity [1].In addition to driving chemical reactions by lowering the energy barrier, it has been reported that mechanical forces can introduce new reaction pathways that are not accessible thermally.Alternative pathways due to mechanical work have been reported for synthesis of dibenzophenazine [2], ring-opening of trans and cis isomers of a 1,2-disubstituted benzocyclobutene [3,4], and degradation of polymers [5].Mechanochemistry is therefore highly relevant for chemical synthesis, for example as applied in ball milling, [6,1] where mechanical stress drives chemical reactions that create products with desired properties.
Another application where mechanically driven reactions are ubiquitous is lubrication of moving mechanical components in the field of tribology [7].During operation, tribological systems frequently run in the boundary lubrication regime.Under such conditions, surface-active additives that are added to lubricant formulations, form protective films, called tribofilms, on the contacting surfaces that facilitate relative motion with controlled friction and reduced wear.[8,9].While tribochemical reactions are present in nearly all lubricated (and some non-lubricated) mechanical systems, their fundamental mechanisms are not fully understood because they occur between two moving surfaces where direct experimental measurement is challenging.However, understanding these mechanisms has the potential for significant impact, specifically for optimizing tribofilm growth and more generally for enhancing our understanding of the pathways by which mechanical stress can drive chemical reactions.
Film formation occurs through adsorption and chemical reactions between additives and surfaces.These chemical reactions are believed to be driven by the local heating and mechanical stress that are inherent to tribological contacts.The study of these reactions is referred to as tribochemistry [10].The amount by which reactions leading to tribofilm formation are accelerated by mechanical stress is determined by its effect on the energy barrier for the reaction to proceed.For tribochemical reactions, shear stress due to forces acting parallel to the sliding surfaces is particularly relevant.[10,11] Experimentally, the effect of shear to accelerate chemical reactions can be quantified by measuring the increase in reaction rate or yield with increasing stress.Examples of this experimental approach are studies of bond cleavage [12], polymerization, [13,14] and film formation [15].Specifically relevant to lubricant additives, recent studies have shown that tribofilm growth rates increase with shear stress for zinc dialkyl dithiophosphate (ZDDP) on steel [16], iron [17], silicon [18], and other non-ferrous surfaces [19], as well as for dimethyl disulfide [20], diethyl disulfide, and dimethyl trisulfide on copper [21].Such studies consistently show that shear stress increases the reaction yield or enables reactions to proceed at lower temperatures than observed for purely thermally driven reactions.
Experimental studies of shear-driven reactions have been complemented by simulation-based investigations that provide atomic-scale information about reactions and reaction pathways.First, density functional theory (DFT) calculations have been carried out to identify the pathways of shear-driven reactions, for example, reactions of propanethiolate on gold [22] and methanethiolate on copper [23].DFT was also used to analyze the effect of shear stress to accelerate the removal of alkoxy species from oxidized iron [24] and the effect of normal stress on molecular dissociation of organophosphorus additives on iron [25].
Reactive molecular dynamics (MD) simulations, i.e., MD simulations based on empirical models that capture chemical reactions, have also been used to study shear-driven reactions [7].Pathways of reactions between lubricants and surfaces were investigated for water on SiO 2 [26], hexadecane-lubricated diamond sliding against W(100) [27], di-tert-butyl disulfide on Fe(100) [28,29], and phosphate esters and inorganic alkali polyphosphates on iron oxide [14,30].Reactive MD simulations have quantified the amount by which a reaction energy barrier is reduced by shear for oligomerization during gas phase lubrication [14], film formation from di-tert-butyl disulfide on Fe(100) [28], as well as chemical mechanical polishing of Cu in hydrogen peroxide and glycine [31].Nudged elastic band calculations with a reactive potential have also been used to show how shear stress can lower energy barriers at the molecular scale [32].
While the goals of experimental work and atomistic simulations are often complementary, there are distinct differences between the material in the models and physical systems.First, in practice, film forming additives are typically used in very low concentrations in a lubricant formulation, and most of the liquid present in an interface is base oil.Previous studies have shown that hydrocarbon-based oil can react chemically with various surfaces to create thin films, even without additive.[33,34] Simulations of a model base oil, hexadecane, have shown thin film formation as well.[35,27] However, previous atomistic simulations have not included both base oil and additives to understand their interactions during film formation.
The surface with which additives react also plays an important role in both determining reaction pathways, as well as reaction rate.Most engineering components comprise ferrous materials, which are expected to have a thin oxide layer.However, these oxides can be removed by wear such that the highly reactive metallic iron surface is exposed.[36] It has been observed that the presence of freshly exposed metal surfaces (due to wear) increases the rate of the chemical reactions.[37] This implies that film formation reactions will occur on both iron and iron oxide during operation.For simplicity, previous simulations have used crystalline iron surfaces [28,29,38,39].Others modeled iron oxide, either amorphous or crystalline.[40,41,42].
However, there has been no comparison of mechanically driven film formation reactions on crystalline iron and iron oxide surfaces.
Previous computational studies have demonstrated the reaction mechanisms for di-tert-butyl disulfide on Fe(100), but if or how these mechanisms are relevant to more representative conditions with base oil and a surface oxide is not known.Here, the reactions of di-tert-butyl disulfide confined and sheared between ferrous surfaces were studied using reactive molecular dynamics simulations at a range of temperature and pressure conditions, with or without base oil molecules, either on Fe(100) or H-passivated Fe 2 O 3 .Di-tert-butyl disulfide is a component of sulfurized isobutylene extreme pressure additive whose number of sulfur atoms can vary between 1 and 6. [43,44] The disulfide compound was selected because similar chemistries, i.e., dimethyl disulfide and diethyl disulfide, have been used in experimental studies on the reaction mechanism with iron surfaces.[45,46] The reaction pathways were analyzed, and then the effects of heat and pressure on each step in the pathway were quantified as the number of observed reaction events in the simulation.
The yield of the rate-limiting step of the reaction was then fit to the classic Bell model [47] to quantitatively compare the three cases.The results provide a better understanding of reactions between di-tert-butyldisulfide and ferrous surfaces, as well as the effects of environment and operating conditions on the film formation reactions of extreme pressure additives.

Methods
Three model systems were created, as shown in Fig. 1.The first model system, referred to subsequently as "Fe(100) + additive", consisted of 54 di-tert-butyl disulfide molecules confined between two Fe(100) slabs.
The second model, called "Fe(100) + additive + base oil", incorporated n-dodecane as a model base oil.This model comprised 25 di-tert-butyl disulfide molecules mixed with 25 n-dodecane molecules (base oil) between two Fe(100) slabs.The number of additive and base oil molecules was chosen to achieve 50% concentration while keeping approximately the same separation between the walls and, therefore, approximately the same shear rate for a given wall speed as the models without base oil.While this concentration does not reflect the low concentration of extreme pressure additives in bulk lubricant formulations, it approximates the expected higher density of adsorption of these additives on surfaces.Lastly, in the third model, called "Fe 2 O 3 + additive", we replaced the ideal crystalline Fe surface of the first model with H-passivated Fe 2 O 3 slabs, with 54 di-tert-butyl disulfide molecules in between.Simulations were performed using LAMMPS [48], Python scripts were used for post-processing, and the visualization of the results was carried out with OVITO [49].
The structures of the di-tert-butyl disulfide and the n-dodecane molecules were obtained from Pub-Chem [50] and duplicated using the Packmol [51] package.The Fe(100) slabs were created in Virtual NanoLab [52].For the Fe(100) + additive and Fe(100) + additive + base oil models, the two Fe(100) slabs each had dimensions of 3.4×3.4×1.0 nm (x × y × z) and were initially positioned 2.0 nm and 3.5 nm apart, respectively.The Fe 2 O 3 slabs were created by heating two cyrstalline Fe 2 O 3 slabs to 4000 K over 2.5 ps and then holding at that temperature for 125 ps.The model was then cooled from 4000 K to 300 K over 500 ps.The H-passivation process was carried out by placing 600 H atoms near the surface of each slab.The simulation was run at 700 K for 500 ps (to speed up the process) followed by cooling to 300 K over 250 ps, where the high temperature was used to accelerate the process, as done in previous simulations [53,54].It was confirmed that the system had reached steady-state at the end of this process by verifying the potential energy and number of H atoms bonded to the surface reached plateaus (Fig. S1).Then, any H atom not bonded to the surface was removed from the model.The dimensions of the H-passivated Fe 2 O 3 slabs were 3.6 × 3.6 × 2.5 nm.The Fe 2 O 3 slabs were thicker than the Fe(100) slabs due to the larger unit cell of the former, but the surface areas of both models were similar, ensuring comparable area for reactions during the simulations.The separation between the two slabs after the initialization process was 2.3 nm.For all models, periodic boundary conditions were applied in the directions parallel to the plane of the surfaces (x and y), and the boundaries were fixed in the direction normal to the surfaces (z).The atoms approximately 0.3 nm from the top and bottom of the models were treated as rigid bodies.
Interatomic interactions, including chemical bonding/de-bonding, were modeled using the ReaxFF potential [55] for Fe/S/C/H/O interactions using a parameter set developed by Shin et al. [56], with a time step of 0.25 fs.This force field has been used to model oxidation of butane on Cr 2 O 3 in the presence of FeS 2 [56], reactions between H 2 SO 4 and Fe 3 C(100) at high temperatures, [57] and carburization of iron nanoparticles in ethylene pyrolysis [58].Further, the accuracy of the force field for di-tert-butyl disulfide on Fe(100) was evaluated previously by comparison of adsorption energies calculated by ReaxFF to those obtained from DFT [28].Thermostating of the model was performed by controlling the temperature of the Fe (and O) atoms in the middle layers of the Fe (H-passivated Fe 2 O 3 ) slabs using a Langevin thermostat with a damping time constant of 25 fs.For simulations with relative motion of the slabs, the component of the kinetic energy in the direction of shear was excluded to remove the contribution of the imposed motion from the calculation of temperature used by the thermostat.
The simulations consisted of three stages -heat, load, and sliding -each 2 ns in duration.First, during the heating stage, the models were subjected only to heat at constant temperatures ranging from 500 K to 700 K.These temperatures were chosen both to accelerate the reaction due to time limitations of a reactive molecular dynamics simulation and to represent the flash temperatures expected in sliding contacts [45].
Next, a constant normal load was applied to the rigid part of the top slab while the positions of the atoms in the bottom slab were fixed, leading to contact pressures ranging from 0.50 GPa to 1.50 GPa.This pressure range is consistent with contact pressures in mechanical components with non-conformal contacts.
During the loading stage, the top slab was allowed to move in the z-direction, which led to a reduction of the gap between the slabs.The distance between the slabs after loading was ≈1 nm, approximating the near-contact conditions of boundary or mixed lubrication at which extreme pressure additives operate in lubricated interfaces.Finally, during the sliding stage, the top and the bottom slabs were moved in opposite directions along the x direction at a speed of ±5 m/s.The positions and bond orders of the atoms were recorded every 2.50 ps, and a covalent bond was identified when the bond order exceeded 0.3.The sliding stage was run at all temperatures and all pressures.The total number of temperature/pressure combinations was 30, 15, and 6 for the Fe(100) + additive, Fe(100) + additive + base oil, and Fe 2 O 3 + additive models, respectively.

Reaction Pathways
Iron sulfide films are known to protect surfaces from friction and wear in lubricated mechanical systems [59].The initiation of these reactions is the reaction between sulfur (S) in additives and surface iron (Fe).In a previous simulation-based investigation of the reaction between di-tert-butyl disulfide and Fe(100), [29] it was shown that this reaction occurs through three steps: S-S bond cleavage, followed by Fe-S bond formation, and lastly dissociation of the S-C bond.This pathway, here referred to as path a and shown in Fig. 2a, was observed in all three model systems studied here.Analysis of the likelihood of chemisorption showed that more tert-butyl radicals chemisorbed at higher pressures and temperatures.Importantly, although the steps of the three pathways differ slightly from one another, they all lead to sulfur bonded to the surface, which can be assumed to be the onset of formation of iron sulfide films that have been observed in tribological experiments.[59]

Reaction yield
The number of bonding/dissociation events for each step in these reactions increased during all three stages (heating, loading, sliding) of the simulation, but the most rapid increase was observed during the heat and sliding stages (Fig. S2).This suggests that the reactions were activated both thermally and mechanically.
However, to compare different temperature and pressure conditions, reaction yield was quantified as the number of bonding/dissociation events calculated at the end of the sliding stage relative to the total number of events possible given the number of di-tert-butyl disulfide molecules in each model system.This analysis was performed separately for each of the three main steps in the reaction: S-S cleavage, S-Fe/S-O formation, and S-C dissociation.
The S-S dissociation was always the first step in the reaction pathway.The yield for this reaction step from all simulations is summarized in Fig. 3a, b, and c, where yield is represented by color (dark corresponding to the highest yield) as a function of temperature on the abscissa and pressure on the ordinate axis.The heat map color indicates the model, where red corresponds to the Fe(100) + additive model, green to the Fe(100) + additive + base oil model, and blue to the Fe 2 O 3 + additive system.
The heat maps show that this step of the reaction depends on temperature and pressure differently for the three model systems.For the Fe(100) + additive case, in Fig. 3a, S-S dissociation increases with both temperature and pressure.As shown in Figs.3d and e, the effect of temperature is more significant at the lower pressures, and the effect of pressure is more significant at the lower temperatures.
In contrast, for the base oil containing model, the yield at any pressure or temperature in Figs.3b, e, and h is above 90%.This indicates that nearly all of the possible bonds were broken by the end of the simulation under any condition, so temperature or pressure dependence could not be observed.
Lastly, for the H-passivated Fe 2 O 3 model, the effect of temperature is more significant than the effect of pressure.As shown in Figs.3c, f, and i, the yield of the S-S dissociation step increases rapidly with the increase in temperature from 500 to 700K, at any pressure.However, the effect of pressure is not significant at any temperature.Fig. 4 shows the reaction yield for the Fe-S or O-S bond formation.The format of the sub-figures is the same as Fig. 3a, d, and g.For the Fe(100) + additive case, in Fig. 4, yield increases with pressure at all temperatures, with the effect of pressure being larger at lower temperatures.The thermal effect depends on pressure.Specifically, at low pressure, yield increases with increasing temperature, but the opposite is observed at high pressure.It has been observed from the simulations that some of the Fe-S bonds that form are subsequently broken by shear force at high pressure, and this process is facilitated by high temperature.
The decrease is attributable to this bond breakage.In the Fe(100) + additive + base oil case, shown in Figs.4b, e, and h, like the S-S dissociation, nearly all possible Fe-S bonds are formed at any pressure or temperature.The maximum yield is reached for all cases except the lowest temperature and pressure conditions.
For the H-passivated Fe 2 O 3 model, shown in Figs.4c, f, and i, the reaction yield does not reach 100%.
There is generally more yield at the higher pressures and temperatures, but yield does not increase significantly with either pressure or temperature.
The final step we investigated was the dissociation of the S-C bond.This was identified as the ratelimiting step of the reaction in our previous study of heating-only simulations, where S-C dissociation commenced at the highest temperature of the three steps in the reaction [28].Our previous studies have also shown that this step of the reaction is accelerated by temperature as well as pressure [29].The results for the three models here are shown in Fig. 5.
For the Fe(100) + additive case (Figs.5a, d, and g), like in the previous two reaction steps discussed, yield increases with both pressure and temperature.However, for the S-C dissociation, the increase in yield with temperature is more significant than pressure.
For the model with base oil, shown in Figs.5b, e, and h, this is the only reaction step that did not reach the maximum yield.Also, increasing the temperature increased yield at all pressures (Fig. 5e), although the reaction step did not appear to be accelerated by pressure (Fig. 5h).
Lastly, for the H-passivated Fe 2 O 3 case, despite the fact that the yield increased slightly with temperature, very few reactions were observed under any set of conditions.The effects of base oil and H-passivated oxide surface on the reactions are discussed next.

Effect of base oil
Previous studies of film forming additives and base oil have shown that a higher concentration of additive relative to the base oil results in faster film formation [60,61,62,63], because base oil that is chemically or physically adsorbed on the surface can limit access of the additives to the surface.[60] This suggests that the simulations with base oil present should have lower yield than those without base oil.However, the analyses in the previous section showed there were only minor differences between the yield with and without base oil.
The first two steps of the reaction achieved nearly 100% yield in both cases (Figs. 3 and 4).The third, rate-limiting step did not reach 100%, but the yield was comparable with and without the base oil (Fig 5).
One reason that lower yield was not observed in the base oil case is that there were fewer di-tert-butyl disulfide molecules in the simulation with base oil; this was done to ensure that the model sizes would be approximately the same.Since yield is reported as the ratio of broken bonds to the total number of available bonds, the smaller number of additive molecules in the model with base oil resulted in a higher yield for the same number of broken bonds.Another difference between these simulations and experiments is the concentration of the additive.In experimental studies, the concentration of film forming additives is very low, typically less than 5 wt.%.[60] Due to the size-scale limitations of a MD simulation, an approximately 50 wt.%ratio was used here.As a result, there was relatively little base oil present in the model, and it had a negligible effect on yield.
Although the yield of the S-C dissociation step was similar for the cases with and without base oil, some difference in the effect of pressure was observed.For the model without base oil, there is an increase in yield with pressure in Figs.5d and 5g.This pressure-dependence was confirmed by running two more simulations at a higher pressure (see Fig. S1).However, when base oil is present, no increase in yield with pressure is observed in Figs.5e and 5h.Previous DFT calculations [23] and reactive MD simulations [29] have shown that the S-C dissociation reaction occurs through lateral displacement of the C atom with respect to the S atom.Shear stress therefore drives this step of the reaction by pushing the C atom along the reaction coordinate.In the simulations, the dodecane molecules were weakly bonded with the Fe(100) surface, but remained near the surface (≈ 0.2 nm) throughout the simulations.Additionally, visualizations of the simulations showed that there were often dodecane molecules in the vicinity of the S-C bonds before dissociation.This suggests that the presence of the base oil impedes the lateral movement of the tert-butyl group, thereby minimizing the effect of mechanical stress to drive the reaction step.
Lastly, our simulations showed that the dodecane molecules became aligned with the direction of sliding during the shear stage (see Fig. S2).This phenomenon has been observed previously for dodecane molecules confined and sheared between mica surfaces, and the alignment was shown to lower friction.[64] Therefore, we infer that the aligned dodecane molecules reduced the shear force transferred from the walls to the reactants and, consequently, decreased the pressure dependence of the reaction yield.

Effect of H-passivated oxide surface
Next, we analyzed the effect of the H-passivated oxide surface on film formation.For all three steps of the reaction, the yield was lower for the model with H-passivated Fe 2 O 3 than either of the models with the Fe(100) surface.This difference can be attributed to several factors.
First, in the simulations with H-passivated Fe 2 O 3 , it was observed that the O and H atoms were released from the surface in the form of water, O 2 , OH due to thermal as well as mechanical effects.A representative snapshot of the model with O-and H-containing species released from the surfaces is shown in Fig. S5.These elemental H and O, as well as OH groups, are highly reactive and bonded with tert-butyl thiyl radicals (C 4 H 9 S), the main decomposition product, to form various oxide species (see Fig. 2).We also observed oxidation of di-tert-butyl disulfide molecules to derivative oxide species.Experimental studies have shown that when ZDDP is oxidized, the oxidation reaction products were ineffective as anti-wear agents.[65] This suggests that oxidation inhibited the ability of the additive to form protective films on surfaces, consistent with the observations in our simulations that the oxides formed from di-tert-butyl disulfide and its moieties were less likely to proceed through any of the reaction pathways.
Second, the chemical stability of the H-passivated Fe 2 O 3 limits reactions with the surface.As seen in Fig. 4, the yield of the chemisorption step of the reaction for the H-passivated Fe 2 O 3 is ≈50%, whereas it reaches nearly 100% for the other two cases.Experimental studies of reactions between disulfides and ferrous surfaces have shown that reaction rates are lower on oxide surfaces.[66] The H-passivated Fe 2 O 3 surface was passivated with H in the simulations, prior to introducing the additive molecules, making it even more chemically stable.[67] Lastly, Fig. 5 shows that the reaction yield of the rate-limiting step was much lower for Fe 2 O 3 + additive than for the other two models.The models with Fe(100) proceed through pathway a, where chemisorption occurs via Fe-S bonding, whereas, in the H-passivated oxide model, S could bond with either Fe or O on the surface.Previous studies of a similar chemical system showed that more S-O bonds (HSO − CH 3 compared to HSO 3 − CH 3 ) increased the S-C bond dissociation energy.[68,69] Therefore, in our simulations, the S-O bonding to the tert-butyl sulfide moieties may have increased in S-C dissociation energy, leading to lower reaction yield.

Bell model
The results above are presented in terms of pressure, since load is the controllable/measurable parameter in an experiment (or a simulation mimicking an experiment).However, it is known that tribochemical reactions are driven by shear stress rather than by pressure alone [8,70].Therefore, the change in yield with pressure observed here can be attributed to an increase in shear stress.The shear stresses in the simulations were calculated as a time average of the lateral force on the rigid part of the bottom and top Fe(100) (or H-passivated Fe 2 O 3 ) slabs divided by the surface area (area of the model system in the xy plane).
The yield of a reaction is exponentially related to the height of the energy barrier that must be overcome for that reaction to proceed.The height of the barrier can be lowered by stress, leading to the classic Bell model [71]: where r y is the reaction yield, A is a pre-exponential factor, k B is the Boltzmann constant, T is the temperature, E 0 is the energy barrier in the absence of stress (thermal activation energy), ∆V * is the activation volume, and τ is the shear stress.By taking the natural logarithm of Eq. ( 1), a linear relationship between the natural logarithm of the reaction yield and (k B T ) −1 is obtained, where the slope corresponds to the reaction energy barrier: This framework is generally applicable to reactions driven by mechanical stress; for tribochemical reactions, it is used to quantify the effect of shear stress [11].
The reaction yield for the rate-limiting step of the reaction for each model system in Fig. 5 can be fit to Eq. ( 2).The multi-parameter fit was optimized by minimizing the mean of the squared differences between ln(r y ) as calculated from the Bell model and ln(r y ) as measured from our simulations (see Fig. S3).
Fit parameters were ln(A), E 0 , and ∆V * .This fitting approach was previously used for shear-driven reactions.[72] The first fit parameter, E 0 , was found to be 11.8 ± 1.5, 12.6 ± 2.8, and 10.2 ± 4.4 kcal/mol for the Fe(100) + additive, Fe(100) + additive + base oil, and Fe 2 O 3 + additive models, respectively.It has been reported that, if the reaction rate or yield is used in Eq. ( 2), as done here, the magnitudes of fit values for E 0 cannot be directly correlated to thermal activation energy [11].Further, it has been shown that the availability of reaction sites is correlated to barrier height [73], so the finite size of the models here is expected to affect the fit value of E 0 .However, for a series of systems that undergo similar chemical reactions, the relative magnitudes of E 0 obtained from linear fits are meaningful [11].Therefore, the fact that the differences between the E 0 values for the three cases are less than the fit error indicates their thermal activation energies are similar.
More differentiation between the three cases was observed in the magnitude of the activation volume, which was found to be ∆V * = 2.10 ± 0.65, 0.21 ± 0.59, and 0.75 ± 0.54 Å3 for Fe(100) + additive, Fe(100) + additive + base oil, and Fe 2 O 3 + additive, respectively.Although the physical meaning of activation volume has been debated, it is a measure of how susceptible a given reaction is to mechanical activation.[11] Therefore, the larger activation volume for the Fe(100) + additive case, with the comparable activation energies for all three cases, indicates that the S-C dissociation is most readily driven by shear for Fe(100) + additive.
As mentioned above, S-C dissociation occurs through lateral displacement of the C atom with respect to the S atom [23,29], a process that is facilitated by lateral force.Shear can be transmitted from the walls to the molecules either directly, as is the case in pathways a and b when the tert-butyl sulfide is bonded to the surface, or indirectly by the strain field within the liquid, as would occur in pathway c.The latter case was demonstrated in a study of ZDDP film formation where the experiment was performed in full film lubrication, i.e., no surface-surface contact.[8].The S-C dissociation is most readily activated by shear on the Fe(100) surface and without base oil because the direction of shear is perfectly aligned with the reaction pathway on the atomically smooth surface (unlike the H-passivated Fe 2 O 3 surface).base oil is not present in the Fe(100) + additive model to impede the transmission of shear force from wall to reactant.
Thus, compared to the Fe(100) + additive case, the activation volume is smaller with H-passivated Fe 2 O 3 due to the lack of alignment of shear force with the reaction pathway and is smaller with base oil due to the interference of the base oil molecules in the transmission of shear force to the reactants.

Conclusions
In summary, the reaction pathways of di-tert-butyl disulfide on ferrous surfaces were investigated using reactive MD simulations run across a range of pressure and temperature conditions.The simulations were designed to evaluate whether mechanisms identified previously for the additive on a crystalline iron surface are relevant to more representative conditions with base oil or a surface oxide.
Comparing the Fe(100) and Fe 2 O 3 surfaces, results showed that the reaction yield was lower and proceeded through different pathways on the oxide due to the availability of O and H in the system.The lower yield on the Fe 2 O 3 was attributed to three factors: the chemical stability of the H-passivated oxide surface, additive oxidation reactions that hindered the additive decomposition and chemisorption steps, and increased S-C dissociation energy that lowered the yield of the rate-limiting step of the reaction.
The effect of dodecane base oil was less significant.The reaction pathway was the same with or without the base oil.However, the pressure-dependence of the yield was larger without the base oil.This was proposed to be due to alignment of the dodecane molecules on the surface, which lowered the shear stress transferred to the reactants at a given pressure.This was corroborated by the fit value of activation volume, which was smallest for the base oil simulations, indicating that the ability of shear stress to drive these reactions was lessened by the presence of base oil.
Generally, the findings of this research provide a better understanding of how thermal and mechanical factors in a sliding interface drive film formation reactions.Further, the study demonstrates an approach to predicting the pressure and temperature dependence of reactions between additives and surfaces, and it can be applied to other additives as well as other relevant surfaces.

Figure 1 :
Figure 1: Side-view snapshots of the model systems during the sliding stage of the simulations for (a) Fe(100) + additive, (b) Fe(100) + additive + base oil, and (c) Fe 2 O 3 + additive.(d) Chemical structures of n-dodecane and di-tert-butyl disulfide.In all figures, the surface Fe and O atoms are shown in brown and blue, respectively.Di-tert-butyl disulfide atoms are shown in yellow (S), dark gray (C), and white (H).Both H and C atoms in the n-dodecane are shown in light gray to distinguish the base oil from the additive.

Figure 2 :
Figure 2: Three distinct reaction pathways were observed in the simulations.In all three pathways, the reaction involved S-S cleavage, S-C dissociation, and S chemisorption to the surface.The detailed steps in each pathway are described in the text.

Figure 3 :
Figure 3: Heat maps (a, b, and c) of reaction yield for S-S bond cleavage (number of bonds broken at the end of the shear stage of the simulation) as a function of temperature and pressure.An increase in yield is shown by a color change from light (lowest) to dark (highest).The maps are created from 30, 15, and 6 data points for the red, green, and blue plots, respectively.Yield for models Fe(100) + additive, Fe(100) + additive + base oil, and Fe 2 O 3 + additive are shown in red, green, and blue, respectively.2D representations show the change in yield as a function of temperature (d, e, and f) at two pressures (0.5 and 1.5 GPa) and as a function of pressure (g, h, and i) at two temperatures (500 and 700 K).

Figure 4 :
Figure 4: Heat maps of reaction yield for Fe-S (a and b) and Fe-S/O-S bond formation (c) as a function of temperature and pressure, and corresponding 2D plots at representative pressures and temperatures.The figure format and number of data points are the same as in Fig. 3.The heat map for the H-passivated Fe 2 O 3 model had a smaller yield range than the plots of the other models.

Figure 5 :
Figure 5: Heat maps of reaction yield for S-C dissociation as a function of temperature and pressure, and corresponding 2D plots at representative pressures and temperatures.Figure format and number of data points are the same as in Fig. 3.The 2D plots have the same range while the heat maps have different color ranges due to variation of the yield.