1 Introduction

Friction properties of lubricated contacts change with pressure and with external electric fields. The ability to electrotune a liquid’s lubricity, i.e. to change its frictional behavior by applying external electric potentials, offers promising applications in diverse fields such as soft robotics, microfluidics or actuation [1,2,3,4]. Over the last years, ionic liquids (IL) have been increasingly studied as candidates for electrotunable lubrication, which can be attributed to their temperature stability and large electrochemical window [4,5,6,7,8].

Aqueous lubricants, on the other hand, experience broadening use due to their environmental friendliness and biocompatibility. As water is relatively volatile and reactive, aqueous lubricants are found mainly in food processing, microfluidics, or in biomimetic systems [4, 9]. Natural water-lubricated systems can exhibit exceptionally low coefficients of friction (COF) of \(\mu <0.002\) [10]. Experiments show that the COF for aqueous solutions on gold electrodes can also be electrotuned [11,12,13,14]. Prior aqueous friction studies have used different experimental techniques, including the surface-force apparatus [11, 15, 16], atomic-force [12] and colloidal-probe [13, 14] microscopy.

With macroscopic friction being governed by nanoscale processes, molecular dynamics (MD) is a valuable tool for insights in electrotunable lubrication. MD simulations from the last decade have, for instance, investigated the influence of molecular structure [7, 17], sliding velocity [18, 19] or surface roughness [20] and charge [8] on the friction properties of some ILs. For simulating voltage effects, such studies mainly use the constant charge (CC) method, where each electrode atom is assigned a fixed charge value over the course of the simulation and image charge effects are neglected. Some equilibrium MD studies suggest only a minor influence of this artifact on ion adsorption structure or ion diffusion [21, 22]. For non-equilibrium phenomena, in particular interfacial friction, other authors argue that the way in which charges are distributed on electrodes has a ’dramatic impact’ [23,24,25]. It is therefore desirable to use a more realistic way for simulating atomic charges on electrodes in MD friction experiments.

We use MD to investigate the electrotunable friction of canonical aqueous systems, namely pure water and a dilute aqueous NaCl electrolyte, at various external pressures. Using a constant potential (CP) method, atomic charge values in the metal are updated at each timestep and, therefore, respond to their environment. This approach yields more physical results than the CC method when modelling electrodes [26,27,28,29,30].

We use a model nanoscale system consisting of two parallel interfaces separated by the lubricant for our friction computer experiments (see Figs. 1, 2). Hydrodynamic lubrication between atomically flat surfaces is governed by fluid viscosity \(\eta\) and interfacial friction coefficient \(\lambda\). Using slip velocity \(v_\mathrm{{s}}\) and shear rate \(\dot{\gamma }\) (see Fig. 1), these quantities can be obtained from an analysis of the velocity profile inside the gap [31] as

$$\begin{aligned} \eta&= \tau / \dot{\gamma } \, \, \text {and} \end{aligned}$$
(1)
$$\begin{aligned} \lambda&= \tau / v_\mathrm {s}\,, \end{aligned}$$
(2)

where \(\tau\) is the shear stress of the system. Based on our simulation data, we isolate and quantify the individual effects of pressure and voltage on \(\eta\) and \(\lambda\). This will lead to a model which describes \(\tau\) as a function of pressure and voltage.

Fig. 1
figure 1

Idealized velocity profile with fluid slip for the present setup. The lower (upper) electrode with total negative (positive) charge has a prescribed velocity of \({u^{-}=0}\,{(u^+=20\,\mathrm {m}\,\mathrm {s}^{-1})}\). A Couette-type velocity profile with shear rate \(\dot{\gamma }\) establishes across the gap. Slip velocities \(v_\mathrm {s}^{\pm }\) correspond to the velocity differences between the (positively/negatively biased) wall and its closest fluid layer

Full size image

2 Methods

