A Molecular Dynamics Study on the Tribological Performance of Imidazolium − Based Ionic Liquids Mixed with Oil in Comparison to Pure Liquids

: The purpose of this work is to propose an advanced lubricant model of ILs used as additives to conventional oil. All-atoms molecular dynamics simulations are used to investigate the structure and tribology of oxidatively stable pure imidazolium-based ionic liquids (ILs), branched alkane low friction oil, and a mixture of ILs and oil conﬁned between iron surfaces. Equilibrium and shear simulations are performed at a temperature of 450 K and undergo different applied loads and shear velocities to mimic engine operations. Density proﬁles reveal the formation of layered structures at the interface. The intensity and number of the density peaks vary according to the composition of the system and the applied pressure. Velocity proﬁles reveal the presence of no-slip conditions in the pure ILs system and very high slip for the oil. The presence of a stable IL layer at the surface of the mixed lubricant fully reduces the slip of oil. Overall, the mixture displays lower friction in comparison to pure ILs. The formed corrosion protective anion layer on the metal surface makes the mixture a potential candidate for a new generation of high-performance lubricants.


Introduction
Lubricants are materials that extend the service life of machinery by reducing the friction of rolling and sliding parts [1].The ultimate objective of the science of tribology is to design high-performance, multi-functional, environmentally friendly, low-cost, and easy-to-process lubricants [2].Conventional lubricants such as natural oils already possess such features as low friction coefficients (<0.15) and non-toxicity that is typically not accessible by other one-component lubricating fluids [3,4].On the contrary, they lack thermal and oxidation stability in addition to solidification and loss of oil under high pressures [5].It is a significant disadvantage as the oils are subject to extreme normal and shear loads in industrial settings (i.e., metal forming and ball-bearing) [6][7][8][9][10].Hydrocarbons or n-alkanes under confined conditions experience an interfacial phenomenon where a solid-like behavior of the confined liquid takes place [11,12].The consequences of this behavior are the ordering of molecules at the solid surface, density oscillations within the oil film, the increase in shear viscosity, and velocity slip [13].It is especially critical when the gap between the confining slabs approaches nanometer thickness under high pressures and shear rates [14][15][16][17][18].The surface roughness [12,13] and the interfacial layer adsorption [19,20] are shown to affect the properties of the lubricant as well.
Ionic liquids (ILs) are molten salts known for their low vapor pressure and promising thermo-rheological behavior.The useful properties of ILs have attracted a substantial amount of attention to be used in various applications.Specifically, ILs have been proposed to become potential candidates for energy conservation solutions, designer solvents, and a new generation of lubricants [21].The latter is of interest for this study of enhancing the tribological performance of conventional oils.The tribological research on ILs often emphasizes the replacement of conventional lubricants.For example, ILs improved friction reduction and antiwear capabilities in comparison to frequently used polyalphaolefins, phosphazene, polyol esters, and perfluoropolyether lubricants [22].Moreover, ILs outperformed many gears, engines, and commercial hydraulic oils [22].ILs demonstrate no significant difference in performance under no pressure and high pressure conditions in steel-steel contacts, showing the ability to withstand extreme loads [23].It is found that the superior tribological performance lies in the deposition of ions on the solid sliding parts.This interfacial confined film does not allow direct contact with metal surfaces and lowers their friction.
The performance of ILs as lubricants relies on the stable film formation mechanism between the interacting surfaces.It has been observed that cations and anions under confinement orient in particular patterns [24,25].Similar to ordered layers in graphite, ILs nanostructures are characterized by a low interlayer strength, making them ideal lubricant candidates [26].In addition, Xiao et al. [27] carried out experiments at very high pressures (up to 3 GPa), comparing the film-forming ability of imidazolium-based ILs against that of conventional silicon oil.It was shown that when the viscosity was equivalent in both systems, ILs formed a thicker film than oil and performed better overall under high loads.Recent molecular dynamics simulations studies showed that ILs tend to bind strongly on the iron surface and can be used as corrosion inhibitors or metal decontamination substrates [28].The intensity of these effects can be tuned by controlling the length of the cation alkyl chain [28].Likewise, ILs films on diamond-like carbon surfaces can reduce friction through the restriction of C-C bond formation between the sliding interfaces [29].
Conventional oils are limited in high-temperature environments due to their tendency to decompose.The thermal stability of ILs thus offers another advantage to available lubricants as they can withstand temperatures up to 670 K without undergoing degradation [30].In fact, some imidazolium ILs display heat resistance to temperatures up to 750 K while maintaining their characteristic of low-temperature fluidity (glass transition temperature , of ILs can vary from 157 K to 433 K) [31,32].This allows ILs to operate in a wide range of temperatures.
Properties such as low volatility and non-flammability are crucial for the design of lubricants, as they can minimize evaporation and fire hazards while transporting and handling the material.In addition, the thermo-oxidative and chemical stability of ILs prevents corrosive processes from taking place, thus protecting and extending the lifetime of machinery [33,34].The low toxicity and high biodegradability of many ILs are also advantageous from an industrial perspective, as they prevent hazardous risks for both workers and the environment [35].Finally, as ILs are composed of at least two ions, a vast database of ILs (>10 6 ) can be explored to design lubricants that can be tuned and adapted according to different sliding contacts or specific requirements of the machinery [36,37].In many research papers, neat ILs are explored to become a new generation of lubricants [26,36,[38][39][40][41].Such studies allow applying all the exceptional physicochemical properties of ionic liquids to lubrication.However, IL chemistry tends to be complex and expensive, leading to financial challenges in industrial settings.Therefore, there are studies done to discover the possibility of using ILs as additives to conventional lubricants [26,36,[38][39][40][42][43][44][45][46][47][48][49].
For example, 1% of [N 12,H,H,H ][Cl] was added to oil to lubricate steel-aluminum contact surfaces.[P 6,6,6,14 ][DEHP] ILs additives achieved steady lubrication and reduction of friction at a high temperature of 373 K.At the same time, a common antiwear agent zinc dialkyl dithiophosphate (ZDDP) only showed good performance under ambient conditions [50].In a steel-steel contact experiment, a substantial wear decrease was seen after adding 1-5 wt% of bisimidazolium IL to poly(ethylene glycol) (PEG) [51].Overall, it was proved that a small amount of IL in a base oil can greatly enhance the tribological performance of the mixture in comparison to neat oil [43,52].Moreover, experimental results concluded that in many cases such mixed lubricants can outperform standalone ILs in wear reduction.In this work, we use atomistically detailed molecular dynamics simulations to predict the structural and tribological properties of imidazolium-based ILs with the linear and branched alkyl chain of cation (1-heptyl-3-methylimidazolium and 1-(5-methylhexyl)-3methylimidazolium) and bis[(trifluoroethane)sulfonyl]imide anion, oil (9-octylheptadecane), and IL-based lubricants under metal confinement.The bis[(trifluoroethane)sulfonyl]imide anion displays superior oxidative stability to other anions [53].Thus, it was selected for our study.Two types of cations with long alkyl chains (branched and linear) are chosen for this study because the length of the alkyl chains of cations is a controlling factor for tribological performance.It has been found that long side chains show improved friction and wear performance [54][55][56][57].The n-alkanes with long molecular chain lengths demonstrate the same advances [16].In addition, the branched n-alkanes yield an increase in shear stresses and stay in the liquid phase under higher applied pressures than linear molecules in confinement [19].That is why 9-octylheptadecane is given a preference to represent a low friction model oil in this research.Lastly, the most common rolling and sliding materials in industry are steel and iron, so we chose iron to confine the liquids [58,59].
Our previous study showed that ions create layered structures near the solid-liquid interface, and these layers are strongly adsorbed to the iron surface [60].Here, we characterize the structure from density profile data after equilibrium simulations to get an insight into molecular ordering under extreme confinement.We then analyze the flow behavior and the stress response of the liquids after shear flow simulations with various flow strengths.To understand the efficiency of the designed lubricant, we calculate and compare the friction coefficient for all the studied systems.Our goal is to keep the low friction coefficient that hydrocarbon lubricants offer while maintaining the stability of the liquid layer on the surface to prevent the metal surface from degradation.Thus, ILs will have enhanced tribological properties for lubrication when used as additives to conventional oils.
This paper is constructed as follows: Section 2 is dedicated to the discussion of detailed simulation methodology.Section 3 uncovers the interfacial and tribological phenomena of the resulting lubricant compared to other studied liquids under different confinement and shear conditions.In Section 4, concluding remarks on the subject and future efforts are presented.

Simulations Preparation and Force Field
The ILs chosen for this study consist of two imidazolium-based cations with long linear and branched alkyl chains (1-heptyl-3-methylimidazolium [C 7 C 1 Im] + and 1-(5methylhexyl)-3-methylimidazolium [5mC 6 C 1 Im] + ), each coupled with an oxidatively stable anion (bis[(trifluoroethane)sulfonyl]imide [NTf 2 ] − ).As a model oil with a low friction coefficient, a branched alkane (9-octylheptadecane [C 25 H 52 ]) is used to simulate pure conventional fuel behavior.The chemical structures of the molecules studied are shown in Figure 1.The liquids are confined between two identical surfaces, each made of 6480 α-iron (Fe) atoms [61].The area of the solid-liquid interface is 5325 Å, and the height of the metal slab is 14 Å (roughly the thickness of 10 atomic layers).Four different confined systems are studied: IL with linear alkyl chain cation (linear IL), IL with branched alkyl cation (branched IL), 9-octylheptadecane, and a mixture of linear IL and oil (i.e., lubricant).Both the linear and branched ILs-Fe systems size is 51,180 atoms while the oil-Fe is 63,472 atoms.The target fraction of ILs in oil is calculated based on the number of ion pairs required to fully cover the metal surfaces.The mixed system is created by placing the IL with linear cation directly on the upper and lower interfaces and filling up the bulk with oil molecules (see Figure 2c).Here, the system size is 65,376 atoms, out of which 11,760 are IL atoms.
required to fully cover the metal surfaces.The mixed system is created by placing t with linear cation directly on the upper and lower interfaces and filling up the bulk oil molecules (see Figure 2c).Here, the system size is 65,376 atoms, out of which 1 are IL atoms.The general AMBER force field is used as a model for the inter-and intramole interactions present in the liquids [62][63][64].Optimization of isolated ion structures of performed using the B3LYP/6-311++g(d,p) level through Gaussian 09 software [65 RESP method is used to assign partial charges with a rescaling factor of 0.8 to accou the polarizability and charge transfer presence in ions [66].Such rescaling has been s to predict ILs bulk properties that agree with experimental observations [6 Additionally, the bulk properties of these ionic liquids were explored in previous res [69].For the 9-octylheptadecane, the atom charges were assigned using the AM1 method [70,71].Intermolecular potentials for ILs and 9-octylheptadecane inter systems are calculated using the 12-6 Lennard-Jones (LJ) potential.The Fe-Fe interac are simulated using the embedded atom method (EAM) [72,73], while cross interac between atoms of liquids and iron are represented by LJ.The remaining interaction computed by mixing rules [74,75].The exact LJ cross parameters used in this stud listed in our previous work [60,61,76,77].Both Coulombic and LJ interaction computed up to a cutoff distance of 12 Å.Long-range interactions are approximated the particle-particle particle-mesh (PPPM) algorithm and tail correction [78].Pe boundary conditions are applied in the x and y directions, while the boundary cond required to fully cover the metal surfaces.The mixed system is created by placing the IL with linear cation directly on the upper and lower interfaces and filling up the bulk with oil molecules (see Figure 2c).Here, the system size is 65,376 atoms, out of which 11,760 are IL atoms.The general AMBER force field is used as a model for the inter-and intramolecular interactions present in the liquids [62][63][64].Optimization of isolated ion structures of ILs is performed using the B3LYP/6-311++g(d,p) level through Gaussian 09 software [65].The RESP method is used to assign partial charges with a rescaling factor of 0.8 to account for the polarizability and charge transfer presence in ions [66].Such rescaling has been shown to predict ILs bulk properties that agree with experimental observations [67,68].Additionally, the bulk properties of these ionic liquids were explored in previous research [69].For the 9-octylheptadecane, the atom charges were assigned using the AM1-BCC method [70,71].Intermolecular potentials for ILs and 9-octylheptadecane interfacial systems are calculated using the 12-6 Lennard-Jones (LJ) potential.The Fe-Fe interactions are simulated using the embedded atom method (EAM) [72,73], while cross interactions between atoms of liquids and iron are represented by LJ.The remaining interactions are computed by mixing rules [74,75].The exact LJ cross parameters used in this study are listed in our previous work [60,61,76,77].Both Coulombic and LJ interactions are computed up to a cutoff distance of 12 Å.Long-range interactions are approximated using the particle-particle particle-mesh (PPPM) algorithm and tail correction [78].Periodic boundary conditions are applied in the x and y directions, while the boundary condition The general AMBER force field is used as a model for the inter-and intramolecular interactions present in the liquids [62][63][64].Optimization of isolated ion structures of ILs is performed using the B3LYP/6-311++g(d,p) level through Gaussian 09 software [65].The RESP method is used to assign partial charges with a rescaling factor of 0.8 to account for the polarizability and charge transfer presence in ions [66].Such rescaling has been shown to predict ILs bulk properties that agree with experimental observations [67,68].Additionally, the bulk properties of these ionic liquids were explored in previous research [69].For the 9-octylheptadecane, the atom charges were assigned using the AM1-BCC method [70,71].Intermolecular potentials for ILs and 9-octylheptadecane interfacial systems are calculated using the 12-6 Lennard-Jones (LJ) potential.The Fe-Fe interactions are simulated using the embedded atom method (EAM) [72,73], while cross interactions between atoms of liquids and iron are represented by LJ.The remaining interactions are computed by mixing rules [74,75].The exact LJ cross parameters used in this study are listed in our previous work [60,61,76,77].Both Coulombic and LJ interactions are computed up to a cutoff distance of 12 Å.Long-range interactions are approximated using the particle-particle particle-mesh (PPPM) algorithm and tail correction [78].Periodic boundary conditions are applied in the x and y directions, while the boundary condition is fixed in the direction of confinement (z-axis).However, to properly account for electrostatic interactions in the z-direction, a volume equivalent to three empty simulation boxes is introduced between the metal slabs, thus eliminating potential slab-slab interactions (e.g., slab correction implementation).All of the simulations are run in LAMMPS using a time step of 1 fs [79].
To confine the liquids between the iron slabs and introduce various degrees of confinement, constant velocities are symmetrically applied towards the center of the simulation boxes to both upper and lower surfaces.In this way, the systems are uniformly compressed.
It is important to note that the velocity for the compression has to be relatively slow to allow the molecules to have enough time to evolve and orient accordingly under the applied pressure.
To study the equilibrium dynamics of the systems, microcanonical ensemble (NVE) simulations are performed at a temperature of 450 K using a Langevin thermostat with a damping factor of 100 [80,81].At this stage of system preparation, the iron atoms do not undergo any dynamics.The degree of confinement is controlled by the final gap size between the metal slabs, which can be varied with simulation run time: the longer the simulation run, the more confinement (smaller gap size) is introduced.Once the compression is complete, the systems are equilibrated in the same Langevin thermostat (T = 450 K) and NVE conditions.The normal force exerted by the liquids on both top and bottom solid surfaces is computed and recorded every 50 ps averaged over each time step.The simulation runtime was in the range of 3-10 ns.The end time was chosen as the time at which the recorded force reaches a steady state and equal values at both slabs, plus an additional 1 ns to ensure equilibration.Typically, the more pressure that was applied, the longer it took for the system to relax.The values of the pressure for all four systems were around 0.2 GPa, 0.5 GPa, and 1 GPa.All the experimental conditions are listed in Table 1.Density profiles along the z-direction are computed by averaging density data over the last 1 ns of the equilibrium simulation.The final relaxed structures are presented in Figure 2.

Shear Simulations Details
After the systems were relaxed, shear simulations were performed by sliding the explicit walls.The walls are set to move in opposite directions along the x-axis with the constant velocity V/2 (refer to Figure 3).The shear velocities chosen are V 1 = 10 m/s, V 2 = 20 m/s, V 3 = 50 m/s, and V 4 = 100 m/s.The various degrees of confinement are maintained throughout the simulation in the systems by applying normal force (or pressure P) on the metal slabs.This is accomplished by dividing the force (calculated from the equilibrium simulations) by the number of atoms (F atom ) in the slab and adding the resulting value uniformly on each Fe atom of the bottom slab.Similarly, -F atom is the applied force at the top slab.For these simulations, only the metal slabs are kept at a constant temperature of 450 K using a Langevin thermostat.Non-equilibrium molecular dynamics are performed by solving Newton's equations of motion in NVE ensemble to shear the fluids: . r i = p i /m i and , where r i is the position of particle i, p i is thermal momentum, m i is mass of the particle, and U is system potential.The periodic boundary conditions are applied in x and y directions, and z is the direction of confinement.The liquids dynamics are computed using an NVE ensemble allowing the temperature of the liquid to evolve freely.In these simulations, the momentum is transferred from the walls into the fluid phase according to the established force field to determine the boundary behavior of the velocity of fluids, and there is no artificial coupling between the velocity of the walls and the confined fluid.The simulation run time is set to 5 ns for all of the systems.Velocity profiles and density profiles are also recorded during the simulations.Data are collected at analyzed after the system reaches a steady state (typically the last 4 ns of the run).Lastly, the time averages of shear force (F x in the x-direction) and normal force (F y in Fluids 2022, 7, 384 6 of 18 the y-direction) of the molecules exerting on the surface are computed and recorded every 100 ps.The schematic of the shear simulation is shown in Figure 3.
transferred from the walls into the fluid phase according to the established force field to determine the boundary behavior of the velocity of fluids, and there is no artificial coupling between the velocity of the walls and the confined fluid.The simulation run time is set to 5 ns for all of the systems.Velocity profiles and density profiles are also recorded during the simulations.Data are collected at analyzed after the system reaches a steady state (typically the last 4 ns of the run).Lastly, the time averages of shear force (Fx in the x-direction) and normal force (Fy in the y-direction) of the molecules exerting on the surface are computed and recorded every 100 ps.The schematic of the shear simulation is shown in Figure 3.All simulations are run under an elevated temperature of 450 K, at a range of external pressures applied to the confined liquids, and at various shear rates to simulate realistic and extreme operational conditions [82-84].

Results and Discussion
In this section, we report the results and carry on discussions about the effect of various degrees of confinement on the molecular structuring, behavior of liquids under shear flow of different strengths, and tribological properties of two types of imidazoliumbased ionic liquid [C7C1Im] + [NTf2] − (linear IL) and [5mC6C1Im] + [NTf2] − (branched IL) as well as branched alkane oil [C25H52] and a mixture of linear IL and oil.The preformed equilibrium simulations reveal complex restructuring of the IL molecules at the metal surface that was previously identified as the main contributor to the stability of the interfacial layer [60].The shear simulations provide insights into tribological capabilities All simulations are run under an elevated temperature of 450 K, at a range of external pressures applied to the confined liquids, and at various shear rates to simulate realistic and extreme operational conditions [82-84].

Results and Discussion
In this section, we report the results and carry on discussions about the effect of various degrees of confinement on the molecular structuring, behavior of liquids under shear flow of different strengths, and tribological properties of two types of imidazolium-based ionic liquid [C 7 C 1 Im] + [NTf 2 ] − (linear IL) and [5mC 6 C 1 Im] + [NTf 2 ] − (branched IL) as well as branched alkane oil [C 25 H 52 ] and a mixture of linear IL and oil.The preformed equilibrium simulations reveal complex restructuring of the IL molecules at the metal surface that was previously identified as the main contributor to the stability of the interfacial layer [60].The shear simulations provide insights into tribological capabilities of the proposed IL-based lubricant to enable the inverse design of such lubricating materials.

Density Profiles
In Figure 4, the density profiles along confinement direction z for the unpressurized case are shown for all four liquids.The introduction of the metal surfaces with no external forces acting on the system gives rise to multilayer formation at both lower and upper walls.This is in agreement with IL behavior reported in previous studies [60].This type of restructuring is due to the interaction between the liquid and solid atoms of the surface [11].In the case of linear and branched ILs (Figure 4a,b), the first interfacial layer is mainly populated by anions, but some cations are present at the surface as well.The second layer mainly consists of cations.Due to the symmetric nature of the [NTf 2 ] − anion, it arranges tightly on the surface.forces acting on the system gives rise to multilayer formation at both lower and u walls.This is in agreement with IL behavior reported in previous studies [60].This of restructuring is due to the interaction between the liquid and solid atoms of the su [11].In the case of linear and branched ILs (Figure 4a,b), the first interfacial layer is ma populated by anions, but some cations are present at the surface as well.The second l mainly consists of cations.Due to the symmetric nature of the [NTf2] − anion, it arra tightly on the surface.The density profile of the oil show at least three distinguishable layers at the surface after which no additional fluctuations are observed.The density profile o mixture of linear IL and oil (Figure 4d) displays two high intensity peaks for anion cation and no vivid peaks at the IL-oil interface.Moreover, the anion's first density has a higher intensity than in the pure IL system.Here, higher ordering is anticip because ions experience confinement not only from the metal slab but also from the liq liquid interface formed with oil.There is a small overlap between liquid-liquid ph but no significant penetration of ions into the oil phase is observed.The phase separa of imidazolium-based ILs and model oil have been previously reported [85].Lastly, to their distance from the surface, no interfacial phenomena are detected in th molecules.The density profile of the oil show at least three distinguishable layers at the iron surface after which no additional fluctuations are observed.The density profile of the mixture of linear IL and oil (Figure 4d) displays two high intensity peaks for anion and cation and no vivid peaks at the IL-oil interface.Moreover, the anion's first density peak has a higher intensity than in the pure IL system.Here, higher ordering is anticipated because ions experience confinement not only from the metal slab but also from the liquid-liquid interface formed with oil.There is a small overlap between liquid-liquid phases, but no significant penetration of ions into the oil phase is observed.The phase separation of imidazolium-based ILs and model oil have been previously reported [85].Lastly, due to their distance from the surface, no interfacial phenomena are detected in the oil molecules.
Interestingly, oil densities in the bulk region in pure (ρ = 0.670 ± 0.003 g/cm 3 ) and lubricant (ρ = 0.648 ± 0.007 g/cm 3 ) cases differ by 3%.This can be due to the discrepancy in the load distribution of the more densely packed ions on the surface than 9-octylheptadecane molecules; hence the oil experiences less pressure.All the densities in bulk at zero pressure and at high pressures are relative to what can be found in literature after extrapolation to the elevated temperature of 450 K (see Table 2 for bulk density values under high degrees of confinement) [69,86,87].

System
ρ (g/cm 3 ) at P 1 ρ (g/cm 3 ) at P 2 ρ (g/cm 3 ) at P 3 At higher pressures, the intensity of density and its behavior change.Both linear (Figure 5a,d,g) and branched (Figure 5b,e,h) ILs show a closely packed first interfacial layer.The density values in all the additional layers grow with increasing degrees of confinement.Generally, in comparison to linear IL, branched IL systems exhibit higher values of density in the initial peak and lower values for any subsequent peak.Most importantly, with the increase of external load on the IL systems, the first density peak of the cation shifts outside of the first anion peak.It means that the layer at the solid interface becomes mostly populated by anions and the next layer becomes mostly populated by cations.This anioncation structure is called a double layer (or electrical double layer) and it usually occurs at charged surfaces [88][89][90].In our case, the asymmetrical bulky cations get "squeezed out" due to the entropic effect.The anion of choice for this study ([NTf 2 ] − ) is known to show low corrosiveness and great thermo-oxidative stability on metal surfaces [34].Additionally, these anions form a stable layer on the iron surface, hence the anticorrosive protective layer is formed on the surface.The bulk density increases with the load for all systems as can be seen in Table 2 and Figure 5.It also reveals higher amplitude density oscillations than in unpressurized systems.In the case of the highest applied pressure on the oil system (Figure 5i), the amplitud of density oscillations in the bulk region is four times higher than for the system with n external pressure on the metal slabs.This suggests that under high normal load, o clusters are present not only at the interface but also in bulk, indicating that the oil star to solidify.This phenomenon is frequently reported in experiments and modeling studi of the tribology of hydrocarbons [91,92].Subjecting such highly confined liquid to she can lead to dry friction due to separation and transverse flow, which can become significant drawback of oil lubricants [5].
The density fluctuations of the oil bulk in the unpressurized system (Figure 4c) ca be attributed to the thermal fluctuations, which are within the statistical error (~1%).the equilibrium simulations, the temperature of the liquid is controlled.As pressu increases, the oil system starts to experience an interesting density behavior.When n load is applied, the first density peak is observed at ρ = 1.7 g/cm 3 and at least thr additional layers are seen (Figure 4c).At P1 = 0.24 GPa (Figure 5c), the first peak recorded at ρ = 2.4 g/cm 3 and an additional three or four layers are observed.At P2 = 0.6  and (c,f,i) 9-octylheptadecane with external pressure of (a-c) P 1 ≈ 0.2 GPa, (d-f) P 2 ≈ 0.5 GPa, and (g-i) P 3 ≈ 1 GPa applied to the walls.The number of layers increases as more pressure is applied to the systems.
In the case of the highest applied pressure on the oil system (Figure 5i), the amplitude of density oscillations in the bulk region is four times higher than for the system with no external pressure on the metal slabs.This suggests that under high normal load, oil clusters are present not only at the interface but also in bulk, indicating that the oil starts to solidify.This phenomenon is frequently reported in experiments and modeling studies of the tribology of hydrocarbons [91,92].Subjecting such highly confined liquid to shear can lead to dry friction due to separation and transverse flow, which can become a significant drawback of oil lubricants [5].
The density fluctuations of the oil bulk in the unpressurized system (Figure 4c) can be attributed to the thermal fluctuations, which are within the statistical error (~1%).In the Fluids 2022, 7, 384 9 of 18 equilibrium simulations, the temperature of the liquid is controlled.As pressure increases, the oil system starts to experience an interesting density behavior.When no load is applied, the first density peak is observed at ρ = 1.7 g/cm 3 and at least three additional layers are seen (Figure 4c).At P 1 = 0.24 GPa (Figure 5c), the first peak is recorded at ρ = 2.4 g/cm 3 and an additional three or four layers are observed.At P 2 = 0.62, the first peak is recorded at ρ = 2.2 g/cm 3 and six additional layers are present.Finally, at P 3 = 1.11GPa (Figure 5i), the density at the surface reaches 2.5 g/cm 3 , and the layers propagate throughout the whole system.This behavior implies that the liquid prefers to rearrange into more ordered clusters over increasing the packing density at the interface.
Figure 6 shows the density profile of the lubricant system (IL linear and oil).Two prominent anion peaks as well as two peaks and a shoulder peak are present for the cation case.In comparison with the unpressurized case, the system with the highest degree of confinement has a more pronounced structuring between cation and oil molecules.The oil also experiences an increase in density at the liquid-liquid interface with the applied load.The number of IL layers does not change with the applied pressure.Like the pure system, the bulk density of oil at high pressure shows the propagation of layering throughout the bulk (Figure 6d).system, the bulk density of oil at high pressure shows the propagation of layering throughout the bulk (Figure 6d).

Velocity Profiles
Velocity profiles (Vx), as a function of the z-coordinate of linear and branched ILs, are plotted in Figure 7a-f.Four different velocities are applied to the metal walls to access the liquid response over a range of shear rates.The data reveal no-slip boundary conditions at the solid-liquid interface for ILs.However, even though the velocity profiles seem to follow linear Couette flow behavior, this is only true for the bulk of the liquids.The phenomenon is called central localization (a specific case of shear banding) when interfacial layers of liquid move with the solid walls as if these layers are a continuation of it [93][94][95][96].Hence, the no-slip condition is extended to the ordered ion structures present on the surface.This phenomenon of central localization is observed due to the high stability of the IL layer on the surface, even under nano-confinement and high share rate conditions [97][98][99][100].

Velocity Profiles
Velocity profiles (V x ), as a function of the z-coordinate of linear and branched ILs, are plotted in Figure 7a-f.Four different velocities are applied to the metal walls to access the liquid response over a range of shear rates.The data reveal no-slip boundary conditions at the solid-liquid interface for ILs.However, even though the velocity profiles seem to follow linear Couette flow behavior, this is only true for the bulk of the liquids.The phenomenon is called central localization (a specific case of shear banding) when interfacial layers of liquid move with the solid walls as if these layers are a continuation of it [93][94][95][96].Hence, the no-slip condition is extended to the ordered ion structures present on the surface.This phenomenon of central localization is observed due to the high stability of the IL layer on the surface, even under nano-confinement and high share rate conditions [97][98][99][100].It is also found that the branching of the cation alkyl chain affects the "no-slip "frozen" layer thickness.In fact, for linear IL the interfacial layer is 5-6 Å thick, comp to the 8 Å layer in the branched IL at P1.However, at the highest pressure P3, the laye thickness of the branched IL is reduced to 6 Å.This difference can be due to the pac size of the branched molecules being slightly larger than the packing size of linear o In other words, the radius of gyration of branched cations is bigger than the radiu gyration of linear cations until the pressure is high enough to compress both typ molecules and reduce the difference in size.Furthermore, these size scales are relati the thickness of ionic liquids' double layers and consistent with other reported data The oil system (Figure 7g,h) shows a large slip, as expected [101].Other observations these shear simulations propose that the bulk density fluctuations smoothen out at shear rates due to the active rearrangement of molecules by the induced flow.T dynamics contribute to the reduction of crystallization of the liquids in bulk u confinement.The density and the number of layers at the surface stay the same equilibrium simulations.
For the lubricant system, due to the high free energy of adsorption of ILs on the surface, the slip of liquid is reduced regardless of the applied load (see Figure 8).phenomenon has relevant industrial applications as it allows for the decreasing of w in the machinery by preventing dry lubrication and the splitting of the oil, and s down the oxidation of the metal surfaces.It was also observed that single cations from the IL cluster into the oil phase when the highest degrees of confinement and la wall velocities were applied.Yet, the first interfacial layer always remains intact.It is also found that the branching of the cation alkyl chain affects the "no-slip" or "frozen" layer thickness.In fact, for linear IL the interfacial layer is 5-6 Å thick, compared to the 8 Å layer in the branched IL at P 1 .However, at the highest pressure P 3 , the layering thickness of the branched IL is reduced to 6 Å.This difference can be due to the packing size of the branched molecules being slightly larger than the packing size of linear ones.In other words, the radius of gyration of branched cations is bigger than the radius of gyration of linear cations until the pressure is high enough to compress both types of molecules and reduce the difference in size.Furthermore, these size scales are relative to the thickness of ionic liquids' double layers and consistent with other reported data [95].The oil system (Figure 7g,h) shows a large slip, as expected [101].Other observations from these shear simulations propose that the bulk density fluctuations smoothen out at high shear rates due to the active rearrangement of molecules by the induced flow.These dynamics contribute to the reduction of crystallization of the liquids in bulk under confinement.The density and the number of layers at the surface stay the same as in equilibrium simulations.
For the lubricant system, due to the high free energy of adsorption of ILs on the iron surface, the slip of liquid is reduced regardless of the applied load (see Figure 8).This phenomenon has relevant industrial applications as it allows for the decreasing of wear in the machinery by preventing dry lubrication and the splitting of the oil, and slows down the oxidation of the metal surfaces.It was also observed that single cations split from the IL cluster into the oil phase when the highest degrees of confinement and largest wall velocities were applied.Yet, the first interfacial layer always remains intact.

Tribological Properties
The next step is to verify whether using ILs additives to conventional lubricants lead to a friction reduction in the system.In Figure 9a,b, as anticipated, the normal fo (Fn) exerted by the walls on the liquid depend on the pressure applied and independent of the velocity.On the other hand, the shear forces (Fs) for ILs and lubri systems (Figure 9c,d,g,h) show dependency on both velocity and pressure.The dat oil are highly noisy, and after averaging it over time, a weak velocity dependence strong pressure dependence is detected.Both normal and shear forces for all syst reach constant values over the simulation run time.

Tribological Properties
The next step is to verify whether using ILs additives to conventional lubricants can lead to a friction reduction in the system.In Figure 9a,b, as anticipated, the normal forces (F n ) exerted by the walls on the liquid depend on the pressure applied and are independent of the velocity.On the other hand, the shear forces (F s ) for ILs and lubricant systems (Figure 9c,d,g,h) show dependency on both velocity and pressure.The data for oil are highly noisy, and after averaging it over time, a weak velocity dependence but strong pressure dependence is detected.Both normal and shear forces for all systems reach constant values over the simulation run time.
Friction coefficients are calculated using µ = F s F n and plotted in Figure 10, where they are presented as a function of pressure (Figure 10a-d) and shear velocity (Figure 10e-h).At low values of pressure, the difference between friction coefficients in the range of shear velocities is more prominent.At higher pressures, almost no variance is observed due to the structuring of the molecules in the liquids.This is consistent with the density and velocity profile data.An increase in the friction coefficient with increasing pressure in a low velocity regime is identified for pure ILs and mixture systems.The opposite happens for a high sliding velocity regime (20-50 m/s), where the friction coefficient decreases with applied load.A general trend of friction coefficient increase with high shearing rates is seen for all the systems.The 9-octylheptadecane as low friction model oil shows the lowest friction coefficients at P 1 = 0.24 GPa and P 2 = 0.62 GPa.A jump in friction coefficient values occurs at the highest loads and velocity due to the oil solidification.Friction coefficients are calculated using  and plotted in Figure 10, where they are presented as a function of pressure (Figure 10a-d) and shear velocity (Figure 10e-h).At low values of pressure, the difference between friction coefficients in the range of shear velocities is more prominent.At higher pressures, almost no variance is observed due to  The mixed IL and oil system shows lower friction coefficient values at all applied loads than that of pure ILs (~0.12 and below for lower shear rates-see Table 3).To the best of our knowledge, this low friction coefficient value has yet to be achieved by pure ILs lubricants in contact with neutral metals at such high pressures and temperatures.The

Conclusions
To better elucidate the lubrication properties of ILs and their mixture with conventional oil under metal confinement, all-atomistic molecular dynamics simulations are performed in this study.A series of microstructural and tribological predictions on pure ILs (1-heptyl-3-methylimidazolium and 1-(5-methylhexyl)3-ethylimidazolium cations, and bis[(trifluoroethane)sulfonyl]imide anion), pure oil (9-octylheptade-cane), and IL-oil mixture are presented.Both equilibrium and shear studies under different applied loads are performed.For all of the systems, the density profiles under equilibrium conditions and no external load reveal the formation of a layered structure at the solid-liquid interface.When external loads are applied to the system, the formation of a double layer of ILs is observed.In addition, the density values of all the layers increase with the degrees of confinement.At high applied pressures, the oil solidified.The mixture of IL and oil shows two layers of anion and cation and a low intensity additional density peak indicated a weak interface between the cations and oil molecules.For all IL-containing systems under high applied pressure, the interfacial layer primarily consists of anions.The formed anion layer protects the iron surface from corrosion due to the oxidatively stable nature of the chosen anion.
Velocity profiles reveal the presence of no-slip boundary conditions at the solid-liquid interface of ILs and through the first double layer of ions.The linear velocity profile in the bulk of the liquids indicates central localization of shear in the ILs systems.As expected, oil exhibited a large slip.When IL is added to the oil, the high free energy of adsorption of ILs on the iron surface helps to reduce the slip in the system.Normal and shear forces are also calculated and used to obtain the friction coefficient of the various configurations studied.Pure ILs systems displayed the highest friction coefficient, making them unideal lubrication candidates.However, when used as additives to conventional oil, a decrease in the friction coefficient is observed in comparison to pure ILs.In addition to the anion protective layer on the metal surface, these findings suggest that ILs as additives to conventional oils are better candidates for the next generation of high-performance and noncorrosive lubricants than pure ionic liquids.
In this work, we explored ILs and IL-based lubricants at elevated temperature conditions.Meanwhile, low-temperature lubrication is just as important in machinery [102].ILs are expected to perform well due to their excellent thermal stability in both high and lowtemperature regimes.Nevertheless, investigating the behavior of ILs at low temperatures is still crucial for understanding the full capability of the proposed lubricant.Moreover, surface irregularity results in a reduction of the ordering of the lubricant at the surface and might increase the overall friction and wear in lubricating systems [12,103,104].The Gibbs free energy of the adsorption of molecules varies with the surface selection [77].Henceforth, surface altercations such as chemical composition and nano roughness can influence the properties of a lubricant and should be studied to give a complete picture of IL tribological performance.

Figure 2 .
Figure 2. Simulation snapshots of (a) IL with linear alkyl chain cation, (b) 9-octylheptadecan (c) mixed linear IL and oil systems after relaxation of the liquids under confinement (P1 ≈ 0.2 The linear IL and branched IL simulation boxes look identical.

Figure 2 .
Figure 2. Simulation snapshots of (a) IL with linear alkyl chain cation, (b) 9-octylheptadecane, and (c) mixed linear IL and oil systems after relaxation of the liquids under confinement (P1 ≈ 0.2 GPa).The linear IL and branched IL simulation boxes look identical.

Figure 2 .
Figure 2. Simulation snapshots of (a) IL with linear alkyl chain cation, (b) 9-octylheptadecane, and (c) mixed linear IL and oil systems after relaxation of the liquids under confinement (P 1 ≈ 0.2 GPa).The linear IL and branched IL simulation boxes look identical.

Figure 3 .
Figure 3. Schematic of shear simulation for linear IL system.

Figure 3 .
Figure 3. Schematic of shear simulation for linear IL system.

Figure 4 .
Figure 4. Density profiles along confinement direction (z) for (a) linear IL, (b) branched IL, octylheptadecane, and (d) lubricant with no external pressure applied (P0 = 0 GPa) to the w Density profiles showed a layered structure for all systems.

Figure 4 .
Figure 4. Density profiles along confinement direction (z) for (a) linear IL, (b) branched IL, (c) 9-octylheptadecane, and (d) lubricant with no external pressure applied (P 0 = 0 GPa) to the walls.Density profiles showed a layered structure for all systems.

Fluids 2022, 7 ,Figure 5 .
Figure 5. Density profiles along confinement direction (z) for (a,d,g) linear IL, (b,e,h) branched I and (c,f,i) 9-octylheptadecane with external pressure of (a-c) P1 ≈ 0.2 GPa,(d-f) P2 ≈ 0.5 GPa, an (g-i) P3 ≈ 1 GPa applied to the walls.The number of layers increases as more pressure is applied the systems.

Figure 5 .
Figure 5. Density profiles along confinement direction (z) for (a,d,g) linear IL, (b,e,h) branched IL,and (c,f,i) 9-octylheptadecane with external pressure of (a-c) P 1 ≈ 0.2 GPa, (d-f) P 2 ≈ 0.5 GPa, and (g-i) P 3 ≈ 1 GPa applied to the walls.The number of layers increases as more pressure is applied to the systems.

Figure 6 .
Figure 6.(a) Simulation snapshot of the lubricant.Density profiles along confinement direction (z) for the lubricant system with applied pressure of (b) P1 = 0.23 GPa, (c) P2 = 0.47 GPa, and (d) P3 = 1.13 GPa to the walls.At high loads, the interfacial layer is mostly populated by anions.

Figure 6 .
Figure 6.(a) Simulation snapshot of the lubricant.Density profiles along confinement direction (z) for the lubricant system with applied pressure of (b) P 1 = 0.23 GPa, (c) P 2 = 0.47 GPa, and (d) P 3 = 1.13 GPa to the walls.At high loads, the interfacial layer is mostly populated by anions.

Figure 7 .
Figure 7. Velocity profiles (Vx) along confinement direction (z) for (a,c,e) linear IL, (b,d,f) branched IL with applied external pressure of (a,b) P1 ≈ 0.2 GPa, (c,d) P2 ≈ 0.5 GPa, and (e,f) P GPa to the walls.The inset of the (b) branched IL velocity profile at pressure P1 is a magnific of the ionic layer with the no-slip boundary condition.The gray dashed line shows the solid-l interface, and the dotted line indicates the no-slip boundary within the ionic layer.For t octylheptadecane system, only (g) unpressurized and (h) under the pressure of P1 case presented as an example of slipping oil at the metal wall.Velocity profiles for ILs systems sho slip and central localization of shear, while oil possesses a large slip even under no external lo

Figure 7 .
Figure 7. Velocity profiles (V x ) along confinement direction (z) for (a,c,e) linear IL, (b,d,f) and branched IL with applied external pressure of (a,b) P 1 ≈ 0.2 GPa, (c,d) P 2 ≈ 0.5 GPa, and (e,f) P 3 ≈ 1 GPa to the walls.The inset of the (b) branched IL velocity profile at pressure P 1 is a magnification of the ionic layer with the no-slip boundary condition.The gray dashed line shows the solid-liquid interface, and the dotted line indicates the no-slip boundary within the ionic layer.For the 9-octylheptadecane system, only (g) unpressurized and (h) under the pressure of P 1 cases are presented as an example of slipping oil at the metal wall.Velocity profiles for ILs systems show no slip and central localization of shear, while oil possesses a large slip even under no external load.

Figure 8 .
Figure 8. Velocity profiles along confinement direction (z) for the lubricant (mixture of linear IL 9-octylheptadecane oil) with the applied external pressure of (a) P1 = 0.23 GPa, (b) P2 = 0.47 GPa (c) P3 = 1.13 GPa.The inset of the (a) at pressure P1 is a magnification of the ionic layer with th slip boundary condition.The gray dashed line shows the solid-liquid interface, and the dotted indicates the no-slip boundary within the ionic layer.The addition of IL fully reduces the slip o

Figure 8 .
Figure 8. Velocity profiles along confinement direction (z) for the lubricant (mixture of linear IL and 9-octylheptadecane oil) with the applied external pressure of (a) P 1 = 0.23 GPa, (b) P 2 = 0.47 GPa, and (c) P 3 = 1.13 GPa.The inset of the (a) at pressure P 1 is a magnification of the ionic layer with the no-slip boundary condition.The gray dashed line shows the solid-liquid interface, and the dotted line indicates the no-slip boundary within the ionic layer.The addition of IL fully reduces the slip of oil.

Figure 9 .
Figure 9. Normal force as a function of time for (a) linear IL and (b) 9-octylheptadecane systems.Shear forces as a function of time for the case of (c,e,g) V = 20 m/s for all pressures and the case of (d,f,h) P = 0.47 GPa for all velocities.The data are presented for (c,d) linear IL, (e,f) 9octylheptadecane, and (g,h) lubricant systems.Normal forces only depend on applied load and shear forces depend both on the external pressure and the sliding velocity of the system.

Figure 9 .
Figure 9. Normal force as a function of time for (a) linear IL and (b) 9-octylheptadecane systems.Shear forces as a function of time for the case of (c,e,g) V = 20 m/s for all pressures and the case of (d,f,h) P = 0.47 GPa for all velocities.The data are presented for (c,d) linear IL, (e,f) 9-octylheptadecane, and (g,h) lubricant systems.Normal forces only depend on applied load and shear forces depend both on the external pressure and the sliding velocity of the system.

Figure 10 .
Figure 10.Friction coefficient as a function of (a-d) pressure and (e-h) velocity for (a,e) linear IL, (b,f) branched IL, (c,g) 9-octylheptadecane,and (d,h) lubricant systems.The friction coefficients of IL-containing systems decrease with applied pressure at high shear rates.

Figure 10 .
Figure 10.Friction coefficient as a function of (a-d) pressure and (e-h) velocity for (a,e) linear IL, (b,f) branched IL, (c,g) 9-octylheptadecane,and (d,h) lubricant systems.The friction coefficients of IL-containing systems decrease with applied pressure at high shear rates.

Table 3 .
Friction coefficients for all the systems.The friction coefficient of the mixed system is consistently lower than that of pure ionic liquids.