Exploring Nanoscale Lubrication Mechanisms of Multilayer MoS2 During Sliding: The Effect of Humidity

Solid lubricants have received substantial attention due to their excellent frictional properties. Among others, molybdenum disulfide (MoS2) is one of the most studied lubricants. Humidity results in a deterioration of the frictional properties of MoS2. The actual mechanism at the nanoscale is still under debate, although there are indications that chemical reactions are not likely to occur in defect-free structures. In this study, we performed nonequilibrium molecular dynamics simulations to study the frictional properties of multilayer MoS2 during sliding in the presence of water. Moreover, we also investigated the effect of sliding speed and normal load. We confirmed earlier results that a thin layer of water organizes as a solidified, ice-like network of hydrogen bonds as a result of being confined in a two-dimensional fashion between MoS2. Moreover, we found that there exists an energy-driven, rotational dependence of the water network atop/beneath MoS2. This orientational anisotropy is directly related to the dissipative character of MoS2 during sliding. Finally, three distinct frictional regimes were identified, two for a thin layer of water and one for bulk water. In the case of a thin layer and low coverage, water represents a solid-like contaminant, causing high energy dissipation. For a thin layer and high coverage, water starts to act as a solid-like lubricant, reducing dissipation during sliding. Finally, a regime where water acts as a liquid lubricant, characterized by a clear velocity dependence was found.


