Nonequilibrium molecular dynamics simulations of stearic acid adsorbed on iron surfaces with nanoscale roughness

Nonequilibrium molecular dynamics (NEMD) simulations have been used to examine the structure and friction of stearic acid films adsorbed on iron surfaces with nanoscale roughness. The effect of pressure, stearic acid coverage, and level of surface roughness were investigated. The direct contact of asperities was prevented under all of the conditions simulated due to strong adsorption, which prevented squeeze-out. An increased coverage generally resulted in lower lateral (friction) forces due to reductions in both the friction coefficient and Derjaguin offset. Rougher surfaces led to more liquidlike, disordered films; however, the friction coefficient and Derjaguin offset were only slightly increased. This suggests that stearic acid films are almost as effective on contact surfaces with nanoscale roughness as those which are atomically-smooth.


Introduction
The requirement to reduce the energy consumption and thus CO 2 emissions from engineering systems has resulted in a general lowering of lubricant viscosity in order to minimise hydrodynamic friction losses. However, this means that an increasing number of engineering components operate under boundary lubrication conditions, where solid asperities come into direct contact, leading to high friction and wear. As a result, lubricant additives which can reduce friction and wear under boundary conditions, such as organic friction modifiers (OFMs), are of increasing importance [1].
OFMs are amphiphilic molecules that contain both a polar head group and a nonpolar hydrocarbon tail group. The most widely studied OFM molecule in both experiments [2][3][4][5][6] and MD simulations [7][8][9][10][11][12][13][14][15] is stearic acid; a carboxylic acid with a saturated, linear C 18 alkyl chain. The acid head group adsorbs to metal or ceramic surfaces and strong, cumulative van der Waals forces between proximal nonpolar tails leads to the formation of incompressible monolayers that prevent contact between solid surfaces and reduce adhesion and friction [1,8].
In order to fully understand the performance of OFM additives, it is necessary to obtain a detailed picture of: (i) the nanoscale structure of their films and (ii) their tribological behavior [7,8]. Detailed structural information on OFM films can be obtained from experimental techniques such as sum frequency spectroscopy (SFS) [16], polarized neutron reflectometry (PNR) [2], the surface force apparatus (SFA) [3], and in situ atomic force microscopy (AFM) [4]. The tribological behavior of OFM films can also be investigated using SFA [3] and AFM [4] as well as dedicated boundary friction experiments [5,6] which employ low sliding velocities (mm s −1 ) and high pressures (GPa) in order to maintain boundary lubrication conditions [8].
Classical nonequilibrium molecular dynamics (NEMD) simulations can be used to simultaneously probe the nanoscale structure and friction of OFM films, making it a valuable complement to experiments [8]. The use of accurate all-atom force fields enables the structure and friction of large molecular systems to be reliably analysed over time [9].
The effect of surface roughness on the performance of OFMs is not fully understood [1]. Most previous NEMD simulations of OFM films on metal surfaces have utilised atomically-smooth slabs, with only atomic corrugation [7][8][9][10][11][12], although the influence of artificial asperities placed on top of atomically smooth slabs have also been considered [13][14][15]. The presence of nanoscale roughness has been shown to significantly influence the adhesion of solid surfaces through both MD simulations and AFM experiments [17]. The effect of 3D nanoscale fractal roughness on the friction between sliding surfaces has also been investigated for both dry sliding [18][19][20] and in the presence of lubricant molecules [21][22][23] in large-scale NEMD simulations. These previous simulations suggested that the friction coefficient depends on how effectively the lubricant is able to prevent the direct contact of asperities on opposing surfaces [22]. Nonpolar lubricant molecules can be easily squeezed-out from between asperities [24,25], but strongly-adsorbed OFMs are expected to be more difficult to remove, making them more effective in reducing asperity contact and thus friction in the boundary lubrication regime [8]. Other NEMD simulations of alkanes between amorphous and crystalline surfaces showed that the layering of lubricant films was strongly suppressed on rougher amorphous surfaces, resulting in a more liquidlike film [26]. This may have significant implications regarding the ability of OFM films to separate opposing solid surfaces with nanoscale roughness.
The kinetic friction coefficient in nanotribological systems can be obtained using the extended Amontons-Coulomb law: F L =μ F N +F 0 , where F L and F N are, respectively, the average total lateral and normal forces acting on the outermost layer of atoms in each of the slabs and F 0 is the Derjaguin offset, the friction force at vanishing load, which is generally attributed to adhesive forces [7,15,27]. Recent SFA experiments of palmitic acid OFMs in PAO base oil highlighted F 0 as a significant contribution to F L under relatively low applied pressures ( < 10 MPa) [27]. Previous NEMD simulations showed that the Derjaguin offset increased with increasing surface roughness, decreased with increasing stearic acid coverage, and became negligible when sufficient lubricant was present to separate the stearic acid films [15]. In our previous NEMD simulations, we demonstrated that the friction coefficient is reduced at high coverage of stearic acid, due to the formation of solidlike monolayers which allow very little interdigitation and yield slip planes between well-defined molecular layers [8]. In this study, we will probe the effect of; (i) higher pressures (0.5-2.0 GPa), in order to test the limits of operability of OFMs, and (ii) surfaces with different levels of 3D root mean square (RMS) roughness (0.2, 0.5, 0.8 nm), in order to assess the effect of nanoscale roughness on the structure and tribological behavior of OFM films. The roughness added in this study is expected to be similar to that found on experimental surfaces, since it is added using a random midpoint displacement (RMD) algorithm, which generates surfaces with a quantifiable RMS roughness [18]. Conversely, in previous confined NEMD simulations of stearic acid, the roughness was varied by changing the size of artificial roughness features [13][14][15].
The 'smooth' steel surfaces used in tribological experiments generally have a RMS roughness of approximately 10 nm [28,29], though other surfaces, such as gold, can be produced with a RMS roughness well below 1 nm [30]. Therefore, all of the RMS roughness values used in this study represent an extremely smooth experimental surface; however, accurate representation of larger RMS roughness values would require prohibitively large simulation areas [26]. Moreover, since the roughness features considered here are on the atomic scale, they have the most potential to weaken OFM films by disrupting intermolecular forces between neighbouring molecules. The current simulations will highlight the effect of changes in the nanoscale structure within the film due to changes in surface roughness on both the friction coefficient and the Derjaguin offset. The results obtained with stearic acid will also be directly compare to results for nhexadecane under similar conditions [22], in order to highlight the reasons for the differences in their friction reduction performance in the boundary lubrication regime [5,6].

