Molecular dynamics simulation of the Stribeck curve: Boundary lubrication, mixed lubrication, and hydrodynamic lubrication on the atomistic level

Lubricated contact processes are studied using classical molecular dynamics simulations for determining the entire range of the Stribeck curve. Therefore, the lateral movement of two solid bodies at different gap height are studied. In each simulation, a rigid asperity is moved at constant height above a flat iron surface in a lubricating fluid. Both methane and decane are considered as lubricants. The three main lubrication regimes of the Stribeck curve and their transition regions are covered by the study: Boundary lubrication (significant elastic and plastic deformation of the substrate), mixed lubrication (adsorbed fluid layer dominates the process), and hydrodynamic lubrication (shear flow is set up between the surface and the asperity). We find the formation of a tribofilm in which lubricant molecules are immersed into the metal surface—not only in the case of scratching, but also for boundary lubrication and mixed lubrication. The formation of a tribofilm is found to have important consequences for the contact process. Moreover, the two fluids are found to show distinctly different behavior in the three lubrication regimes: For hydrodynamic lubrication (large gap height), decane yields a better tribological performance; for boundary lubrication (small gap height), decane shows a larger friction coefficient than methane, which is due to the different mechanisms observed for the formation of the tribofilm; the mixed lubrication regime can be considered as a transition regime between the two other regimes. Moreover, it is found that the nature of the tribofilm depends on the lubricant: While methane particles substitute substrate atoms sustaining mostly the crystalline structure, the decane molecules distort the substrate surface and an amorphous tribofilm is formed.


Introduction
Understanding the fundamental processes of lubrication is critical for many technical applications. In many technical applications, e.g. machining in production processes, gears, and bearings, extreme conditions (pressure and temperature) prevail in the contact. Moreover, the actual contact process and dominating mechanisms occur on very small scales. Therefore, experimental in situ methods for investigating the mechanisms in tribocontacts-in particular on the atomistic level-are presently not feasible. Molecular dynamics (MD) simulations [1,2] using classical force fields [3][4][5] are an attractive tool for studying the fundamental mechanisms in tribological processes, e.g. the mechanisms of the different lubrication regimes of the Stribeck curve. MD provides a great level of insights for a large number of observables. Also, the strong physical basis of the simulation method requires no parametrization on experimental data of the actually studied process. However, the accuracy of the results strongly depends on the correctness of the applied force field describing the molecular interactions and the simulation setup.
Lubrication regimes are usually classified using the Stribeck curve [6]. Thereby, depending on the gap height between two solid bodies, different lubrication conditions and mechanisms dominate the contact process [7][8][9]. Figure 1 schematically shows the lubrication regimes applicable when the two solid bodies are in close proximity. In the hydrodynamic lubrication (HL) regime, the two solid bodies are separated by the fluid in a way that a bulk fluid phase exists in the gap. Upon relative motion between the solid bodies, a shear flow is established in the fluid. However, the fluid pressure in the gap may be high enough such that elastic deformation occurs in the solid bodies. Viscous friction dominates the system in the HL regime. In the mixed lubrication (ML) regime, the two solid bodies are in close proximity and the gap height has approximately the thickness of the adsorption layers of the fluid on the solid surfaces. Both elastic and plastic deformation may occur. The friction is influenced by a large number of mechanisms in the ML regime, e.g. the distortion of the adsorbed fluid layer, viscous friction, and plastic deformation. In the boundary lubrication (BL) regime, the solid bodies are in contact and plastic deformation is caused by the asperity penetrating the substrate. The friction is dominated by the elastic and plastic deformation. For studying the interrelations between bulk phase properties and mechanisms in the direct vicinity of the contact zone, e.g. the interplay of the squeeze-out of the fluid and the formation of a tribofilm, large system sizes have to be considered. In this work, we used large-sized molecular dynamics (MD) simulations to study different lubrication regimes along the Stribeck curve in a model system to elucidate the dominating mechanisms in the different regimes.
MD simulation has been extensively used in recent years for studying friction on the atomic scale. In particular, the hydrodynamic lubrication regime has been studied many times in the literature using MD simulations, e.g. Refs.  for a recent review. Also, contact processes comprising deformation of the solid bodies of dry contacts have been extensively studied in the past [32][33][34][35][36][37][38]. Lubricated systems with elastic and/ or plastic deformation on the other hand have been less frequently studied using classical MD simulations [39][40][41][42][43][44][45][46][47][48][49][50][51][52]. Moreover, lubricated contact processes comprising plastic deformation, where the solid bodies are fully immersed by the fluid such that also a fluid bulk phase exists that exhibits heat conduction, load transfer, back pressure against the squeeze-out etc. has been only rarely investigated [39, 40, 42-46, 53, 54]. Also, to the best of our knowledge, the transition of the different lubrication regimes has not yet been studied using MD.
In this work, we study large-sized systems where the contact process is fully immersed in a fluid to elucidate the interplay of different friction mechanisms (e.g. the formation of tribofilm, break-up of fluid adsorption layers, squeeze-out, elastic and plastic deformation, dissipation, and heat conductivity) in the different lubrication regimes. In the system, a spherical indenter (modeling the tip of an asperity) carries out a lateral movement on an atomistically flat substrate. The studied setup is a simplified yet a representative one. The substrate is an iron single crystal and the indenter is a single crystal diamond. The effect of two lubricants is compared: Methane and n-decane (for brevity, referred to as 'decane' in the following). The influence of the lubrication gap height is systematically studied to investigate the three different lubrication regimes BL, ML, and HL along the Stribeck curve. Hence, simulations were carried out in this work at constant gap height, which is different to the usual approach, where the load is prescribed. We consider both a conceptual model system using the simple Lennard-Jones model as lubricant (modeling methane) as well as studying a practically relevant system using decane as lubricant.
The kinematics of the indenter are fully prescribed. Hence, the forces on the indenter at a given gap height are a response to the kinematics. A large number of mechanical and thermodynamic properties of the system were studied in detail. For the mechanical properties, the total forces on the indenter, the coefficient of friction, the groove and chip formed on the substrate, the behavior of the fluid in the adsorption layer, and the formation of the tribofilm in the substrate surface are discussed. For the thermodynamic properties, the dissipated energy, the temperature field, and the pressure field of the fluid are studied. This comprehensive concept study enables a systematic assessment of the fundamental mechanisms on the different lubrication regimes of the Stribeck curve, as shown in Fig. 1.
The present work is outlined as follows: Section 2 introduces the simulation setup and the observables used to characterize the contact process. The results for both the mechanical and thermodynamic properties are discussed in section 3.1 and 3.2, respectively. www.Springer.com/journal/40544 | Friction solid substrate, and a fluid lubricant. Two types of simulations were carried out: One with methane as fluid and one with decane as fluid. For each type, simulations at different gap heights were carried out. The size of the substrate block was 401, 763, and 31 Å in x-, y-, and z-direction, respectively. The surface of the substrate was in the xy-plane. The indenter was cap-shaped with a radius of R = 500 Å. The height of the indenter cap was 70 Å. It carried out two sequential movements: (i) positioning/ indentation, during which the indenter moves in -z-direction normal to the surface; and (ii) lateral movement, during which the indenter moves in +x-direction. To systematically study the three lubrication regimes BL, ML, and HL, the indentation depth was varied systematically. A given indentation depth h  yields a certain lubrication gap h between the indenter and the substrate, as shown in Fig. 2. The indentation gap h  was defined as the distance between the center of the particles of the bottom of the spherical cap and the top of the substrate (Fig. 2). However, the lubrication gap height h was used throughout the study for characterizing the simulations, which was defined as