INTRODUCTION
The presence of friction and wear has been linked to a loss of almost a quarter of the world's total energy consumption (Holmberg and Erdemir, 2017). The current demand of our industry has surpassed the application of "classical" liquid lubricants in peculiar conditions, such as ultrahigh vacuum, extreme contact pressures, super low/high temperatures, and extremely small dimensions. A new generation of lubricants, the so-called solid lubricants, offers new possibilities. Moreover, solid lubricants have been linked to superlubricity, a state in which the frictional forces vanish (Baykara et al., 2018;Martin and Erdemir, 2018). Among others, transition metal dichalcogenides (TMDs), and in particular molybdenum disulfide (MoS 2 ), are prototypical examples. TMDs are layered materials, where the monolayers have the so-called two-dimensional (2D) crystal structure with a transition metal covalently bound between chalcogens. The layers are held together by weak electrostatic and van der Waals interactions. As a result of this unique binding character, easy shearing between the layers is accommodated, creating an enormous potential for achieving low coefficients of friction (COF) (Martin et al., 1993;Scharf and Prasad, 2013). This is confirmed by numerous studies on their remarkable tribological characteristics (Martin et al., 1993;Watanabe et al., 2004;Cho et al., 2006;Evaristo et al., 2008;Pimentel et al., 2011;Mutafov et al., 2015). This, in combination with their low toxicity (Teo et al., 2014;Chng and Pumera, 2015), versatile chemistry (Chhowalla et al., 2013), and peculiar electronic properties, has resulted both in macroscale applications, such as the automotive, aerospace, and space industries (Cho et al., 2006;Nian et al., 2017;Vazirisereshk et al., 2019) and nanoscale applications, such as photovoltaic and optoelectronic devices, transistors, and solar cells (Li and Zhu, 2015;Choi et al., 2017).
In a significant number of (potential) applications, ambient conditions, such as humidity, are unavoidable, resulting in contaminated surfaces. Humidity however has shown to induce a significant alteration in the tribological properties of TMDs and other 2D materials, such as graphite/graphene (Evaristo et al., 2008;Scharf and Prasad, 2013;Levita et al., 2016;Levita and Righi, 2017;Hasz et al., 2018). In the case of MoS 2 , humidity not only leads to an increase in wear but also to a significant degradation in frictional properties (Panitz et al., 1988;Chhowalla and Amaratunga, 2000;Zhao and Perry, 2010;Khare and Burris, 2013;Levita and Righi, 2017;Serpini et al., 2017;Arif et al., 2019;Serpini et al., 2019). For example, it has been found that clusters of contaminating adsorbates can induce a static friction, seemingly "locking" surfaces in relative motion that would otherwise slide smoothly (Ouyang et al., 2018). As a result, MoS 2 solid lubrication currently has a relatively small effective working range, limited to outer space or vacuum (Hasz et al., 2018). To extend its applicability, one should fully understand how humidity affects the lubricating mechanisms at the nanoscale during sliding.
Although the scientific community agrees that the beneficial frictional properties of MoS 2 deteriorate in case of increasing humidity (Zhao and Perry, 2010;Dudder et al., 2011;Khare and Burris, 2014;Serpini et al., 2017), the actual mechanism at the nanoscale is still under debate (Hasz et al., 2018). Some ascribe it to oxidation, whereby MoS 2 oxidizes to MoO 2 (Panitz et al., 1988;Tagawa et al., 2007;Curry et al., 2017). Others showed that humidity does not necessarily lead to oxidation but that water physically interacts with the surface, thereby increasing interlayer friction (Khare and Burris, 2013;Levita and Righi, 2017;Serpini et al., 2017;Arif et al., 2019). Some use an argument based on adhesion, when explaining the changing frictional properties as a result of humidity (Hasz et al., 2018). In an experimental study, Lee et al., (2020) showed that when a layer of water is intercalated between mica and exfoliated MoS 2 , the frictional properties deteriorated. Interestingly, similar results were found for graphene and mica , contradicting previous reports that humidity leaves the frictional behavior unchanged (Berman et al., 2014). These results hint at a universal frictional mechanism for intercalated water layers that is independent of the type of solid lubricant itself (Lee et al., 2020).
The frictional response of a material is highly complex and characterized by nonlinear, out-of-equilibrium dissipative processes (Vanossi et al., 2013). Simulation techniques offer possibilities to gain insights by disentangling the phenomena and using simplified models. Among others, Levita et al. made several computational contributions (Levita et al., 2016;Levita and Righi, 2017). In one study (Levita et al., 2016), the type of interaction between MoS 2 and water was determined using firstprinciple calculations (FPCs). The literature was discordant on whether this should be depicted as hydrophobic or hydrophilic (Late et al., 2012;Gaur et al., 2014). In addition, the tendency of MoS 2 toward "adsorption of" vs. "being oxidized by" water was investigated. Due to the stable basal plane of MoS 2 (lacking any dangling bond and characterized by low chemical reactivity), it was found that also in case of defects, such as S-vacancies, only physisorption occurs. This observation was confirmed by others (Arif et al., 2019). Upon introducing water within the bilayer interface, the interlayer distance was found to increase, weakening the interlayer binding energy and destabilizing the structure. It was confirmed that MoS 2 does not favor oxidation (e.g., to MoO 2 ) and acts rather hydrophobically (Levita et al., 2016). This conclusion has been shared by others (Luan and Zhou, 2016) using molecular dynamics (MD) simulations.
To study the frictional properties, sliding should be considered; after all, kinetic friction is a dynamical response. In another study, based on ab initio MD, the chemical and physical interactions taking place during sliding of a single asperity of MoS 2 , in the presence of water, were studied (Levita and Righi, 2017). At high temperature and load, it was found that increasing water coverage present within the interface reduces the distance traversed after an initial kick. Again, and as expected from previous studies (Khare and Burris, 2013;Levita et al., 2016), no oxidation occurred. Furthermore, it was found that the water molecules formed a water network of H bonds. Finally, S-vacancies did not induce significant deviations in the results. These observations were confirmed mainly by an experimental study (Arif et al., 2019), indicating that the main mechanism affecting the frictional properties of MoS 2 , in the presence of water, is physisorption.
The earlier mentioned results are relevant, but some remarks are in place. Levita and Righi (2017) used a high load (11 GPa) and temperature (320 K) to provide enough energy to the system to overcome potential energy barriers, which might affect the outcomes. Additionally, the simulation time was limited to the picosecond range, affecting the statistical precision of the results. Moreover, a structural analysis of the water network was not properly considered. When liquids, and water in particular, are confined to a small space, their properties might alter significantly (Bampoulis et al., 2016); this is in fact one of the problems of classical lubricants. In the case of water confined between two hydrophobic layers, ice-like water films have been found (Bampoulis et al., 2016), raising the question whether water acts as a (solid) lubricant or as a friction inducer (Lee et al., 2020). In the case of MoS 2 , Kwac et al. (2017) found that using MD, a characteristic tile-like pattern of squares and diamond forms, resulting from the MoS 2 surface charge. However, the sliding of the MoS 2 layers was not considered in this study.
The main aim of this study is to present a qualitative and quantitative analysis of the nanoscale frictional properties when sliding multilayer MoS 2 in the presence of water. More generically, our goal is to reveal broader features of the frictional behavior of solid lubricants in the presence of contaminants. Based on previous results (Khare and Burris, 2013;Levita and Righi, 2017), tribochemical reactions were not expected to occur. Therefore, classical nonequilibrium MD simulations were performed. This was done by considering the intercalation of various number of water molecules, varying from dry sliding to bulk water, located at the center interface. Furthermore, the effects of normal load and of sliding speed were investigated. Apart from gaining insights into the conformational landscape of water, this setup allows for insights into the mechanisms of the additional channel of dissipation, present when considering intercalated water, in the case of bulk MoS 2 and that of solid lubricants in general.

