The influence of water polarization on slip friction at charged interfaces

Albert Szent-Györgyi called water “the matter and matrix, mother and medium of life,”¹ echoing water’s active role in various biochemical and physiological processes.² Water exhibits unusual and essential physical and chemical properties, viz., the ability to form strong hydrogen bonds, large dipole moment, and high dielectric strength. Water is the natural lubrication medium in living systems, and the friction in biological environments is particularly relevant.^3,4 The peculiar property of water, which is special to friction, is its ability to retain a bulk-like fluidity when confined to ultra-thin films.^5,6 Most non-associating liquids, including organic solvents, tend to solidify, and the viscosity diverges when confined to films of thickness 5–8 molecular layers.^7,8 This is because the translational entropy available to the molecules of non-associating liquids decreases with the increase in confinement.⁹ Meanwhile, water is unique, and the solidification of water is associated with forming a highly directional hydrogen-bonded network, which is suppressed with thin-film confinement.⁶ This enables water to retain its viscosity when confined to films of thickness as narrow as 2–3 nm. Moreover, the viscosity of water remains roughly three times the bulk value, even confined to a mono-molecular layer thickness.⁶ Recent molecular dynamics (MD) studies have revealed a spatially varying viscosity of confined water due to changes in its structural and dynamical properties caused by the interaction forces exerted by the confining surfaces.^10–13 While the spatial inhomogeneities of the transport properties occur in confined systems, calculating meaningful transport coefficients, especially at extreme confinements, remains a challenge.^14–16 In such situations, viscosity can be approximated as an effective viscosity, while in reality, it is a strongly nonlocal function, which depends on both strain rate and density and their associated couplings.^17–21 By “effective viscosity,” we mean a pore-averaged viscosity, rather than a position-dependent viscosity incorrectly defined by Newton’s local law of viscosity. To date, exact calculations of the nonlocal viscosity for highly confined fluids remain elusive due to their complexity,¹⁷ although they have been performed for idealized systems.^20,21

The peculiar nature of water to retain its bulk property when confined to ultra-thin films spurs increased interest in the field of nanofluidics—the study of fluid transport at nanoscale dimensions (0–100 nm). The wall–fluid interaction dominates in the nanofluidic system due to the high surface-to-volume ratio, making the wall effects increasingly important. One notable wall effect is the occurrence of a non-zero slip velocity near the wall, a phenomenon rarely observed for macroscale fluid flow.^22,23 For the macroscopic length scale, slip velocity, the difference in the velocity between the fluid and the wall, is negligible compared to the maximum velocity of the fluid within the channel. Hence, it can be assumed to be zero for simplicity. However, the slip velocity can be a significant fraction of the fluid’s maximum velocity in narrow channels at the nanometer scale. Therefore, a proper understanding of the fluid behavior adjacent to solid boundaries is important when dealing with nanoscale fluid flows. The fluid slip at the surface is described using a material property called the slip friction coefficient (ξ). Navier²⁴ defined the slip friction coefficient as a proportionality coefficient that relates the shear stress at the wall (σ_xz) and slip velocity (u_s), given by σ_xz = ξu_s. Here, the fluid flows in the x direction over the solid boundary with the unit normal pointing in the z direction (i.e., into the liquid).

A fundamental notion in nanoscale fluid flow is that the properties of the nanofilm and the surface material entirely determine the interfacial friction. However, the new avenues to control friction rely on reversibly and remotely controlling the fluid–solid interaction through non-invasive mechanisms. One such non-invasive mechanism, controlling the friction between shearing surfaces using an electrical potential, has been known since the pioneering work of Edison.²⁵ He observed a change in the frictional character of the lubricated contacts when a potential difference was applied across them. Historically, lubricants based on aqueous salt solution received the most research attention due to their high electrical conductivity.^26,27 The electrostatic fields emanating from the surfaces influence the electrolyte ions, thereby controlling the structure and dynamics of the lubricant. A more detailed discussion regarding various approaches to the effect of electrical potential on friction is reviewed.²⁸ The electrical control of liquid–solid friction has been studied extensively for hydrophobic surfaces.^29,30 The initial attempts to model the electrical control of friction neglected the pure electrostatic contribution to friction, assuming that the electric contribution is small compared to the non-electric contribution.²⁹ Later, Xie et al.³¹ reworked this model and included the effects of pure electrostatic contribution to friction. The improved model overcomes the previous model’s²⁹ limitations by allowing for a broader range of surface charge densities and clarifying the dependency of liquid–solid friction on the wall and liquid atomic sizes. Recently, Siedl et al.³² studied the wall friction dependency on the electrostatic potential difference between water-lubricated gold electrodes and found a different behavior of wall friction at the anode and cathode. According to their study, the increase in the voltage between two planar surfaces results in two separate interfacial friction coefficients. However, they employed the non-equilibrium molecular dynamics (NEMD) method to calculate the interfacial friction coefficient, which is known to have limitations in high slip systems.³³ This emphasizes the need for further improvement in this area of study.

