Effects of intercalated water on the lubricity of sliding layers under load: a theoretical investigation on MoS2

Two-dimensional (2D) materials, such as graphene and molybdenum disulfide (MoS2) have recently become some of the most studied nano-materials due to their wide array of technological and industrial applications. Among these, they display great potential as solid lubricants. Friction properties of 2D-materials, however, are very sensitive to environmental conditions, e.g. humidity. In MoS2, for instance, humidity can hamper its tribologic performances. Past experiments and recent ab-initio molecular dynamics simulations have highlighted that, at ordinary temperatures, a possible reason for lower lubricity is the physical interaction of water with the layers. It is, therefore, crucial to better understand the microscopic mechanisms underlying this behaviour, in order to optimise the lubrication performance of these materials. In this paper we used density functional theory simulations and classical molecular dynamics simulations to provide a multi-scale description of how external load affects the energetic, structural and dynamic properties of intercalated water between MoS2 layers. As a result of combining these two different approaches, we provide an atomistic description of the role of intercalated water in modifying the frictional behaviour of physically interacting layers, e.g. MoS2. The identified interlocking mechanism, which is enhanced under load, is relevant for understanding the frictional effects observed for water confined in slit nanopores, and for nanofluidics applications.


Introduction
Gaining more control over energy and material losses in mechanical processes due to friction and wear is one of the greatest and most fundamental challenges of our time [1,2].As a cooperative effort worldwide, the scientific community is committed to finding new materials, coatings and lubricants that can potentially provide low friction and wear especially in moving mechanical processes operating under severe conditions (large pressures and fast shear).However, finding the most effective lubricants that are also environmentally friendly and cost-competitive [3] remains something of an open question.Layered materials, such as molybdenum disulfide (MoS 2 ) and graphene, have recently obtained exceptional attention amongst nanomaterials, due to the wide range of technological applications they can be employed in, e.g. as optical devices, catalytic interfaces for hydrogen storage and as solid lubricants [4][5][6][7][8][9][10][11][12][13][14][15].Solid lubricants are, in fact, gaining more and more interest for applications, such as space applications and nanodevices with moving components, where liquid lubrication is not feasible [14,16,17].
The versatility of two-dimensional (2D) materials is inherent in the lamellar structure such materials exhibit: the layers, held together by weak forces [18][19][20], display large surface areas and large mechanical resistivity.Furthermore, thanks to their well-known low shear strength [14,[21][22][23][24][25][26][27][28][29][30], the lubrication performance of 2D materials, e.g.graphene, MoS 2 and hexagonal-BN, is superior to that of other nanomaterials [31].The link between the outstanding tribological properties of these materials and their mechanical properties (e.g.elastic modulus, strength and fracture) has also been widely explored [32,33].All of these features make them very suitable for applications in the lubrication field, such as lubricant additives, space lubrication materials, and nanoscale lubricating films.
The different effects of humidity on the lubricating properties of graphite and MoS 2 may be related to the different nature of the water-layer interactions.It has, in fact, been suggested that the chemical interaction of water with the reactive defects, such as edges, produced during graphite rubbing, may be the key mechanism to enable graphite lubricity [42], while the physical attraction that dominates the water-MoS 2 interaction may produce an antilubrication effect [43].Indeed, ab-initio molecular dynamics (MD) simulations showed that intercalated water molecules impede interlayer sliding in MoS 2 [44].Such an antilubricant effect has been observed for a single water layer intercalated between both graphene [45] and MoS 2 layers deposited on mica [46].The different hydrophilic character of the investigated interfaces points at the existence of a universal frictional mechanism underlying the antilubricant behaviour of water intercalated within planar layers.Despite the vast experimental advances in the field, a detailed and unique understanding, at the molecular and atomistic level, of how water interacts with layered solid lubricants and what the effect of these interactions are on their lubricity is still lacking.In this paper we investigate such a mechanism at an atomistic level and explain the reason why water molecules are able to interlock MoS 2 layers, thus increasing friction, especially under load.Our results might be relevant to the emerging field of nanofluidics [45], where water and aqueous solutions are squeezed into slit nanopores [47].

MoS 2 : a case study
As a solid lubricant, MoS 2 is vastly employed in applications ranging from the aerospace to automotive industries, thanks to a number of desirable properties such as its ultralow friction coefficient [8,33,35,48,49] and low toxicity [50].However, the aforementioned loss of lubricity in humid environments can pose severe limits on the range and efficiency of MoS 2 applicability.Whereas some evidence points to the importance of chemical oxidation of MoS 2 layers activated by water vapour [51,52], which leads to the formation of MoO 3 with worse frictional profiles [53][54][55]; other studies are bringing the focus on the molecular interactions of water with the layers as the fundamental mechanism behind the observed loss of lubricant properties [43].Theoretical and computational methods can be uniquely placed to shed light on these kind of problems, especially if one is interested in a variety of external conditions that include large pressures and fast shear rates.In particular, atomisitc simulations can provide a detailed insight into the mechanism of interaction of water with solid lubricants.Although not a lot of theoretical literature is currently available on this subject, some recent work reporting the results of static and dynamic (Car-Parinello molecular dyanamics) ab-initio density functional theory (DFT) simulations has been published [44,56] where the authors investigate the nature of the interaction of water molecules intercalated on a single sheet and between slabs of MoS 2 .In these studies they concluded that at ordinary temperatures and in the absence of molecular oxygen, the main mechanism behind the loss of lubricity in humid environments is the physical interaction of molecular water with the MoS 2 layers.Despite being only mildly physisorbed on the basal plane, water molecules seem to disrupt the sliding motion without involving chemical reactions.These results confirm that, at ordinary temperatures in anaerobic conditions, water can increase MoS 2 friction even in the absence of any chemical oxidation.Such findings, even if not entirely conclusive, represent a crucial starting point for carrying out further atomistic investigations into the role of water in MoS 2 lubricity.
While ab-initio calculations are essential in order to accurately describe the energetics of a system, their applicability is limited by the size of the system and the computational costs associated with the methods.In this paper, we will take the findings from recent ab-initio studies of MoS 2 as our starting point to carry out a classical MD investigation.Classical MD simulations enable us to explore various severe external conditions and model much larger systems than the ones investigated with DFT.This kind of investigation can be particularly informative for lubricants like MoS 2 , considering the extreme conditions under which they are expected to perform.
In particular, motivated by the evidence that no chemical reaction is detected between water molecules and the slabs [44,56], we will use the MoS 2 case study (a force field for the wettability of MoS 2 has been recently published [57]) to investigate, at a larger scale, what the bulk effects of intercalated water molecules on the relative motion of two layers in the presence of an applied load are.In this particular study we consider only infinite layers and therefore we have not explicitly investigated possible edge effects such as those described in reference [52].Furthermore, we are aware defects, doping [16,58] and chemical modifications can influence the lubrication properties of MoS 2 and represent very active fields of Both cells are optimised at PBE level in Quantum Espresso.We replicated the cell in x-and y-direction in order to build the larger systems for our MD simulations.
research, however, these mechanisms are out of the scope of the present investigation.

Methods
We performed classical MD simulations with the LAMMPS MD simulation package [59].We investigated the behaviour of water molecules confined between two slabs of MoS 2 under shear in the x-direction and external normal load.In particular, following the work by Levita and co-workers in reference [56], we looked at two water coverages (25% and 50%, where a full layer coverage is defined as one water molecule per MoS 2 unit) which were chosen to the end of investigating the effect of different submonolayer water concentrations.To explore the effects of the normal load on the lubricity of MoS 2 in various simulation conditions, including harsh external conditions, we applied loads ranging from 0.2 to 11 GPa.Indeed the absence of any chemical event in such harsh conditions allowed us to rule out the role of water tribochemistry in affecting the frictional properties of MoS 2 layers.
The final system, with water molecules uniformly distributed between the slabs, was built starting from the geometry of two super-cells (one for each covarage), which were optimised using DFT simulations with the Perdew-Burke-Ernzerhof (PBE) functional [60] in the Quantum Espresso code [61]; the optimised cells are reported in figure 1.We build the extended systems for use in our classical MD simulations by replicating the cells in x-and y-directions such that the systems are 100 times larger than the optimised cells.After replicating the systems, the total number of atoms for the 25% and 50% cases is 10 800 and 12 800 respectively.We consider a system with a normal load applied on the top slab only; in the SI an example of the full set-up used in our simulations is shown.For each external load investigated we perform 5 ns equilibration simulations in the microcanonical ensamble (NVE).In this study we are considering the slabs to be frozen during the whole length of the simulation and, hence, they are treated as rigid bodies.A delicate step of the aformentioned equilibration set up is the control of the target external load.We employed a particularly stable barostat approach to ensure the external load is maintained at the desired value.Following the equilibration step, we then perform shear simulations for 10 ns for each load configuration.A velocity in the xdirection is applied to all of the atoms in both slabs keeping constant the external load on the top slab.In particular, we use v(x) = 100 m s −1 on top slab and v(x) =−100 m s −1 on bottom slab.This choice is further motivated by recent computational works which find MoS 2 to be more suitable for solid lubrication under heavy load and large shear [62,63].Furthermore, within the restriction in which microscopic simulations are done, the sliding velocities we used  were one order of magnitude lower than the 'thermal velocities' of the particles.
For modelling the interactions between MoS 2 and water, we used the force field parameters provided in [57] whereas for modelling waterwater interactions we employed the standard TIP3P model for consistency with the MoS 2 -water force field [64].

Macroscopic friction coefficients
In order to gain a deeper understanding of the mechanism of interaction of water with the slabs, we start by looking at the macroscopic coefficients of friction for both the coverages studied.The coefficient of friction µ is obtained by considering that friction depends upon both the load and the area, i.e.F = αA + µF N , where F is the shear stress (Shearing Force), α is a constant, A is the contact area and F N is the normal load applied [65] (Normal Force).Therefore by plotting F against F N , the slope of the plot gives the friction coefficient µ.
The Normal and the Shearing force for each external load are obtained by taking the cumulative average (every 100 points to control the height of the error bars) of these components during the whole length of the shear simulation and plotting them as reported in figure 2.
We identify two distinct friction regimes.In the low normal load regime (0.2-3 GPa, dotted lines in figure 2) we see larger friction for the lower water coverage (µ l 25% = 0.16 and µ l 50% = 0.11).At around 5 GPa there is a crossing point and the larger coverage system shows a larger friction coefficient (5-11 GPa, solid lines in figure 2, µ h 25% = 0.05 and µ h 50% = 0.10).In this second regime it was more difficult to perform simulations with a stable force behaviour and a larger standard deviation of measured forces emerges for each point.To understand how such macroscopic behaviour correlates with the atomistic mechanism we then look more in detail at the results of our shear simulations.

Sheared water structures
The length of the shear simulations is determined upon checking the convergence of the total energy of the system with respect to time.For all the systems investigated a simulation time of 10 ns corresponded to converged energy values.By looking at the configurations of our MD simulations at t = 10 ns, we can see that, upon thermal equilibration followed by load and shear, the water molecules undergo some significant structural rearrangements.When the systems are compressed to the desired loads the water molecules began to arrange themselves into clusters (during the equilibration step).As a result of the normal load and shear, these clusters then turn into wellordered stripes.
Figure 3 illustrates this behaviour for the 25% and 50% case at three different external loads to represent the low, medium and large pressure regimes.In the 25% coverage system, the stripes of water molecules are oriented in the same direction as the shear velocity only for the very low load cases (0.2 and 0.5 GPa in the SI) but starting from 1 GPa the stripes tend to sit diagonally with respect to the shear direction.As the external load increases the water stripes tend to be less ordered and less aligned with the direction of shear.In particular, the stripes of water molecules stop aligning with the direction of shear at 1 GPa for the low coverage case.In the 50% coverage case, this transition occurs at 5 GPa.Then as the load increases, the stripes of water molecules continue to become thinner and remain oriented such that they are at a diagonal with the direction of the shear velocity.At the lower loads (0.2, 0.5 and 1 GPa) the stripes of water molecules are oriented in the same direction as the shear velocity.An example of these structures can be observed in panel d in figure 3 for the 1 GPa case; 0.2 and 0.5 GPa cases are reported in the SI.

Cluster analysis
The built-in cluster analysis program in the Ovito simulation analysis package [66] was used to monitor the evolution of the aggregation of the water molecules throughout the shear simulations.In doing so, if any two water molecules were within 3.2 Å of each other in a given timestep then they were said to be in the same aggregate or cluster.Therefore we measure the size of hydrogen bonded clusters of water molecules based on nearest neighbour distances at the different loads, but we do not measure the structure of the water molecules within these clusters.What we observe can be put into direct context with the snapshots from the final frames reported in figure 3.

25% case
As shown in figure 4 the oxygen cluster analysis shows that at larger pressures we observe a larger number of clusters, which is in line with the observation that the stripes are less ordered in the large load regime (panel a in figure 4).The size of the largest oxygen cluster in each system was also determined (panel b in figure 4).At low load (1 GPa), the system quickly settles into a state where the size of the largest cluster is nearly constant in size, containing 400 water molecules which corresponds exactly to number of water molecules in the system.This is in line with the visual observation in this system of a single stable stripe for the whole duration of the shear simulation.Upon increasing the load the system displays wider fluctuations in the maximum cluster size.At 5 GPa, the size fluctuates significantly between 100 and 400 water molecules throughout the simulation.This suggests the formation of multiple but smaller stripes, as seen in the snapshots.At 11 GPa, we observe smaller fluctuations and smaller clusters, again in line with what we found by looking at the snapshots from the simulations.As a result, the effect of the external load on the largest oxygen clusters is to break the big water stripes into smaller and equally sized clusters.
It is important to point out that, the drastic fluctuations displayed in the data are also partly due to the normal thermal movements of water molecules which may cause for clusters to break/reform.However, the average behaviour that emerges is in line with the observations we gather from simulations.

50% case
The results of the cluster analysis for the 50% case display significantly less fluctuations.The number of clusters increases with the normal load, again displaying a less structured system at larger pressures (panel c in figure 4).Interestingly, the number of clusters in all of the cases is ten times smaller than the 25% case.By looking at the largest cluster size plot (panel d in figure 4), they also show that at 1 GPa all of the water molecules in the system (800) are involved in forming a single cluster (stripe).The maximum size decreases with load suggesting the formation of multiple stripes and the maximum size shows significant fluctuations of the measured value at 11 GPa.Nevertheless, compared to the 25% coverage systems, even at 11 GPa the cluster size remains large at around 700 water molecules.
In summary, it appears that one of the effects of the external load on the systems is to prevent the water molecules from forming ordered stripes; this means that at larger loads water molecules form smaller networks with each other.

Discussion and conclusion
Previous theoretical and studies of xenon on a graphene monolayer revealed that changes in the friction behaviour of the substrate can occur above a critical coverage of the adsorbate.In reference [67], the authors have investigated the atomistic mechanisms underlying the observed phenomenon by means of MD simulations.The results uncover the existence of a close correlation between the island size and commensurability: small islands are in register with the substrate, while larger islands are less commensurate.Simulations carried out under equivalent conditions on gold unveils a similar trend.However, a much lower commensurability is found for each investigated temperature and coverage due to a larger mismatch between the adlayer and the substrate [68].On the premises of the relation between static friction and interfacial commensurability, such results can, thus, support the existence of a critical coverage for the depinning transition and can explain its dependence on temperature as observed in quartz crystal microbalance (QCM) experiments [68].Furthermore, these observations confirm the theory on the size dependence of static friction, according to which nominal incommensurate interfaces commensurate below a critical size [69].The critical size, which depends on the relation between interparticle interaction strength and the potential corrugation, is dominated by the lattice mismatch.The aforementioned studies were the first, to our knowledge, to demonstrate that the size of adsorbate islands at a material's surface depends not only on the interparticle interaction strength of the adsorbate and the corrugation in the potential energy surface (PES) but also the lattice mismatch between the substrate's crystal and the structure of the adsorbate island.These are important findings when combining the results of our MD and DFT simulations in order to provide a more wholistic view of the results.As a result, in order to gain a deeper understanding of such relations for the system investigated in this paper, in figure 5 we report the computed PES for a water molecule on a MoS 2 layer obtained within the DFT-D scheme [70] including Van der Waals interactions to correct the exchange correlation energy, here described with the PBE parametrization [60].
The adsorption energy (E ads ) is expressed in mJ m −2 and is calculated by keeping fixed the lateral coordinates of the oxygen atom at the high symmetry sites indicated in figure 5 and relaxing the adsorption distance to its equilibrium value for every lateral position.In agreement with previous calculations [56], the most energetically favourable adsorption location is obtained when the water molecule occupies the hollow site at the centre of three sulfur atoms of the outermost atomic plane (in green).The less energetically favourable position corresponds, instead, to the on-top configuration, where the O atom is located above a S atom.The energy difference between the PES maxima and minima, i.e. the PES corrugation corresponding to the maximum energy barrier that the molecule has to overcome when diffusing along the layer, is ∆E ads = 6 mJ m −2 .We investigated how this picture is affected by an applied load by calculating the PES in the presence of a vertical force on the O atom and keeping the MoS 2 layer frozen at its equilibrium configuration.The results obtained for 1 and 5 GPa loads, reported as supplementary information (is available online at stacks.iop.org/2DM/8/035052/mmedia), show that the PES does not change its shape with load, the most (un)favourable adsorption configurations remain the (top) hollow sites.However, the PES corrugation increases by orders of magnitudes with load, becoming equal to ∆E ads = 50 mJ m −2 at 1 GPa and ∆E ads = 280 mJ m −2 at 5 GPa pressures.As a result, we expect that the molecular diffusion on the layer becomes more and more difficult as the load increases, with the waterlayer interaction prevailing over the water-water interaction.
The behaviour of the friction coefficients we computed from our MD simulations can be correlated to the water structures we observe in the final simulation snapshots.In the low pressure regime (0.2-1 GPa) we expect the water-water interactions to prevail over the water-slab interactions.As a result, we observe the formation of large and ordered stripes as water will more easily diffuse on the PES and occupy positions that optimise the inter-molecular interactions rather than those that optimise the moleculesubstrate interaction.This can be directly related to the fact that surface corrugation increases with the external load (so lower pressure means less corrugation).To further explain why in this regime we observe a larger friction coefficient with the lower coverage we need to consider that, as discussed earlier in this paragraph and in [69], smaller islands become more commensurate with the surface which results in larger friction coefficients.Therefore, in the low pressure regime, we find that the water forms smaller clusters (that are more commensurate with the MoS 2 surface) in the 25% coverage system than in the 50% coverage system, which in turn results in the larger friction coefficients within this regime for the 25% coverage systems.
At large pressure we need to take into account the effect of the external load.In this case the water-slab interactions start to compete with the water-water interactions.We expect that, at larger pressures, the corrugation of the surface also becomes larger (the barriers for water to move and form networks become larger) and water tends to remain trapped in the minima of the PES for longer time forming less ordered stripes and smaller clusters.This behaviour emerges for both of the coverages for which we can observe the disruptions of the aligned and ordered patterns and formation of clusters.However, as a result of having more water molecules in the 50% case, the formation of small water clusters is even more evident at large pressures, which therefore leads to the larger measured friction coefficient for this system.
In conclusion, we believe this atomistic investigation can shed light on the bulk effects of intercalating water molecules within layered structures under load and shear.Such effects can be related to the observed loss of lubricity in some 2D materials, e.g.MoS 2 , when in presence of humidity.We show this can be directly correlated to the external working conditions.In particular, in the presence of large sheer velocity and large external load on the slabs, the friction behaviour of MoS 2 is influenced by the way water interacts with the slabs.The effect of external load is, hence, to prevent water molecules from forming large and hydrogen-bonded clusters between slabs.The larger external loads favour the formation of smaller water islands that commensurate with the surfaces; this ultimately leads to an increase of the friction coefficients.
Despite having used MoS 2 as a case study, we think that the investigation proposed in this paper, consisting of the analysis of the effects of the loaddependent PES corrugation (obtained with DFT) on the water structures and diffusion properties (observed with classical MD simulations) can be applied to any intercalated species at 2D materials in order gain a better understanding of the correlation between microscopic and macroscopic friction properties.

Figure 1 .
Figure 1.The optimised geometries employed to build our systems of water molecules between MoS2 (yellow atoms are S and pink atoms Mo) slabs are here reported: (a) and (b) show side and top view of the system where the number of water molecules placed between the slabs correspond to the 25% coverage of the slab surface; (c) and (d) show side and top view of the 50% case.Both cells are optimised at PBE level in Quantum Espresso.We replicated the cell in x-and y-direction in order to build the larger systems for our MD simulations.

Figure 2 .
Figure2.Average Normal Force vs. average Shearing Force for the two coverages investigated (blue circles are for 25% case and turquoise triangles for 50% case).Black solid lines refer to full curves whereas coloured lines are linear regressions.We identify two trends in the curves: a low pressure regime (dotted lines) where the slopes of the lines, which represents our friction coefficients, are µ l 25% = 0.16 and µ l 50% = 0.11; and a large pressure regime (solid lines) where we have µ h 25% = 0.05 and µ h 50% = 0.10.

Figure 3 .
Figure 3. Snapshots taken at the end of 10 ns shear simulations with three different normal loads are shown.

Figure 4 .
Figure 4. (Left column) Number of clusters of water molecules in the simulated systems as a function of time.(Right column) Number of water molecules in the largest cluster as a function of time.The top row contains plots representing the 25% coverage systems; while the bottom row contains plots from the 50% coverage systems.

Figure 5 .
Figure 5. Potential energy surfaces at zero applied load for a water molecule (not shown) on a MoS2 layer (Mo atoms are in grey and sulfur atoms in green).The PES is constructed by calculating the molecule adsorption energy at the different high-symmetry positions indicated by the black dots in the unit cell.The data are computed within the DFT-D scheme, including the Van der Waals interactions.