Simulation setup
A representative example of the systems simulated in this study is shown in Fig. 1a. It consists of two stearic acid monolayers adsorbed on α-iron slabs with 3D nanoscale RMS roughness (Fig. 1b). The α-iron slabs are chosen as a model for steel, which is of significant academic and industrial interest. Though α-iron oxide would be a more accurate representation of a steel surface [8], no classical MD force-field is currently available which can accurately model its deformation under the high pressures applied in these simulations. Stearic acid was chosen as a model OFM, an important class of boundary lubricant additive, which has been used in numerous previous experiments [2][3][4][5][6] and MD simulations [7][8][9][10][11][12][13][14][15]. Our previous NEMD study showed that whilst there were significant variations in the structure and friction of stearic acid films and those formed by other types of OFM (amides and glycerides), all gave the same general trends [8]. Unlike in our previous study [8], no lubricant molecules were added between the stearic acid films. Preliminary squeeze-out simulations showed that only two molecular layers of n-hexadecane remained between stearic acid films at 0.5 GPa [8], and since higher pressures (0.5-2.0 GPa) were applied in the current simulations, very few n-hexadecane molecules were expected to remain between the asperities. A more complete understanding of lubricant squeeze-out in the presence of OFM films would certainly be an interesting target of future MD simulations and AFM experiments, but is beyond the scope of this current study. All structures were constructed using the Materials and Processes Simulations (MAPS) platform from Scienomics SARL.
Most surfaces have roughness on several length scales, including the nanoscale, that can be described by a self-affine fractal scaling law [18,32]. Here, the Hurst exponent and the root mean square (RMS) roughness can be used to quantify the amount of roughness [19]. Using the random midpoint displacement (RMD) algorithm, rough surfaces can be generated which are periodic across their boundaries [19]. The RMD algorithm, with a Hurst exponent of 0.8, was used to independently generate the same RMS nanoscale roughness (0.2, 0.5 or 0.8 nm) in the top and bottom slabs [18,22]. In order to avoid generating only a few large asperities, the RMD algorithm did not start from the centre of only one square (entire slab), but rather four smaller, equally sized squares [19]. The slabs themselves had approximate x, y, z dimensions of 11, 11, 5 nm. Periodic boundary conditions were applied in the x and y directions. The dimensions of the solid surfaces confining the stearic acid films are much larger in the current simulations than those used in many previous NEMD studies [7][8][9][13][14][15], in order to provide a more faithful representation of the statistical distribution of the heterogeneous surface morphology [26]. On experimental surfaces, the typical height of a roughness feature is expected to be approximately 2-3 orders of magnitude smaller than its lateral dimensions [33]. To replicate this in NEMD simulations would require a prohibitively large system size and hence the 'steepness' of the roughness features in these and previous NEMD simulations [18,33] is expected to be somewhat exaggerated. Nonetheless, the use of an RMD algorithm provides a more realistic representation of nanoscale surface roughness than harmonic [34] or other artificially introduced [15] roughness features.
Stearic acid molecules were oriented perpendicular to, and initially 3 Å from, the interior surfaces of the two slabs (Fig. 1a). The surface coverage, Γ, can be defined as the average number of stearic acid molecules present in a given surface area (nm −2 ). Three coverages of stearic acid were considered; a high surface coverage (Γ=4.56 nm −2 ) close to the maximum theoretical value [35]; a medium coverage (Γ=3.04 nm −2 ) approximately 2/3 of the maximum coverage; and a low coverage (Γ=1.52 nm −2 ) around 1/3 of the maximum coverage (Fig. 1b). Note that the Γ values assume an atomically smooth surface. The high, medium and low coverages correspond to 600, 400, and 200 stearic acid molecules adsorbed on each of the 131.6 nm 2 slabs respectively. The highest coverage simulated has also been observed experimentally [35]. Simulations with no stearic acid molecules between the slabs were also conducted for comparison.
Classical MD simulations were performed using LAMMPS [36]. The MD equations of motion were integrated using the velocity-Verlet algorithm with an integration time-step of 1.0 fs. Fast-moving bonds involving hydrogen atoms were constrained with the SHAKE algorithm [37]. The stearic acid molecules were represented by the L-OPLS-AA force-field [38]. This is an updated form of the OPLS-AA force-field [39] which was explicitly parameterised for long-chain molecules. We have shown previously that the use of accurate all-atom force-fields is critical to obtain accurate tribological behavior in confined systems which include long-chain molecules [9]. Lennard-Jones interactions were cut-off at 10 Å and 'unlike' interactions were evaluated using the geometric mean mixing rules [39]. Electrostatic interactions were evaluated using a slab implementation of the particle-particle-particle-mesh (PPPM) algorithm [40] with a relative accuracy in the forces of 1×10 −5 . Iron-stearic acid interactions were represented by the Lennard-Jones potential; the parameters for iron were developed by Savio et al. [34] for the adsorption of alkanes. The αiron slabs contain no partial charges so there is no electrostatic interaction between the stearic acid molecules and the α-iron slab; however, there is still preferential adsorption of the head groups due to the stronger iron-oxygen interactions (ε≈0.06 eV) relative to the ironcarbon interactions (ε≈0.01 eV) [15]. Iron-iron interactions within the slabs are represented by the Embedded Atom Model (EAM) [41], which can accurately model deformation of the slabs under the high pressures applied [22]. Iron and steel surfaces quickly become oxidised when exposed to air, which significantly reduces the adhesion force between contact surfaces in experimental systems [22]. Therefore, a Lennard-Jones potential was used for iron-iron interactions between atoms in opposing slabs ( Fig. 1a) in order to mimic the reduced adhesion between oxidised surfaces [18]. The Lennard-Jones parameters, ε=0.02045 eV and σ=0.321 nm, were used for this interaction, as have been successfully utilised in previous NEMD simulations of tribological systems [22,42].