Simulation scenario
where IS  is the size parameter of the interaction potential between the substrate and the indenter. Hence, for 0, h  the substrate and the indenter would directly interact with each other, even if no fluid was present.
For both fluids, 14 gap height values h were considered. The following values for the lubrication gap were used h = -7, -6, -5, -4, -3, -2, -1, 0, 1, 2, 4, 6, 8, 12 Å. Hence, for h < 0, elastic and/ or plastic deformation will certainly occur, whereas for h > 0, de-formation of the substrate may in general occurdepending on the behavior of the fluid in the gap. For a more detailed analysis (see results section), the focus is put on the simulations with h = -6, -4, -2, 0, 4, 8 Å. Hence, in total 28 simulations were carried out. The indenter velocity during the vertical positioning was 10 m/s and during the lateral movement 20 m/s. It is a long-standing and often discussed problem that velocities used non-equilibrium MD (NEMD) friction simulations do not comply with those used in classical friction experiments [77,78]. Accordingly, the speed in most classical friction experiments is significantly smaller than our scratching velocity [79,80]. However, interestingly, 20 m/s is a realistic cutting speed used in manufacturing processes such as grinding and turning [81].
The simulation time 0   indicates the beginning of the indentation phase, which positions the indenter at a given gap height h. Depending on the gap height h, the positioning phase was between 0 440ps    and 0 630ps.

  
The lateral movement was carried out for 1,400 ps. Therefore, the scratching length was 280 Å.
The origin of the Cartesian coordinate system was defined as: 0, x  0 y  in the corner of the simulation domain and 0 z  for the initially unpenetrated surface of the substrate, as shown in Fig. 2. The substrate consisted of 5 8 10  particles. The substrate was fixed in the simulation box by fixing 3 atom layers at the bottom of the substrate (Fig. 2). Periodic boundary conditions were applied in x-, and y-direction. It should be noted that this might in general affect the fluid flow in the system. However, as the simulation box is very large, artifacts from the periodic boundary conditions are probably small.
The motion of the indenter was imposed by prescribing the velocity of the inner core atoms of the indenter, as shown in the dark green region in Fig. 2. Hence, the outer layer atoms of the indenter moved according to Newton's equation of motion, as shown in the light green region in Fig. 2. The indenter consisted of 5 7 10  particles. The radius of the base of the indenter cap was 500 Å. To limit the total number of particles in the box, the cap indenter was trimmed in x-direction to a total width of 350 x   Å (Fig. 2). Hence, the cap indenter had a front and rear face side. However, due to the overall large indenter size, no significant effect of the sideways trimming on the contact process at the bottom of the indenter are expected.
At the top, the simulation domain was confined by a soft repulsive wall. For the initialization, the fluid particles were placed on a regular grid. To establish a liquid state, the soft repulsive wall was lowered such that the fluid particles were compressed. Before the compression, the entire simulation box was equilibrated for 1 ns. This compression was carried out before the movement of the indenter. At the beginning and during the compression, the entire simulation box was 1 44 10   interaction sites in the case of decane. Each methane molecule consists of a single interaction site; each decane molecule consists of 10 interaction sites. The simulation box had in total (indenter, fluid, and substrate) 6 2 37 10   interaction sites in the methane simulation cases and 6 2 94 10   interaction sites in the decane simulation cases.
The time step was 1 fs and 0.5 fs for the methane and decane simulations, respectively. The simulations were carried out with the open-source simulation code LAMMPS [84]. The visualizations were performed using the open access tools ParaView [85] and OVITO [86]. Table 1 gives an overview of the symmetric interaction matrix present in the simulations. The diagonal entries of the matrix indicate the interactions between interaction sites of the same substance. The off-diagonal entries in Table 1 specify the cross-interactions between sites of different substances, i.e. indenter-substrate, fluid-substrate, and indenter-fluid. The substrate and the indenter force field was the same in all simulations. As two different substances (methane and decane) were considered for the fluid, two different force fields were used for their modeling.