COMPUTATIONAL METHODS
The simulations have been carried out by means of the LAMMPS package (Plimpton, 1995), employing a force field developed by Sresht et al. (2017) for MoS 2 . Its foundation is based upon harmonic bonds and angles, with partial charges set to Mo +0.50e and S −0.25e, in combination with the Lennard-Jones (LJ) potential (Sresht et al., 2017). Water was modeled through a nonpolarizable q-SPC/Fw model (Paesani et al., 2006), characterized by a three-site water molecule, where the charges were set to O −0.84e and H +0.42e. Using Lorentz-Berthelot mixing rules, σ ij (σ ii + σ jj )/2 and ϵ ij ϵ ii ϵ jj √ , cross-term parameters were obtained. A cutoff of 12 Å was used for the LJ potential. To account for electrostatic interactions, the particle-particle particle-mesh (PPPM) method was used, with a threshold of 1.0 × 10 −6 for the relative error in the forces. In Figure 1 and in the Supplementary Material, a benchmark of this force field is presented against FPC results (Levita et al., 2014) and experimental data (Schönfeld et al., 1983), by means of the sliding potential energy surface (PES) and structural parameters, respectively. The PES was obtained by performing static calculations on a partially rigid MoS 2 bilayer. Here, only the top layer was allowed to breathe in the z-axis and was translated along the x-axis and the y-axis in twelve steps, for a total distance equal to the size of the period in that direction. (For more details, we refer to Nicolini and Polcar (2016).). The employed force field shows good agreement with both ab initio and experimental data.
First, the six-layered 2H c -MoS 2 system was built using the initial coordinates from Schönfeld et al. (1983). The resulting structure consists of 4,608 atoms and a surface of approximately 44 Å by 51 Å. An energy minimization was performed using a sequence of conjugate gradient (CG) algorithm, with a relative energy tolerance of 1 × 10 −15 , and damped dynamics (DD) algorithm, with a relative energy tolerance of 1 × 10 −10 ; the simulation box was also allowed to relax. Second, an isolated water molecule was optimized according to the force field using CG. The resulting structure was replicated in three dimensions to create a large cube, with a density equivalent to that of liquid water at room temperature, and optimized using both CG and DD, while relaxing the box. Third, the cube of water was put atop three layers of MoS 2 and subsequently optimized using both CG and DD. Thereafter, the water cube was heated up to 300 K using the Nose-Hoover thermostat (Martyna et al., 1994), with a temperature damping parameter of 100 time steps. In the fourth step, the actual structures were created. Random clusters of water molecules of the melted cube, located within FIGURE 1 | Static potential energy surfaces calculated per area for the long (inset: short) diagonal of the hexagonal unit cell. Four different cases are considered: i) in black, pure MoS 2 obtained via DFT (Levita et al., 2014), ii) in red, pure MoS 2 obtained with the force field employed in this study (Sresht et al., 2017), iii) in green, MoS 2 and an ideal, parallel water network via MD, and iv) in blue, MoS 2 and an ideal, antiparallel water network via MD. The dashed lines are obtained via akima spline interpolation and serve merely as a guide to the eyes.
Frontiers in Chemistry | www.frontiersin.org June 2021 | Volume 9 | Article 684441 a 5 Å cutoff above the nearest sulfur layer, were selected. The reason that we chose clusters, instead of randomly spread water molecules, had to do with the fact that in preliminary simulations, we found that the sliding of MoS 2 drives the formation of clusters and therefore are the best representation. Finally, the missing three top layers of MoS 2 were pasted above. By considering the "tiled" structure of the water network mentioned in Kwac et al. (2017), full coverage requires 304 H 2 O molecules. In total, 11 different systems were created consisting of 0 (dry), 1, 2, 5, 10 (3% coverage), 26 (8.5%), 52 (17%), 103 (33%), 206 (67%), 304 (100% full coverage), or 1,236 (bulk) water molecules. Finally, a vacuum of 20 Å was added above each structure in order to avoid interaction with the repeated images while applying periodic boundary conditions. After obtaining the structures, all were once more optimized using CG and DD, while keeping the box dimensions fixed. This time, the bottommost sulfur layer was held fixed in the x-, y-, and z-direction. Moreover, the topmost sulfur layer was made rigid in the xy-plane and the forces in the z-direction were averaged over all sulfur atoms in this layer. This restriction on the top sulfur layer was added to mimic adhesion to the substrate and to be able to control the sliding. After the minimization, a 2-step equilibration was performed, in which seven configurations of every system were created, each with a different normal load applied to the topmost sulfur layer. The normal load was ranging from 0 to 3.0 GPa, with increments of 0.5 GPa. In the first step (100 ps), the load was slowly increased from 0 GPa to the target value, whereas in the second step (100 ps), it was held constant. During the equilibration, the system was thermalized by initializing the temperatures for every atom using the Boltzmann distribution. Moreover, the bottom and top layers of MoS 2 were coupled to a Nose-Hoover thermostat (Martyna et al., 1994). On the one hand, this allows for maintaining the temperature during equilibration, while on the other hand, this avoids overheating from Joule heat production during the actual sliding runs (Vanossi et al., 2013).
Three independent production runs for every configuration were created. In all runs, the bottommost layer of sulfur atoms was held fixed in the x-, y-, and z-directions. Next, the topmost sulfur layer was moved rigidly in the x-direction, along the Mo-S bonds, with three different velocities (20.0, 2.0, and 0.2 m/s) and for a total of ∼27.5Å. For the analysis, the first part of the MD trajectories was discarded and considered as the run-in phase, and only the last ∼22Å was used, which is the distance based on four periods of the cell in that direction (∼5.5Å). The velocities were maintained constant along the trajectories, which was achieved by setting an initial velocity and canceling the force in this direction. On the one hand, this setup allows for the study of the effect of velocity in the sliding direction, and on the other hand, the system keeps the freedom to slide the MoS 2 layers in the y-direction and to breathe in the z-direction. The reason to slide along this direction is because this direction represents the highest friction (Claerbout et al., 2019), and therefore, the most interesting dynamics was expected to be observed. Along all runs, a time integration step of 0.1 fs was used. In Figure 2, the computational setup is depicted.