Our model system consists of two parallel interfaces separated by the lubricant, corresponding to an infinite electrochemical parallel plate capacitor with moving interfaces. A snapshot of this system is shown in Fig. 2. Two Au(111) electrodes are stacked in z-direction, with fluid confined in the gap between upper and lower electrode. The simulation cell is periodically replicated in the xy-plane. When shearing the capacitor plates relative to each other with a velocity u, interfacial friction (quantified by \(\lambda\)) and the fluid’s shear resistance (viscosity \(\eta\)) cause a lateral force \(F_{x}\) in the electrode slabs. Each electrode has three regions: the interfacial atom layer that is in direct contact with the electrolyte, a rigid outer layer that is either at rest (bottom) or moved at a constant velocity (top) and three thermalized atomic layers in between the former two. An electric potential difference \({\Delta \Phi }\) between the gold electrodes causes water molecules to change their (preferred) orientation and ions to move towards their respective counter-electrode.

Fig. 2
figure 2

Model capacitor system with five atomic layers of gold for each electrode. A normal force acting on the top electrode allows pressurization. A voltage difference is achieved by offsetting the electronegativity values of the electrode atoms within our charge optimization approach. The gap is filled with SPC/E water, Na and Cl ions. Temperature is maintained via the three inner layers of each electrode. Prescribing a relative electrode velocity \(u=u^+-u^-\) introduces a (voltage- and pressure-dependent) shear force into the system. Simulated systems are larger than the one depicted here

Full size image

We note that contrary to “grand-canonical” schemes with an explicit reservoir for the lubricant (see e.g. Refs. [7, 17, 20, 32]), the film thickness or number of lubricant molecules is a parameter to these calculations. However, most prior calculations with a reservoir consider parallel surfaces that cannot sustain hydrodynamic lift. Such calculations are therefore suitable to understand squeeze-out phenomena and are unsuitable to compute the film thickness in a macroscopic tribocontact. A full elastohydrodynamic [33] treatment of the contact is necessary for a computation of the film thickness, that could then serve as an input to molecular calculations of the type presented here.

Interactions between gold atoms are modelled using the embedded atom method (EAM) potential by Grochola et al. [34]. Water-gold interactions follow a Morse-type potential model parameterized by Berg et al. [35]. All remaining atoms interact via Coulomb and Lennard-Jones (LJ) interactions using the parameters given in Table 1. LJ cross-interactions are found using Lorentz-Berthelot mixing [36]. Note that LJ parameters for gold are used only for calculating cross-interactions of sodium and chloride ions with gold. Ion charges are scaled by 0.75 to account for electronic polarization effects in water [37]. Scaled ionic charges necessitate the use of adapted LJ parameters [38].

Table 1 Lennard-Jones and Coulomb parameters used
Full size table

We use the SPC/E water model [40], which we chose over the conceptually identical TIP3P model [45] because the latter strongly underestimates solubility of NaCl [46, 47]. Water molecules are kept rigid using the RATTLE algorithm [48].

At each timestep, we recalculate charges on the gold atoms by minimizing a global energy functional. The particular flavor of the CP method used by us is known as the electronegativity equalization method (EEM) [49,50,51,52]. We solve the governing system of linear equations [53] with a relative residual of \(10^{-6}\) using the conjugate gradient solver described by Aktulga et al. [54]. The atom-specific values needed for the EEM algorithm, i.e. electronegativity, Coulomb self-energy and Coulomb shielding constant, along with EEM-related Coulomb cutoff radii, were taken from Monti et al. [55]. A voltage \({\Delta \Phi }\) between the capacitor plates is simulated by offsetting electronegativity values of the gold atoms for the EEM by \(\pm \, e\,{\Delta \Phi }/2\) [56, 57]. We note that this scheme treats the electrodes as perfectly metallic, i.e. the density of states is infinite and an arbitrary number of electrons can be injected into the electrode at the Fermi-level at no cost, while real materials have a finite density of states. This effect may, however, only matter for low-dimensional systems [26, 58, 59]. We note that this scheme does not capture redox reactions within the solvent at the electrodes. However, all our simulations are carried out for bias voltages within the window of electrochemical stability of water.