Molecular model
The substrate was an iron bcc single crystal with an initially atomistically flat surface and a lattice constant of 2.855 Å. It should be noted that these are strong simplifications compared to a real system exhibiting a surface roughness [87,88] and dislocations [89]. The substrate-substrate interactions were prescribed by the iron EAM potential proposed by Mendelev et al. [90]. The crystal had a (100) surface orientation and the lateral movement was performed in [01 _ 1] direction (it should be noted that the lattice orientation with respect to the scratching direction has in general an influence on the process properties [91]).
The diamond indenter was modeled by a Tersoff potential [92]. The potential was developed to describe the intermolecular interactions in carbon diamond. The atoms of the indenter were initialized on a diamond lattice with a lattice constant of 3.5667 Å.
Two different fluids were used as lubricants: Methane and decane. Methane was modeled using a BZS force field proposed by Vrabec et al. [4,93]; decane was modeled using a TraPPE force field proposed by Martin and Siepmann [3,94]. Both the BZS and TraPPE force field are known to describe properties Table 1 Interaction potentials employed in the present study.
www.Springer.com/journal/40544 | Friction of fluids very well, e.g. for the vapor pressure [93,95], the saturated densities [93,95] as well as transport properties such as the viscosity [96] and interfacial properties [93,95]. Also, this force field class is known to usually yield reasonable extrapolation behavior to state points that were not considered for the parametrization [97,98]. The indenter-substrate interactions were modeled using the purely repulsive Weaks-Chendler-Anderson (WCA) potential [99]. Hence, the classical 12-6 Lennard-Jones potential with a cut-off at the minimum of the potential was used. Fluid particles interact with both substrate and indenter atoms. In this work, the solid-fluid interactions were modeled with the (full) 9-6 Lennard-Jones potential. In general, the parameters  and  for the solid-fluid interactions can be determined by different routes, i.e. a top-down or a bottom-up approach. In the top-down approach, the interaction parameters are chosen such that the system described by the force field model captures a certain macroscopic behavior, e.g. an experimentally determined contact angle or adsorption isotherm. In the bottom-up approach, the interaction parameters are derived from quantum chemical calculations. In the present work, the bottom-up approach was used. We used the force matching method outlined by Ercolessi and Adams [100]. The quantum chemical data was taken from Ref. [101]. Details are given in the Electronic Supplementary Material (ESM). It should be noted that the solid-fluid cross-interactions-in particular the solid-fluid interaction energy  -have an important influence on interfacial and tribological processes [42,43,64,65,70,71], such as adsorption layers on the surface [64,65,71], static and dynamic wetting effects [65,70], the heat transfer across the interface [64], and accordingly an influence on the thermal balance of tribological processes [43]. Also, the squeeze-out of particles from the lubrication gap is known to depend on the solid-fluid interactions [42,43].
The chosen molecular models and interaction potentials provide a good trade-off between accuracy and computational efficiency, which is crucial for the large-scale simulation study carried out here.

Definition of observables
A large number of primary simulation data was sampled and then further post-processed such that multiple aggregated observables became available for the evaluation of the process. Both mechanical and thermodynamic properties were considered. The focus is on the mechanisms occurring in the direct vicinity of the contact zone.
The total tangential force T F and normal force N F on the indenter were calculated as the sum of all pair interactions (indenter-substrate (IS) and indenter-fluid (IF)) on the indenter in yand z-direction.
where m is the total number of interaction sites of the indenter. The tangential force is defined here in positive y-direction (scratching direction) and thereby indicates the reaction force of the scratching process. The coefficient of friction was computed from the absolute value of the tangential force and the normal force on the indenter During the process, the forces on the indenter were written to the hard-drive every 1 ps. The error bars of the forces and the coefficient of friction were calculated as their standard deviation during the quasi-stationary phase of the lateral movement (see below). The work done by the indenter W was computed by numerical integration of the total force on the indenter, where the forces N F and T F were computed according to Eqs. (1) and (2). The heat of the thermostat Thermostat Q is defined as the cumulative energy that was added or removed from the particles in the thermostat layer of the substrate (Fig. 2) where Scratch 0   indicates the beginning of the lateral movement.
The temperature and the pressure of the fluid were computed locally during the contact process. The fluid temperature and pressure were computed as per-atom values by averaging the respective property over the atoms within the cut-off radius of each atom. The average aligned velocity from each per-atom collective was subtracted from the total average kinetic energy for the computation of the local temperature. The fluid region was divided into bins in the y-x plane (Fig. 3) over the entire simulation domain. The edge length of the bins was 4 properties were computed bin wise, i.e. averaged over all molecules in a bin. Only fluid particles below the indenter were considered, i.e. below the respective z-coordinate of the upper edge of the indenter. These temperature and pressure profiles were computed in block averages every 20 ps. The maximum temperature and pressure difference max T  and max p  were derived from these temperature and pressure fields, respectively, as the difference of the maximum value of all bins and the mean bulk value derived in the area unaffected by any adsorption layers. This was done for all block N The black rectangle with 4 6      x y Å and the volume dV indicates the shape of an averaging bin for the computation of the temperature, pressure, and density. The bin represents a grid that is distributed over the entire simulation box. averages in the quasi-stationary regime. The final value of a given observable and its error bar was calculated as the average and the standard deviation of all values sampled during the quasi-stationary phase.
The substrate and the indenter surface was analyzed using an alpha shape algorithm [102]. Thereby, the bounding hull surface of the indenter and the substrate is obtained. These surfaces are defined by sets of points. Using the alpha shape method, the obtained surfaces capture both convex and concave parts of the solid body surfaces. For determining the particles that constitute the surface, an alpha radius of Alpha,c 1 05R  was used, where Alpha,c R indicates the critical radius [102]. Here, Alpha,c R was increases by 5% such that holes and gaps in the deformed surface were avoided. The resulting geometrical bodies were used for further process characterizations. The volume of the chip Chip V was calculated from the volume bounding hull obtained from the alpha shape method. The chip was defined here as the geometric body formed by substrate particles above the initial height of the surface and ahead of the center of the indenter, as shown in Fig. 3. The depth of the groove Groove z  was defined as the difference of lowest z-coordinate of the upper surface of the substrate to its initial surface height. Both observables were evaluated every 20 ps in the quasi-stationary phase. From these block average values, the mean value and its standard deviation was computed. The number of fluid particles (i.e. interaction sites in the case of decane) in the gap Gap N and in the tribofilm Tribofilm N were calculated based the alpha shape surface analysis. The evaluation of the number of fluid particles in the gap Gap N and in the tribofilm Tribofilm N was carried out at the very end of the contact process. The number of fluid particles in the gap Gap N counts all fluid interaction sites below the surface of the indenter and above the surface substrate. Only particles within a radius of R yx = 100 Å from the center of the indenter and in front of the indenter are considered for Gap N (Fig. 3). and Tribofilm , N no statistical uncertainty is reported as they were only calculated from the atomistic configuration at the end of the contact process. Moreover, the density profiles of the adsorption layers ( ) z  were calculated from data taken after the equilibration in areas that were unaffected by the indenter.