Sliding Dynamics
We present a full analysis of the results obtained for 2.0 m/s and describe some of the differences and similarities to sliding with FIGURE 2 | Schematic overview of the computational setup. The side view of the six-layered MoS 2 structure, with a cluster of 26 H 2 O molecules intercalated in the center interface, is displayed. Molybdenum, sulfur, oxygen, and hydrogen atoms are displayed in silver, yellow, red, and white, respectively. This specific case resembles the system after equilibration at 0 GPa load and 300 K. From this picture, there is a clear structural deformation visible in the MoS 2 layers directly encapsulating the H 2 O cluster. Figure was obtained with VMD (Humphrey et al., 1996).
Frontiers in Chemistry | www.frontiersin.org June 2021 | Volume 9 | Article 684441 5 20.0 and 0.2 m/s. First, a qualitative analysis on the water dynamics is presented, followed by MoS 2 . As stated earlier, the water molecules in the clusters were randomly selected from a melted ice cube. This allows to observe the naturally resulting water network after being transferred from a 3D "free" space to an encapsulated 2D confinement. In general, a flat, single layer of water molecules was formed (see Figures 3A-J). For bulk water, the formation of two water layers was observed, facing MoS 2 on either side, between which water molecules in a more disordered state were found.
In Figures 4A-C, a top view of the water network for three systems is displayed. For up to 8.5% coverage, we found single clusters, kept together by a network of H bonds. In the present study, a hydrogen bond (H bond) is defined as two oxygen atoms, between which there is a maximum distance of 3.5 Å, and where the O-H-O angle is 150°< α < 210° (Kwac et al., 2017). Through this H bonds, rings of either three or four molecules were formed. Similar observations were made for increasing water coverage; however, the majority of the rings now consist of four molecules. In fact, a similar type of water network, as described in Kwac et al. (2017), can now be distinguished. This ice-like, the so-called perpendicular crossing diamond network (pcd) consists of rings with four molecules, which are either diamond-like, where the four H bonds are "donated" by two of the water molecules, or square-like, where each water molecule contributes with one H bond. As an example, we put the diamond shapes atop the water network in case of full coverage and 67% coverage. From these observations, it can be concluded that water acts as a solidified film when strictly confined, despite of temperatures above the melting point. Changing the load or the sliding velocity did not alter significantly the formation of the structure.
In an "ideal" water network, the pcd configuration should result in four different H bond angles, two per diamond orientation, which sum up to 180°. We found two main orientations of the water network with respect to the simulation box. We define the first orientation, which was the most common in our sliding MD simulations, as the antiparallel one, where none of the H bonds align along the sliding direction. We define the second orientation, which was only present in a few cases, as the parallel one, where half of the H bonds are oriented practically along the sliding direction. More specifically, in the antiparallel case, these angles are ∼ 37.5°, ∼ 52.5°, ∼ 127.5°, and ∼ 142.5°(obtained by doing a separate simulation using an ideal "defect-free" water structure without MoS 2 , see the Supplementary Material), whereas in the parallel case, these angles are ∼ 7.5°, ∼ 82.5°, ∼ 97.5°, and ∼ 172.5°. In Figures  4D-F, the angles of the H bonds with respect to the x-axis, present at every time step, are displayed, including the corresponding histograms. In the case of 8.5% coverage or less, three main angles (the middle one with a wider distribution) could be distinguished, which are ∼ 30°, ∼ 90°, and ∼ 150°. This is consistent with the fact that in this case, both FIGURE 4 | Top view representation of the water network for 8.5%, full, and 67% coverage (A-C). The black arrow represents the sliding direction. The blue diamonds serve as a guide to the eyes in distinguishing the pcd configuration. The bottom row (D-F) represents the angles of all H bonds, at every time step, with respect to the x-axis; in red, the corresponding histograms are depicted. Figure was obtained with VMD (Humphrey et al., 1996).
Frontiers in Chemistry | www.frontiersin.org June 2021 | Volume 9 | Article 684441 three-and four-molecule rings are present. Overall, the H bonds were not aligned along the sliding direction. By increasing the number of water molecules, the pcd network started to form, showing up in six different angles for a water coverage up to 67%, which are ∼ 30.0°, ∼ 50.0°, ∼ 70.0°, ∼ 110.0°, ∼ 130.0°, and ∼ 150.0°(data not shown). Here, the summation of one of the three lower angles combined with one of the three higher ones gives 180°, resulting in a diamond shape. Once more, none of the H bonds were aligned along the sliding direction. For full coverage, four sets of angles were found, again whose sum totals 180°(see Figures 4B,E). On the contrary, in the parallel case as depicted in Figures 4C,F as an example, three angles were found which are ∼ 10°, ∼ 90°, and ∼ 170°, where the central one is the superposition of the two distinct peaks found in the ideal case. The parallel orientation is energetically less favorable (by about 0.04 kcal mol −1 Å −2 , with respect to the static calculations presented in Figure 1) than the parallel one, and generally, along a sliding MD trajectory parallel to antiparallel transitions (but not vice versa) were observed. The deviation between the actual and the expected angles is due to the sliding energy provided to the system, which allows several metastable states where the H bond network experiences a breathing motion. Moreover, the discrete character of the relative orientation of the clusters above/beneath MoS 2 hints at a rotational disorder energy dependence, which was also revealed by the energy difference between the two ideal cases, referred to above. At any rate, to have a consistent quantitative comparison in terms of the frictional properties, only antiparallel configurations were considered. Finally, considering bulk water, we note that this is the only case where the water actually behaves as a liquid. Next, we analyze the MoS 2 dynamics, first from a top view and thereafter from a side view perspective. From the top view, we found a typical characteristic observed in nanoscale friction for 2D solids, the so-called stick-slip dynamics (Lee et al., 2010). In the case of stick-slip, two velocity regimes can be distinguished, a high speed one during which a sudden jump in position occurs and a stick phase. Stick-slip happens as a result of the underlying interactions between two surfaces relative in motion, represented by the potential energy surface (PES), and can be explained using the principles of the Prandtl-Tomlinson (PT) or the more generalized Frenkel-Kontorova model (Socoliuc et al., 2004;Vanossi et al., 2013). Consider a ball being dragged over a PES by a moving spring. When encountering an energy barrier, the ball will stay trapped in the minimum and the spring will be stretched, that is, the stick phase. At a certain point, the energy stored in the spring will be sufficient to pass the energy barrier, the ball "slips" to the next local energy minimum, and as a result, the energy stored in the spring is dissipated. The higher the barrier, the higher will be the dissipation. In general, it is believed that in the temperature-extended PT model and in friction force microscopy experiments, stick-slip is a thermally activated process, whereas in MD simulations, this is only believed to be true for very low sliding speeds (Vanossi et al., 2013).
In our case, the moving spring is indirectly represented by the topmost sulfur layer, which was moved with a constant velocity in the x-direction and the van der Waals interactions between the layers. From the above, it becomes clear that whenever the PES has no effective corrugation, stick-slip will not occur. It is found that the corrugation of the PES goes hand in hand with the commensurability of the system (Irving et al., 2017;Wang et al., 2017), as in the case considered in this study. Whenever there is a mismatch, be it a rotational, structural, or a combination of both, the PES flattens out due to the cancellation of forces. In this case, one will observe smooth dynamics, characterized by low energy dissipation.
We observed the stick-slip dynamics for the low water coverage regime (up to 17%), as can be observed in the third column of Figures 3K-O, where the x and y components of the layer's center of mass position are reported. From the trajectories, it becomes clear that the slip events can take place in both x-and y-directions for layers 2-5 (albeit slips are observed more frequently in the x-direction). For layer 6, the slip can only occur in the y-direction since the velocity in the x-direction is held constant. In the Supplementary Material, the reader can find a set of graphs where the profiles of both components of the center of mass position as a function of time are reported. In addition, also the time profiles of the lateral force acting on the top layer can be found in the Supplementary Material. From these graphs, a clear sawtooth-like pattern is present for low water coverages, once more confirming the presence of stick-slip dynamics. All these results are in line with the PES of commensurate MoS 2 , displayed in Figure 1. When the load was increased, stick-slip is enhanced and present for higher water coverage. This observation is in line with earlier results. In 2014, Levita et al. (2014) studied the sliding properties of dry MoS 2 using first-principle calculations. With an increasing load, the corrugation of the PES increased as a result of changing van der Waals, electrostatic, and Pauli interactions. Moreover, purely by considering the functional form of the force field, especially the LJ part describing the interactions between MoS 2 layers and between MoS 2 and water, one could also have predicted this outcome. Increasing the load reduces the interatomic distances and therefore results in a regime with more repulsive LJ interactions. When the velocity was increased, we found that the stick-slip dynamics disappear, even for dry sliding. A possible explanation for this observation is that the adhesion between the layers is reduced due to the increased velocity. This lowers the probability for the atoms to relax to their local minimum energy positions (Onodera et al., 2010). In other words, the layers tend to "fly" over each other, without experiencing a strong interlayer interaction. Although not presented in Figure 3, we also found that the water clusters themselves display stick-slip dynamics up to this low coverage regime.
In general, when the water coverage was increased, we observed a transition to smooth dynamics, both for the MoS 2 layers and the water clusters. However, in the special case of a parallel water network orientation, stick-slip dynamics persisted again, both in the MoS 2 and the water clusters. This behavior is also suggested by the static PES profiles (Figure 1), where the PES is displayed for the "ideal" parallel and antiparallel water network compared to "pure" MoS 2 . The antiparallel PES is flat due to lattice incommensurability with MoS 2 . The parallel PES however shows that there is in fact a commensurability present along the Frontiers in Chemistry | www.frontiersin.org June 2021 | Volume 9 | Article 684441 sliding direction, where twice the length of the H bond gives the lattice period of MoS 2 . Next, we analyze the MoS 2 dynamics from the side view perspective. For low water coverage up to 3%, we observed a "card-deck" shearing behavior, where the bottommost layer remains in place, while the topmost layer moves according to the driving velocity. The intermediate layers present increasing displacements, depending on their relative position in the layered system. At the value of 8.5% water coverage, there is a transition toward a new type of dynamics, where we found that MoS 2 starts to divide into two blocks, one above and one beneath the water molecules. For moderate coverage (17% and above), the layers in the bottom block, where the bottommost sulfur atoms are fixed in the x-, y-, and z-directions, remain together and do not reposition. Additionally, the layers in the top block also adhere to each other; however, this block slides as a whole according to the driving velocity given to the topmost sulfur layer. This is confirmed by Figures 3A-J, where snapshots at the beginning and end of the trajectory are displayed for different amounts of water coverage. In these specific examples, the water network is oriented in an antiparallel fashion for full coverage, whereas for 67%, this is parallel. We found that for the low-water coverage regime, increasing the velocity slightly shifts the transition from the card-deck shearing to higher water coverage; the same accounts for increasing the load. In other words, card-deck shearing is more pronounced for higher loads and velocities.
Using the PES profiles once more, we can also explain the card-deck shearing vs. block sliding. Let us consider the six layers of MoS 2 in a minimum energy configuration and start sliding solely the top layer. At a certain moment, an energy barrier will be reached according to the PES, resulting from the interaction with the fifth layer. In this case, instead of passing the barrier, it is energetically more favorable to slide both the top and the fifth layers until the fifth layer reaches the energy barrier from the interaction with the fourth layer. This process continues until the bottom layer is reached, which is static, resulting in the card-deck shearing. We saw that upon introducing water at the interface, the card-deck shearing diminishes until it completely disappears. This behavior can be explained by means of the actual values of the LJ parameters. The interlayer interaction between MoS 2 is stronger than that between a layer of MoS 2 and water. Thus, when the top layer slides, the two layers beneath will follow since their adhesion is stronger than that with the water layer. For the low-water coverage regime, the third and fourth layers are still able to interact, albeit less than dry sliding, by means of local bending of the layers around the water cluster. Once the threshold is reached, the layers increase the interlayer distance, instead of local bending, at which point the interaction between water and MoS 2 dominates.
Finally, we comment on the structural behavior of MoS 2 . For the low-coverage regime, we found a structural distortion of MoS 2 , representing a puckered shape (out-of-plane elastic deformation), most strongly pronounced in the layers directly facing the water. The MoS 2 layer bends outward from the center interface, creating a locally enhanced interlayer distance. A similar observation was already visible in the computational setup ( Figure 2). Increasing the number of water molecules increases the intralayer structural distortion up to a threshold, after which MoS 2 no longer bends around the water cluster but increases more significantly the interlayer distance of the center interface. This reveals the delicate interplay between the energy cost of the structural distortion and maintaining the interlayer van der Waals interactions. In other words, the LJ interactions between the layers are replaced by LJ contributions between water molecules and the MoS 2 layers.