Simulation procedure
The systems were initially energy minimized and then compressed (F N =0.5, 1.0, 1.5, 2.0 GPa), thermostatted (300 K) in the directions perpendicular to the compression (x and y), and allowed to equilibrate. During the compression stage, the temperature is controlled using a global Langevin thermostat [43], with a time relaxation constant of 0.1 ps. The pressure was controlled by applying a constant normal force to the outermost layer of atoms in the upper slab, keeping the zcoordinates of the outermost layer of atoms in the lower slab fixed (Fig. 1a), as is common in confined NEMD simulations [13][14][15][21][22][23]. The slab separation initially varied in a damped harmonic manner, so sliding was not applied until a constant average slab separation was obtained and the hydrostatic pressure within the stearic acid film was close to its target value. The compression simulations were approximately 0.5 ns in duration.
A velocity of +5 m s −1 was then added in the x direction to the outermost layer of atoms in the top slab and −5 m s −1 to the bottom slab (total sliding velocity=10 m s −1 ), as shown in Fig. 1a, and sliding simulations were conducted for 2.0 ns. All simulations were run for long enough to yield a sufficient sliding distance (20 nm) to obtain representative values for the lateral (friction) forces (uncertainty < 10%). While lower sliding velocities are desirable to match those used in boundary friction experiments (typically mm s −1 ), they are not yet accessible using NEMD simulations of this scale [7,8].
During the sliding simulations, any heat generated by shearing of the molecules was dissipated using a Langevin thermostat [43], with a time relaxation constant of 0.1 ps, which acted only on the 10 Å of αiron slabs closest to the fixed layers (Fig. 1a), and was applied in the direction perpendicular to both the sliding and compression (y). This is known to be advantageous over direct thermostatting of the fluid, which has been shown to significantly alter the behavior of confined fluids under sliding conditions [44]. This boundary thermostatting method is common in confined NEMD simulations [13][14][15][21][22][23] and has been shown to be effective in controlling the temperature of similar systems at the sliding velocity applied applied here (10 m s −1 ) [45]. At the onset of sliding, an expansion due to the increase of temperature by shear heating was expected, so it was ensured that steady state sliding [7] had been attained before sampling of the lateral and normal forces began (See Appendix A). The time taken to achieve steady state sliding equated to approximately 0.5 ns (5 nm of sliding).