Results and discussion
The results are structured as follows: First, the mechanical properties of the process are discussed, then, the thermodynamic properties. The mechanical properties comprise the kinetics of the process, the mechanical work done by the indenter, the formation of the chip and the groove in the substrate as well as the characterization of the lubrication gap, the adsorption of the fluid on the substrate, and the formation of the tribofilm in the substrate surface. The thermodynamic properties comprise an analysis of the overall dissipated energy as well as the temperature and pressure field of the fluid in the lubrication gap. From the 14 gap heights that were considered in this work, six characteristic gap heights that capture the main characteristics of the Stribeck curve lubrication regimes are discussed in detail in the following. Figure 4 shows the results for the total forces on the indenter during the entire contact process for the exemplaricly chosen gap height 4 h   Å for illustration of the overall contact process kinetics. The normal force N F and the tangential force T F on the indenter show the characteristic shape for a subsequent indentation and lateral movement process [42,43,103]. During the indentation, the normal force exhibits a steep increase as the indenter causes (at the exemplaricly chosen h) elastic and plastic deformation on the substrate. Moreover, the bulk fluid phase is squeezed out of the lubrication gap, which requires work done by the indenter [43]. During the lateral movement, the normal force quickly drops and reaches an approximately constant value. The tangential force fluctuates around zero during the indentation. At the beginning of the lateral movement, the tangential force exhibits a jump and then fluctuates around an approximately constant value, which is due to the formation of the chip. During the lateral movement, the process reaches a quasi-stationary state at approximately 1,000 ps.

 
During that quasi-stationary state, observables were averaged to yield heightspecific process parameters (see methods section for details). Figure 5 shows mechanical process parameters sampled and averaged during the quasi-stationary phase as a function of the gap height: Fig. 5(a) shows the normal and tangential total force on the indenter;   Figure 6 shows screenshots of the visualization of the contact process during the quasi-stationary phase taken at the end of the lateral movement Scratch 400 ps,   which provides direct insights into the atomistic mechanisms of the process. From Fig. 5, it can be seen that the normal total force on the indenter N F decreases with increasing gap height. This is as expected: As the process kinematics are prescribed, the normal force on the indenter is large at small gap height. Also for the tangential force, the (absolute) value T F   obtained during the quasi- is no direct contact between the indenter and the substrate and only the fluid shear flow acts as a mediator between the relative motion of the bodies. For small (negative) gap heights, i.e. the indenter penetrates the substrate, the normal force is about an order of magnitude larger than the corresponding tangential force. This is mostly due to the indenter geometry, which is relatively flat and therefore has a high contact area with the substrate. The comparison of the results for the methane and the decane simulations shows that, qualitatively, the results for the two fluid cases are similar. Yet, important quantitative differences are observed. For small gap heights h < 1 Å, the absolute values of the forces N F   and T F   in the decane simulation cases exceed those in the corresponding methane cases. Vise versa, for large gap heights h > 1 Å, the absolute values of the forces are larger in the methane cases. Interestingly, the cross-over of the methane and decane force occurs for smaller gap height for the tangential force than for the normal force. This transition can also be clearly seen in the behavior of the coefficient of friction, as shown in Fig. 5(b).
From the results for the coefficient of friction, the three friction regimes and the typical shape of the Stribeck curve can be clearly identified: For small gap heights values h < -1 Å, BL occurs. For intermediate gap heights -1 Å < h < 4 Å, ML prevails. For large gap heights 4 Å < h, HL occurs. The identification of the three lubrication regimes is also supported by the simulation screenshots depicted in Fig. 6. For 4 Å < h, the coefficient of friction is very small, which is due to the fact that the contact process is fully mediated by the hydrodynamic shear flow of the fluid and only very little elastic and plastic deformation is done on the substrate. Hence, HL prevails. In that regime the fluid bulk phase properties dominate the process behavior. In the BL regime h < -1 Å, significant elastic and plastic deformation of the substrate occurs as a result of the penetration of the indenter into the substrate. Moreover, a tribofilm ① is formed on the ① The term 'tribofilm' is usually associated with chemical reactions between the lubricant and the substrate surface. The term tribofilm is adopted here despite the fact that nonreactive force fields were used since the heterogeneous phase formed from substrate and fluid particles has similar characteristics as tribofilms.
substrate surface in the BL regime (Fig. 6). In the BL regime, the coefficient of friction increases with increasing gap height. Moreover, the coefficient of friction in the BL regime is about an order of magnitude larger than that in the HL regime. In the ML regime -1 Å < h < 4 Å, a smooth transition occurs from the BL regime to the HL regime. The adsorption layer of the fluid on both solids plays a dominant role in the ML regime as most fluid particles in the gap are part of the adsorption layer (Fig. 6). The coefficient of friction exhibits a maximum at the lower bound of the ML regime and strongly decreases with increasing gap height towards the HL regime. The results for the Stribeck curve obtained in this work from MD simulations ( Fig. 5(b)) are in good qualitative agreement with experimental data, e.g. Leong et al. [104]. Experimental Stribeck curves usually show the evolution of the coefficient of friction with a dimensionless parameter, such as the suitably scaled speed of contact. In our case, the contact separation, i.e. gap height h appears, which is usually scaled with the surface roughness. Since the latter vanishes in our case, the choice of such a dimensionless parameter is non-trivial here and we study the Stribeck curve as a function of h.
The three lubrication regimes can moreover be characterized based on the slip in the system. In the HL regime, the shear slip occurs in the bulk fluid phase. In the ML regime, the slip plane lies within the adsorption layer of the fluid molecules on the solid surfaces. The slip in this solid-fluid interfacial region is energetically less favorable and yields an increased coefficient of friction. In the BL regime, the slip lies in parts below the substrate surface, which leads to the formation of the chip and the groove. This slip plane is energetically less favorable and more work has to be done by the indenter.
A detailed evaluation of the force contribution is given in the ESM. Therein, it is shown that the substrate-indenter interactions contribute significantly to the total force on the indenter in the BL regimefor both the normal and the tangential force. This is despite the fact that at all studied gap heights, a closed fluid monolayer remains between the two solid bodies (Fig. 6). However, as the indentersubstrate interactions have a larger range than the size of the fluid particle beads, there is a direct interaction The comparison of the methane and decane results provides further insights. Interestingly, in the direct vicinity of the cross-over between the methane and the decane results for the coefficient of friction, the results for ( ) h  for both fluid cases exhibit a maximum. It is not clear if this is a coincidence. In the ML and HL regime, the methane simulations yield a higher coefficient of friction than the decane simulations. This is as expected: In the ML regime, the adsorption behavior of the fluid is crucial. Methane, as a simple practically spherical molecule, forms very structured adsorption layers (see details below), which dampens the slip-stick motion. Decane, on the other hand, forms a more unraveled adsorption layer, which enhances the slip-stick motion on the surface. In the HL regime, the fluid bulk phase properties prevail and decane has more favorable tribological properties than methane.
The work done by the indenter is depicted in Fig. 5(c). As the tangential force absolute value in the decane simulations exceeds that in the methane simulations at given h, also the work done by the indenter is larger in the decane simulations. This holds for the entire BL regime and for the lower bound of the ML regime. However, in the range 1 Å < h < 6 Å, more work is done by the indenter in the methane cases. This is due to the fact that the strongly structured adsorbed methane molecules interact energetically unfavorable in the bulk fluid inlet. This also results in a higher pressure increase in the lubrication gap, e.g. at h = 4 Å, as shown in Fig. 15 below.