The evaluation of electrostatic potentials within the EEM uses a switching scheme [60, 61] with an inner and outer cutoff of 0.0 and \(1.2\,\mathrm {nm}\), respectively, that cuts off the Coulomb contribution to the functional. This significantly speeds up the evaluation of the Coloumb interaction but introduces errors. We here use this tapered Coloumb interaction only during charge optimization but revert to the full long-ranged interaction evaluated with the particle-particle particle-mesh (PPPM) method [62] when evaluating the forces. We need to use this mixed scheme because a full PPPM evaluation of the electrostatic potential in EEM is computationally prohibitive. To account for the slab geometry of the system, the Yeh/Berkowitz correction for the Coulomb forces is used within PPPM [63].

For optimized charges, forces can be obtained straightforwardly from the spatial gradient of the total energy by virtue of the Hellmann–Feynman theorem [64]. Since our EEM does not optimize the PPPM total energy expression, we generally expect that forces and energies are not consistent, leading to heating of the system. We only see heating during the initial stages of our simulation where ionic positions have not yet equilibrated. This indicates that electrostatic interactions are effectively screened once the electrochemical system is near thermodynamic equilibrium and the error imposed by our mixed cutoff procedure is small.

We would like to point out that our choice of force-field starts from well-established parametrizations of as few simple components as possible together with one particular modification that has not been explored in this context yet, namely the EEM. Drawing upon commonly accepted water, ion, and gold models that offer decades of research into their various properties provides us with a certain confidence in their specific realms of applicability. Morse interaction between fluid and substrate better reproduces interfacial behavior (e.g. water dipole orientation) [35]. Charge scaling on the ions is a simple but effective method to improve ion solubility and dynamic behavior [37, 65].

A periodic simulation cell with approximate dimensions 18 × 18 × 7.5 nm3 is filled with an fcc lattice of Au atoms equilibrated at standard conditions. Removing the inner Au atoms exposes atomically flat Au(111) surfaces on the upper (positive total charge) and lower (negative total charge) electrodes. Each electrode has 22,320 gold atoms in five atomic layers. The gap is filled with SPC/E water molecules and equal numbers of sodium and chloride ions. At a concentration of 0.1 mol kg−1, this translates to 55,428 water molecules and 102 ion pairs. A short energy minimization is then performed to avoid atom overlap at this stage.

The top shear layer has a prescribed velocity of u = 20 m s−1, while the bottom shear layer remains static. We impose pressure by adding an external force to the topmost layer, but additionally increase its total mass and add a damping force as described in Ref. [66]. Temperature control is implemented by a Langevin thermostat [67, 68] applied to the three thermalized layers of each electrode with a target temperature of \(300\,\mathrm {K}\) and a time constant of \(200\,\mathrm {fs}\). Atoms are thermalized only in the y-direction to avoid contamination of the measured friction force with the viscous damping intrinsic to this thermostat.

3 Results

We now describe the results from our sliding simulations, run parametrically for three different voltages (\({{\Delta \Phi }=0\text {,}\,0.6~\text {and}~ 1.2\,\mathrm {V}}\)), 6 normal stresses (\({\sigma =0.1}\) ... \({1.0\,\mathrm {GPa}}\)) and 2 ion concentrations (\({\text {[NaCl]}=0}\) and \({0.1\,\mathrm {mol}\,\mathrm {kg}^{-1}}\)). We first equilibrated each configuration for \(5\,\)ns at the given voltage and pressure. Subsequent production runs lasted another \(2\,\mathrm {ns}\).

Figure 3 shows (time-averaged) data points for shear stresses \(\tau\), along with the model functions for \(\tau (\sigma , {\Delta \Phi })\) developed in Sect. 4. With the present system reflecting the hydrodynamic lubrication regime, the COF (defined as the slope in a \(\sigma (\tau )\)-plot [69]) is quite low, i.e. \({\mu \approx 10^{-3}}\). Although (bulk) SPC/E water has been shown to remain liquid for the temperature and pressure values used in this work [70], we suspect the onset of crystallization to be responsible for outliers at the highest pressure and voltage values probed, most visible in Figs. 3 and 7.

Fig. 3
figure 3