Frictional Behavior
Using the same analysis as described in Claerbout et al. (2019), we can make an estimation of the dissipated energy as a result of friction coming from the PES corrugation and other channels such as structural distortions. By considering the instantaneous force that is required to slide all the atoms in the uppermost sulfur layer rigidly with a constant velocity along the sliding direction, we can calculate the dissipated work as follows: where F ext − N i 1 F i,y is the sum of the atomic contributions of the top sulfur layer along the sliding direction. The integral ranges along time τ, the total time of the trajectory, and v is the sliding velocity. In Figure 5A, the work profiles at 0 and 3.0 GPa load are displayed for dry sliding and 8.5% water coverage. Furthermore, in this panel, the dissipation is displayed for the parallel water network in the case of 67% coverage. In Figure 5B, the full coverage with an antiparallel orientation and bulk water, again at 0 and 3.0 GPa load, are displayed. From these graphs, several observations can be made. First of all, an increasing load leads to higher dissipation. This observation was expected, and it can be explained with the same argument used earlier. Increasing the load reduces the interlayer distance, resulting in a more repulsive regime of the LJ interactions and a more corrugated PES. A second observation is the dependence of the energy dissipation on the number of water molecules. Starting from dry sliding, increasing the water coverage results in an increasing dissipation. In other words, water acts as a contaminant that increases the frictional forces. This continues up to 3% coverage, after which the dissipation starts to decrease. Reaching moderate coverage (17%), the dissipation falls below dry sliding where the water starts to act as a friction-weakening lubricant. It should be noted that the dissipation behavior follows the same trends as the card-deck shearing vs. block sliding, with respect to the amount of water coverage. These observations can be explained as follows. For dry sliding and low water coverage, a PES corrugation is present, resulting in a significant dissipation during sliding, as explained above. When the number of water molecules is increased in the low coverage regime, an additional dissipation occurs due to the fact that the number of degrees of freedom and the possible channels for energy dissipation increases. This can, for example, be seen by the puckering of the layers in the low water coverage regime. The water causes a puckered shape, resulting in an enhanced strain in the structure, which comes with an energy cost. Next, the puckered shape promotes a deformation of the neighboring layers during sliding. The PES profiles can also explain the frictional behavior in the high coverage regime. There is no effective corrugation in the PES in the case of an antiparallel water network, compared to that of dry MoS 2 ; thus, when increasing the water coverage, the dissipated energy decreases. However, in the case of a parallel oriented network, the PES corrugation is of the same order as that of dry MoS 2 , which is translated in a similar dissipation behavior.
Considering different velocities, an interesting behavior was observed. In the case of dry sliding and the low-water coverage regime (below 17%), increasing the velocity decreased the dissipation. Above this coverage threshold, increasing the velocity increased the dissipation, with the most pronounced signature for bulk water. To explain this phenomenon, we divided the coverage into three regimes. The first regime, which we termed the contamination regime, is characterized by low water coverage and is governed by stick-slip dynamics, where, as stated above, a higher velocity results in a possible reduced adhesion and a lower probability for the atoms to move to local minimum energy positions. These factors reduce the energy dissipated during the slip phase. An explanation provided by Zeng et al. (2018) on the frictional velocity dependence of graphene studied by using atomic force microscopy, entailed an aging argument, resulting in a stronger interlayer contact state for lower velocities. Although this case might not directly relate to our computational study (Vanossi et al., 2013), it corroborates the fact that for an enhanced slip, more energy is irreversibly and instantaneously transformed into vibrations and heat. We defined the second regime (ranging roughly between 17% and full coverage) as the solidified water lubrication regime. Here, the lubrication type changes: stick-slip is no longer present and the block sliding mechanism is dominant. Despite the smooth dynamics and the seemingly flat PES compared to dry MoS 2 , some corrugation still remains. Here, thermal activation can be exploited to overcome the low energy barriers, reducing the dissipation. However, with increasing velocity, the system might not explore the local energy minima long enough to benefit from thermal fluctuations toward overcoming the barriers (Szlufarska et al., 2008). Moreover, high velocity regimes are in general further away from equilibrium and characterized by more complex and nonlinear processes, and thereby more dissipative. The third regime (which we termed the liquid water lubrication regime) considers bulk water lubrication, where an enhanced signature of increased dissipation was observed. Apart from the discussion above, here, one should take into account the fact that we are no longer only dealing with solids but also with liquids. Whenever the hydrodynamic fluid film creates enough MoS 2 interlayer separation, the MoS 2 -water interaction dominates. Liquids are characterized by viscosity, which is known to be strongly affected by the sliding velocity (Vanossi et al., 2013;Ni et al., 2019). In our case, we found out a Newton-like behavior of the fluid characterized by an almost linear increase in the average velocity of the water molecules when moving from the static bottom three layers to the mobile top three layers of MoS 2 (see the Supplementary Material).
In general, two types of friction can be distinguished, namely, static and kinetic friction. Both types of friction can be calculated via the Amontons' law of friction (Amontons, 1699;Popova and Popov, 2015): where F is the average over the external force in the sliding direction, α 0 is the frictional force at zero load, μ is the coefficient of friction, and F ⊥ is the normal load applied to the system. Although this law is in general valid only at the macroscopic scale, Mo et al. (2009) have shown that in the case of a single asperity contact, it can work at the nanoscale. In the final part of the analysis, we computed the coefficients of friction present during sliding. Since we calculated the coefficients of friction through a linear fit of the frictional forces experienced by the top layer FIGURE 5 | Profiles of the dissipated work for MD trajectories sliding in the x-direction. The left panel (A) represents the case for dry sliding (black and red) and 8.5% coverage (green and blue) at zero and 3.0 GPa load, respectively, and one case with 67% coverage (orange) at 0.5 GPa load which is characterized by a parallel water network. The right panel (B) represents the dissipation at zero and 3.0 GPa load for full coverage with an antiparallel water network (black and red) and bulk water (green and blue).
Frontiers in Chemistry | www.frontiersin.org June 2021 | Volume 9 | Article 684441 (which is in constant motion), we are only dealing here with kinetic friction. To obtain the average lateral force, we used a procedure applied before in one of our other studies (Claerbout et al., 2019), the so-called bootstrap method. More details on this averaging technique can be found in the Supplementary Material.
In Figure 6, we present the kinetic coefficients of friction as a function of the water coverage for all three velocities. In this analysis, we selected only the trajectories in which the water network for the high coverage regime was characterized by an antiparallel orientation. This approach is justified since this orientation is more energetically favorable. Again, three frictional regimes can be distinguished. In the first contaminated regime, the coefficient of friction increases with the water coverage, reaching a maximum at a coverage of 3%, after which the coefficient of friction decreases and drops below the dry sliding value for a coverage of 17%. When the water coverage increases further (solidified water lubrication regime), the coefficient of friction monotonically decreases, eventually falling below the superlubricity threshold (COFs below 0.01) (Baykara et al., 2018;Martin and Erdemir, 2018). After 70% water coverage (liquid water lubrication regime), the friction remains seemingly stable till bulk water; at this point, the water molecules are no longer solely in a solid form, but they are also present as a liquid, bearing the usual velocity dependence (Popov, 2017). We found similar trends for the other velocities, albeit that for the contamination regime; the lower velocity has higher coefficients of friction (and vice versa), in agreement with the earlier discussion about energy dissipation.