Characterization of chip and groove
The formation of the chip and the groove in the substrate are characteristic properties of the contact process. Figure 7 shows the substrate surface at the end of the contact process (view direction is in negative z-direction) for the cases h = -6, -4, -2, 0, 4, 8 Å; indenter and fluid particles are not shown. Figure 8 shows the chip volume Chip V and the groove depth  the case of decane. Also the groove depth decreases with increasing gap height.
In the HL regime h > 4 Å, the magnitude of the groove depth and the chip volume is very small, as shown in Fig. 7. Hence, both the elastic and plastic deformation is very small in that regime-as expected. In the ML regime -1 Å < h < 4 Å, the adsorbed fluid layers strongly interact with the substrate surface. This causes elastic and plastic deformation and thereby the formation of both a faint chip and a groove in the substrate, as shown in Fig. 7. In the BL regime h < -1 Å, significant elastic and plastic deformation of the substrate www.Springer.com/journal/40544 | Friction occurs (Fig. 7). As the indenter has a very flat shape, the chip forms primarily sideways.
For both the chip volume and the groove, important differences are observed for the methane and decane simulations. In the BL regime, the Chip V results for the two fluid cases are essentially the same, whereas significant differences are observed in the ML regime. There, at a given gap height, the methane case yields a larger chip formation compared to the decane case. For all three lubrication regimes, the presence of decane molecules as a lubricant causes a deeper groove compared to methane molecules, as shown in Fig. 8.
In the case of methane, the groove depth is surprisingly approximately constant in the range h = -3…4 Å (mostly ML), as shown in Fig. 8. This is due to the fact that the adsorbed methane particles (especially the ones adsorbed on the substrate) strongly interact with the substrate surface (Fig. 6).
The spherical particles substitute lattice atoms in the substrate, which effectively increases the volume of the substrate. This compensates the volume loss due to chip formation (Figs. 7 and 8). This contributes also to the fact that the methane simulation results yield a larger coefficient of friction in the ML regime than the decane simulations.
Due to the presence of adsorbed fluid layers on the surfaces in the case of decane, the indenter size is effectively increased, which increases also the groove depth Groove z  at a given gap height h. Therefore, the groove depth Groove z  exceeds the gap height h. This holds for both the BL and the ML regime. On the contrary, for the methane simulations, the groove depth Groove z  is smaller than the gap height in the BL regime. This is due to the fact that a significant amount of the methane particles is imprinted into the substrate surface forming a tribofilm. In parts, the methane particles substitute lattice atoms of the substrate (see details below). This effectively increases the volume of the substrate. The fact that we obtain A detailed analysis of the machined surface (Fig. 7) provides further interesting insights. The groove and the chip formed by the indentation (left part of the trace on the substrate surface) in both the HL and ML regime is more pronounced in the methane cases compared to the decane cases. The fact that the imprint caused by the indentation is significantly more prominent than the groove and chip formed during the lateral movement in the HL and ML regime is due to the squeeze-out behavior of the fluid. For gap heights h > 0 Å, a significant amount of fluid particles remains in the gap at the end of the indentation. At the beginning of the lateral movement, more fluid particles are squeezed out of the gap. This decreases the effective size of the indenter and thereby also decreases the groove and chip formed during the lateral movement compared to the indentation.
A further interesting aspect lies in the squeeze-out behavior of the fluid during the lateral movement in the BL regime. For h = -4…-2 Å, the width of the groove clearly decreases during the lateral movement. This is also due to the fact that more fluid particles are squeezed-out of the gap, which effectively decreases the indenter size. At very small gap heights, e.g. h = -6 Å, no such transient squeeze-out at the beginning of the lateral movement is observed, i.e. the formed groove and chip pile has a constant width x  during the entire contact process. Overall, in the BL regime, the squeeze-out of the fluid has important consequences for the formed groove and chip-depending on the gap height. The comparison of the chip formation in the methane and the decane case at h = -6 Å (Fig. 7) reveals that the chip in the decane case is larger (more pronounced and darker red coloring of the chip). This is due to the fact that the decane molecules form a more stable adsorption layer on both surfaces sideways of the indenter. This effectively increases the indenter size sideways. This is consistent with the finding that the chip removal in the decane case requires more work than in the methane case, as shown in Fig. 5-bottom. Interestingly, the fluid adsorption layer on the indenter surface does not break up (Fig. 6), whereas the adsorption layer on the substrate surface is distorted in the BL regime and in parts of the ML regime. Evidently, the (diamond) indenter has a higher hardness than the (iron) substrate. Therefore, no significant elastic deformation is observed for the indenter, whereas important elastic and plastic deformation is observed for the substrate. These deformations in the substrate surface are probably related to the break-up of the adsorption layer on the substrate surface: When a dislocation occurs on the substrate surface, also the adsorption layer in the direct vicinity of the dislocation experiences a perturbation that can initiate a break-up of the highly structured adsorption layer. A detailed investigation of this interesting phenomenon is out of the scope of this work.