Hansen et al.³⁴ developed a theoretical model to estimate the slip friction coefficient using equilibrium molecular dynamics (EMD) simulations. They calculated the friction coefficient utilizing the autocorrelation function of the center of mass velocity of the slab of particles near the wall (slab velocity) and the cross correlation function of the slab velocity and the force due to wall–slab interactions. This approach differs from the previous EMD method,³⁵ which determines the friction coefficient using the time integral of the autocorrelation function of the tangential force exerted on the confined fluid by the wall. The advantage of using the Hansen et al.³⁴ method is that it does not suffer from the so-called plateau problem encountered in the previous EMD method.^35–37 However, their method sometimes results in poor statistical averaging due to the limited number of particles in the slab region and becomes statistically less reliable at higher correlation lag times.

In this context, we carried out molecular dynamics simulations to study the electrical control of friction of aqueous electrolytes confined in a graphene nanochannel. We adopted a fixed charge method for simulating the electrostatic potential difference between the nanochannel walls, resulting in an asymmetrically charged system. We then evaluated the slip friction coefficient at both interfaces of the nanochannel, utilizing the EMD method.³⁸ 

Figure 1 shows the simulation domain, where water and alkali chloride solutions are confined within the graphene nanochannel. The solution consists of 1662 water molecules and 64 ions, with an equal number of cations and anions. Each wall of the nanochannel consists of three layers of graphene sheets stacked together. The surface area of each graphene sheet is 3.5 × 3.5 nm². The width (L) of the nanochannel refers to the distance between the centers of mass of the inner graphene layers (or first layer). The outermost graphene layers (or third layer) are frozen throughout the simulation and were added to confine the system. The intralayer carbon–carbon interactions were modeled using the optimized Tersoff potential.⁴⁹ The interlayer carbon–carbon interactions were modeled using the LJ parameters of Hummer et al.⁴⁶ 

The boundary conditions employed were periodic in the x and y directions, while they were fixed in the z direction. The EMD simulations were performed in the NAP_zzT ensemble. This means that the number of particles N, surface area A, temperature T, and normal pressure P_zz were all held constant. The normal pressure was applied by exerting a distributed normal force on the wall atoms. The truncation of LJ forces affects different liquid properties to varying extents. For instance, the density of water at the center of the channel was slightly lower (987 kg/m³) compared to the density of a homogeneous system under a pressure of 1 atm, which includes long-range corrections (998 kg/m³). In the pre-equilibration runs, the (uncorrected) normal pressure was adjusted to 1 atm by utilizing one of the walls as a piston and monitoring its average equilibrium position. We refer to this pressure as “uncorrected” to distinguish it from the pressure obtained by adding long range corrections to the direct pressure due to truncated LJ forces. This correction is only suitable for spatially homogeneous systems at equilibrium. Following this, we fixed the fluctuating wall at its average equilibrium position for the subsequent equilibration and production run. The mean width (L) of the nanochannel varies from 4.0 to 4.1 nm across different systems examined, which is ∼14 molecular diameters of water. To ensure isothermal conditions, the inner graphene layers (first and second layer) were thermostated using a Langevin thermostat⁵⁰ at a temperature of 300 K. The liquid is unthermostated since it has been found that thermostating confined fluids can lead to simulation artifacts.⁵¹ While not shown, we checked that the temperature profile of water indeed remained constant at 300 K. Throughout the thermostating process, the center of mass position of the inner graphene layers was consistently constrained to its initial position at each time step, maintaining a constant channel width. The duration of the simulations was 4.0 ns equilibration run followed by 2.0 ns production run. During the production runs, the data for post-processing were extracted. In order to obtain better statistics, a total of 50 independent simulations were run, with each simulation starting from different initial velocities.