System shear stress \(\tau\) vs. normal stress \(\sigma\) for a pure SPC/E water and b the dilute electrolyte. Data points are given as averages over time and both electrode slabs. Error bars, only given for the extreme pressure values, show the standard error of the mean. Lines represent the model functions for \(\tau (\sigma , {\Delta \Phi })\) developed below

Full size image

Shear stresses are primarily affected by the applied normal stress, with the influence of voltage being in the range of the measurements’ standard deviations. However, particularly in the electrolyte case, a small influence of voltage on friction can be observed. This influence will become clearer once the contributions from fluid viscosity and interfacial friction are separated.

The average spatial distribution of atom species Cl, Na, O are shown in Fig. 4 for selected pressure values, along with the velocity profiles \(v_{x}(z)\) of the confined fluid. With no voltage applied, ions tend to be found away from the surface. Some Na atoms (partly) shed their solvation shell and adsorb on the gold electrodes. As voltage increases, many Na atoms squeeze in between the negatively polarized electrode and the first water layer, forming the so-called inner Helmholtz plane (IHP). In the high-voltage case, local Na concentration in the IHP reaches approximately 6 mol kg−1. Cl ions, on the other hand, remain fully solvated and are found mainly beyond the first two water layers. The anion distribution is therefore more spread out inside the gap, with the peak Cl concentration being roughly 10 times its average value. For both ion types, distributions away from the surface resemble the exponential decay expected from Gouy-Chapman theory [71].

Fig. 4
figure 4

Average distributions for ions and water oxygen atoms at a–d the lowest and e–h the highest pressure values probed. Shown are a, e the \(\hbox {Cl}^-\) concentration, b, f the \(\hbox {Na}^+\) concentration, c, g the oxygen density and d, h the fluid velocity. Except for panels d, h, curves for different voltage values are shifted for better visibility and normalized to the average atom concentrations. Insets in panels b and f show same data with different axes

Full size image

Beyond the first dense fluid layers, a Couette velocity profile establishes inside the gap. Nonlinearities in the velocity profiles arise at higher voltages, reflecting the slanted distribution of ions. In the no-voltage case, slip at both interfaces is roughly equal. With increasing voltage, the sum of both slip velocities (therefore mean fluid velocity) remains more or less unchanged. However, the ratio of \(v_{s}^{+}\) and \(v_{s}^{-}\) changes, with slip at the Na/Au interface reducing and slip at the Cl/Au interface growing.

The slip velocity \(v_\mathrm {s}\) corresponds to the difference in speed between a solid wall (\(u^{-}=0,\,u^{+}=20\,\mathrm{m}\,\mathrm{s}^{-1}\)) and the closest fluid layer. Our protocol for determining slip velocities from time-averaged velocity profiles of the water molecules, as shown in Fig. 4 (d, h), is as follows: We neglected contributions from the fringes of the solid-liquid interface, i.e. we started measuring \(0.05\,\mathrm{nm}\) after the first non-zero value for \(v_{x}(z)\) on either side. Slip velocities were then calculated from spatial average over the next \({\Delta } z = 0.1\,\mathrm {nm}\). \(v_{s}^{+}\) for the positively polarized top interface and \(v_{s}^{-}\) for the negatively polarized bottom interface are shown in Fig. 5. In general, slip decreases as pressure increases. While slip velocities of both interfaces are similar at zero voltage, increasing \({\Delta \Phi }\) separates the slip velocities for both interfaces. This effect is stronger in the electrolyte case.

Fig. 5
figure 5

Slip velocities \(v_\mathrm {s}\) extracted from velocity profiles for all pressure and voltage values for a pure water and b the dilute electrolyte. Markers of the same color indicate top and bottom electrode at the same voltage difference. Slip velocities decrease as pressure increases. With \({\Delta \Phi }=0\), slip velocities \(v_\mathrm {s}^{-}\) for the bottom and \(v_\mathrm {s}^{+}\) for the top interface are nearly identical. As voltage increases, the difference between these slip velocities grows

Full size image