Characterization of tribofilm and lubrication gap
A distinct feature observed in the simulations is the formation of the tribofilm caused by fluid particles being imprinted into the substrate surface. This can be seen from Fig. 7, where the groove shape is not colored smoothly, but with a color scattering. Also, the formation of the tribofilm can be clearly seen in the simulation screenshots depicted in Fig. 6. The mechanisms forming the tribofilm in the methane and decane case show important differences. Moreover, the tribofilm properties strongly depend on the gap height h. As the tribofilm formation is directly related to the behavior of the fluid particles in the lubrication gap, the two are discussed in the following corporately. Figure 9 shows the results of the analysis of number of fluid particles, i.e. interaction sites, in the tribofilm as well as those in the lubrication gap; Fig. 10 provides a detailed analysis of the unperturbed fluid adsorption layer; Fig. 12 provides a detailed analysis of the formation of the tribofilm; and Fig. 11 provides additional detailed screenshots illustrating the findings. Figure 9 shows the results for the total number of fluid particle interaction sites in the lubrication gap Gap N and in the tribofilm Tribofilm N (Section 2.3 for the     definition of the observables) as a function of the gap height h. With increasing gap height, the number of fluid particle interaction sites in the lubrication gap between the indenter and the substrate surface Gap N increases and the number of fluid particle interaction sites in the tribofilm in the substrate surface Tribofilm N decreases. These trends are obtained for both methane and decane. However, important differences are observed for the two fluids: At a given gap height, more decane particle interaction sites than methane particle interaction sites remain in the lubrication gap, whereas more methane particle interaction sites are imprinted into the substrate surface forming the tribofilm. Yet, the differences between the methane and decane results are significantly larger for Gap N compared to Tribofilm . N A second effect that increases the number of decane particle interaction sites in the lubrication gap in comparison to the methane case lies in the different squeeze-out behavior of the two fluids. The combination of the high pressure in the fluid and the molecular structure of the decane molecules results in a roughening of the substrate surface on the atomistic level in the gap, as shown in Figs. 6 and 11. This acts as a hindrance for the squeeze-out of the decane molecules and more molecules remain in the lubrication gap during the contact process.
The three lubrication regimes can be clearly identified from the behavior of Tribofilm ( ) N h (Fig. 9): In the BL regime, a significant amount of fluid particle interaction sites is imprinted into the substrate surface due to the close contact of the two solid bodies and the high contact forces. The trend of Tribofilm ( ) N h is approximately linear in the BL regime. In the ML regime, the adsorbed fluid layers are disrupted (Fig. 6) and some fluid molecules are imprinted into the substrate surface; the results for Tribofilm ( ) N h converge to zero in that regime. In the HL regime, practically no fluid particles are imprinted into the substrate surface and no tribofilm is formed-as expected. Figure 10 shows the density profiles of the adsorption layer of the fluid at the unpenetrated substrate surface, i.e. the initial state of the fluid in the gap before the starting of the indentation. In the case of methane, molecules form a strongly structured adsorption layer at the substrate surface. The first adsorption layer shows a distinct peak followed by a groove, where the local density decreases essentially to zero. Four more distinct adsorption peaks can be observed in the methane density profiles. The distance between these peaks is approximately δ z = 2.5 Å, which is in good agreement to 2   with  being the size parameter of the methane particles [93]. Hence, the second layer particles take body centered lattice positions with respect to the first layer. Opposite to the adsorption behavior of methane, a more complex behavior is observed for decane. This is simply due to the fact that different decane interaction sites can adsorb on the surface. Here, the CH 3 sites and CH 2 sites are distinguished. Due to the force field settings (see ESM for details), the CH 2 sites exhibit a stronger adsorption on the substrate surface. Therefore, the first adsorption peak stems from CH 2 sites. Due to the fact that the bond length between CH 2 and CH 3 sites is smaller than the size parameter of the two sites, an overlap of the first CH 2 and CH 3 adsorption peak is observed. With increasing z, only two more distinct peaks are observed in both the CH 2 and CH 3 profile. Also, due to the molecular structure of decane, the depth of the first groove does not decay to zero as observed for methane. Overall, the total fluid density profile in the case of decane shows a smoother transition from the surface into the bulk phase in comparison to methane due to the differences in the molecular structure. This has important consequences for the tribological behavior of the contact processin particular the stick-slip motion. Thereby, the more disordered adsorption layers in the case of decane contributes to a lower coefficient of friction in the ML and the HL regime compared to the methane case. Figure 11 shows detailed screenshots of the tribofilm and the lubrication gap from the atomistic visualization of the rear side of the indenter. The gap heights h = 8 Å and -6 Å are exemplarily shown. In the HL regime (Figs. 11(a) and 11(c)), the methane particles form a strongly structured solid-like adsorption layer, whereas the decane molecules form a more disordered adsorption layer. It can be clearly seen that, at h = -6 Å (Figs. 11(b) and 11(d)), the fluid adsorption layer is broken up and fluid particles are imprinted into the surface. The screenshots shown in Fig. 11-bottom moreover reveal important differences www.Springer.com/journal/40544 | Friction regarding the formation of the tribofilm during the contact process: During the indentation (left part of the groove), the number of methane particles imprinted into the substrate is high enough such that the substrate lattice structure is broken up and the tribofilm became amorphous. With ongoing lateral movement (middle and right part of the tribofilm depicted in Fig. 11(b)), less methane fluid particles are imprinted into the substrate and the tribofilm remains a crystalline structure. In the case of decane (Fig. 11(d)), the imprinted relatively large and elongated decane molecules cause a semi-amorphous-crystalline structure during the entire contact. Figure 12 shows the histograms of the fluid particle interaction sites in the tribofilm as a function of the coordinate z. The origin 0 z  lies on the unperturbed substrate surface, as shown in Fig. 2. In the HL regime, practically no fluid particles are imprinted into the substrate surface. In the ML regime, the break-up of the adsorption layer in combination with the high contact forces initializes the formation of the tribofilm. In the BL regime, the close contact of the two solid bodies in combination with elastic and plastic deformation of the substrate and the very high contact forces cause the formation of a prominent tribofilm. The formation of dislocations at the substrate surface in conjunction with deformations (formation of groove and chip) makes space for fluid particles. While it is energetically unfavorable for fluid particles to penetrate the atomistically flat surface, it becomes easier to penetrate the surface in the vicinity of a dislocation, i.e. a lattice defect. Hence, the deformation of the substrate increases the amount of energetically favorable locations, where fluid particles can be imprinted into the surface. Moreover, the imprinting of methane particles is energetically evidently less costly compared to the imprinting of elongated decane particles. Thus, methane particles are imprinted at lower gap height h, as shown in Fig. 12-top, compared to decane particles.
In the BL regime, the tribofilm has a thickness of approximately 10 Å (Fig. 12). Important differences are observed for the methane and the decane cases. At a given gap height, the thickness of the tribofilm is larger in the case of methane. More importantly, the structure of the tribofilm formed by the two fluids differs significantly: In the case of decane, a smooth profile of the number of imprinted fluid interaction sites is observed over the entire thickness of the tribofilm, which exhibits a single maximum slightly below the substrate surface 0. z  In the case of methane on the other hand, a profile with distinct oscillations is observed indicating that the methane particles substitute lattice atoms of the substrate. This is due to the fact that the methane particles have a similar size as the lattice atoms. This again shows that the substrate surface in the BL regime is strongly distorted in the decane cases, whereas the lattice structure remains more intact in the methane casesdespite the fact that more fluid particles are imprinted into the substrate surface.
The tribofilm comprises particles also at z > 0 Å (Fig. 12). This is due to the fact that fluid particles are also imprinted into the chip surface sideways of the indenter (Fig. 7). Moreover, the maximum peak of the tribofilm histogram slightly shifts deeper into the substrate body with decreasing gap height h. This is a result of the formation of the groove in the substrate (Fig. 7). Hence, effectively the substrate surface is pushed to lower z values, but the sampling of Tribofilm N  and Tribofilm N remains in the fixed coordinate system. Therefore, this shift is mostly an artifact of the definition of the observable.
Summarizing, for the tribofilm, methane particles substitute lattice atoms in the substrate, whereas decane particles cause a severe distortion of the surface upon being imprinted into the substrate surface. For the lubrication: While the methane adsorption layer is more structured than that of decane, the methane adsorption layer is broken up more easily. The tribofilm behavior is more relevant in the BL regime; the adsorption layer behavior is more relevant in the ML regime. Figure 13 shows the heat removed by the thermostat Thermostat Q during the contact process. At the beginning of the lateral movement, some transient processes occur (Fig. 13-top), which is in accordance with the course obtained for the forces on the indenter (Fig. 4). Then, at approximately 900 fs  