In an asymmetrically charged system, the graphene nanochannel walls are charged oppositely and homogeneously such that a uniform electric field is induced between them. Each carbon atom of the graphene walls was assigned a fixed partial charge of ±0.003, ±0.006, and ±0.009 e corresponding to the surface charge density (Σ) of ±0.02, ±0.04, and ±0.06 Cm⁻², respectively. In a more realistic system, the potential difference between the surfaces is maintained constant, and in contrast to the fixed charged system, the former considers the fluctuation of surface charges induced by the local density fluctuations of electrolyte solutions.⁵² However, at lower potential differences, the ion and solvent density profiles were nearly comparable for both approaches.^52,53 Moreover, the approach where the potential in the wall is fixed, which is valid for metallic surfaces, has recently been called into question for semimetallic surfaces such as graphene.⁵⁴ 

The number density distribution of pure water oxygen and alkali chloride (LiCl, NaCl, and KCl) solutions normalized with respect to the bulk (center channel) value confined in a graphene nanochannel for Σ = 0.0 Cm⁻² and Σ = ± 0.06 Cm⁻² is shown in Figs. 2 and 3, respectively. When Σ = 0, the density distributions are symmetric across the nanochannel for all atomic species, with each species adopting a three-peak structure near the wall and a uniform region in the center. Among the different atomic species, the oxygen and the Cl⁻ ion density distributions do not change appreciably. Water molecule density peaks are located at ∼0.32 and 0.61 nm from the graphene electrode, consistent with previous observations.^55,56 Cl⁻ ions undergo repulsion from the neutral graphene surface, indicating the general repulsion behavior of strongly hydrated anions from graphene.⁵⁷ 

Regarding cations, Li⁺ ions exhibit major and minor adsorption peaks at distances of 0.40 and 0.74 nm from the wall, respectively, while Na⁺ ions display two minor adsorption peaks at 0.53 and 0.78 nm. A layer of water molecules is observed between the graphene electrodes and the Li⁺ and Na⁺ ionic layers, suggesting full solvation of both ions and residing in the outer Helmholtz plane (plane passing through the locus of centers of fully solvated ions). Li⁺ shows a higher concentration at the interface region with a peak intensity of 1.7 times the bulk value. The greater enhancement of Li⁺ is due to its ability to retain its primary hydration shell, thereby maintaining favorable cation–water electrostatic interactions in the interface region.^58,59 K⁺ ions form a major adsorption peak in the outer Helmholtz plane, with some of the K⁺ ions partly shedding the water molecules in their primary hydration shell and entering into the inner Helmholtz plane (the plane passing through the locus of centers of partially hydrated ions)—the partial dehydration is due to its weak interaction with the neighboring water molecules compared to other cations.

The asymmetric charges in the nanochannel walls influence the density distribution of water molecules. There is a noticeable change in the first density peak of water molecules at both interfaces, where a higher concentration of water molecules is observed near the positively charged surface [see Fig. 3(a)]. The asymmetry in the water density peak is attributed to the free energy change associated with the adsorption of each state of the water dipole.⁶⁰ According to the three-state water model,⁶¹ the water molecules near a surface exist in three possible forms: two types of monomers—flip up (oxygen atom facing the surface and hydrogen atom pointing toward the solution) and flop down (reverse of the flip up state), and one dimer—two associated monomers with no net dipole. The percentage of these water species depends upon the charge of the surface. For example, most water molecules in the interfacial region exist as dimers at a potential of zero charges.⁶⁰ As the surface charge increases from the potential of zero charges, more dimers convert to monomers. Figure 4 illustrates this, showing very limited interfacial water molecules oriented toward the wall at zero surface charge density. As the surface charge density increases, the water dipole orientation toward the wall increases for the negatively charged interface and reverses for the opposite side. When the surface becomes more positively charged, the preferential adsorption of flipped dipole increases as the free energy change due to adsorption associated with the flipped dipole is greater than its counterpart. That means that the number of water molecules with oxygen atoms pointing to the surface increases, leading to an increased density of water molecules near the positively charged surface.