CONCLUSION
In this study, we investigated the frictional properties of bulk commensurate MoS 2 in the presence of water. By employing nonequilibrium MD simulations, the effects of load, sliding velocity, and the number of water molecules were elucidated. Several conclusions can be drawn from our simulations. First, when placed within a 2D confinement, water acts as a solidified lubricant. The resulting ice-like structure, where the H bonds lead to a characteristic tile-like network, has a discrete (e.g., parallel or antiparallel) orientation above/beneath MoS 2 . Second, two classes of sliding dynamics were distinguished, namely, card-deck shearing vs. block sliding and stick-slip dynamics vs. smooth dynamics. Both show a dependence on the amount of water coverage and the sliding velocity. Within the low water coverage regime, card-deck shearing and stick-slip dynamics are dominant. These types of dynamics are characterized by high dissipative signature and strong load dependence. For moderate to high water coverage, block sliding and smooth dynamics were observed with low energy dissipation and low load dependence. One exception is the metastable case, where the water network aligns along the sliding direction. In this case, a behavior comparable to that of low water coverage is observed. Combining the results above, we found that three lubrication regimes could be identified. In the first one (contamination regime and low water coverage), water acts as a contaminant and friction strengthening was observed with increasing coverage. In the second one (solidified water lubrication regime and high water coverage), water acts as a solid-like lubricant, and superlubricity was observed. In the third one (liquid water lubrication regime and bulk water), water acts as a liquid Frontiers in Chemistry | www.frontiersin.org June 2021 | Volume 9 | Article 684441 lubricant, and low coefficients of friction were found, also showing a more pronounced viscosity-like velocity dependence. Finally, we showed and quantified that due to the solidified structure of water, there was a rotational dependence present for the frictional behavior of water between a multilayer system of MoS 2 . Our results directly contribute to strengthen the general understanding of the nanoscale frictional mechanisms at play when considering MoS 2 (and 2D materials in general) in the presence of humidity. Once more, it was confirmed that the processes defining friction at this scale are anisotropic and complex in nature. Our results shed light upon the design of nanoscale devices, such as nano-and micro-electromechanical systems (NEMS/MEMS), in which humidity might play an important role, to ramp up the efficiency of the devices. In future work, the scalability of our results to the multi-asperity macroscale will be studied.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
VC carried out the conceptualization, methodology, validation, calculation, formal analysis, investigation, visualization, and writing of the original draft. PN supervised the study and contributed to the conceptualization, methodology and validation, and reviewed and edited the manuscript. TP guided the project and reviewed and edited the manuscript. The funding was secured by PN and TP.