Results and discussion
Variations in the nanoscale structure within the films with nanoscale RMS roughness and stearic acid coverage were monitored through visualised trajectories, atomic mass density profiles, velocity profiles, and radial distribution functions. Thermal and mechanical equilibration were confirmed before sampling began (See Appendix A). The tribological behavior of the stearic acid films were probed by examining the change in lateral forces on the outer layer of atoms in the slabs, F L , with normal force, F N . The forces are presented as pressures (force divided by the contact area) in order to aid experimental comparisons. A linear plot of F L against F N was then used to establish the Derjaguin offset, F 0 , and the friction coefficient, µ, for each of the surface roughness and coverages combinations simulated. The variation in F 0 and µ with RMS surface roughness and coverage was then plotted in order to establish the relative influence of each. 3.1. Structure   Fig. 2 shows the variation in the interface between the stearic acid films on the top and bottom slabs with surface coverage and RMS roughness. Even with the exaggerated steepness of the nanoscale roughness features [33], and high pressures applied in these simulations, the stearic acid films maintain separation of the asperities at all coverages, which demonstrates the robustness of the films.
These results are in contrast to those obtained for similar coverages of n-hexadecane, which allowed direct contact of asperities under the same conditions [22]. At increasing stearic acid coverage, there is a larger separation between the opposing slabs. Moreover, the terminal carbon atoms, highlighted in yellow, form more ordered layers, resulting in a smoother interface between the films adsorbed on the upper and lower slabs. Increased RMS surface roughness seems to supress this layering somewhat and leads to more disordered interfaces; however, these changes appear to be less significant than the dependence on surface coverage. Fig. 3 shows the change in the atomic mass density profiles for all atoms in the stearic acid molecules in z, ρ(z), and x-velocity profile in z, v x (z), with coverage at 0.2, 0.5 and 0.8 nm RMS roughness. The velocity profiles at all levels of coverage and roughness show that there is no slip at the surface, as expected for the strongly absorbed carboxylic acid head groups [7,8]. The tail groups also move at equal velocity to the slabs to which they are absorbed ( ± 5 m s −1 ), but only up to the region where they become interdigitated with the other film, at which point the atoms move slower as the interface is sheared [8]. In common with results from atomically smooth slabs [8], there is a decrease in the interdigitation between films adsorbed on opposing slabs as the coverage of stearic acid is increased, as indicated by reduced overlap of the mass density profiles. Therefore, the interface is more ordered in the 4.56 nm −2 coverage systems, resulting in steeper velocity profiles compared to those at lower coverage. This suggests the presence of slip planes between the films adsorbed to the top and bottom slabs, which leads to less shearing of the film and thus a lower viscous contribution to the friction coefficient [8]. On slabs with greater RMS surface roughness, there is more interdigitation between films adsorbed on opposing slabs; although this effect is less significant than the change with coverage. The velocity profiles are similar at all levels of RMS surface roughness, though they are generally steeper at 0.2 nm, suggesting a clearer slip plane due to a smoother interface. At 0.2 nm RMS roughness, the mass density profiles are very similar to those observed on atomically smooth α-iron oxide surfaces [8]; however, as the surface roughness increases, the outermost head group peaks (closest to the iron surface) become considerably less intense. This is similar to what has been highlighted previously for linear alkanes, where less ordered, more liquidlike films were observed on surfaces with greater nanoscale roughness [26]. This is explored further for stearic acid through the film radial distribution function (RDF) in Fig. 4. Fig. 4 shows the change in the radial distribution function (RDF) for the carbonyl carbon and terminal carbon atoms with coverage at 0.2, 0.5 and 0.8 nm RMS roughness. In Fig. 4, the carbonyl carbon shows long-range order for all levels of RMS roughness and coverages, with the major peaks occurring at multiples of r=4 Å. This slightly exceeds the unit-cell dimension of the α-iron surface (2.86 Å) [41], suggesting that, while the surface structure determines the head group packing of the stearic acid films on α-iron oxide [8], on α-iron, the limiting head group area [35] is too large to occupy every lattice site. The carbonyl peak at 4 Å is of similar intensity at all coverages; however, at 4.56 nm −2 , there are sharper, more intense peaks at 8 and 12 Å, indicating long-range ordering of the head groups and more solidlike films. The terminal peaks at 4, 12 and 16 Å are most intense at 4.56 nm −2 coverage, suggesting that the tail groups remain solidlike at the interface [8]. The change in the RDFs with coverage is similar at 0.2 nm RMS roughness as on atomically smooth α-iron oxide surfaces [8]; however, there is much less variation on the RDFs with coverage at 0.8 nm RMS roughness. Films on surfaces with larger RMS roughness generally have weaker carbonyl and terminal peaks, particularly at 8 and 12 Å, indicating less long-range order and more liquidlike films.  atoms in the slabs (Fig. 1a) are both time-averaged over the data acquisition period (1.5 ns). This approach is commonplace in confined NEMD simulations of tribological systems [7][8][9][13][14][15] and is similar how these forces are measured experimentally [5,6]. All of the combinations in Fig. 5 yield a linear increase in F L with F N and a non-zero intercept, indicating that all of the systems follow the extended Amontons-Coulomb friction law [15]. A steeper gradient in Fig. 5 indicates a greater friction coefficient while a larger intercept shows a larger Derjaguin offset, F 0 , which represents load-independent adhesive forces [7]. Generally, as the surface coverage of stearic acid increases, there is a reduction in lateral force through a reduction in both the gradient (µ) and intercept (F 0 ). Conversely, as the RMS surface roughness increases, there is an increase in lateral forces due to an increase in both gradient (µ) and intercept (F 0 ). The variation in the F 0 and µ values (determined from Fig. 5) with surface coverage and RMS roughness are presented in Figs. 6 and 7 respectively. Fig. 6a shows that the variation of the friction coefficient, µ, with  coverage depends on the RMS surface roughness. At 0.2 nm RMS roughness, the friction coefficient increases by 3% between 1.52 nm −2 and 3.04 nm −2 coverage and then decreases by 12% between 3.04 nm −2 and 4.56 nm −2 coverage. This is the same pattern observed in NEMD simulations of stearic acid adsorbed on atomically smooth α-iron oxide surfaces at high sliding velocity (10 m s −1 ) [8]. However, even with only 0.2 nm RMS roughness, the friction coefficient is higher, and is reduced by less between low and high coverage, than on atomically smooth surfaces [8]. The friction-coverage behavior changes at 0.5 nm RMS roughness, where the friction coefficient at 1.52 nm −2 and 3.04 nm −2 are virtually identical, but there is still a decrease of 11% between 3.04 nm −2 and 4.56 nm −2 coverage. At 0.8 nm RMS roughness, there is a continuous decrease in friction coefficient with increasing coverage, by 5% between 1.52 nm −2 and 3.04 nm −2 and by 12% between 3.04 nm −2 and 4.56 nm −2 . The general reduction in in friction coefficient with increasing coverage can be explained through the reduction in interdigitation between stearic acid films on the upper and lower slabs (Fig. 3). The reduction in interdigitation means that the stearic acid molecules do not need to rearrange as much in order to facilitate sliding of the surfaces, resulting in a reduction in lateral force [46]. The increase in friction coefficient between 1.52 nm −2 and 3.04 nm −2 at 0.2 nm roughness is due to slower rearrangement of the more closelypacked stearic acid molecules, which leads to a higher friction coefficient at high sliding velocity (10 m s −1 ) [8]. At higher roughness, the films generally become more liquidlike (Figs. 3 and 4) and hence molecular rearrangement is faster, leading to a reduction in the friction coefficient at 3.04 nm −2 . Fig. 6b shows that F 0 decreases linearly with increasing surface coverage to a similar degree at all levels of RMS surface roughness. There is a decrease in F 0 of approximately 40% between 1.52 nm −2 and 3.04 nm −2 as well as between 3.04 nm −2 and 4.56 nm −2 . Previous NEMD simulations [15] also observed a similar decrease in F 0 with increasing coverage on atomically smooth slabs; however, there was much less change in F 0 on slabs which contained nanoscale roughness features. Analysis using the smooth particle method (SPM) showed that F 0 was exponentially related to the interface contact area between the upper and lower films [15]. In the current simulations, the interdigitation between the upper and lower films is reduced at increasing coverage regardless of the RMS surface roughness (Fig. 3). This is because the smoother interface (Fig. 2) results in a smaller contact area and thus a lower F 0 . The discrepancy between the changes in F 0 with coverage at different levels of roughness may arise from differences in the method used to impose the surface roughness, or perhaps even the force-field applied; this is discussed further below. Fig. 7 shows that both µ and F 0 increase linearly with increasing RMS surface roughness. From Fig. 2 and Fig. 3, as well as previous SPM analysis [15], this can be attributed to an increase in interdigitation, leading to increased interface contact area and thus increased adhesion and resistance to sliding. It is not due to local breakdown of the OFM film or solid-solid contact (Fig. 2). The gradients of the increase in µ and F 0 with RMS roughness are similar at all coverages. The gradient for the change in F 0 with RMS roughness suggests a value of just 0.03 GPa at 4.56 nm −2 and 0.08 at 1.52 nm −2 coverage on atomically smooth surfaces, which agrees well with previous NEMD simulations [15]. However, reference [15] found that, at low coverage, the change in µ and F 0 were less dependent on roughness than at high coverage. This was explained through the fact that, in solidlike, high coverage films, roughness features below the film are almost perfectly reproduced on top of it; whereas in liquidlike, low coverage films, the stearic acid molecules have more freedom to rearrange themselves according to the surface features of the slab, filling the troughs between asperities and ensuring a more even sliding interface [15]. However, Figs. 2 and 3 demonstrate that on these larger slabs, which provide a more faithful representation of the statistical distribution of the heterogeneous surface morphology [26], the interface is much smoother at higher coverage, even at 0.8 nm RMS roughness. Another factor to the differences between the current and previous results may be the force-field used to represent the stearic acid molecules. Reference [15] represented stearic acid with the OPLS-AA force-field [39], which has been shown to yield elevated melting points for long-chain molecules [9]. The use of OPLS-AA [39] may lead to more solidlike films than when a more accurate force-fields for long-chain molecules, such as L-  OPLS-AA [38] are used (as they are in the current simulations). Therefore, the use of OPLS-AA [39] may result in high coverage films which are less deformable, leading to rougher interfaces and thus more adhesion and resistance to interfacial sliding at high coverage than if L-OPLS-AA [38] were employed. From these simulations, it is clear that µ and F 0 are generally decreased with increasing stearic acid coverage and decreased by increasing nanoscale RMS roughness. Generally, the gradients in Fig. 7 are shallower than those in Fig. 6, suggesting that µ and F 0 are influenced more heavily by stearic acid coverage than nanoscale RMS surface roughness. Moreover, comparing Fig. 6a to Fig. 7b and Fig. 6a to Fig. 7b, it seems that F 0 is more susceptible to changes in stearic acid coverage and nanoscale RMS roughness than µ (note that the µ y-axis is truncated).