The increase in the surface charge density results in the redistribution of the electrolyte particles in the interphase region, with the positively charged ions concentrating closer to the negatively charged surface (i.e., at z = 0) and vice versa. The position of the ion density peak remains the same during charging, while the respective peak intensities vary drastically. The peak position for the ion density remains consistent because the LJ interactions dominate during the trade-off between the Coulomb potential and the repulsive part of the LJ potential.⁶² The peak intensity of Cl⁻ increases to 3.8 times its bulk value near the positively charged electrode and remains consistent across all electrolyte solutions. Near the negatively charged electrode, the first density peak of K⁺ increases by fivefold, while Li⁺ and Na⁺ show only a twofold rise. This increase in K⁺ density is attributed to its tendency to partially dehydrate, leading to enhanced electrostatic interaction with the negatively charged graphene electrode and, consequently, higher ion concentration at the interface. The concentration of counterions near their respective electrodes drastically decreases, nearly reaching zero for K⁺ and Na⁺. However, in the LiCl solution, a notable concentration of Li⁺ remains near the positively charged graphene electrode due to its favorable interactions with water molecules in the interface region.^58,59

Here, we turn our attention to the influence of surface charge density on the slip friction coefficient in an asymmetrically charged graphene nanochannel. The slip friction coefficient is calculated by following the EMD method³⁸ based on the Hansen et al.³⁴ model. The slab width is a crucial parameter in determining the slip friction coefficient using this method. Finding the optimal slab width requires striking a balance between a size large enough to include the fluid layer experiencing slip over the surface and a size small enough to exclude the fluid that exhibits normal shear flow.³⁴ Niavarani and Priezjev⁶³ have highlighted the significance of the first density peak near the wall in determining the friction coefficient. Moreover, previous studies^64,65 also reported a strong correlation between the surface-induced structure of the first fluid layer and the friction coefficient. Hence, we chose a slab width of 0.5 nm from the wall, the distance over which the first density peak of fluid extends.

Figure 5 depicts the dependence of the slip friction coefficient on surface charge density for water and LiCl, NaCl, and KCl solutions confined in a graphene nanochannel with asymmetrically charged walls. Our calculated value of slip friction coefficient (1.74 ± 0.08 × 10⁴ kg m⁻² s⁻¹) for a water–graphene system at zero surface charge density closely matches the findings reported by Varghese et al.³⁸ As illustrated in Fig. 5, the slip friction gradually increases with increasing surface charge density for all systems investigated in our study. Furthermore, this increase in slip friction is similar for both walls of the asymmetrically charged graphene nanochannel, indicating that the differential adsorption of ions and water dipoles at both interfaces (as discussed in Sec. III A) has a negligible impact on slip.

Our findings on slip friction behavior for electrolyte solutions exhibit a trend akin to that observed in the work of Xie et al.³¹ Their simulated system involves aqueous NaCl confined in a symmetrically charged graphene nanochannel with homogeneously charged walls. They noted a gradual increase in friction with the increase in surface charge density and obtained identical results for positively and negatively charged nanopores. Furthermore, they developed a theoretical model that rationalizes their observed friction-charge behavior. However, the electric contribution to friction in their model was based on the Gouy–Chapman theory, which focuses on the effect of ions on charged interfaces and neglects the physical properties of the solvent despite its significant presence on the electrode surface. Consequently, an alternative approach is needed to elucidate the friction-charge behavior observed for pure water in our study.

Water molecules possess permanent dipole moments due to the differences in electronegativity among their constituent atoms. When a charged electrode is present, the generated electric field exerts a torque on the nearby water dipoles, causing them to align with the field direction. This field-induced ordering is counteracted by the thermally driven disorder, resulting in the partial orientation of the dipole along the field direction, called orientation polarization. The orientation polarization significantly affects the physical properties of water molecules,^66–68 with our focus here being on its influence on hydrodynamic properties. One such property is the viscosity of the liquid at the charged interface, which deviates from its bulk value due to the energy change associated with the orientation polarization of polar molecules.^69,70 Since slip friction results from the interaction between the first adsorbed water layer and the electrode surface, an important question arises whether the energy change due to the orientation polarization of water dipoles in the primary hydration sheath of the charged electrode modifies the slip friction coefficient. This question is addressed through theoretical reasoning.