For determining fluid viscosities \(\eta\) according to Eq. (1), we time-averaged the z-dependent lateral velocities \(v_{x}(z)\) of all fluid atoms for each parameter set. We then fitted linear functions to the velocity profiles (using only the inner third of z-values, thus neglecting nonlinearities in the velocity profiles that are due to interfacial processes) to determine the shear rates \(\dot{\gamma }\). These results are shown in Fig. 6 and indicate that \(\eta\) is mainly influenced by external pressure \(\sigma\), and not by the applied voltage \({\Delta \Phi }\). The relationships \(\eta (\sigma )\) for both ion concentrations were (heuristically) fitted with a quadratic law.

Fig. 6
figure 6

Viscosities determined for all configurations using Eq. (1). No systematic trend is visible for different voltages and ion concentrations, only for external pressure \(\sigma\). Viscosity was therefore (heuristically) fitted to a quadratic function of \(\sigma\). Error bars for fluid viscosity represent the standard error of the measurement. The main uncertainty comes from fluctuations of the shear stress. Error bars for normal stresses are smaller than the symbol size. We only show the error bars exemplarily for NaCl at 1.2 V. Errors for the other data points are comparable

Full size image

As normal stress \(\sigma\) increases, so does the system’s shear stress \(\tau\), while slip velocities \(v_\mathrm {s}^{\pm }\) decrease. Therefore, the interfacial friction coefficient \(\lambda\) defined in Eq. (2) increases with pressure. A potential difference \({\Delta \Phi }\) across the interfaces barely changes the sum of both slip velocities, which, through Eq. (1), is consistent with the apparent lack of voltage influence on viscosity as depicted in Fig. 6. \({\Delta \Phi }\) does, however, change the ratio of the two slip velocities. This means that the interfacial friction coefficients \(\lambda ^{\pm }\) are voltage-dependent. Data shown in Fig. 5 suggests that the effects of \(\sigma\) and \({\Delta \Phi }\) on \(\lambda\) are independent of each other. We therefore assume that two contributions can be factorized, such that

$$\begin{aligned} \lambda (\sigma , {\Delta \Phi }) \approx \lambda _\mathrm {ref} \, \varLambda _{\sigma }(\sigma) \, \varLambda _{\Phi }({\Delta \Phi })\, . \end{aligned}$$
(3)

Thus, measured values for \(\lambda\) are shown in two different ways. First, data series for each voltage value are normalized with respect to the values at some fixed \(\sigma _\mathrm {ref}\). The influence of \(\sigma\) on \(\lambda\), termed \(\varLambda _{{\sigma }}\), is then found from interpolation between these normalized data series (Fig. 7). Similarly, normalizing \(\lambda (\sigma , {\Delta \Phi })\) with respect to the values at some fixed \({\Phi _\mathrm {ref}}\) will lead to \(\varLambda _{\Phi }\) (Fig. 8). We chose \(\sigma _\mathrm {ref}=0.6\,\mathrm{GPa}\,\) and \({\Phi }_\mathrm {ref}=0\,\mathrm{V}\) as reference values.

Fig. 7
figure 7

The effect of pressure \(\sigma\) on interfacial friction \(\lambda\) for all configurations. Each data point is normalized to \({\lambda (\sigma _{\mathrm {ref}}=0.6\,\mathrm {GPa})}\) at the corresponding voltage and interface. \(\varLambda (\sigma )\) appears linear in \(\sigma\) below 1 GPa. Crystallization is suspected for highest pressure and voltage values

Full size image
Fig. 8
figure 8

The effect of potential \({\Phi }\) on interfacial friction for all configurations. Data is shown separately for a pure water and b the dilute electrolyte. Data points are normalized to \({\lambda ({\Phi }=0\,\mathrm {V})}\) at the corresponding pressure and interface. The data is consistent with a quadratic relationship between \({\Phi }\) and \(\varLambda\). The difference of \(\varLambda (\Phi )\) between the pure water and the electrolyte systems are more prominent at the negatively polarized electrode, suggesting a strong influence of the Na ions adsorbed to the interface

Full size image