Friction
Relating these simulations to experiments, the lateral (friction) force measured in OFM-lubricated systems is likely to be heavily influenced by F 0 at lower F N , i.e. lower applied loads and thus contact pressures ( < 0.5 GPa) [27]. However, in this case, there is likely to be a significant number of lubricant molecules present which separate the OFM films; and this has been shown to eliminate F 0 [8,15]. Conversely, when F N is high, as in the boundary lubrication regime, the lateral force experienced is more heavily influenced by changes in µ, and any variation in F 0 becomes less significant (Fig. 5). Therefore, the extended Amonton-Coulomb law under the high load approximation [7,8], µ≈F L /F N , should provide an accurate estimation of the friction coefficient in OFM-lubricated systems in the boundary regime.
The µ value observed with no stearic acid molecules present between the α-iron slabs decreases with increasing RMS roughness; from 0.2 nm (µ=0.28), to 0.5 nm (µ=0.25) and 0.8 nm (µ=0.24), which agrees well with previous NEMD simulations using the same α-iron force-field and similar conditions [22]. This trend is due to the reduction in solid-solid contact area at higher roughness. In these previous simulations, slabs with 0.5 nm RMS roughness lubricated with n-hexadecane reduced the friction coefficients by 10% at 1.48 nm −2 (µ=0.22) and 30% at 2.96 nm −2 (µ=0.17) compared to the unlubricated case [22]. In the current simulations, at 0.5 nm RMS roughness, a similar coverage of stearic acid yielded friction coefficients reduced by 40% at 1.52 nm −2 (µ=0.149) and 47% at 4.56 nm −2 (µ=0.132) compared to the unlubricated case. The µ value obtained with n-hexadecane is both higher and more sensitive to changing coverage than when stearic acid is used. This is because when the rough surface is lubricated with n-hexadecane, direct contact of asperities on opposing surfaces occurs; however, this is avoided when the rough surfaces are lubricated by stearic acid, even at low coverage. This can be attributed to the strong adsorption of the carboxylic acid head group to the α-iron surface, which prevents molecules being squeezed out from between asperities [8,27].
In High Frequency Reciprocating Rig (HFRR) experiments [6], µ has also been observed to decrease with an increasing concentration of stearic acid. The concentration of OFMs in a lubricant is generally used as the variable in tribology experiments rather than surface coverage, because it is far easier to measure and control [8]; however, recent depletion isotherm experiments have shown that stearic acid adsorption follows the Langmuir isotherm model [2]. In HFRR experiments on steel surfaces, a reduction in µ with an increasing concentration of stearic acid in n-hexadecane has been observed; from 0 mmol dm −3 (µ=0.25), to 0.1 mmol dm −3 (µ=0.18), to 1 mmol dm −3 (µ=0.12), and 10 mmol dm −3 (µ =0.10) [6]. These results agree well with the decrease in µ with increasing coverage of stearic acid observed in these NEMD simulations. The µ values from the simulations are somewhat higher than those from the experiments due to the higher sliding velocities employed (10 vs. 0.05 m s −1 ) [5], as predicted by stress-promoted thermal activation theory [6,9].
There is much interest in the different friction behavior of saturated OFMs and those with Z-unsaturation in their tail groups [2][3][4][5]. Mixtures of these are generally used collectively in commercial additive formulations [1], which may affect their friction reduction performance [5,8]. Our previous NEMD simulations [8] showed that, at a given coverage, films of Z-unsaturated OFMs, such as oleic acid, actually had very similar nanoscale structures as those formed by saturated OFMs, such as stearic acid. However, by comparing the friction-velocity behavior from NEMD simulations and experiments, it was suggested that OFMs with Z-unsaturated tail groups yield lower coverage films than saturated OFMs on steel surfaces [8]. Specifically, in the NEMD simulations low coverage films gave higher friction coefficients which were independent of sliding velocity [8], as experimentally observed for oleic acid [5], whereas high coverage films yielded a low friction coefficient which increased linearly with the logarithm of sliding velocity [8], as experimentally observed for stearic acid [5]. This postulation has recently been confirmed using depletion isotherm experiments [2] which showed a much lower plateau coverage for oleic acid (Γ≈2 nm −2 ) than stearic acid (Γ≈4 nm −2 ) on iron oxide from nhexadecane. Therefore, assuming a high experimental concentration, the low coverage films in these simulations are representative of those observed experimentally for oleic acid, whist the high coverage films are comparable to stearic acid. The structure of the stearic acid films varied less with coverage in the presence of nanoscale RMS roughness (Figs. 3 and 4); however, the high coverage films gave lower lateral forces (Fig. 5), through reductions in both µ and F 0 (Figs. 6 and 7). This suggests that OFMs with Z-unsaturation will be significantly less effective in reducing friction than those with completely saturated tails on surfaces with nanoscale RMS roughness as well as atomically smooth ones.