Thus, we developed an expression that relates the slip friction coefficient to the electric field strength. According to Eq. (17), the slip friction coefficient increases proportionally to the square of the uniform electric field emanating from the charged electrode. This formulation offers an alternative approach to understanding slip friction behavior at charged interfaces by incorporating the physical properties of the solvent, differing from previous theoretical models.^29,31 The derived expression is analogous to the one for the viscoelectric effect,^69,70 which describes the variation in liquid viscosity perpendicular to the surface under an external electric field. This similarity is expected, as both phenomena stem from the energy change due to the orientation polarization of polar molecules at charged interfaces. In the viscoelectric effect, a strong electric field within the double layer interacts with dipolar molecules, increasing the viscosity of polar liquids near the charged surface. This phenomenon differs from the electroviscous effect, which refers to the viscosity change of a suspension of particles due to the presence of surface charge and the surrounding ionic double layer.^79,80 When a particle in suspension acquires a charge, an electric double layer forms around it, and under shear (from bulk fluid movement relative to the particles), this double layer deforms, leading to increased energy dissipation and viscosity.^81,82

To validate our simulation results and determine the proportionality constant, we plotted the normalized friction coefficient (normalized w.r.t. the unperturbed slip friction coefficient) against the square of the electric field strength, given by E_⊥ = $Σε0$⁠. We utilized the average value of the slip friction coefficient to calculate the normalized friction coefficient since it is similar for both interfaces. As shown in the inset of Fig. 5, the normalized friction coefficient increases linearly with the square of the electric field strength, as described in Eq. (17). The values of the proportionality constant, k, obtained from a linear fit to the simulation data for pure water and aqueous LiCl, NaCl, and KCl are 1.36, 1.28, 1.14, and 1.48 × 10⁻²⁰ m²/V², respectively. The presence of ions causes only minimal changes in the proportionality constant values, indicating that water polarization primarily determines the slip friction at charged interfaces. To fully understand the minimal variation in the proportionality constants requires a deeper investigation into their contributing factors.

In summary, our study utilized molecular dynamics simulations to evaluate the slip friction coefficient for pure water and aqueous electrolytes confined in a graphene nanochannel with asymmetrically charged walls. The slip friction coefficient is determined using the EMD method³⁸ based on the Hansen et al.³⁴ model. This method calculates the slip friction coefficient by considering the statistical fluctuations of the slab of particles near the surface. Our results reveal a gradual increase in the slip friction coefficient with an increase in the surface charge density of the nanochannel walls. Moreover, the slip friction coefficient is nearly identical for both walls of the nanochannel.

To elucidate the observed slip friction behavior, we developed an expression describing the effect of the electric field on the slip friction coefficient. According to this expression, the slip friction coefficient increases proportionally to the square of the uniform electric field emanating from the charged electrode. This increase in the slip friction coefficient stems from the energy change resulting from the orientation polarization of polar molecules, similar to the viscoelectric effect.^69,70 Unlike previous theoretical models,^29,31 our formulation incorporates the physical properties of the solvent, addressing a significant limitation. The proportionality constant in our expression, determined empirically for pure water and aqueous electrolytes, exhibits only minimal variations, suggesting that water polarization primarily governs the slip friction behavior. These findings align with a recent study on the electrolyte Seebeck effect, demonstrating almost the same thermovoltages for pure water and aqueous alkali halides confined between graphene electrodes.⁸³ Our study underscores the importance of water polarization in determining slip friction at charged interfaces and paves the way for exploring a similar behavior in other polar solvents.

The authors acknowledge the Swinburne supercomputing OzSTAR facility for providing computational resources for this work. The authors want to thank Dr. Sobin Alosious, Postdoctoral Research Associate at the University of Notre Dame, and Mr. Binu Varghese, Research Scholar at the Indian Institute of Technology—Madras, for the insightful discussions during the revision of the manuscript.

The authors have no conflicts to disclose.

Amith Kunhunni: Conceptualization (equal); Data curation (lead); Formal analysis (lead); Investigation (lead); Methodology (equal); Software (equal); Validation (equal); Visualization (lead); Writing – original draft (lead); Writing – review & editing (equal). Sleeba Varghese: Data curation (supporting); Formal analysis (supporting); Investigation (supporting); Methodology (supporting); Software (equal); Writing – review & editing (supporting). Sridhar Kumar Kannam: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Supervision (equal); Validation (equal); Writing – review & editing (equal). Sarith P. Sathian: Conceptualization (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Resources (equal); Supervision (equal); Validation (equal); Writing – review & editing (equal). Peter J. Daivis: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Project administration (equal); Supervision (equal); Validation (equal); Writing – review & editing (equal). B. D. Todd: Conceptualization (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Resources (equal); Supervision (equal); Validation (equal); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.