Based on the available data, \(\varLambda _{\sigma }\) is linear in \(\sigma\)—at least until the onset of water crystallization at very high pressure [70]. With regards to \(\varLambda _{\Phi }\), the data is less conclusive. Theoretical considerations suggest three contributions to \(\lambda\): density and structure factor in the first fluid layer, acting linearly on \(\lambda\), as well as a lateral energy barrier term, which reflects fluid-substrate corrugation and influences \(\lambda\) quadratically [72,73,74,75]. In a first-order approximation, applied electric potentials are here assumed to linearly influence only this lateral energy barrier. Therefore, \(\varLambda _{\Phi }\) is modelled as being quadratic in \({\Phi }\). Consequently, the grey lines in Fig. 8 are quadratic fits to our data.

For both pure water and the NaCl electrolyte, \(\varLambda _{\Phi }\) increases for \({\Phi }<0\) and decreases with \({\Phi }>0\), meaning that the voltage effect on interfacial friction is higher for the negatively polarized interface which attracts sodium ions and hydrogen atoms of the water molecules.

4 Discussion

For the parameter space probed here, pressure \(\sigma\)—through its effect on bulk fluid viscosity \(\eta\)—clearly has a larger influence on lubrication than applied voltages (cf. Fig. 3). This is because our gap widths were around \(4\,\)nm throughout, and therefore well above surface separation values for which interfacial phenomena (e.g. hydration lubrication) in such systems dominate [15, 16]. Even the most carefully crafted smooth gold surfaces will show roughness comparable to our gap width [11]. In line with previous studies, we would expect different results (in particular a stronger influence of voltage on friction) as less smooth surfaces are introduced [17, 76, 77].

Nevertheless, this idealized setup has enabled an analysis of wall slip, which generally fades with increasing surface roughness [78]. From analyzing slip velocities across all parameters, we have found that the effects of pressure \(\sigma\) and voltage \({\Phi }\) on interfacial friction coefficient \(\lambda\) are in a similar range. Separating \(\varLambda _{\sigma }\) and \(\varLambda _{{\Phi }}\) shows a (roughly) linear influence of \(\sigma\) and quadratic influence of \({\Phi }\) on \(\lambda\), which is in agreement with theory [73, 74] and other MD results [8, 25].

For pure SPC/E water, the difference between \(\lambda ^{+}\) and \(\lambda ^{-}\) for the two interfaces with applied voltage can be traced back to the preferential orientation of water molecules at the interface (see Fig. 9). With no potential bias, hydrogen atoms in our simulations tend to face towards the electrodes. The creators of our Au-water potential reported the same behavior [35], which is also backed by first principles calculations [79]. When voltage is applied, this trend increases for the negatively polarized interface (\(\lambda ^{-}\uparrow\)) and reverses at the opposite electrode (\(\lambda ^{+}\downarrow\)). We hypothesize that rotation of the water molecules affects the commensurability between Au substrate and contact layer, thus increasing or decreasing the lateral interaction energy barrier and thus interfacial friction [74, 80]. For the NaCl electrolyte, these effects are even larger and likely due to the different size of the two ions.

Fig. 9
figure 9

Mean dipole orientation \(\langle \cos {\Theta }\rangle\) across the nanogap for the pure water case. See inset for definition of the angle \({\Theta }\). Spatial structure in \(\langle \cos \,{\Theta }\rangle\) disappears towards the center of the gap. With no applied potential, hydrogen molecules tend to face towards both Au(111) surfaces. As voltage increases, this trend increases (decreases) at the negatively (positively) charged electrode

Full size image

The effect of \({\Delta \Phi }\) on \(\lambda\) is larger for the electrolyte case, which must be attributed to the ions at the interface. With electrostatic bias enabled, sodium ions tend to (partly) leave their solvation shell behind and leave the first water layer towards the electrode, see Figs. 2 and 4. For chloride ions, on the other hand, this event occurs less frequent. Momentum transfer between electrode and bulk fluid is facilitated through this first layer of sodium atoms.

A friction model for the present setup can be formulated as follows. Shear stress \(\tau\) is related to \(\eta\) and \(\lambda ^{\pm }\) via Eqs. (1) and (2). A fourth equation between these quantities is gained from the (ideal) relationship between shear velocity u, slip lengths \(v_\mathrm {s}^{\pm }\) and slope of the Couette velocity profile \(\dot{\gamma }\) inside a gap of height h (cf. Fig. 1). This yields