Thermodynamic properties
, the contact process and the heat removed via the thermostat reach a  quasi-stationary state, as shown in Fig. 13-top. This holds for all studied gap heights for both the methane and the decane simulations.  Fig. 13-bottom shows that 80% ... 95% of the work done by the indenter dissipates into heat that is then removed by the thermostat. In the BL regime, the plastic deformation of the substrate for the formation of the groove and the chip as well as the formation of the tribofilm causes large dissipation of heat in comparison to the ML and HL regime. In the HL regime, the heat flux Thermostat Q  removing the dissipated heat from the system is about two orders of magnitude smaller compared to the BL regime. In the HL regime, the dissipated heat hardly depends on the gap height. This is as expected, as viscous friction dominates the system in that regime. The comparison of the methane and the decane cases regarding the heat removed from the system provides further interesting insights. In the BL regime, more energy dissipates in the decane cases compared to the methane cases. This is due to the fact that the imprinting of decane molecules and the formation of the amorphous tribofilm requires more energy compared to the methane case. In the ML regime, more dissipated energy is removed from the system in the methane cases compared to the decane cases. For h ≈ 2 Å, the dissipated heat in the case of methane is about twice that in the case of decane. This is due to the dominance of the adsorbed fluid layers and their break-up, which is more energy intense for the strongly structured methane particles (Fig. 10). Moreover, for the methane simulations in the thin film regime, a tribofilm formation already takes place (Fig. 6), which also requires energy. Figures 14 and 15 show the fluid temperature and pressure field, respectively, averaged over the entire quasi-stationary part of the simulations. Results for h = -6, -4, -2, 0, 4, 8 Å are shown for both methane and decane. The temperature and the pressure field were computed with respect to the initial bulk values. The temperature field shows an elliptic shape, whereas the pressure field is more complex and unsymmetric. This confirms predictions from mesoscopic modeling approaches [105,106]. The pressure profile is strongly influenced by the differences between the fluid inlet and outlet flow. Both the temperature and pressure fields reflect the shape of the indenter, i.e. the fluid flow passing through the constriction underneath the indenter. Overall, both the fluid temperature and pressure in the gap increases with decreasing gap height h. The pressure profiles are primarily influenced by the obtained normal force N F at a given gap height, as shown in Fig. 5(a). The temperature profiles are primarily influenced by the dissipation behavior of the system and the work done by the indenter and thereby by the tangential force T F at a given gap height, as shown in Figs. 5(a) and 5(c). For both the temperature and pressure fields, the center of the pattern is shifted slightly in positive y-direction, i.e. is located slightly towards the front side of the indenter.
In the HL regime, practically no temperature increase T  is observed for the fluid in the gap (Fig. 14), which is in accordance with the dissipation behavior (Fig. 13). Only a faint pressure increase p  is observed in that regime-for the methane case (Fig. 15). In the ML regime, a temperature increase of about 25 K T   is detected-which is similar for the two fluids. This is probably due to the fact that the thermodynamic properties in adsorbed fluid layers are very similar to the corresponding bulk phase behavior [107]. In the BL regime, both the temperature and the pressure in the gap are strongly increased compared to the bulk phase values, as shown in Figs. 14 and 15. Fig. 15 Pressure field of the fluid in the lubrication gap with respect to the moving coordinate system (Fig. 2). The lubrication gap was determined by the slit between the -shape surface of the indenter and the substrate. Results shown for h = 8, 4, 0, -2, -4, -6 Å (top to bottom). The left column shows the results for the methane simulations; the right column those for the decane simulations. Data was sampled during the quasi-stationary phase ( Scratch 400 ps   ). The color code is the same for all shown plots.
The pressure fields for the methane cases (Fig. 15) in the BL regime exhibit a significant pressure drop, i.e. the pressure at the rear side of the indenter (out flow of the fluid) is smaller than the bulk phase pressure. This is due to the expansion of the fluid in the outlet of the rear side of the indenter.
The comparison between the methane and the decane cases provides additional insights in the contact process. Figure 16-top and bottom shows the maximum temperature and maximum pressure increase of the fluid in the gap as a function of the gap height during   in the fluid decreases linearly with increasing gap height h. Moreover, in the BL regime, both the temperature and pressure of the fluid in the gap increase more strongly during the lateral movement in the decane case. This is in line with the findings for the contact forces (Fig. 5), the formation of the tribofilm (Fig. 12), and the dissipated energy (Fig. 13). For h = -6 Å, the temperature increases about 50 K more in contact zone in the case of decane in comparison to the methane case. Also in line with the dissipated heat removed by the thermostat, a cross-over between the methane and decane simulations is observed for max T  in the ML regime. In the HL regime, the temperature increase of the two fluids is essentially the same and hardly depends on the gap height.
The maximum pressure of the fluid in the gap max p  (Fig. 16-bottom) overall shows a similar behavior as the maximum temperature (Fig. 16-top). The observable max ( ) p h  has a large gradient in the BL regime, then flattens in the ML regime, and finally converges in the HL regime. However, some differences in the max T  and max p  behavior can be pointed out: The cross-over of the methane and decane results for max T  at h = -1 Å is in excellent agreement with the behavior of the dissipated heat removed from the system Thermostat Q  (Fig. 13). The cross-over of the methane and decane results for max p  on the other hand at h = -2 Å is in excellent agreement with the behavior of the total normal force on the indenter z F (Fig. 5). This shows that for the thermodynamic properties of the fluid, the fluid pressure is dominated by the contact pressure imposed by the indenter, whereas the fluid temperature is a result of multiple simultaneous phenomena, e.g. the distortion of the adsorption layer, the formation of the tribofilm, and the elastic and plastic deformation of the substrate.