Conclusions
In this study, we have used large-scale nonequilibrium molecular dynamics (NEMD) simulations to examine the nanoscale structure and friction of stearic acid, a model organic friction modifier (OFM), adsorbed on iron surfaces with 3D nanoscale RMS roughness. Three different coverages of stearic acid (1.52, 3.04, 4.56 nm −2 ) were adsorbed between iron slabs with three different levels of RMS roughness (0.2, 0.5, 0.8 nm). High (0.5-2.0 GPa) pressures, as experienced between asperities in the boundary lubrication regime, were simulated in order to investigate the robustness of stearic acid monolayers under these rather extreme conditions. The stearic acid films were able to maintain separation of asperities on opposing surfaces under even the highest RMS roughness (0.8 nm) and lowest coverage (1.52 nm −2 ) systems simulated, due to strong adsorption of the head groups which prevented squeeze-out. As a result, the friction coefficient was reduced by more than 40% compared to when no stearic acid molecules were present between the slabs. Moreover, comparing the results of the current study to previous NEMD simulations [22], the stearic acid films were found to be significantly more effective than n-hexadecane in reducing friction on surfaces with 0.5 nm RMS roughness. The decrease in friction coefficient moving from n-hexadecane to stearic acid is consistent with that observed in boundary friction experiments.
The results of this study suggest that the lateral (friction) force measured experimentally is likely to be heavily influenced by the Derjaguin offset at lower pressure ( < 0.5 GPa). However, in this case, there is likely to be lubricant present to separate the stearic acid films, which eliminates the Derjaguin offset. Conversely, when the pressure is high, as in the boundary lubrication regime, the lateral force experienced is more heavily influenced by changes in the friction coefficient, and any variation in the Derjaguin offset becomes less significant.
These simulations reiterate the key role of OFM coverage in the film structure and boundary friction reduction. On surfaces with nanoscale roughness, systems with a higher coverage of stearic acid generally yielded lower lateral (friction) forces due to reductions in both the friction coefficient and Derjaguin offset. This is due to the formation of solidlike films, which allow little interdigitation with the opposing film, resulting in smooth interfaces and low resistance to interfacial sliding.
Surfaces that have greater levels of nanoscale RMS roughness have more disordered, liquidike stearic acid films; however, the friction coefficients and Derjaguin offsets is only slightly increased, despite the rather extreme steepness of the roughness features. Therefore, these results suggest that OFMs are only slightly less effective in reducing friction on surfaces with nanoscale roughness as those which are atomically smooth.
This study showed no evidence of any collapse of stearic acid films at asperity tips even at the highest contact pressures; the changes in friction observed with pressure and coverage originated solely from changes in the structure of the OFM films. Since the surface roughness simulated was considerably steeper than the level expected in real engineering components, this suggests that such OFM film collapse between asperities is unlikely on engineering surfaces under the conditions studied. The temperature profiles at all coverages are consistent with what is expected for well-equilibrated boundary-thermostatted confined systems subjected to relatively low sliding velocities [45,47]. The temperature is lower in the centre of the stearic acid films than in the slabs, suggesting that the systems are in quasi-equilibrium at the sliding velocity employed here (10 m s −1 ) [45].