$$\begin{aligned} u=v_\mathrm {s}^{-} + \dot{\gamma }\, h + v_\mathrm {s}^{+} = \tau \left( \frac{1}{\lambda ^{-}(\sigma , {\Phi })} + \frac{h(\sigma )}{\eta (\sigma )} + \frac{1}{\lambda ^{+}(\sigma , {\Phi })} \right) \,. \end{aligned}$$
(4)

We now also need constitutive expressions for the gap height h, viscosity \(\eta\) and interfacial friction \(\lambda ^{\pm }\). Assuming that h only depends on normal stress \(\sigma\) through elastic deformation (compressibility of the water) yields

$$\begin{aligned} h(\sigma )=h_\mathrm {0}+\beta \, \sigma \,. \end{aligned}$$
(5)

We further make the ansatz that viscosity \(\eta\) is a quadratic function in \(\sigma\) [81]:

$$\begin{aligned} \eta (\sigma ) = \eta _\mathrm {0} + \alpha \, \sigma ^2\,. \end{aligned}$$
(6)

Finally, we need an expression for the interfacial friction \(\lambda\) that depends linearly on normal stress \(\sigma\) and quadratically on applied voltage \({\Delta \Phi }\) (see Ref. [74]), yielding

$$\begin{aligned} \lambda ^{\pm }(\sigma , {\Delta \Phi }) \approx \lambda _\mathrm {ref}^{\pm } \, \underbrace{ \left( 1+\zeta \, {\Delta } \sigma \right) }_{\varLambda _{\sigma }} \, \underbrace{\left[ a\left( \pm {\Delta \Phi }/2-x_0\right) ^2+y_0\right] }_{\varLambda _{\Phi }}\,. \end{aligned}$$
(7)

In these expressions, \(\eta _{0}\) and \(h_{0}\) are viscosity and gap height values, respectively, extrapolated to zero pressure. \(\lambda _\mathrm {ref}\) is the interfacial friction coefficient at \({{\Delta \Phi _{\mathrm{ref}}}=0\,\mathrm {V}}\) and \({\sigma _\mathrm {ref}=0.6\,\mathrm {GPa}}\). \({\Delta }\sigma =\sigma - \sigma _\mathrm {ref}\) denotes the difference of pressure to our reference pressure \(\sigma _\mathrm {ref}\). \(\alpha ,~\beta ,~ \zeta ,~a,~x_0,~y_0\) are fitting parameters, of which \(\beta\) can be related to the fluid’s compressibility. Solving Eq. (4) for \(\tau\) and inserting these independently obtained parameters leads to the functions for \(\tau (\sigma ,{\Phi })\) shown in Fig. 3.

We note that this model is of an empirical quality and should not be relied on outside the range of parameters investigated here. While some of the assumptions are backed up by theoretical considerations (e.g. the quadratic dependence on voltage), others (such as the variation of height with pressure) simply describe the linear response of the system that will surely break down for larger perturbations.

5 Summary and Conclusions

At molecular scales, the way charges are distributed on an electrode has a large effect on interfacial friction [24, 25]. This emphasizes the need for realistic methods to modelling  electrode charges beyond a simple constant charge approach. We here allow electron charges to react to the dynamics of the electrolyte [27,28,29,30], yielding a constant potential approach. This work has explored how such constant potential methods can be used for molecular simulations of electrotunable lubrication.

Using an atomically flat plate capacitor with nanoconfined lubricant in its gap, we probed the effect of pressure and voltage on hydrodynamic aqueous lubrication with different salt concentrations. We found that the interfacial friction coefficient \(\lambda\) behaves differently at anode and cathode. In agreement with theoretical analyses [73, 74] and other MD results [8, 25], \(\lambda\) is influenced linearly by external pressure and approximately quadratically by external electrostatic bias. We finally developed and tested an expression for shear stress \(\tau\) that incorporates the effects of pressure and voltage on viscosity and wall slip.