Conclusions
The atomistic mechanisms of the three main lubrication regimes of the Stribeck curve were studied using a simplified yet representative large-sized molecular dynamics (MD) simulation setup. Thereby, the transition from boundary lubrication (BL), to mixed lubrication (ML), and hydrodynamic lubrication (HL) was elucidated. The behavior of two fluids in the system were studied: Methane (as a simple model fluid) and decane (as a model lubricant).
It is found that different mechanisms on the atomistic level dominate the different lubrication regimes. In the BL regime, the formation of a tribofilm, the squeeze-out, and the holding up of the lubrication film are crucial. In the ML regime, the behavior of the adsorption layer in the gap is crucial-in particular a possible break-up of the adsorption layer structure. In the HL regime, the influence of the fluid bulk phase properties prevail-as expected. Yet, the gap height of the transition between the BL, ML, and HL regime depends on the studied system.
The simulation scenario used in this work can be favorably used to study different lubrication regimes. For future work, it would be interesting to study the effect of different simulation parameters, e.g. to assess the influence of the periodic boundary conditions applied in x-and y-direction, which may in general affect the fluid flow and fluid pressure field. Moreover, www.Springer.com/journal/40544 | Friction it would be interesting to modify the scenario such that the normal load instead of the gap height is prescribed. Also, it would be interesting to study the relation of the formation of the tribofilm and the dislocation behavior as well as to consider the thermal effects of the workpiece, i.e. the heat transfer in the bulk material, in more detail. Yet, a larger workpiece size would have to be considered in the simulations. Moreover, it would be interesting to consider a reactive force field for the fluid to obtain more realistic predictions for the formation of the tribofilm.
Important insight was gained into the formation and influence of the tribofilm. If long chain molecules such as decane form a tribofilm, this formation requires a significant amount of energy as the surface structure of the substrate is distorted and the relatively large molecules are imprinted into the surface. Opposite to that, if small molecules such as methane form a tribofilm, this is found to require less energy as the lattice structure remains in most parts intact and the fluid particles substitute lattice atoms. Nevertheless, these energetic relations also result in the fact that the formation of a tribofilm occurs already at larger gap height for methane used as a fluid. The formation of a tribofilm with a decane fluid on the other hand requires larger forces, contact pressure, etc.